我一直在尝试将信号拟合到 2D 高斯函数，虽然我能够使用 sciKit-image 的curve_fit函数来找到参数的协方差矩阵，但我不知道如何估计测量的高斯函数的误差，以及总积分强度的误差。

我的二维高斯的当前参数是：

• amplitude- 比例因子，
• mu_x和mu_y，沿 x 轴和 y 轴的高斯中心，
• sig_u和sig_v，高斯的宽度，定义为沿和垂直于其对角线，
• theta, 旋转度数, 和
• offset，从 0 开始的偏移量。

根据我的阅读，我尝试了几种方法来量化高斯中的误差（适合 64 x 64 像素网格，供参考），并想知道哪种方法更合适（或者如果它们都不是):

1. 原始信号和拟合信号之间的平方差

首先，我将高斯方差作为每个像素的原始信号和拟合信号之间的平方差。使用的方程如下：

# Variance for each pixel
var_g = (signal - model)**2  # returns a 64 x 64 array with variances at each point

然后，为了找到总积分强度，我对模型数组求和：

# Calculating the total integrated intensity
int_sum = model.sum()  # This is summed over 64 x 64 = 4096 pixels
int_sum = g[0][0] + g[0][1] + g[0][2] + ... + g[63][63]

据我了解，错误组合的公式为：

var_int = diff(I, g)**2 * var_g[0][0] + diff(I, g)**2 * var_g[0][1] + ...

其中diff(I, g)意味着对沿高斯函数的每个单独点的强度函数进行微分。由于它是一个线性函数，我的理解是这diff(I, g)只是一个，var_int因此是：

var_int = var_g.sum()

2. 错误组合

我尝试的另一种方法涉及从 2D 高斯获取参数并使用它们来计算每个像素的高斯方差，再次使用误差组合定律：

# Variance at each pixel
var_g = diff(g, amp)**2 * var_amp + diff(g, mu_x)**2 * mu_x + ... + diff(g, offset)**2 * offset

同样，因为总积分信号是高斯函数的线性方程，所以总积分信号的方差为：

var_int = var_g.sum()

在这种情况下哪种方法更合适，如果不是，您会建议我做出哪些建议/更改？

谢谢！

更新：我用来生成二维高斯的代码如下。

import numpy as np

def Gaussian2D_v5(coords=None,  # x and y coordinates for each image.
                  amplitude=1,  # Highest intensity in image.
                  xo=0,  # x-coordinate of peak centre.
                  yo=0,  # y-coordinate of peak centre.
                  sig_u=1,  # Standard deviation in u.
                  sig_v=1,  # Standard deviation in v.
                  theta=0,  # Orientation of the diagonal matrix, in degrees
                  offset=0):  # Offset from zero (background radiation).

    # Set x and y coordinates
    x, y = coords
    # Convert numbers to floats
    xo = float(xo)
    yo = float(yo)

    # Convert theta to radians
    theta = np.deg2rad(theta)

    # Define diagonal matrix
    mat_diag = np.asarray([[sig_u**2, 0],
                           [0, sig_v**2]])
    # mat_diag = np.asarray(mat_diag)
    # Define rotation matrix
    mat_rot = np.asarray([[np.cos(theta), -np.sin(theta)],
                          [np.sin(theta), np.cos(theta)]])
    # mat_rot = np.asarray(mat_rot)
    # Define covariance matrix
    mat_cov = np.matmul(np.matmul(mat_rot, mat_diag),
                        np.linalg.inv(mat_rot))

    # Set the coordinates and stack them along the last axis
    mat_coords = np.stack((x - xo, y - yo), axis=-1)
    mat_coords_row = mat_coords[:, :, np.newaxis, :]
    mat_coords_col = mat_coords[:, :, :, np.newaxis]

    # Perform calculation across the meshgrid
    G = amplitude * np.exp(-0.5 * np.matmul(np.matmul(mat_coords_row,
                                                      np.linalg.inv(mat_cov)),
                                            mat_coords_col)) + offset
    # Reduce the dimensions of the answer
    return G.squeeze().ravel()

# Generating the noisy Gaussian
coords = np.meshgrid(np.arange(0, 64), np.arange(0, 64))
model = Gaussian2D_v5(coords,
                      amplitude=20,
                      xo=32,
                      yo=32,
                      sig_u=6,
                      sig_v=3,
                      theta=20,
                      offset=20).reshape(64, 64)
noise = np.random.normal(0,         # mean
                         10,        # standard deviation
                         (64, 64))  # array shape
signal = model + noise

# Curve fitting
# Set initial guess
p0 = [20,  # amplitude
      32,  # xo
      32,  # yo
      2,   # sig_u
      2,   # sig_v
      0,   # theta
      0]   # offset
# Set upper and lower bounds
bound_lo = [0, 0, 0, 0, 0, 0, 0]
bound_hi = [100, 64, 64, np.inf, np.inf, 90, np.inf]

bounds = (bound_lo, bound_hi)

# Fit to 2D Gaussian equation
popt, pcov = curve_fit(Gaussian2D_v5,
                       coords,
                       signal.ravel(),
                       p0=p0,
                       bounds=bounds)

# Generated fitted curve
model = Gaussian2D_v5(coords, *popt).reshape(64, 64)