function [ params,perrors,chisqmin,exitflag,errormatrix ] = fitpchisq( x,y,model,guess,options )
%FITCHISQ Perform a Pearson's chi-squared fit of a statistical model to data.
%   PARAMS = FITCHISQ( x,y,model,guess,options) returns the values of the
%   parameters PARAMS at the global minimum of the chi-squared. A
%   statistical model MODEL, which is a function of PARAMS and the
%   control variable X (implemented as a so-called Anonymous function), 
%   is fitted to measurements Y(X) with uncertainties give by the square
%   root of the prediction, i.e. it assumes that we are fitting a binned
%   distribution of data by minimizing Pearson's chi-squared. 
%
%   In order for the fit to converge rapidly a
%   reasonable GUESS of the values of the parameters must be given.
%   OPTIONS to the minimization function that is used, FMINUNC, may be
%   passed. They should be set up with function OPTIMSET.
%
%   [PARAMS,PERRORS,CHISQMIN,EXITFLAG,OUTPUT,ERRORMATRIX] = fitchisq( x,y,dy,model,guess,options )
%   also returns the uncertainties on the parameters PERRORS, the chi-squared
%   at the global minimim CHISQMIN, the exit flag EXITFLAG from fminunc
%   (everything positive is good), some OUTPUT from fminunc, and the NxN 
%   error matrix for the parameters (N is the number of parameters).
% 
%   Example:
%         % Define the measurements x,y,error on y (call it dy)
%         x  = [0 1 2 3 4 5];
%         y  = [0 1 2 3 4 5];
%         dy = [1 1 1 1 1 1];
% 
%         % Define the model to perform a chi-squared fit of
%         guess = [1,1];
%         pname = cellstr(['Intercept';
%                          'Slope    ']);
%                         % 'Parabolic']);
%         model = @(p,x) p(1)+x*p(2); %+ x.^2*p(3);
% 
%         % Fit options (for fminunc)
%         options = optimset(options,'Display','off'); % Supress most output of minimization.
% 
%         % Do the fit and print/display some return values
%         [params,perrors,chisqmin,exitflag,errormatrix ] = fitchisq( x,y,dy,model,guess,options )
% 
%         % Plot the the data
%         errorbar(x,y,e,'r*')
% 
%         % Plot the fitted model
%         xfine = linspace(min(x),max(x),100);
%         yfit  = model(params,xfine);
%         hold on
%         plot(xfine,yfit,'b-')
%         legend('data','fit')
% 
%         % Print the values and errors of the fitted parameters
%         for i=1:length(guess)
%             fprintf('%9s: %f +- %f\n',char(pname(i)),params(i),perrors(i))
%         end


% Define the chi-squared function to minimize, taking the data and
% fit model from the arguments.

chisq = @(params) sum((y-model(params,x)).^2./model(params,x));

% Perform the minimization with fminunc.
options = optimset(options,'LargeScale','off'); % We aren't providing an external gradient function
[params,chisqmin,exitflag,OUTPUT,GRAD,HESSIAN] = fminunc(chisq,guess,options);

% The inverse of 1/2 of the hessian is the error matrix. The
% square roots of the diagonal elements are the error estimates
% for the parameters. In general there is a non-zero correlation...
%
% For whatever reason, a warning is printed about using the line-search
% algorithm because the gradient is not supplied - this can be ignored in
% general.

if exitflag >= 1
    errormatrix = inv(0.5*HESSIAN);
    perrors=sqrt(diag(errormatrix));
else
    params = guess;
    chisqmin = 999.;
    perrors = 999.;
    errormatrix = 999.;
    error('fitpchisq - fit using fminunc did not converge')
end

end