Main Content

nlinfit

Nonlinear regression

collapse all in page

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.

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.

example

Examples

collapse all

Open Live Script

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
Open Live Script

Generate sample data from the nonlinear regression model y=b1+b2⋅exp{-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 = 3×1

    1.0041
    3.0997
    2.1483
Open Live Script

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
Open Live Script

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
Open Live Script

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: @(x)mse*(errorparam(1)+errorparam(2)*abs(model(beta,x))).^2
                     MSE: 1.6245
          ScheffeSimPred: 6
          WeightFunction: 0
            FixedWeights: 0
    RobustWeightFunction: 0

Input Arguments

collapse all

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

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

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

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

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

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.

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.

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

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.

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 FunctionEquationDefault Tuning Constant
'' (default) No robust fitting—
'andrews'

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

1.339
'bisquare'

w=I[|r|<1]×(1−r2)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.

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.

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.

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 pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

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

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+(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'.

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

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

collapse all

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[yi−f(xi,b)]  2.

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

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.

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.

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.

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.

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 nlinfit

More About

collapse 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=yi−f(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(yi−f(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 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.

Version History

Introduced before R2006a

See Also

| | |

Topics