Code covered by the BSD License  

Highlights from
Generalized Nonlinear Non-analytic Chi-Square Fitting

image thumbnail

Generalized Nonlinear Non-analytic Chi-Square Fitting

by

 

11 Jan 2006 (Updated )

Performs chi-square fit with uncertainty estimation when measurement errors are known.

[params,dParams,gof,stddev]...
function [params,dParams,gof,stddev]...
    = fitChiSquare(xData,yData,modelFun,guess,dxData,dyData,options)
%FITCHISQUARE Performs a chi-square fit returning parameter uncertainty.
%
%   Usage:
%
%   [params,dParams,gof,stddev] =
%       fitChiSquare(xData,yData,modelFun,guess,dxData,dyData)
%
%       Performs a least-squares fit of the parameters of modelFun,
%       minimizing the chi-square between the model and the data, using the
%       specified uncertainties in x and y data.
%
%   ... = fitChiSquare(xData,yData,modelFun,guess,dxData,dyData,options)
%
%       As above, with the specified options.
%
%   ... = fitChiSquare(xData,yData,modelFun,guess)
%
%       Uses empty data uncertainty arrays, sets ErrorsUnknown to true.
%
%   Author:     N. Brahms
%   Contact:    Contact via Matlab File Exchange website
%   Version:    2.5
%   Updated:    5/26/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.  May be [] if
%                   Options.ErrorsUnknown = 1.
%       dyData  - A length m vector of uncertainties corresponding to each
%                   measurement in yData.  May be [] blank if
%                   Options.ErrorsUnknown = 1.
%       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 struct to pass to the fitting routine;
%                 this is the same as described in optimset, but with the
%                 following additional fields.  If any of the following
%                 fields are absent from the struct, the default value is
%                 used.  This parameter is optional.
%
%           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
%               Set this option to 1 if data uncertainties (dxData or dyData)
%               are unknown. fitChiSquare produces an error if
%               ErrorsUnknown is 0 and either dxData or dyData is empty.
%               If 0, then
%                   chi2 may be interpretred as a statistic indicating the
%                   likelihood that the model is appropriate for describing
%                   the observed data.  See definitions of chi2 and
%                   probChi2 above.  The delta chi2 = 1 bounds (dParams.dl
%                   and dParams.du) may be interpreted as 68% confidence
%                   intervals for the parameters (given that the parameters
%                   are gaussian-distributed).  stddev is given by
%                       stddev = sqrt(dy^2 + sum( (df/dxi*dxi)^2 ) )
%                   where dy and dxi are the data uncertainties and df/dxi
%                   is the finite-difference derivative of the model versus xi.
%               If 1, then
%                   chi2 is set to 1, and neither chi2 nor probChi2 carry
%                   meaning.  The delta chi2 = 1 bounds may be interpreted
%                   as "68% confidence intervals for the parameters given
%                   the model is correct."
%
%               If 1, assumes the model is correct (that is, chi2 is set to 1)
%               and uses the model to calculate the data standard deviation.
%               Assumes standard deviation is uniform by datum.
%               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.  Chi2 is a meaningful estimate of model
%               accuracy.  The delta chi2 = 1 bounds may be interpreted as
%               the "68% confidence interval" (if the parameters are
%               gaussian-distributed).
%           FitUncertainty = 0 | {1} | k-vector
%               If 1 (default), fits uncertainties
%               If 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];
%           FitIndex = {ones(size(guess))} | k-vector
%               fitChiSquare only fits each parameter if the corresponding
%               index in FitIndex is true.
%           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
%           MaxUncIter = {10}
%               maximum number of iterations in uncertainty solver
%           MaxUncCount = {20}
%               maximum number of function calls in uncertainty solver
%           UncOptions  = optimset options |
%                            {optimset('TolFun',1e-2,'Display','off')}
%               optimset options to pass to uncertainty solver (fzero)
%           Scale = {0} | 1
%               If 1, uses parameter scaling to make optimization radix
%               equally sensitive to each parameter.  This is useful if
%               parameters take on values different by more than an order
%               of magnitude (i.e. p = [1 1e10]).  Note that the scaling
%               factor is determined from the guess, so zero guesses are
%               not scaled.
%
%   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 of the indexed
%                 parameter
%           du  - absolute upper deviation at dChi^2=1 of the indexed
%                 parameter
%           lVal- a length k vector of the value of each parameter at the
%                 dChi^2=1 projection point for the lower deviation of the
%                 indexed paramter
%           uVal- a length k vector of the value of each parameter at the
%                 dChi^2=1 projection point for the upper deviation of the
%                 indexed parameter
%       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]
%           probChi2 - the probability of obtaining a chi-square equal to or
%                      greater than chi2 given the data is drawn from the
%                      model using the best-fit parameters
%       stddev  - the expected deviation from the fit for each observation
%
%   Example:
%   Run and edit fitChiSquareExample for example use.
%
%   Notes:
%   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).  Alternatively, it can be used to
%   determine measurement errors when the theoretical function is known.
%
%   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 and/or speed of
%   the calculation.
%
%   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

%--------------------------------------------------------------------------
%   (c) 2006 N. Brahms
%
%   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
%       2.4.2   - Added probchi2 field to gof struct
%       2.4.3   - M-Lint optimization
%       2.4.4   - Changed options handling to make code cleaner, added
%                 ability to fit only a variable subspace
%       2.5     - Added Scale option
%       2.5.1   - Help updates, options is now optional :)
%       2.6     - Added 4 parameter syntax
%
%   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

if ~exist('options','var')
    options = optimset;
end
options = initializeOptions(options,guess);
opEPar = initializeUncOptions(options);

if nargin==4
    options.ErrorsUnknown = true;
    dxData = [];
    dyData = [];
end

% 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);

% Scale input for unit radix
scale = ones(size(guess));
if options.Scale
    scale(guess~=0) = abs(guess(guess~=0));
    guess = guess./scale;
    if ~isempty(options.OutputFcn)
        unscaledOutputFcn = options.OutputFcn;
        options.OutputFcn = @(p,ov,s) feval(unscaledOutputFcn,p*scale,ov,s);
    end
end

% 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)*scale(ia),[options.LowerBound(ia) options.UpperBound(ia)],1)
        error('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)*ones(length(yData),1);
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;
gof.probChi2 = chi2prob(chiSquare,dof);

% Fits parameter uncertainties
% Be careful not to use ia in any nested functions!
dParams = initializeDParams;
dParams(1:length(params)) = dParams;
for ia=1:length(params)
    if options.FitUncertainty(ia)
        dParams(ia)=delta(ia);
    end
end

% Rescale data back from unity
params = params.*scale;
for ia=1:length(dParams)
    dParams(ia).du = dParams(ia).du*scale(ia);
    dParams(ia).dl = dParams(ia).dl*scale(ia);
    dParams(ia).d  = dParams(ia).d*scale(ia);
    dParams(ia).lVal = dParams(ia).lVal*scale(ia);
    dParams(ia).uVal = dParams(ia).uVal*scale(ia);
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 = %g + %g - %g',...
            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./scale).*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./scale),'.');
    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'));
    end
end

if options.Plot || options.PlotResiduals
    drawnow
end
% Last executing line of fitChiSquare    

% --- Internal functions follow ------------------------------------------

    % ----------- model --------------------------------------------------

    % This calls the model function.  It returns the residual at each
    % datum.
    function F = model(x)
        
        % If we scaled to unit radix, we must scale back before calling the
        % model.  Since scale was initialized to ones, we just multiply in
        % any case.
        x = x.*scale;
        
        % 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)   
        % See model
        x = x.*scale;
        
        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

    %-------------------- fzeroOutput ---------------------------
    % Used with fzero to stop when UncIter is exceeded
    function stop = fzeroOutput(x,o)
        
        if o.iteration > options.MaxUncIter || o.funccount > options.MaxUncCount
            stop = 1;
        else
            stop = 0;
        end
            
        % Display iterative function information
        if options.DisplayUncVal
            str1 = sprintf('Iteration = %d\tVarA = %6.4g\tVarB = %6.4g\tFValA = %6.4g\tFValB = %6.4g\tEVal = ',...
                o.iteration,o.intervala,o.intervalb,o.fvala,o.fvalb);
            str2 = sprintf('%g\t',pEVarVal);
            disp([str1 str2]);
        end
        
    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);
                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
            else
                F = (chiSquare+1) - sum( model(var).^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
    end % modelFix

    % ------------ reconstruct -------------------------------------------
    % Builds the parameter array from pEVar and pVar

    function x = reconstruct(in)
        x = zeros(size(params));
        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)
        xu = zeros(size(params));
        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)
            
            % 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);
            opParDev.OutputFcn = @(x,o,s) fzeroOutput(x,o);
            % Find the parameter deviation to get delta chi2 = 1
%             [bnd,fval,exitflag] = fzero(@modelFix,pFix,opParDev);
            bnd = fzero(@modelFix,pFix,opParDev);
            
            % Assign the result
            val = reconstruct(bnd);
                
            if options.DisplayUncVal
                disp(sprintf('Fix = %g\tVar = %g\tChi2 = %g',...
                    pFix,pVar,chiSquare));
            end
        end % uncMin

    end % delta

end % fitChiSquare

% ----- initializeUncOptions ------------------------------------------

function uncOp = initializeUncOptions(options)

% Set a relatively high finishing tolerance
uncOp = options.UncOptions;

end %initializeUncOptions
    
% ----- initializeOptions ---------------------------------------------

function options = initializeOptions(optionsIn,guess)

% Manufacture the fit options defaults
options = optimset;
options.Plot = false;
options.PlotResiduals = false;
options.PlotVariable = true;
options.FitUncertainty = true(1,length(guess));
options.FitIndex = true(1,length(guess));
options.ErrorsUnknown = false(1,length(guess));
options.UncOptions=optimset('TolFun',1e-2,'Display','off');
options.DisplayUncVal=false;
options.LowerBound=-Inf*ones(1,length(guess));
options.UpperBound=Inf*ones(1,length(guess));
options.Scale = false;
options.MaxUncIter = 10;
options.MaxUncCount = 20;

names = fieldnames(options);

for i = 1:length(names)
    if isfield(optionsIn,names{i})
        if length(optionsIn.(names{i})) == 1 &&...
                length(options.(names{i})) == length(guess) &&...
                length(guess)>1
            optionsIn.(names{i}) = optionsIn.(names{i})*ones(1,length(guess));
        end
        if islogical(options.(names{i}))
            options.(names{i}) = convertBoolean(optionsIn.(names{i}));
        else
            options.(names{i}) = optionsIn.(names{i});
        end %if islogical
    end %if isfield
end %for

end
    
% ----- chi2fun ----------------------------------------------------------

% Returns the probability of chi-square equal to or greater than x, with r
% degrees of freedom.  This is 1-D(x,r), where D is the cumulative
% distribution function.  See
% http://mathworld.wolfram.com/Chi-SquaredDistribution.html
% Note that MATLAB's gammainc is the "regularized gamma function" used in
% MathWorld.
function p = chi2prob(x,r)

p = gammainc(r/2,x/2);

end

% ----- trylsq ------------------------------------------------------------

% Uses lsqnonlin if it exists
function [sol,exitflag] = trylsq(fun,x0,lb,ub,options)

    if(exist('lsqnonlin','file'))
        [sol,trash,trash,exitflag] = lsqnonlin(fun,x0,lb,ub,options);
    else
        [sol,trash,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');
end
yData = yData(:);
siy = size(yData);

if doError
    sidy = size(dyData);
    if sidy(1)>1 && sidy(2)>1
        error('dY is not a vector');
    end
    dyData = dyData(:);
    sidy = size(yData);
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
    % Grow error vectors if the user passed scalar errors
    if sidx(1) == 1 && sidx(2) == six(2)
        for i=1:six(2)
            dxData(2:length(xData(:,1)),i)=dxData(1,i);
        end
        sidx = size(dxData);
    end
    if sidy(1) == 1 && sidy(2) == 1
        dyData(2:length(yData(:,1)),1)=dyData(1);
        sidy = size(dyData);
    end
    if any((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 MATLAB:sprintf:InputForPercentSIsNotOfClassChar;
errstr = ['Input vector' sprintf(' %s',in)...
            ' can not be converted to a boolean value'];
warning on MATLAB:sprintf:InputForPercentSIsNotOfClassChar;

% Handles cell arrays
if iscell(in)
    % All cells are numbers
    if all(isnumeric(cell2mat(in)))
        value = cell2mat(in) > 0;
    elseif islogical(cell2mat(in))
        value = cell2mat(in);
    % All cells are 'on' or not 'on'
    elseif iscellstr('in');
        in = lower(in);
        for i=1:length(in)
            switch in{i}
                case 'on'
                    value(i)=true;
                case 'off'
                    value(i)=false;
                otherwise
                    error(errstr);
            end %switch
        end %for
    else
        error(errstr);
    end %if number or string
% Array of numbers or booleans, or single string
else
    if all(isnumeric(in))
        value = in > 0;
    elseif all(islogical(in))
        value = in;
    elseif isa(in,'char') && all(isstrprop(in,'alpha'));
        switch lower(in)
            case 'on'
                value=true;
            case 'off'
                value=false;
            case 'true'
                value=true;
            case 'false'
                value=false;
            otherwise
                error(errstr);
        end %switch
    else
        error(errstr);
    end %if number or string
end %if cell

end %convertBoolean

Contact us