Finding the error in the total integrated intensity of a fitted 2D Gaussian

Question

I have been trying to fit signals to a 2D Gaussian function, and while I have bene able to use sciKit-image's curve_fit function to find the covariance matrix for the parameters, I am at a loss as to how to estimate the error of the measured Gaussian, as well as the error in the total integrated intensity.

The current parameters I have for my 2D Gaussian are:

• amplitude - The scaling factor,
• mu_x and mu_y, the centres of the Gaussian along the x- and y-axes,
• sig_u and sig_v, the width of the Gaussian as defined along and perpendicular to its diagonal,
• theta, the degree of rotation, and
• offset, the offset from 0.

Based on my reading, I have tried a couple of methods to quantify the error in the Gaussian (fitted over a 64 x 64 pixel grid, for reference), and would like advice as to which is more appopriate (or if none of them are):

1. Squared Difference Between Raw and Fitted Signals

To start with, I have taken the variance of the Gaussian to be the squared difference between the raw signal and the fitted signal for each pixel. The equation used is as follows:

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

Then, in order to find the total integrated intensity, I sum the model array:

# 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]

From my understanding, the formula for the combination of errors would then be:

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

Where diff(I, g) means differentiating the intensity function with respect to each individual point along the Gaussian function. Since it's a linear function, my understanding is that diff(I, g) is just one, and var_int is therefore:

var_int = var_g.sum()

2. Combination of Errors

The other method I tried involved taking the parameters from the 2D Gaussian and using them to calculate the variance the Gaussian for each pixel, again using the law of combination of errors:

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

Again, because the total integrated signal is a linear equation in terms of the Gaussian function, the variance of the total integrated signal is then just:

var_int = var_g.sum()

Which method would be more appropriate in this case, and if not, what suggestions/changes would you recommend I make?

Thanks!

UPDATE: The code I used to generate the 2D Gaussian is as follows.

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)