How to fit a gaussian to data in matlab/octave?

Question

I have a set of frequency data with peaks to which I need to fit a Gaussian curve and then get the full width half maximum from. The FWHM part I can do, I already have a code for that but I'm having trouble writing code to fit the Gaussian.

Does anyone know of any functions that'll do this for me or would be able to point me in the right direction? (I can do least squares fitting for lines and polynomials but I can't get it to work for gaussians)

Also it would be helpful if it was compatible with both Octave and Matlab as I have Octave at the moment but don't get access to Matlab until next week.

Any help would be greatly appreciated!

Imbecilic

Do you have a single peak (only 1 Gaussian)? Or multiple peaks (multiple, overlapping Guassians)? – Inca

It's just a single peak per file. – Imbecilic

If it's just one peak, take the mean and standard-dev of the numbers and that defines your sample normal distribution. Have you tried that? Otherwise, if you have the stats toolbox, use normfit(). – Hild

@Justin: Your first statement is wrong. If I have a bunch of data points around x=-100, with y-values corresponding to the standard normal there, and a bunch of similar values around x=-2, the mean of all those points will obviously not be zero, and the standard deviation will obviously not be unity. – Inca

Answer 1

Fitting a single 1D Gaussian directly is a non-linear fitting problem. You'll find ready-made implementations here, or here, or here for 2D, or here (if you have the statistics toolbox) (have you heard of Google? :)

Anyway, there might be a simpler solution. If you know for sure your data y will be well-described by a Gaussian, and is reasonably well-distributed over your entire x-range, you can linearize the problem (these are equations, not statements):

   y = 1/(σ·√(2π)) · exp( -½ ( (x-μ)/σ )² )
ln y = ln( 1/(σ·√(2π)) ) - ½ ( (x-μ)/σ )²
     = Px² + Qx + R         

where the substitutions

P = -1/(2σ²)
Q = +2μ/(2σ²)    
R = ln( 1/(σ·√(2π)) ) - ½(μ/σ)²

have been made. Now, solve for the linear system Ax=b with (these are Matlab statements):

% design matrix for least squares fit
xdata = xdata(:);
A = [xdata.^2,  xdata,  ones(size(xdata))]; 

% log of your data 
b = log(y(:));                  

% least-squares solution for x
x = A\b;                    

The vector x you found this way will equal

x == [P Q R]

which you then have to reverse-engineer to find the mean μ and the standard-deviation σ:

mu    = -x(2)/x(1)/2;
sigma = sqrt( -1/2/x(1) );

Which you can cross-check with x(3) == R (there should only be small differences).

Inca

Thanks very much. I had only been able to find the first of the links via google and that wasn't working with my data, the second one works a treat though. Also thanks for the explanation/equations. :D – Imbecilic

@user1806676: I haven't tried the linearized approach, but at least the math is correct. You should do some experimenting and validating there. – Inca

Tried the linearized approach. Works well. – Cella

+1. For linearizing and taking the anti-log after to find your coefficients. Good least squared error solution! – Konstantine

Answer 2

Perhaps this has the thing you are looking for? Not sure about compatability: http://www.mathworks.com/matlabcentral/fileexchange/11733-gaussian-curve-fit

From its documentation:

[sigma,mu,A]=mygaussfit(x,y) 
[sigma,mu,A]=mygaussfit(x,y,h)

this function is doing fit to the function 
y=A * exp( -(x-mu)^2 / (2*sigma^2) )

the fitting is been done by a polyfit 
the lan of the data.

h is the threshold which is the fraction 
from the maximum y height that the data 
is been taken from. 
h should be a number between 0-1. 
if h have not been taken it is set to be 0.2 
as default.

Festive

Better suited as a comment, no? – Imogeneimojean

While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. - From Review – Imogeneimojean

@ScottHoltzman Thanks for the heads up, I have included the relevant description. – Festive

Answer 3

i had similar problem. this was the first result on google, and some of the scripts linked here made my matlab crash.

finally i found here that matlab has built in fit function, that can fit Gaussians too.

it look like that:

>> v=-30:30;
>> fit(v', exp(-v.^2)', 'gauss1')

ans = 

   General model Gauss1:
   ans(x) =  a1*exp(-((x-b1)/c1)^2)
   Coefficients (with 95% confidence bounds):
      a1 =           1  (1, 1)
      b1 =  -8.489e-17  (-3.638e-12, 3.638e-12)
      c1 =           1  (1, 1)

Basifixed

Note that fit is not built-in; it's part of the curve fitting toolbox – Inca

Answer 4

I found that the MATLAB "fit" function was slow, and used "lsqcurvefit" with an inline Gaussian function. This is for fitting a Gaussian FUNCTION, if you just want to fit data to a Normal distribution, use "normfit."

Check it

% % Generate synthetic data (for example) % % %

    nPoints = 200;  binSize = 1/nPoints ; 
    fauxMean = 47 ;fauxStd = 8;
    faux = fauxStd.*randn(1,nPoints) + fauxMean; % REPLACE WITH YOUR ACTUAL DATA
    xaxis = 1:length(faux) ;fauxData = histc(faux,xaxis);

    yourData = fauxData; % replace with your actual distribution
    xAxis = 1:length(yourData) ; 

    gausFun = @(hms,x) hms(1) .* exp (-(x-hms(2)).^2 ./ (2*hms(3)^2)) ; % Gaussian FUNCTION

% % Provide estimates for initial conditions (for lsqcurvefit) % % 

    height_est = max(fauxData)*rand ; mean_est = fauxMean*rand; std_est=fauxStd*rand;
    x0 = [height_est;mean_est; std_est]; % parameters need to be in a single variable

    options=optimset('Display','off'); % avoid pesky messages from lsqcurvefit (optional)
    [params]=lsqcurvefit(gausFun,x0,xAxis,yourData,[],[],options); % meat and potatoes

    lsq_mean = params(2); lsq_std = params(3) ; % what you want

% % % Plot data with fit % % % 
    myFit = gausFun(params,xAxis);
    figure;hold on;plot(xAxis,yourData./sum(yourData),'k');
    plot(xAxis,myFit./sum(myFit),'r','linewidth',3) % normalization optional
    xlabel('Value');ylabel('Probability');legend('Data','Fit')

Kampmeier

Caveat: Requires the optimization toolbox since the code uses lsqcurvefit. – Konstantine

Answer 5

Allow me to add a simple function that does what you ask for - it seems to do the job for reasonably well-behaved data, but you want to make sure your select an appropriate range of X value to use (default = one point either side of the maximum, but you should probably go down to about 0.1 max. Should be easy to see how you can do that).

function [m,s,h]=gaussEstimate(x,y,n)
% fit gaussian to curve defined by x, y
% by taking log(y) and fitting a parabola to the max and points on either
% side (or optionally n points on each side, where n >= 1)
% the resulting fit values can be used with normpdf:
% fit = h * normpdf(xValues, m, s);

if ~exist('n','var')
    n = 1;
end

[~,mi] = max(y);
xi = mi+(-n:n);
ly = log(y(xi));
pf = polyfit(x(xi),ly,2);
m = x(mi)- pf(2)/(2*pf(1));
s = sqrt(-1/(2*pf(1)));
% fit to the bin at mi to get the amplitude "right"
% other methods could be imagined, obviously
h = y(mi)*s*sqrt(2*pi)/exp(-(m).^2/(2*s*s));

end

Moyer