Tutorial

Assume that, for the th observation, the relationship between independent variable and dependent variable follows:

where is a non-linear model function depends on the independent variable and the parameter vector . In order to find the parameter that "best" fit our data, we choose the parameter which minimizes the sum of squared residuals from our data, i.e. solves the problem:

Given that the function is non-linear, there’s no analytical solution for the best . We have to use computational tools, which is LsqFit.jl in this tutorial, to find the least squares solution.

One example of non-linear model is the exponential model, which takes a one-element predictor variable . The model function is:

and the model becomes:

To fit data using LsqFit.jl, pass the defined model function (m), data (tdata and ydata) and the initial parameter value (p0) to curve_fit(). For now, LsqFit.jl only supports the Levenberg Marquardt algorithm.

It will return a composite type LsqFitResult, with some interesting values:

The Jacobian of a vector function is defined as the matrix with elements:

The matrix is therefore:

The Jacobian of the exponential model function with respect to is:

By default, the finite differences is used (see NLSolversBase.jl for more information), is used to approximate the Jacobian for the data fitting algorithm and covariance computation. Alternatively, a function which calculates the Jacobian can be supplied to curve_fit() for faster and/or more accurate results.

The non-linear function can be approximated as a linear function by Talor expansion:

where is a fixed vector, is a very small-valued vector and is the gradient at .

Consider the residual vector functon

with entries:

Each entry’s linear approximation can hence be written as:

Since the th row of equals the transpose of the gradient of , the vector function can be approximated as:

which is a linear function on since is a fixed vector.

The linear approximation of the non-linear least squares problem leads to the approximation of the covariance matrix of each parameter, from which we can perform regression analysis.

Consider a least squares solution , which is a local minimizer of the non-linear problem:

Set as the fixed point in linear approximation, and . A parameter vector near can be expressed as . The local approximation for the least squares problem is:

which is essentially the linear least squares problem:

where , and . Solve the equation where the partial derivatives equal to , the analytical solution is:

The covariance matrix for the analytical solution is:

Note that is the residual vector at the best fit point , with entries . is very close to and therefore can be replaced by .

Assume the errors in each sample are independent, normal distributed with zero mean and same variance, i.e. , the covariance matrix from the linear approximation is therefore:

where could be estimated as residual sum of squares devided by degrees of freedom:

In LsqFit.jl, the covariance matrix calculation uses QR decomposition to be more computationally stable, which has the form:

vcov() computes the covariance matrix of fit:

The standard error is then the square root of each diagonal elements of the covariance matrix. stderror() returns the standard error of each parameter:

margin_error() computes the product of standard error and the critical value of each parameter at a certain significance level (default is 5%) from t-distribution. The margin of error at 10% significance level can be computed by:

confint() returns the confidence interval of each parameter at certain significance level, which is essentially the estimate value ± margin of error. To get the confidence interval at 10% significance level, run:

curve_fit() also accepts weight parameter (wt) to perform Weighted Least Squares and General Least Squares, where the parameter minimizes the weighted residual sum of squares.

Weight parameter (wt) is an array or a matrix of weights for each sample. To perform Weighted Least Squares, pass the weight array [w_1, w_2, ..., w_n] or the weight matrix W:

The weighted least squares problem becomes:

in matrix form:

where

is a residual vector function with entries:

The algorithm in LsqFit.jl will then provide a least squares solution .

In LsqFit.jl, the residual function passed to levenberg_marquardt() is in different format, if the weight is a vector:

Cholesky decomposition, which is effectively a sqrt of a matrix, will be performed if the weight is a matrix:

The solution will be the same as the least squares problem mentioned in the tutorial.

Set and , the linear approximation of the weighted least squares problem is then:

The analytical solution to the linear approximation is:

Assume the errors in each sample are independent, normal distributed with zero mean and different variances (heteroskedastic error), i.e. , where:

We know the error variance and we set the weight as the inverse of the variance (the optimal weight), i.e. :

The covariance matrix is now:

If we only know the relative ratio of different variances, i.e. , the covariance matrix will be:

where is estimated. In this case, if we set , the result will be the same as the unweighted version. However, curve_fit() currently does not support this implementation. curve_fit() assumes the weight as the inverse of the error covariance matrix rather than the ratio of error covariance matrix, i.e. the covariance of the estimated parameter is calculated as covar = inv(J'*fit.wt*J).

Pass the vector of 1 ./ var(ε) or the matrix inv(covar(ε)) as the weight parameter (wt) to the function curve_fit():

Assume the errors in each sample are correlated, normal distributed with zero mean and different variances (heteroskedastic and autocorrelated error), i.e. .

Set the weight matrix as the inverse of the error covariance matrix (the optimal weight), i.e. , we will get the parameter covariance matrix:

Pass the matrix inv(covar(ε)) as the weight parameter (wt) to the function curve_fit():

In most cases, the variances of errors are unknown. To perform Weighted Least Square, we need estimate the variances of errors first, which is the squared residual of th sample:

Unweighted fitting (OLS) will return the residuals we need, since the estimator of OLS is unbiased. Then pass the reciprocal of the residuals as the estimated optimal weight to perform Weighted Least Squares:

Hansen, P. C., Pereyra, V. and Scherer, G. (2013) Least squares data fitting with applications. Baltimore, Md: Johns Hopkins University Press, p. 147-155.

Kutner, M. H. et al. (2005) Applied Linear statistical models.

Weisberg, S. (2014) Applied linear regression. Fourth edition. Hoboken, NJ: Wiley (Wiley series in probability and statistics).