hfit

 [p,pe,chisq,stat]=hfit(x,n,func,inpar,lims,opt,freepar)

 A driver for matmin for easy fitting of histograms
 using correct Poisson bin contents probability.

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

 p,pe are parmeters and errors
 redchi is not the quantity minimized and is for rough
 guidance only.
 stat is the return status of the fit.

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

 See also hfitgauss

0001 function [p,pe,redchi,stat]=hfit(x,n,func,inpar,lims,opt,freepar)
0002 % [p,pe,chisq,stat]=hfit(x,n,func,inpar,lims,opt,freepar)
0003 %
0004 % A driver for matmin for easy fitting of histograms
0005 % using correct Poisson bin contents probability.
0006 %
0007 % x,n is histogram data
0008 % func is function to fit
0009 % inpar are starting parameters of func
0010 % lims is optional x range to fit over
0011 % opt is optional 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 % redchi is not the quantity minimized and is for rough
0018 % guidance only.
0019 % stat is the return status of the fit.
0020 %
0021 % eg: [x,n]=hfill(randn(1,1000),100,-3,3);
0022 %     [p,pe]=hfit(x,n,'gauss',[30,0,1])
0023 %     [p,pe]=hfit(x,n,'gauss',[30,0,1.5],[-2,2],'L',[1,1,0])
0024 %
0025 % See also hfitgauss
0026 
0027 if(~exist('lims'))
0028   lims=[];
0029 end
0030 if(~exist('opt'))
0031   opt=' ';
0032 end
0033 if(~exist('freepar'))
0034   freepar=ones(size(inpar));
0035 end
0036 
0037 % vectorize data (if not already)
0038 x=x(:)';
0039 n=n(:)';
0040 
0041 if(isempty(lims))
0042   lims(1)=x(1);
0043   lims(2)=x(end);
0044 end
0045 
0046 ind=find(x>=lims(1)&x<=lims(2));
0047 xl=x(ind);
0048 nl=n(ind);
0049 
0050 % Do the fit
0051 [p,pe,gof,stat]=matmin('logl',inpar,freepar,func,nl,xl);
0052 e=sqrt(nl);
0053 chi2 = chisq(p,func,nl,e,xl);
0054 npt=length(e(e~=0));
0055 redchi=chi2/npt;
0056 disp(sprintf('\nChiSq/npnt = %.2f / %.0f = %.2f\n\n',chi2,npt,redchi));
0057 
0058 if(~any(opt=='0'))
0059   hplot(x,n);
0060   
0061   span=x(end)-x(1); bw=x(2)-x(1);
0062   
0063   hold on;
0064   
0065   % Plot the fit function over the range used in red
0066   x2=xl(1)-bw/2:span/100:xl(end)+bw/2;
0067   plot(x2,feval(func,p,x2),'r');
0068   
0069   % Plot the fit function over the full range in green
0070   x2=x(1)-bw/2:span/100:xl(1)-bw/2;
0071   hold on; plot(x2,feval(func,p,x2),'g');
0072   x2=xl(end)+bw/2:span/100:x(end)+bw/2;
0073   hold on; plot(x2,feval(func,p,x2),'g');
0074   
0075   hold off;
0076   
0077   if(any(opt=='L'))
0078     ylim([3e-1,max(n)*2]);
0079     set(gca,'YScale','log');
0080   else
0081     ylim([0,max(n)*1.1]);
0082   end
0083 end