# 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

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

collapse all

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

Generate sample data from the nonlinear regression model $y={b}_{1}+{b}_{2}\cdot exp\left\{-{b}_{3}x\right\}+ϵ$, 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 ```

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

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

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

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\left[|r|<\pi \right]×\mathrm{sin}\left(r\right)/r$`

1.339
`'bisquare'`

`$w=I\left[|r|<1\right]×{\left(1-{r}^{2}\right)}^{2}$`

4.685
`'cauchy'`

`$w=1}{\left(1+{r}^{2}\right)}$`

2.385
`'fair'`

`$w=1}{\left(1+|r|\right)}$`

1.400
`'huber'`

`$w=1}{\mathrm{max}\left(1,|r|\right)}$`

1.345
`'logistic'`

`$w=\mathrm{tanh}\left(r\right)}{r}$`

1.205
`'talwar'`

`$w=I\left[|r|<1\right]$`

2.795
`'welsch'`

`$w=\mathrm{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.

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+\left(a+b|f|\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'`.

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'` 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`

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\left({X}_{i},b\right)$ 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,

`${\sum }_{i=1}^{N}{w}_{i}\left[{y}_{i}-f\left({x}_{i},b\right)\right]\text{\hspace{0.17em}}{\text{\hspace{0.17em}}}^{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.

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:

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

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

`${r}_{i}={y}_{i}-f\left({x}_{i},b\right),$`

where yi is the observed response and $f\left({x}_{i},b\right)$ is the fitted response at predictors ${x}_{i}.$

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

`${r}_{i}^{*}=\sqrt{{w}_{i}}\left({y}_{i}-f\left({x}_{i},b\right)\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\left({x}_{i},b\right)$ has elements

`${J}_{ij}=\frac{\partial f\left({x}_{i},b\right)}{\partial {b}_{j}},$`

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}^{*}={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.

## Version History

Introduced before R2006a