# Weird result of fitting a 2D Gauss to data

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.

• Similar question posted here. Dec 11 '20 at 21:57

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]


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.

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


• Thank you for your answer. Your and @corey979 's posts were both illuminating, but I accepted the latter one since it answered my question more directly. Nevertheless, thanks. Dec 13 '20 at 11:46