Gaussian fit - Singular matrix in Invert_RtR

[Andrew Finlay]

Hello,

I'm trying to write code to automatically fit binned simulation data with a Gaussian function, but I'm getting a "Singular matrix in Invert_RtR" error when performing the fit.

I'm currently using gnuplot version 5.0.3 on Windows 8. Code and data are attached below.

Any help / guidance would be much appreciated.

Thanks.


gnuplot code:

reset
unset key
set xlabel "Time of Flight"
set ylabel "Counts"
A = 218.0
Mean = 3561.978563
StdDev = 0.008026
f(x) = A*exp(-(x-Mean)**2/(2*StdDev**2))
set fit errorvariables
fit [3561.851:3562.008] f(x) "histoTest.txt" u 1:2:3 yerror via A,Mean,StdDev
plot "histoTest.txt" u 1:2:3 w yerrorbars, f(x)


Contents of my data file (please note, the error values are the square root of the counts, but when this would be 0, I chose a value of 0.0001 to try to avoid issues):

# ToF Count Error
3561.851 1 1.000000
3561.855 0 0.000100
3561.859 0 0.000100
3561.863 0 0.000100
3561.867 0 0.000100
3561.871 0 0.000100
3561.875 0 0.000100
3561.878 0 0.000100
3561.882 0 0.000100
3561.886 0 0.000100
3561.890 1 1.000000
3561.894 2 1.414214
3561.898 0 0.000100
3561.902 2 1.414214
3561.906 3 1.732051
3561.910 2 1.414214
3561.914 2 1.414214
3561.918 3 1.732051
3561.922 6 2.449490
3561.926 6 2.449490
3561.930 4 2.000000
3561.933 10 3.162278
3561.937 8 2.828427
3561.941 13 3.605551
3561.945 10 3.162278
3561.949 11 3.316625
3561.953 26 5.099020
3561.957 24 4.898979
3561.961 41 6.403124
3561.965 55 7.416198
3561.969 58 7.615773
3561.973 110 10.488088
3561.977 218 14.764823
3561.981 170 13.038405
3561.984 49 7.000000
3561.988 17 4.123106
3561.992 3 1.732051
3561.996 1 1.000000
3562.000 3 1.732051
3562.004 3 1.732051
3562.008 0 0.000100

[Hans-Bernhard Bröker]

Am 27.04.2016 um 19:38 schrieb Andrew Finlay:

> I'm trying to write code to automatically fit binned simulation data

Bad idea. 'fit' is not a tool to be used blindly, or automatically.

Non-linear least-squares fitting needs guidance and consideration of the
condition of the input, or you'll run into numerical instabilities much
of the time.

So let's see why this particular fit failed. Well, for starters the
data is quite clearly not Gaussian at all. It's highly asymmetrical,
and the tails are way to wide for compared to the FWHM.

Trying to fit a Gaussian to that data is an exercise in futility.

[Andrew Finlay]

On Wednesday, April 27, 2016 at 1:29:37 PM UTC-7, Hans-Bernhard Bröker wrote:
> Am 27.04.2016 um 19:38 schrieb Andrew Finlay:
>
> > I'm trying to write code to automatically fit binned simulation data
>
> Bad idea. 'fit' is not a tool to be used blindly, or automatically.
>
> Non-linear least-squares fitting needs guidance and consideration of the
> condition of the input, or you'll run into numerical instabilities much
> of the time.

I understand where you're coming from here, but I thought my initial parameters were close enough that the algorithm should have been able to converge, though maybe I was wrong. And they weren't chosen blindly, the starting parameters are based on the mode and the standard deviation of the unbinned data points. I played around with those parameters some, I got rid of the "Singular matrix in Invert_RtR" error, but now it just doesn't modify my parameters at all. I don't know if I've simply hit the limits of what fit can do. Here's a link to an image of what I came up with anyway: https://drive.google.com/open?id=0B2-M7wNoiJi2bUF3MTkwMmN4bFE

> So let's see why this particular fit failed. Well, for starters the
> data is quite clearly not Gaussian at all. It's highly asymmetrical,
> and the tails are way to wide for compared to the FWHM.
>
> Trying to fit a Gaussian to that data is an exercise in futility.

I would argue that the peak is Gaussian to first order approximation. It's not truly Gaussian obviously, but you should still be able to produce an approximate fit. Though perhaps not in gnuplot...

Thanks anyway.

[Hans-Bernhard Bröker]

It's not. And if it were, the first order of business would have been
to limit the fit to the region of the peak. The dataset has 40 bins,
but at most 6 of them could possibly be said to look Gaussian. Fitting
3 parameters to those 6 inputs might appear to work, and gnuplot will
eventually do it. It still doesn't make a useful fit, though.

Oh, and two more things: The starting value of StdDev appears about a
factor of 2 off.

And for a StdDev around 0.005, a bin width of 0.004 can't really work,
either.

[Stephan Böttcher]

Andrew Finlay <andrew...@gmail.com> writes:

> Contents of my data file (please note, the error values are the square
> root of the counts, but when this would be 0, I chose a value of
> 0.0001 to try to avoid issues):

0.0001 is too low, you weigh those bins to much. I usually do

using 1:2:(sqrt($2+1))

Zero is quite probably for bins with expectation value 1

Stephan

--
Stephan Böttcher Tel: +49-431-880-2508
Extraterrestrische Physik, IEAP, Leibnizstr. 11, 24118 Kiel