5
$\begingroup$

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

Fit result

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.

raw data plot interpolated plot

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.

$\endgroup$
1
  • 1
    $\begingroup$ Similar question posted here. $\endgroup$ Dec 11 '20 at 21:57
10
$\begingroup$

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]

enter image description here

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

fit["ParameterTable"]

enter image description here

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

enter image description here

$\endgroup$
1
  • $\begingroup$ Thank you for your answer $\endgroup$ Dec 13 '20 at 11:43
7
$\begingroup$

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

enter image description here

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Dec 13 '20 at 11:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.