Yohanan Sivan

I disagree with Matheca. The function is intended to fit a general gaussian, not necessarily a probability distribution function. The equation is correct.

However, the user should be aware that removing data points in a deterministic manner (i.e. by thresholding) definitely skews the resulting fit.

Rather than fitting to the whole series with negatives removed, try finding the largest contiguous positive subset of the original data series and fitting to that. This method won't work when the noise amplitude is greater than the distribution amplitude, but in most cases it will give you a better fit.

Even better yet: if accuracy is more important than computation speed, use fmincon with a least-squares difference cost function:

p = fmincon(@(p) sum((y-(p(1)*exp(-(x-p(2)).^2/2/p(3)^2))).^2),...
[max(y) mean(x) 1],[],[],[],[],[0 -inf 0],[2*max(y) inf inf])
A=p(1);
mu=p(2);
sigma=p(3);

"modifying the title bar icon constitutes a breach of the license agreement"

=> http://www.mathworks.com/support/solutions/en/data/1-16N4J/?solution=1-16N4J

Not terribly useful. This is an example, not a fitting function and it is not very well written. It is, however, slightly more useful than the poor documentation on the function nlinfit.

The code does not seem to be performing a fit to a Lorentzian, though, as the function fun.m (why not give it a more useful name?) is multiplied by 2 for some reason and the (x-x_0) term in the denominator is multiplied by 4 for some reason. Possibly this is a standard form in some field I have never come across, but this is certainly not the Lorentzian function I am interested in fitting. It is easy enough to fix, but worth pointing out to others who might use this.