nlinfit

Nonlinear regression

Syntax

  • beta = nlinfit(X,Y,modelfun,beta0) example
  • beta = nlinfit(X,Y,modelfun,beta0,options) example
  • beta = nlinfit(___,Name,Value) example
  • [beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(___) example

Description

example

beta = nlinfit(X,Y,modelfun,beta0) returns a vector of estimated coefficients for the nonlinear regression of the responses in Y on the predictors in X using the model specified by modelfun. The coefficients are estimated using iterative least squares estimation, with initial values specified by beta0.

example

beta = nlinfit(X,Y,modelfun,beta0,options) fits the nonlinear regression using the algorithm control parameters in the structure options. You can return any of the output arguments in the previous syntaxes.

example

beta = nlinfit(___,Name,Value) uses additional options specified by one or more name-value pair arguments. For example, you can specify observation weights or a nonconstant error model. You can use any of the input arguments in the previous syntaxes.

example

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(___) additionally returns the residuals, R, the Jacobian of modelfun, J, the estimated variance-covariance matrix for the estimated coefficients, CovB, an estimate of the variance of the error term, MSE, and a structure containing details about the error model, ErrorModelInfo.

Examples

expand all

Nonlinear Regression Model Using Default Options

Load sample data.

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

Fit the Hougen-Watson model to the rate data using the initial values in beta0.

beta = nlinfit(X,y,@hougen,beta0)
beta =

    1.2526
    0.0628
    0.0400
    0.1124
    1.1914

Nonlinear Regression Using Robust Options

Generate sample data from the nonlinear regression model

y=b1+b2exp{b3x}+ε,

where b1, b2, and b3 are coefficients, and the error term is normally distributed with mean 0 and standard deviation 0.1.

modelfun = @(b,x)(b(1)+b(2)*exp(-b(3)*x));

rng('default') % for reproducibility
b = [1;3;2];
x = exprnd(2,100,1);
y = modelfun(b,x) + normrnd(0,0.1,100,1);

Set robust fitting options.

opts = statset('nlinfit');
opts.RobustWgtFun = 'bisquare';

Fit the nonlinear model using the robust fitting options.

beta0 = [2;2;2];
beta = nlinfit(x,y,modelfun,beta0,opts)
beta =

    1.0041
    3.0997
    2.1483

Nonlinear Regression Using Observation Weights

Load sample data.

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

Specify a vector of known observation weights.

W = [8 2 1 6 12 9 12 10 10 12 2 10 8]';

Fit the Hougen-Watson model to the rate data using the specified observation weights.

[beta,R,J,CovB] = nlinfit(X,y,@hougen,beta0,'Weights',W);
beta
beta =

    2.2068
    0.1077
    0.0766
    0.1818
    0.6516

Display the coefficient standard errors.

sqrt(diag(CovB))
ans =

    2.5721
    0.1251
    0.0950
    0.2043
    0.7735

Nonlinear Regression Using Weights Function Handle

Load sample data.

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

Specify a function handle for observation weights. The function accepts the model fitted values as input, and returns a vector of weights.

 a = 1; b = 1;
 weights = @(yhat) 1./((a + b*abs(yhat)).^2);

Fit the Hougen-Watson model to the rate data using the specified observation weights function.

[beta,R,J,CovB] = nlinfit(X,y,@hougen,beta0,'Weights',weights);
beta
beta =

    0.8308
    0.0409
    0.0251
    0.0801
    1.8261

Display the coefficient standard errors.

sqrt(diag(CovB))
ans =

    0.5822
    0.0297
    0.0197
    0.0578
    1.2810

Nonlinear Regression Using Nonconstant Error Model

Load sample data.

S = load('reaction');
X = S.reactants;
y = S.rate;
beta0 = S.beta;

Fit the Hougen-Watson model to the rate data using the combined error model.

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(X,y,@hougen,beta0,'ErrorModel','combined');
beta
beta =

    1.2526
    0.0628
    0.0400
    0.1124
    1.1914

Display the error model information.

ErrorModelInfo
ErrorModelInfo = 

              ErrorModel: 'combined'
         ErrorParameters: [0.1517 5.6783e-08]
           ErrorVariance: @(x)mse*(errorparam(1)+errorparam(2)*abs(model(beta,x))).^2
                     MSE: 1.6245
          ScheffeSimPred: 6
          WeightFunction: 0
            FixedWeights: 0
    RobustWeightFunction: 0

Input Arguments

expand all

X — Predictor variablesmatrix

Predictor variables for the nonlinear regression function, specified as a matrix. Typically, X is a design matrix of predictor (independent variable) values, with one row for each value in Y, and one column for each coefficient. However, X can be any array that modelfun can accept.

Data Types: single | double

Y — Response valuesvector

Response values (dependent variable) for fitting the nonlinear regression function, specified as a vector with the same number of rows as X.

Data Types: single | double

modelfun — Nonlinear regression model functionfunction handle

Nonlinear regression model function, specified as a function handle. modelfun must accept two input arguments, a coefficient vector and an array X—in that order—and return a vector of fitted response values.

For example, to specify the hougen nonlinear regression function, use the function handle @hougen.

Data Types: function_handle

beta0 — Initial coefficient valuesvector

Initial coefficient values for the least squares estimation algorithm, specified as a vector.

    Note:   Poor starting values can lead to a solution with large residual error.

Data Types: single | double

options — Estimation algorithm optionsstructure created using statset

Estimation algorithm options, specified as a structure you create using statset. The following statset parameters are applicable to nlinfit.

DerivStep — Relative difference for finite difference gradienteps^(1/3) (default) | positive scalar value | vector

Relative difference for the finite difference gradient calculation, specified as a positive scalar value, or a vector the same size as beta. Use a vector to specify a different relative difference for each coefficient.

Display — Level of output display'off' (default) | 'iter' | 'final'

Level of output display during estimation, specified as one of 'off', 'iter', or 'final'. If you specify 'iter', output is displayed at each iteration. If you specify 'final', output is displayed after the final iteration.

FunValCheck — Indicator for whether to check for invalid values'on' (default) | 'off'

Indicator for whether to check for invalid values such as NaN or Inf from the objective function, specified as 'on' or 'off'.

MaxIter — Maximum number of iterations100 (default) | positive integer

Maximum number of iterations for the estimation algorithm, specified as a positive integer. Iterations continue until estimates are within the convergence tolerance, or the maximum number of iterations specified by MaxIter is reached.

RobustWgtFun — Weight functionstring | function handle | []

Weight function for robust fitting, specified as a valid string or function handle.

    Note:   RobustWgtFun must have value [] when you use observation weights, W.

The following table describes the possible string values. Let r denote normalized residuals and w denote robust weights. The indicator function I[x] is equal to 1 if the expression x is true, and 0 otherwise.

Weight FunctionEquationDefault Tuning Constant
'' (default) No robust fitting
'andrews'

w=I[|r|<π]×sin(r)/r

1.339
'bisquare'

w=I[|r|<1]×(1r2)2

4.685
'cauchy'

w=1(1+r2)

2.385
'fair'

w=1(1+|r|)

1.400
'huber'

w=1max(1,|r|)

1.345
'logistic'

w=tanh(r)r

1.205
'talwar'

w=I[|r|<1]

2.795
'welsch'

w=exp{r2}

2.985

You can alternatively specify a function handle that accepts a vector of normalized residuals as input, and returns a vector of robust weights as output. If you use a function handle, you must provide a Tune constant.

Tune — Tuning constantpositive scalar value

Tuning constant for robust fitting, specified as a positive scalar value. The tuning constant is used to normalize residuals before applying a robust weight function. The default tuning constant depends on the function specified by RobustWgtFun.

If you use a function handle to specify RobustWgtFun, then you must specify a value for Tune.

TolFun — Termination tolerance on residual sum of squares1e-8 (default) | positive scalar value

Termination tolerance for the residual sum of squares, specified as a positive scalar value. Iterations continue until estimates are within the convergence tolerance, or the maximum number of iterations specified by MaxIter is reached.

TolX — Termination tolerance on estimated coefficients1e-8 (default) | positive scalar value

Termination tolerance on the estimated coefficients, beta, specified as a positive scalar value. Iterations continue until estimates are within the convergence tolerance, or the maximum number of iterations specified by MaxIter is reached.

Robust — Indicator for robust fitting'off' (default) | 'on'

Indicator for robust fitting, specified as 'off' or 'on'.

    Note:   Robust will be removed in a future software release. Use RobustWgtFun for robust fitting.

WgtFun — Weight function for robust fittingstring | function handle

Weight function for robust fitting, specified as a string indicating a weight function, or a function handle. WgtFun is valid only when Robust has value 'on'.

    Note:   WgtFun will be removed in a future software release. Use RobustWgtFun instead.

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example: 'ErrorModel','proportional','ErrorParameters',0.5 specifies a proportional error model, with initial value 0.5 for the error parameter estimation

'ErrorModel' — Form of error term'constant' (default) | 'proportional' | 'combined'

Form of the error term, specified as the comma-separated pair consisting of 'ErrorModel' and one of the following strings indicating the error model. Each model defines the error using a standard mean-zero and unit-variance variable e in combination with independent components: the function value f, and one or two parameters a and b.

'constant' (default)y=f+ae
'proportional'y=f+bfe
'combined'y=f+(a+b|f|)e

The only allowed error model when using Weights is 'constant'.

    Note:   options.RobustWgtFun must have value [] when using an error model other than 'constant'.

'ErrorParameters' — Initial estimates for error model parameters1 or [1,1] (default) | scalar value | two-element vector

Initial estimates for the error model parameters in the chosen ErrorModel, specified as the comma-separated pair consisting of 'ErrorParameters' and a scalar value or two-element vector.

Error ModelParametersDefault Values
'constant'a1
'proportional'b1
'combined'a, b[1,1]

For example, if 'ErrorModel' has the value 'combined', you can specify the starting value 1 for a and the starting value 2 for b as follows.

Example: 'ErrorParameters',[1,2]

You can only use the 'constant' error model when using Weights.

    Note:   options.RobustWgtFun must have value [] when using an error model other than 'constant'.

Data Types: double | single

'Weights' — Observation weightsvector | function handle

Observation weights, specified as the comma-separated pair consisting of 'Weights' and a vector of real positive weights or a function handle. You can use observation weights to down-weight the observations that you want to have less influence on the fitted model.

  • If W is a vector, then it must be the same size as Y.

  • If W is a function handle, then it must accept a vector of predicted response values as input, and return a vector of real positive weights as output.

    Note:   options.RobustWgtFun must have value [] when you use observation weights.

Data Types: double | single | function_handle

Output Arguments

expand all

beta — Estimated regression coefficientsvector

Estimated regression coefficients, returned as a vector. The number of elements in beta equals the number of elements in beta0.

Let f(Xi,b) denote the nonlinear function specified by modelfun, where xi are the predictors for observation i, i = 1,...,N, and b are the regression coefficients. The vector of coefficients returned in beta minimizes the weighted least squares equation,

i=1Nwi[yif(xi,b)]2.

For unweighted nonlinear regression, all of the weight terms are equal to 1.

R — Residualsvector

Residuals for the fitted model, returned as a vector.

  • If you specify observation weights using the name-value pair argument Weights, then R contains weighted residuals.

  • If you specify an error model other than 'constant' using the name-value pair argument ErrorModel, then you can no longer interpret R as model fit residuals.

J — Jacobianmatrix

Jacobian of the nonlinear regression model, modelfun, returned as an N-by-p matrix, where N is the number of observations and p is the number of estimated coefficients.

  • If you specify observation weights using the name-value pair argument Weights, then J is the weighted model function Jacobian.

  • If you specify an error model other than 'constant' using the name-value pair argument ErrorModel, then you can no longer interpret J as the model function Jacobian.

CovB — Estimated variance-covariance matrixmatrix

Estimated variance-covariance matrix for the fitted coefficients, beta, returned as a p-by-p matrix, where p is the number of estimated coefficients. If the model Jacobian, J, has full column rank, then CovB = inv(J'*J)*MSE, where MSE is the mean squared error.

MSE — Mean squared errorscalar value

Mean squared error (MSE) of the fitted model, returned as a scalar value. MSE is an estimate of the variance of the error term. If the model Jacobian, J, has full column rank, then MSE = (R'*R)/(N-p), where N is the number of observations, and p is the number of estimated coefficients.

ErrorModelInfo — Information about error model fitstructure

Information about the error model fit, returned as a structure with the following fields:

ErrorModelChosen error model
ErrorParametersEstimated error parameters
ErrorVarianceFunction handle that accepts an N-by-p matrix, X, and returns an N-by-1 vector of error variances using the estimated error model
MSEMean squared error
ScheffeSimPredScheffé parameter for simultaneous prediction intervals when using the estimated error model
WeightFunctionLogical with value true if you used a custom weight function previously in nlinfit
FixedWeightsLogical with value true if you used fixed weights previously in nlinfit
RobustWeightFunctionLogical with value true if you used robust fitting previously in

More About

expand all

Weighted Residuals

A weighted residual is a residual multiplied by the square root of the corresponding observation weight.

Given estimated regression coefficients, b, the residual for observation i is

ri=yif(xi,b),

where yi is the observed response and f(xi,b) is the fitted response at predictors xi.

When you fit a weighted nonlinear regression with weights wi, i = 1,...,N, nlinfit returns the weighted residuals,

ri*=wi(yif(xi,b)).

Weighted Model Function Jacobian

The weighted model function Jacobian is the nonlinear model Jacobian multiplied by the square root of the observation weight matrix.

Given estimated regression coefficients, b, the estimated model Jacobian, J,for the nonlinear function f(xi,b) has elements

Jij=f(xi,b)bj,

where bj is the jth element of b.

When you fit a weighted nonlinear regression with diagonal weights matrix W,nlinfit returns the weighted Jacobian matrix,

J*=W1/2J.

Tips

  • To produce error estimates on predictions, use the optional output arguments R, J, CovB, or MSE as inputs to nlpredci.

  • To produce error estimates on the estimated coefficients, beta, use the optional output arguments R, J, CovB, or MSE as inputs to nlparci.

  • If you use the robust fitting option, RobustWgtFun, you must use CovB—and might need MSE—as inputs to nlpredci or nlparci to ensure that the confidence intervals take the robust fit properly into account.

Algorithms

  • nlinfit treats NaN values in Y or modelfun(beta0,X) as missing data, and ignores the corresponding observations.

  • For nonrobust estimation, nlinfit uses the Levenberg-Marquardt nonlinear least squares algorithm [1].

  • For robust estimation, nlinfit uses an iterative reweighted least squares algorithm ([2], [3]). At each iteration, the robust weights are recalculated based on each observation's residual from the previous iteration. These weights downweight outliers, so that their influence on the fit is decreased. Iterations continue until the weights converge.

  • When you specify a function handle for observation weights, the weights depend on the fitted model. In this case, nlinfit uses an iterative generalized least squares algorithm to fit the nonlinear regression model.

References

[1] Seber, G. A. F., and C. J. Wild. Nonlinear Regression. Hoboken, NJ: Wiley-Interscience, 2003.

[2] DuMouchel, W. H., and F. L. O'Brien. "Integrating a Robust Option into a Multiple Regression Computing Environment." Computer Science and Statistics: Proceedings of the 21st Symposium on the Interface. Alexandria, VA: American Statistical Association, 1989.

[3] Holland, P. W., and R. E. Welsch. "Robust Regression Using Iteratively Reweighted Least-Squares." Communications in Statistics: Theory and Methods, A6, 1977, pp. 813–827.

Was this topic helpful?