nlinfit

Nonlinear regression

Syntax

beta = nlinfit(X,Y,modelfun,beta0)

beta = nlinfit(X,Y,modelfun,beta0,options)

beta = nlinfit(___,Name,Value)

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(___)

Description

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.

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.

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.

[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

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 = 5×1

    1.2526
    0.0628
    0.0400
    0.1124
    1.1914

Nonlinear Regression Using Robust Options

Generate sample data from the nonlinear regression model $$y=b_1+b_2\cdot exp\{-b_{3}x\}+ϵ$$, where $$b_1$$, $$b_2$$, and $$b_3$$ 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 = 3×1

    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 = 5×1

    2.2068
    0.1077
    0.0766
    0.1818
    0.6516

Display the coefficient standard errors.

sqrt(diag(CovB))

ans = 5×1

    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 = 5×1

    0.8308
    0.0409
    0.0251
    0.0801
    1.8261

Display the coefficient standard errors.

sqrt(diag(CovB))

ans = 5×1

    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 = 5×1

    1.2526
    0.0628
    0.0400
    0.1124
    1.1914

Display the error model information.

ErrorModelInfo

ErrorModelInfo = struct with fields:
              ErrorModel: 'combined'
         ErrorParameters: [0.1517 5.6783e-08]
           ErrorVariance: [function_handle]
                     MSE: 1.6245
          ScheffeSimPred: 6
          WeightFunction: 0
            FixedWeights: 0
    RobustWeightFunction: 0

Input Arguments

X — Predictor variables

matrix

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 predictor. However, X can be any array that modelfun can accept.

Data Types: single | double

Y — Response values

vector

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 function

function 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 values

vector

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 options

structure 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 gradient

eps^(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 iterations

100 (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 function

character vector | string scalar | function handle | []

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

Note

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

The following table describes the possible character vectors and string scalars. 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 Function	Equation	Default Tuning Constant
'' (default)	No robust fitting	—
'andrews'	$$w=I\left[\left|r\right|<\pi \right]\times \sin (r)/r$$	1.339
'bisquare'	$$w=I\left[\left|r\right|<1\right]\times {\left(1-r^2\right)}^2$$	4.685
'cauchy'	$$w=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left(1+r^2\right)$}\right.$$	2.385
'fair'	$$w=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left(1+\left|r\right|\right)$}\right.$$	1.400
'huber'	$$w=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\max \left(1,\left|r\right|\right)$}\right.$$	1.345
'logistic'	$$w=\raisebox{1ex}{$\tanh (r)$}\!\left/ \!\raisebox{-1ex}{$r$}\right.$$	1.205
'talwar'	$$w=I\left[\left|r\right|<1\right]$$	2.795
'welsch'	$$w=\exp \left\{-r^2\right\}$$	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 constant

positive 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 squares

1e-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 coefficients

1e-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.

Name-Value 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 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 'constant', 'proportional', or 'combined' 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+\left(a+b\left|f\right|\right)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 parameters

1 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 Model	Parameters	Default Values
'constant'	a	1
'proportional'	b	1
'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 weights

vector | 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

beta — Estimated regression coefficients

vector

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

Let $$f(X_i,b)$$ denote the nonlinear function specified by modelfun, where $$x_i$$ 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,

$${\displaystyle {\sum }_{i=1}^N{w}_i\left[y_i-f(x_i,b)\right]}{}^2.$$

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

R — Residuals

vector

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 — Jacobian

matrix

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 matrix

matrix

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 error

scalar 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 fit

structure

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

ErrorModel	Chosen error model
ErrorParameters	Estimated error parameters
ErrorVariance	Function handle that accepts an N-by-p matrix, X, and returns an N-by-1 vector of error variances using the estimated error model
MSE	Mean squared error
ScheffeSimPred	Scheffé parameter for simultaneous prediction intervals when using the estimated error model
WeightFunction	Logical with value true if you used a custom weight function previously in nlinfit
FixedWeights	Logical with value true if you used fixed weights previously in nlinfit
RobustWeightFunction	Logical with value true if you used robust fitting previously in nlinfit

More About

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

$$r_i=y_i-f(x_i,b),$$

where y_i is the observed response and $$f(x_i,b)$$ is the fitted response at predictors $$x_i.$$

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

$$r_i^{*}=\sqrt{w_i}\left(y_i-f(x_i,b)\right).$$

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(x_i,b)$$ has elements

$$J_{ij}=\frac{\partial f(x_i,b)}{\partial b_j},$$

where b_j 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^{*}=W^{1/2}J.$$

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 the algorithm of Iteratively Reweighted Least Squares ([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.