function [params,dParams,gof,stddev]...
= fitChiSquare(xData,yData,modelFun,guess,dxData,dyData,options)
%--------------------------------------------------------------------------
% [params,dParams,gof,stddev] =
% fitChiSquare(xData,yData,modelFun,guess,dxData,dyData,options)
%
% fitChiSquare is a generalized chi-square fitting routine for any
% model function when data measurement errors are known; it returns the
% model parameters and their uncertainties at the delta chi-squared = 1
% boundary (68% confidence interval).
%
% This function uses the optimization toolbox function lsqnonlin to first
% perform a nonlinear least squares fit of the data to the model. It
% then calculates the standard deviation of each point to the fit value,
% then calculates the chi squared. Finally, the function fzero is
% used on each parameter to find the value at which delta chi-squared is
% equal to 1 while minimizing chi^2 with respect to the other parameters.
% In the case that the parameter uncertainties are normally distributed,
% the delta chi^2 = 1 method gives the 68% confidence limit for the
% parameters. Monte Carlo or investigations of many data sets should be
% used to confirm the parameter uncertainties are normally distributed.
% This method gives upper and lower bounds for each parameter.[1,2]
%
% If one encounters the following error message:
%
% Unexpected termination flag 0 in non-estimating variable
% minimization during uncertainty estimation
%
% this is because the non-varying parameter minimization routine has
% encountered its iteration or evaluation limit. Raise
% options.UncOptions.MaxFunEvals or options.UncOptions.MaxIter.
%
% Note: If the user can not use lsqnonlin (i.e., the optimization
% toolbox is not installed), the program will use the built-in function
% fminsearch instead. This may reduce the robustness of the calculation.
%
% Author: N. Brahms
% Contact: Contact via Matlab File Exchange website
% Version: 2.1
% Updated: 1/25/06
%
% Inputs:
% xData - a m x n matrix of independent variable data, where
% columns represent independent variables, and rows
% represent successive observations thereof
% yData - a length m vector of dependent variable data
% dxData - a m x n matrix of uncertainties corresponding to each
% measurement in xData
% dxData - a length m vector of uncertainties corresponding to each
% measurement in yData
% modelFun- a pointer to a vectorized modeling function:
%
% [result] = modelFun(params,xData)
%
% which returns a length m vector of f(x) using the parameters
% passed
%
% guess - a length k vector of initial guesses for the parameters
% options - an options object to pass to the fitting routine;
% this is the same as described in optimset, but with the
% following additional fields (values in braces are defaults):
%
% Display = 'iter' | {'on'} | 'off'
% passed to the fitting routine. If set to 'on' or 'iter',
% also causes fitChiSquare to display the fitting results in
% the command window.
% Plot = {0} | 1
% plots the fit result. If PlotResiduals is also on, the fit
% result is plotted in subplot [1 2] and the residuals are
% plotted in subplot 3.
% PlotVariable = {1} | ... | k
% which independent variable to plot against.
% PlotResiduals = {0} | 1
% plots the weighted fit residuals
% ErrorsUnknown = {0} | 1
% if 1, assumes the model is correct and uses the model to
% calculate the data standard deviation. Assumes deviation
% is uniform.
% If using ErrorsUnknown, the delta chi2 = 1 bounds for the
% parameters can only be interpreted as the "68%
% confidence interval given the model is correct" (and then
% only if the parameters are gaussian-distributed).
% if 0, uses measurement errors to calculate the data
% standard deviation, and fits parameter uncertainties
% FitUncertainty = 0 | {1} | k-vector
% if 1 (default), fits uncertainties
% of 0, does not fit uncertainties (returns empty dparams)
% if length k boolean vector, fits uncertainties only to 'on'
% indices. Example for a three parameter fit:
% FitUncertianty = [0 1 1];
% LowerBound = vector | {[-Inf ... -Inf]}
% a length k row vector of parameter lower bounds
% UpperBound = vector | {[Inf ... Inf]}
% a length k row vector of parameter upper bounds
% DisplayUncVal = {0} | 1
% shows additional information when fitting the uncertainties
% UncOptions = optimset options |
% {optimset('TolFun',1e-2,'Display','off')}
% optimset options to pass to uncertainty solver (fzero)
%
% Outputs:
% params - a length k vector of fit model parameters
% dParams - a length k cell array of structures with the following
% fields:
% d - equals (dl+du)/2
% dl - absolute lower deviation at dChi^2=1
% du - absolute upper deviation at dChi^2=1
% lVal- a length k vector of the value of each parameter at the
% dChi^2=1 projection point for the lower deviation
% uVal- a length k vector of the value of each parameter at the
% dChi^2=1 projection point for the upper deviation
% gof - a goodness-of-fit struct with the following fields:
% chi2- the chi-square of the fit
% dof - degrees-of-freedom of the fit
% eta2- the eta-square of the fit (a.k.a. correlation index) -
% this quantity is analogous to the r-square (correlation
% coefficient) of a linear least-squares fitting [3]
% stddev - the expected deviation from the fit for each observation
%
% References:
% 1. W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling.
% Numerical Recipes; The Art of Scientific Computing. (Cambridge
% University Press: Cambridge). 1986.
% 2. P.R. Bevington, D.K. Robinson. Data Reduction and Error Analysis
% for the Physical Sciences. (McGraw-Hill: New York). 1992.
% 3. http://mathworld.wolfram.com
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%
% History
% 1.0 - Initial version
% 1.1 - Changed syntax and added the doError switch
% 1.1.1 - Added more detailed reporting
% 1.2 - Changed parameter uncertainty minimization function from
% fminbnd to fminsearch with coded constraints
% 1.3 - Changed non-estimating parameter minimization to occur
% inside fminsearch in uncertainty minimization. Also,
% fixed a bug causing improper weighting of initial chi2
% minimization. Removed doError switch.
% 1.3.1 - Minor code rearrangement
% 1.3.2 - Can run without the optimization toolbox
% 1.3.3 - Added UncMaxFunEvals & UncMaxIter options
% 1.3.4 - Fixed a bug whereby the uncertainty estimator would
% return bogus values instead of an error if the included
% non-varying parameter minimization routine halted
% unexpectedly.
% 2.0 - File name change to fitChiSquare. Added example use
% file and included between function in release. Added
% stddev output.
% 2.1 - Added bounds, uncertainty options. Updated
% documentation.
% 2.2 - Added ErrorsUnknown option, cleaned up documentation.
% 2.3 - Added plotting options, changed to allow FitUncertainty
% with ErrorsUnknown. Changed dParams to normal array.
% 2.4 - Added gof struct - note syntax change
% 2.4.1 - Fixed bug where struct was misassigned if
% FitUncertainty(1) was equal to 0
%
% Known problems / suggested features
% - Add flags indicating that the fit has stopped at bounds
%
%--------------------------------------------------------------------------
% Uncertainty loop constants and defaults
LARGE_NUM = realmax; % The number used in the minimization routine when
% the minimizing variable is out of bounds
UNC_OPT = optimset('TolFun',1e-2,'Display','off');
initializeOptions;
opEPar = initializeUncOptions(options);
% Check degrees of freedom
if (length(yData) <= length(guess))
error('Nonpositive degrees of freedom');
end
% Condition input
[xData, dxData, yData, dyData] = conditionInput(xData, dxData,...
yData, dyData, ~options.ErrorsUnknown);
% Force bound vectors to be row vectors
options.LowerBound = conditionBound(options.LowerBound);
options.UpperBound = conditionBound(options.UpperBound);
% Check guess is in bounds
for ia=1:length(guess)
if any(options.LowerBound>options.UpperBound)
error('Lower bound exceeds upper bound')
end
if ~between(guess(ia),[options.LowerBound(ia) options.UpperBound(ia)],1)
error(sprintf('Guess %d is not in bound',ia));
end
end
% Nested function variable initialization
pVarIndex = 0;
pVar = 0;
pFix = 0;
pEVarVal = zeros(length(guess)-1,1);
pEVarInd = zeros(length(guess)-1,1);
pVarBnd = [-Inf, Inf];
pVarMult = 1;
if ~strcmp(options.Display,'off')
disp('Fitting function parameters...');
end
% Calculate the degrees of freedom of the problem
dof = length(yData)-length(guess);
% Calculate the fit with uniform uncertainty
stddev = 1;
% Average deviation from the guess
stddev = sqrt(sum(model(guess).^2)/dof);
if ~options.ErrorsUnknown
startOptions = options; startOptions.TolFun = 1e-2; startOptions.TolX = 1e-8;
params = trylsq(@model,guess,[],[],startOptions);
% Calculate the standard deviation for each datum
stddev = calcStdDev;
% Calculate the least-squares fit
params = trylsq(@model,params,[],[],options);
else
% Calculate the least-squares fit
params = trylsq(@model,guess,[],[],options);
stddev = 1;
stddev = sqrt(sum(model(params).^2)/dof);
end
% Calculate the chi-square function
chiSquare = sum( model(params).^2 );
% Fill goodness-of-fit struct
gof.chi2 = chiSquare;
gof.dof = dof;
gof.eta2 = (std( (yData-model(params).*stddev) ) / std( yData) )^2;
% Fits parameter uncertainties
% Be careful not to use ia in any nested functions!
for ia=1:length(params)
if options.FitUncertainty(ia)
dParams(ia)=delta(ia);
else
% As of R14SP3, the fields below MUST be in the same order as they
% are defined in the delta function
dParams(ia) = initializeDParams;
end
end
% If a verbose mode is on, display the answer and reduced chi-squared
if ~strcmp(options.Display,'off')
for id=1:length(params)
disp(sprintf('Parameter %d = %f + %f - %f',...
id,params(id),dParams(id).du,dParams(id).dl));
end
disp(sprintf('Reduced chi-square = %f',chiSquare/dof));
end
% Plot the fit
if options.Plot
% Subplot if also plotting residuals
if options.PlotResiduals
subplot(3,1,[1 2]);
end
holdstate = ishold;
errorbar(xData(:,options.PlotVariable),yData,stddev,'.');
plotHandle = gca;
hold on
plot(xData(:,options.PlotVariable),...
yData - model(params).*stddev,'-r');
switch holdstate
case 0
hold off
end
xlabel('x1');
ylabel('y');
title('Fit plot');
lstring = sprintf('Fit: ~\\chi^2 = %0.2g\n',chiSquare/dof);
for ia = 1:length(params);
lstring = [lstring ...
sprintf('Par. %d = %0.2g + %0.1g - %0.1g\n',...
ia,params(ia),dParams(ia).du,dParams(ia).dl)];
end
legend({sprintf('Data \\pm \\sigma'),lstring(1:end-1)});
end
% Plots the weighted fit residual
if options.PlotResiduals
if options.Plot
subplot(3,1,3);
end
holdstate = ishold;
plot(xData(:,options.PlotVariable),model(params),'.');
hold on
line([min(xData(:,options.PlotVariable)) max(xData(:,options.PlotVariable))],...
[0 0],'Color',[1 0 0]);
switch holdstate
case 0
hold off
end
xlabel('x1');
ylabel('Residual');
title('Weighted residual plot');
if options.Plot
set(gca,'XLim',get(plotHandle,'XLim'));
% set(gca,'YLim',get(plotHandle,'YLim'));
end
end
drawnow
% --- Internal functions follow ------------------------------------------
% ----------- model --------------------------------------------------
% This calls the model function. It returns the residual at each
% datum.
function F = model(x)
% Handle boundary
testInBound = between(x,[options.LowerBound; options.UpperBound],1);
if ~all(testInBound)
% Return ridiculously large number
F = sqrt(LARGE_NUM)/(length(x)+1)*ones(length(xData),1);
else
% Return the chi-value
F = (yData - feval(modelFun,x,xData))./(stddev+eps);
%eps to avoid problems with zero weights
end
end
% ----------- modelError ---------------------------------------------
% Used in finding the point-by-point standard deviation
function F = modelError(x,xbndData)
F = feval(modelFun,x,xbndData);
end
% ---------- calcStdDev ------------------------------------------------
% Calculates the std. dev. of (f(x_i)-y_i) using finite methods
function sigma = calcStdDev
dim = size(xData);
dm2 = zeros(dim(1),1);
for ib=1:dim(2)
xtestu = xData;
xtestl = xData;
xtestu(:,ib) = xtestu(:,ib)+dxData(:,ib);
xtestl(:,ib) = xtestl(:,ib)-dxData(:,ib);
dmodelu = modelError(params,xtestu);
dmodell = modelError(params,xtestl);
dm2 = dm2 + (dmodelu-dmodell).^2/4;
end
sigma = sqrt(dyData.^2+dm2);
end
%----------------------- modelFix ---------------------------
% This function is used to fit only one variable to delta chi2 = 1. It
% is used with fzero, not lsqnonlin.
function F = modelFix(var)
if between(var,pVarBnd,1)
pVar = var;
if length(params>1)
[pEVarVal,eflag] = trylsq(@modelNFix,pEVarVal,[],[],opEPar);
end
if(eflag<=0)
disp(sprintf(['\tUnexpected termination flag %d in'...
' non-estimating variable minimization\n\tduring'...
' uncertainty estimation'],eflag));
F=NaN;
else
x = reconstruct(pVar);
F = (chiSquare+1) - sum( model(x).^2 );
end
% If past the variable-side bound, make F large. If past the
% extent-side bound, make F infinitely negative.
elseif var<min(pVarBnd)
F = pVarMult*LARGE_NUM;
eflag = 1;
elseif var>max(pVarBnd)
F = -pVarMult*LARGE_NUM;
eflag = 1;
end
% Display iterative function information
if options.DisplayUncVal
str1 = sprintf('Flag = %d\tVar = %6.4g\tFVal = %6.4g\tEVal = ',...
eflag,var,F);
str2 = sprintf('%g\t',pEVarVal);
disp([str1 str2]);
end
end
% ------------ reconstruct -------------------------------------------
% Builds the parameter array from pEVar and pVar
function x = reconstruct(in)
for ic=1:length(pEVarInd)
x(pEVarInd(ic)) = pEVarVal(ic);
end
x(pVarIndex) = in;
end
% ------------ modelNFix ---------------------------------------------
% This function is used to fit all but one variable
function F = modelNFix(x)
for ig=1:length(pEVarInd)
xu(pEVarInd(ig)) = x(ig);
end
xu(pVarIndex) = pVar;
F = model(xu);
end
% ------------ delta -------------------------------------------------
% This function finds the variation to delta chi-squared = 1
% One variable minimizes chi-square+1, while the other variables minimize
% chi-squared
function bnd = delta(varyingIndex)
bnd = initializeDParams;
% These are the shared variables for the fitting functions
pVarIndex = varyingIndex; % The index of the variable being solved
j=1; % Load all the other variables in an array
for ie=1:length(params)
if ie~=varyingIndex
pEVarVal(j)=params(ie);
pEVarInd(j)=ie;
j=j+1;
end
end
pEVarValStatic = pEVarVal; % The non-varying initial array value
pVar = params(pVarIndex); % The solution is stored here, but must
% start at the minimized value
% This is the original variable value
pFix = params(pVarIndex);
if ~strcmp(options.Display,'off')
disp(sprintf('Finding parameter %d lower bound',varyingIndex));
end
% Calculate the lower bound
[bnd.dl,bnd.lVal] = uncMin('lb');
bnd.dl = pFix-bnd.dl;
% Here we have to reinitiate the variables
pEVarVal = pEVarValStatic;
pVar = params(pVarIndex);
% Calculate the upper bound
if ~strcmp(options.Display,'off')
disp(sprintf('Finding parameter %d upper bound',varyingIndex));
end
[bnd.du,bnd.uVal] = uncMin('ub');
bnd.du = bnd.du-pFix;
bnd.d = mean([bnd.du bnd.dl]);
% This function finds the parameter deviation for delta chi-squared
% = 1, and then reminimizes the other variables.
function [bnd,val] = uncMin(lbub)
% Initialize loop variables
bnd = 0;
% Set the search bounds on the uncertainty
switch(lbub)
case 'ub'
pVarBnd = [pFix, options.UpperBound(varyingIndex)]; A=1;
case 'lb'
pVarBnd = [options.LowerBound(varyingIndex), pFix]; A=-1;
end
pVarMult = A;
opParDev = initializeUncOptions(options);
% Find the parameter deviation to get delta chi2 = 1
% [newDelta,fval] = fminsearch(@modelFix,pFix,opParDev);
[result,fval,exitflag] = fzero(@modelFix,pFix,opParDev);
% disp(sprintf('Exitflag = %d\tFval = %g',exitflag, fval));
% Assign the result
bnd = result;
val = reconstruct(bnd);
if options.DisplayUncVal
disp(sprintf('Fix = %g\tVar = %g\tChi2 = %g',...
pFix,pVar,chiSquare));
end
end
end
% ----- initializeUncOptions ------------------------------------------
function uncOp = initializeUncOptions(options)
% Set a relatively high finishing tolerance
uncOp = options.UncOptions;
end
% ----- initializeOptions ---------------------------------------------
function initializeOptions
% Manufacture the fit options if they don't exist yet
if ~exist('options','var') | isempty(options)
options=optimset;
end
% Plot options
try
options.Plot = convertBoolean(options.Plot);
catch
options.Plot = 0;
end
try
options.PlotResiduals = convertBoolean(options.PlotResiduals);
catch
options.PlotResiduals = 0;
end
try
options.PlotVariable;
catch
options.PlotVariable = 1;
end
% FitUncertainty
try
options.FitUncertainty;
catch
options.FitUncertainty=1;
fitUncertaintySpec = 0;
end
options.FitUncertainty = convertBoolean(options.FitUncertainty);
if(length(options.FitUncertainty) ~=1 &...
length(options.FitUncertainty) ~= length(guess))
disp(['Badly formed options.FitUncertainty vector, skipping '...
'uncertainty fit']);
options.FitUncertainty = zeros(1,length(guess));
elseif(length(options.FitUncertainty) == 1)
options.FitUncertainty = ones(1,length(guess))*...
options.FitUncertainty;
end
% ErrorsUnknown
try
options.ErrorsUnknown;
catch
options.ErrorsUnknown=0;
end
%UncOptions
try
options.UncOptions;
catch
options.UncOptions=optimset('TolFun',1e-2,...
'Display','off');
end
%DisplayUncVal
try
options.DisplayUncVal = convertBoolean(options.DisplayUncVal);
catch
options.DisplayUncVal=0;
end
%Bounds
try
options.LowerBound;
catch
options.LowerBound=-Inf*ones(1,length(guess));
end
try
options.UpperBound;
catch
options.UpperBound=Inf*ones(1,length(guess));
end
end
end
% ----- trylsq ------------------------------------------------------------
% Uses lsqnonlin if it exists
function [sol,exitflag] = trylsq(fun,x0,lb,ub,options)
if(exist('lsqnonlin','file'))
[sol,resnorm,residual,exitflag] = lsqnonlin(fun,x0,lb,ub,options);
else
[sol,fval,exitflag] = fminsearch(@modelChi2,x0,options);
end
% For use when substituting fminsearch for lsqnonlin
function F = modelChi2(x)
F = sum(feval(fun,x).^2);
end
end
% ----- between -----------------------------------------------------------
% Returns the portion of vector in and out of the range ( as determined by
% the inclusive tag )
function [inner,outer] = between(vector,range,inclusive);
if ~exist('inclusive','var')
inclusive = 0;
end
if inclusive
inner = vector<=max(range) & vector>=min(range);
outer = vector<=min(range) & vector>=max(range);
else
inner = vector<max(range) & vector>min(range);
outer = vector<min(range) & vector>max(range);
end
end
% ----- conditionInput ----------------------------------------------------
function [xData, dxData, yData, dyData] = conditionInput(xin, dxin, yin, dyin,...
doError)
xData = xin; dxData = dxin; yData = yin; dyData = dyin;
% Force yData and dyData to be column vectors
siy = size(yData);
if siy(1)>1 && siy(2)>1
error('Y is not a vector');
elseif siy(2)>1
yData = yData';
siy = size(yData);
end
if doError
sidy = size(dyData);
if sidy(1)>1 && sidy(2)>1
error('dY is not a vector');
elseif sidy(2)>1
dyData = dyData';
sidy = size(yData);
end
end
% Transpose xData if it is n x m, not m x n
six = size(xData);
if (siy(1) == six(2) & siy(1) ~= six(1))
xData = xData';
six = size(xData);
end
if doError
sidx = size(dxData);
if (siy(1) == sidx(2) & siy(1) ~= sidx(1))
dxData = dxData';
sidx = size(dxData);
end
end
% Check that vectors match
if (siy(1) ~= six(1))
error('X and Y data vectors do not match.')
end
if doError
if (sidx ~= six) | (sidy ~= siy)
error('Data vectors do not match error vectors');
end
end
end
% ----- initializeDParams -------------------------------------------------
function dParamsEmpty = initializeDParams()
dParamsEmpty = struct('dl', NaN, 'lVal', [], 'du', NaN, 'uVal', [],...
'd',NaN);
end
% ----- conditionBound ----------------------------------------------------
function boundVecOut = conditionBound(boundVecIn)
sib = size(boundVecIn);
if (sib(1)~=1 & sib(2)==1)
boundVecOut = boundVecIn';
elseif (sib(1)~=1)
error('Badly formed bound vector');
else
boundVecOut = boundVecIn;
end
end
% ----- convertBoolean ----------------------------------------------------
% Converts vectors of {>0, <=0} or {'on', 'off'} to
% vectors of {1, 0}
function value = convertBoolean(in)
warning off;
errstr = ['Input vector' sprintf(' %s',in)...
' can not be converted to a boolean value'];
warning on;
if iscell(in)
if all(isnumeric(cell2mat(in)))
value = cell2mat(in) > 0;
elseif all(isstrprop(cell2mat(in),'alpha'))
value = strcmp(in,'on');
else
error(errstr);
end
else
if all(isnumeric(in))
value = in > 0;
elseif isstrprop(in,'alpha');
switch in
case 'on'
value=1;
case 'off'
value=0;
otherwise
error(errstr);
end
else
error(errstr);
end
end
end
