hfitgauss

 [p,pe,prob,stat]=hfitgauss(x,n,lims,opt,freepar)

 A driver for matmin to make it as easy as possible
 to fit a Gaussian.
 Very similar to hfit except pre-calcs the starting
 values for Gaussian fit

 x,n is histogram data
 lims is optional x range to fit over
 opt is plotting options:
   'L' = log scale
   '0' = don't draw
 optional freepar specifies which parameters are free

 p,pe are parmeters and errors
 stat is the return status of the fit.
 The fit is done log likelihood style considering all bins.
 To allow a rough goodness of fit assessment chi-squared is
 calculated using non-empty bins only, and converted to the
 corresponding point prob on the chi2cdf.

 eg: [x,n]=hfill(randn(1,1000),100,-3,3);
     [p,pe]=hfitgauss(x,n)
     [p,pe]=hfitgauss(x,n,[-2,2],'L',[1,1,0])

 See also hfit

0001 function [p,pe,prob,stat]=hfitgauss(x,n,lims,opt,freepar)
0002 % [p,pe,prob,stat]=hfitgauss(x,n,lims,opt,freepar)
0003 %
0004 % A driver for matmin to make it as easy as possible
0005 % to fit a Gaussian.
0006 % Very similar to hfit except pre-calcs the starting
0007 % values for Gaussian fit
0008 %
0009 % x,n is histogram data
0010 % lims is optional x range to fit over
0011 % opt is plotting options:
0012 %   'L' = log scale
0013 %   '0' = don't draw
0014 % optional freepar specifies which parameters are free
0015 %
0016 % p,pe are parmeters and errors
0017 % stat is the return status of the fit.
0018 % The fit is done log likelihood style considering all bins.
0019 % To allow a rough goodness of fit assessment chi-squared is
0020 % calculated using non-empty bins only, and converted to the
0021 % corresponding point prob on the chi2cdf.
0022 %
0023 % eg: [x,n]=hfill(randn(1,1000),100,-3,3);
0024 %     [p,pe]=hfitgauss(x,n)
0025 %     [p,pe]=hfitgauss(x,n,[-2,2],'L',[1,1,0])
0026 %
0027 % See also hfit
0028 
0029 if(~exist('lims'))
0030   lims=[];
0031 end
0032 if(~exist('opt'))
0033   opt=' ';
0034 end
0035 if(~exist('freepar'))
0036   freepar=ones(1,3);
0037 end
0038 
0039 % vectorize data (if not already)
0040 x=x(:)';
0041 n=n(:)';
0042 
0043 if(isempty(lims))
0044   lims(1)=x(1);
0045   lims(2)=x(end);
0046 end
0047 
0048 ind=find(x>=lims(1)&x<=lims(2));
0049 xl=x(ind);
0050 nl=n(ind);
0051 
0052 % Calc starting values for fit
0053 inpar(1)=max(nl);
0054 N=sum(nl);
0055 if(freepar(2)==0)
0056   inpar(2)=0;
0057 else
0058   inpar(2)=(1/N)*sum(xl.*nl);
0059 end
0060 inpar(3)=sqrt((1/N)*sum((xl-inpar(2)).^2.*nl));
0061 
0062 % Do the fit
0063 [p,pe,gof,stat]=matmin('logl',inpar,freepar,'gauss',nl,xl);
0064 e=sqrt(nl);
0065 chi2 = chisq(p,'gauss',nl,e,xl);
0066 npt=length(e(e~=0));
0067 redchi=chi2/npt;
0068 prob=chi2cdf(chi2,npt);
0069 disp(sprintf('\nChiSq/npnt = %.2f / %.0f = %.2f (%.1f%s)\n\n',...
0070     chi2,npt,redchi,prob*100,'%'));
0071 
0072 if(~any(opt=='0'))
0073   hplot(x,n);
0074   
0075   span=x(end)-x(1); bw=x(2)-x(1);
0076   
0077   hold on;
0078   
0079   % Plot the fit function over the range used
0080   x2=xl(1)-bw/2:span/100:xl(end)+bw/2;
0081   plot(x2,gauss(p,x2),'r');
0082   
0083   % Plot the fit function over the full range dashed
0084   x2=x(1)-bw/2:span/100:xl(1)-bw/2;
0085   hold on; plot(x2,gauss(p,x2),'g');
0086   x2=xl(end)+bw/2:span/100:x(end)+bw/2;
0087   hold on; plot(x2,gauss(p,x2),'g');
0088   
0089   hold off;
0090   
0091   if(any(opt=='L'))
0092     ylim([3e-1,max(n)*2]);
0093     set(gca,'YScale','log');
0094   else
0095     ylim([0,max(n)*1.1]);
0096   end
0097 end