Weird result of fitting a 2D Gauss to data

Question

I wanted to fit a 2D Gaussian function $$f(x,y) = \frac{1}{2 \pi \sigma^2} \exp \left(-\frac{(x-\mu_1)^2 + (y-\mu_2)^2}{2\sigma^2}\right),$$ to a set of points prepared from an image file (triplets $\{x_i, y_j, f_{ij}\}$), and NonlinearModelFit returns quite a confusing fit.

I used the following code

v = Import["https://pastebin.com/raw/6JuygfJ1", "Table"];
f[x_, y_] := 1/(2 Pi σ^2) Exp[
        -((x - Subscript[μ, 1])^2 + (y - Subscript[μ, 2])^2)/(2 σ^2)
];
fit = NonlinearModelFit[
        v, f[x, y], { Subscript[μ, 1], Subscript[μ, 2], σ}, {x,y}
];

which results in something like this

My problem here is that it even though the $\mu_1, \mu_2$ values look promising, the $\sigma$ is way off. For comparison, here are plots of the data itself, and the data interpolated.

The data comes from an image of a laser beam coming from an optic fiber. There is some visible interference pattern, but I have no idea if that is the problem.

Any help in finding out the strange behaviour will be greatly appreciated.

Answer 1

The function you are fitting can never account properly for your data: your function is a probability density function (PDF) of a bivariate Gaussian distribution – as such, the volume under it (the integral of $x$ and $y$ over the whole plane) is by definition of the PDF equal to unity. Without going into any interpolation, integration, or any fancy estimation – the volume covered by your data can be estimated as follows: the height is of order of $h=1$, and making a very crude approximation of the profile as a triangle with base $2r=100$, you get a cone with volume $\frac{1}{3}\pi r^2 h = \frac{1}{3}\pi\cdot 50^2\cdot 1 \approx 2618$ – that is not even close to unity.

Therefore, you need to fit the function $A\cdot f(x,y)$, not just $f(x,y)$. And providing reasonable starting values is generally a good habit, too (if not provided explicitly, Mathematica assumes all starting values are 1):

f[x_, y_] := Exp[-(((x - μ1)^2 + (y - μ2)^2)/(2*σ^2))]/(2*Pi*σ^2)

ListPointPlot3D[v]

fit = NonlinearModelFit[v, A f[x, y], {{A, 1000}, {μ1, 50}, {μ2, 100}, {σ, 10}}, {x, y}]

fit["ParameterTable"]

Plot3D[fit[x, y], {x, 0, 150}, {y, 0, 150}, PlotRange -> All]

Answer 2

Maybe you do not need to fit, just integrate. Analytically we have

Integrate[x f[x, y], {x, -∞, ∞}, {y, -∞, -∞}, Assumptions -> σ > 0]
(*μ1*)
Integrate[y f[x, y], {x, -∞, ∞}, {y, -∞, -∞}, Assumptions -> σ > 0]
(*μ2*)
Integrate[((x - μ1)^2 + (y - μ2)^2)f[x, y],{x, -∞, ∞},{y, -∞, -∞},Assumptions -> σ > 0]
(*2 σ^2*)

This is sufficient to get the parameters you need by numerical integration of your data.

Using your data we obtain

 n = Total[v[[All, 3]]];
μ1 = Total[v[[All, 1]] v[[All, 3]]/n];
μ2 = Total[v[[All, 2]] v[[All, 3]]/n]
mxy = Total[((v[[All, 1]] - mx)^2 + (v[[All, 2]] - my)^2) v[[All, 3]]/n];
σ = Sqrt[mxy/2];

{n, mx, my, σ}
(*{2284.83, 43.3037, 97.2928, 23.0803}*)

f[x_, y_] := n/(2 Pi σ^2) Exp[-((x - μ1)^2 + (y - μ2)^2)/(2 σ^2)];
GraphicsRow[{ListDensityPlot[v], DensityPlot[f[x, y], {x, 1, 128}, {y, 1, 160}]}]