Copyright (c) 2007, Yohanan Sivan
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
Yu-jiao Shen (view profile)
I am badly in need of this file.Thank you very much for your help.
Dan (view profile)
Felipe Santibañez (view profile)
Achut (view profile)
Georges (view profile)
This file poorly fits a gaussian (and has a much higher average rating) as compared to others here...
The test main file is wrong: the "rand" noise on the gaussian biases the measurements by 0.5, randn should be used instead.
Same main file used for comparison with another gaussian fit (gaussfit.m):
%% data
x=1:100;
sigma=15; mu=40; A=3;
plot(normpdf(x,mu,sigma))
y=A*exp(-(x-mu).^2/(2*sigma^2))+randn(size(x))*0.5;
hold all;
plot(x,y/(sum(y)),'.');
%% fitting
[sigmaNew,muNew,Anew]=mygaussfit(x,y);
y2=Anew*exp(-(x-muNew).^2/(2*sigmaNew^2));
plot(x,y2/(sum(y2)));
%% bis
[sigmaNew2,muNew2]=gaussfit(x,y/sum(y));
plot(normpdf(x,muNew2,sigmaNew2));
legend('truth','measurements','mygaussfit','gaussianfit')
Qing (view profile)
chile1987 (view profile)
Great idea,
though I wonder when you do the retransformation from second order polynomial (a0 + a1*x + a2*x^2)to logarithmic gaussian (log(A) - x^2/2*sigma^2 + x*mu/2*sigma^2 - mu^2/2*sigma^2)(line 46 in mygaussfit) shouldn't you calculate the mean from p(2) via mu = A1/-A2 or equivalently mu = A1*2*sqrt(-1/2*A2) instead of mu = A1*2*sigma^2 ?
I might be totally wrong or missed sth, just a quick idea.
Manuel A. Diaz (view profile)
Omer (view profile)
Namra Aftab (view profile)
hi...i have a sequence of 40 frames and i want to plot gaussian distribution of, say, first pixel of all 40 frames.the result does not look like a gaussian at all.what should i correct?
Muhammad khan (view profile)
hello every one any one can tell me about gussian curve fitting back groung why we use this instead of other curve fitting method what is the benift of this from other curve fitting method.If any one have some good data regarding gussian curve fitting kindly inform me.
Advance thx to all
Avi (view profile)
GREAT
Kiran (view profile)
hello...could someone please explain to me how the approximation from polynomial back to gaussian is done in the code? why are the coefficients equated in the way that they are below :
sigma=sqrt(-1/(2*A2));
mu=A1*sigma^2;
A=exp(A0+mu^2/(2*sigma^2));
Thank you very much
Jan (view profile)
Where do you get he formulas from? and where do you use h for? Can it also without h?
Thanks in advance
David Holz (view profile)
Bob Marley is absolutely correct. The use of a for loop is unnecessary and slows the function down by ~20%.
If you are confident in your inputs you can also remove the warning lines #74-85 of polyfit to knock out 87% of the total previous computation time.
Ryan (view profile)
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);
Ryan (view profile)
Michael Jordan (view profile)
an ‰Ê (view profile)
thanks a lot
The formula used for a Gaussian pdf is wrong. pdf(x)=(A/sqrt(2*sigma)) * exp( -(x-mu)^2 / (2*sigma^2) )...should be used.
The value of parameter h severely influence the result, in last comment, I use the default value, the fitting result is not correct, it looks more better when h = 0.1, how to solve this automatically?
there is a problem,
When I fit a data, for example,
y = [0.0651 0.0548 0.0461 0.0686 0.1268 0.2266 0.2292 0.1187 0.0299 0.0146 0.0092 0.0048 0.0032 0.0024];
it gives out a result like this:
yout = [0.0470 0.0594 0.0743 0.0918 0.1120 0.1352 0.1611 0.1897 0.2208 0.2538 0.2884 0.3238 0.3592 0.3937];
it is not a good fitting, how to solve this problem?
great, thanks!
Thank you very much!
This is exactly what I was looking for!
The code could be written in a more efficient manner -- i.e., using matlab syntax instead of the 'for' loop. Something like:
%% threshold
if nargin==2, h=0.2; end
%% cutting
indx = (y > max(y)*h);
%% fitting
p=polyfit(x(indx), log(y(indx)), 2);
sigma=sqrt(-1/(2*p(1)));
mu=p(2)*sigma^2;
A=exp(p(3)+mu^2/(2*sigma^2));
Thanks, this was a great help.
very helpful but I don't think you should use logx for calculation and you should just do the polyfit with ylog and x.
Very smart idea
works as advertised
working fine. how can i rewrite it for a lorentzian function of the form,
y=y0+(2*A/pi)*(w./(4*(x-x0).^2+w.^2))