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

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 + g + g + ... + g


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

var_int = diff(I, g)**2 * var_g + diff(I, g)**2 * var_g + ...


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)

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

• Thanks for responding! Rather than solving it, I would like to understand the theory behind how it would work, as I need to do so for work. If you could share and explain your working to me, that would be much appreciated. I can upload the Python test script that I was working with to the end of this question, if that might be of interest to you? – TheEponymousProgrammer Nov 29 '19 at 15:14
• If you have the time, I would appreciate that, thanks! I've added the Python code I used to generate the test data at the bottom of the question, if you need to refer to it. – TheEponymousProgrammer Nov 29 '19 at 19:03