Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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

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!

like image 714
user1806676 Avatar asked Nov 08 '12 14:11

user1806676


1 Answers

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')
like image 72
MJRunfeldt Avatar answered Oct 03 '22 06:10

MJRunfeldt