二维高斯拟合20240815

二维高斯拟合

高斯函数表达式

二维高斯函数是一个在二维空间中用来表示高斯分布的函数，常用于统计学、图像处理和机器学习等领域。其数学表达式通常为：

$$f(x,y)=\frac{1}{2\pi {\sigma }_x{\sigma }_y\sqrt{1-{\rho }^2}}\exp \left(-\frac{1}{2(1-{\rho }^2)}\left(\frac{(x-{\mu }_x{)}^2}{{\sigma }_x^2}-\frac{2\rho (x-{\mu }_x)(y-{\mu }_y)}{{\sigma }_x{\sigma }_y}+\frac{(y-{\mu }_y{)}^2}{{\sigma }_y^2}\right)\right)$$

其中，( (x,y) ) 是二维平面上的点，( ${\mu }_x$ ) 和 ( ${\mu }_y$ ) 是分布的均值，( ${\sigma }_x$ ) 和 ( ${\sigma }_y$ ) 是两个方向的标准差，( $\rho$ ) 是两个维度之间的相关系数。这个表达式描述了在二维平面上的高斯分布[函数]

方案一 Matlab：

MATLAB 提供多种方法来计算二维高斯函数的积分。

方案二 Python

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import pandas as pd

# 读取数据
# df = pd.read_excel(path)

# 生成数据
x = np.linspace(-2, 2, 30)
y = np.linspace(-2, 2, 30)
X, Y = np.meshgrid(x, y)
Z = np.exp(-X**2 - Y**2)

# 添加随机噪音
Z_noisy = Z + 0.05*np.random.randn(*Z.shape)

# 二维高斯模型
def TwoD_gaussian(xy, a, x0, y0, sigma_x, sigma_y):
    """
    
    """
    x, y = xy
    z=a*np.exp(-((x-x0)**2/(2*sigma_x**2) + (y-y0)**2/(2*sigma_y**2)))
    return z.ravel()

# 初始参数
a0 = [1, 0, 0, 1, 1]

# 拟合
popt, pcov = curve_fit(TwoD_gaussian, (X, Y), Z_noisy.ravel(), p0=a0)

# 绘图
fig = plt.figure()
ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.plot_surface(X, Y, Z)
ax.set_title('Original')
ax = fig.add_subplot(1, 2, 2, projection='3d')
ax.plot_surface(X, Y, TwoD_gaussian((X, Y), *popt).reshape(X.shape))
ax.set_title('Fitted')
plt.show()

程序运行后图像如下

如下代码有写问题还需要改进。

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

def gaussian_2d(x, y, amplitude, xo, yo, sigma_x, sigma_y, theta):
    a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
    b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
    c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
    return amplitude * np.exp(-(a*(x - xo)**2 + 2*b*(x - xo)*(y - yo) + c*(y - yo)**2))

def generate_data():
    # 生成示例数据
    x = np.linspace(0, 10, 100)
    y = np.linspace(0, 10, 100)
    x, y = np.meshgrid(x, y)
    amplitude = 3.0
    xo = 5.0
    yo = 5.0
    sigma_x = 2.0
    sigma_y = 2.0
    theta = np.pi / 4
    z = gaussian_2d(x, y, amplitude, xo, yo, sigma_x, sigma_y, theta)
    z += 0.2 * np.random.randn(x.shape[0], x.shape[1])  # 添加一些噪声
    return x, y, z

def fit_gaussian(x, y, z):
    # 初始猜测参数
    initial_guess = [3.0, 5.0, 5.0, 2.0, 2.0, np.pi / 4]
    popt, pcov = curve_fit(gaussian_2d, np.vstack((x.ravel(), y.ravel())), z.ravel(), p0=initial_guess)
    return popt

def plot_fit(x, y, z, popt):
    # 绘制原始数据和拟合结果
    fig, ax = plt.subplots(1, 1)
    ax.imshow(z, origin='lower', extent=[x.min(), x.max(), y.min(), y.max()])
    ax.contour(x, y, gaussian_2d(x, y, *popt), colors='r')
    plt.show()

if __name__ == "__main__":
    x, y, z = generate_data()
    popt = fit_gaussian(x, y, z)
    plot_fit(x, y, z, popt)