Simple Least Squares fitting

The CurveFit module provides functions that implement a few least squares approximations.

It is simple in an engineering sense. It simply returns the coefficients and does not do any error analysis.

It does, however, provide a simple and common interface to the routines.

The package also includes nonlinear least squares fitting using a Newton type algorithm

The fitting algorithms include

• Straight lines
• Polynomial fitting
• Power laws
• Log laws
• Exp laws
• Rational polynomial fitting
• A generic non-linear fitting algorithm
• King's law (used in hotwire anemometry)

High Level interface for fitting exp laws

(sol::ExpSumFit)(x)

Calculate the sum of exponentials using solution sol at points x.

Type that represents the modified King's law

High Level interface for fitting straight lines

Type that represents a Linear (original) King's law

High Level interface for fitting log laws

High Level interface for fitting power laws

Type defining a rational polynomial

A rational polynomial is the ratio of two polynomials and it is often useful in approximating functions.

call overload for calling directly ratval

#Uses the object created by curve_fit to estimate values

The call method is overloaded so that the fit object can be used as a function:

Example:

x = 1.0:10.0 @. y = 2*x + 1 + randn()

fit = curve_fit(LinearFit, x, y)

y1 = fit(5.1) y2 = apply_fit(fit, 5.1)

calc_integral_rules(n; m = 2)

Determine coefficients of the rules cumulative integrals of order n using method of undetermined coefficients. Interpolation order is m.

• n=1, m=1 Trapezoidal rule
• n=1, m=2 Simpson's 1/3 rule
• n=1, m=3 Simpson's second (3/8) rule

cumints!(init::ExpSumFitInit, x::AbstractVector{T}, y::AbstractVector{T}) where {T <: Real}

Calculates in-place init.n cumulative integrals using coeffients init.coeff.

Generic interface for curve fitting.

The same function curve_fit can be used to fit the data depending on fit type, shich is specified in the first parameter. This function returns an object that can be used to estimate the value of the fitting model using function apply_fit.

A few examples:

• f = curve_fit(LinearFit, x, y)
• f = curve_fit(Polynomial, x, y, n)

Fits an exp through a set of points: y = a₁*exp(a₂*x)

expsum_fill_X!(init::ExpSumFitInit, x, λ, X)

Fills matrix init.X

\[ [e^{x\lambda_1} e^{x\lambda_2} \ldots e^{x\lambda_n}]\]

or if fitting is done with constant

\[ [e^{x\lambda_1} e^{x\lambda_2} \ldots e^{x\lambda_n} 1]\]

expsum_fill_Y!(init::ExpSumFitInit, x, y)

Fills matrix Y with

\[ [\int^1 \int^2 \ldots \int^n x^0 x^1 \ldots x^{m-n-1}]\]

where ∫ⁱ means the ith cumulative integral of y vs x, and m is the number of columns in matrix Y.

    expsum_fit(x::Union{T,AbstractVector{T}}, y::AbstractVector{T}, n::Int; m::Int = 1, withconst::Bool = true, init::Union{Nothing,ExpSumFitInit} = nothing) where {T <: Real}

Fits the sum of n exponentials and a constant.

\[ y = k + p_1 e^{\lambda_1 t} + p_2 e^{\lambda_2 t} + \ldots + p_n e^{\lambda_n t}\]

If the keyword withconst is set to false, the constant is not fitted but set k=0.

Uses numerical integration with m strips, where the default m=1 uses linear interpolation. m=2 and higher require uniform interval and usually lead to better accuracy.

Passing init preallocates most needed memory, and can be initialized with expsum_init.

Returns a ExpSumFit struct containing a constant k and vectors p, λ, τ, where τ = -1/λ.

The algorithm is from Matlab code of Juan Gonzales Burgos.

expsum_init(x::AbstractVector{T}, n::Int; m::Int = 1, withconst::Bool = true) where {T <: Real}

Initialize most of the memory needed for expsum_fit.

Calculates the coefficients of a linear model using Least Squares

Given a Vandermonde matrix, this function uses the QR factorization to compute the Least Squares fit of a linear model. The code probably could be optimized.

a = gaussnewtonfit(x, y, fun, ∇fun!, a0[[, eps,] maxiter])

Gauss-Newton nonlinear least squares. Given vectors x and y, the tries to fit parameters a to a function f using least squares approximation:

$y = f(x, a₁, ..., aₙ)$

For more general approximations, see gauss_newton_fit.

Arguments:

• x Vector with x values
• y Vector with y values
• fun a function that is called as fun(x, a) where a is a vector of parameters.
• ∇fun! A function that calculares the derivatives with respect to parameters a
• a0 Vector with the initial guesses of the parameters
• eps Maximum residual for convergence
• maxiter Maximum number of iterations for convergence

Return value

A vector with the convrged array. If no convergence is achieved, the function throws an error.

Specification of the fitting function

The function that should be fitted shoud be specified by Julia funcion with the following signature:

fun(x::T, a::AbstractVector{T}) where {T<:Number}

The derivatives with respect to each fitting parameter a[i] should have the following signature:

∇fun!(x::T, a::AbstractVector{T}, df::AbstractVector{T}) where {T<:Number}

No return value is expected and the derivatives are returned in argument df.

Initial approximation (guess)

If the initial approximation is not good enough, divergence is possible.

Careful with parameters close to 0. The initial guess should never be 0.0 because the initial value of the parameter is used as reference value for computing resiuduals.

Convergence criteria

The argumento maxiter specifies the maximum number of iterations that should be carried out. At each iteration,

$aₖⁿ⁺¹ = aₖⁿ + δₖ$

Convergence is achieved when

$|δᵏ / aₖ⁰| < ε$

Example

x = 1.0:10.0
a = [3.0, 2.0, 1.0]
y = a[1] + a[2]*x + a[3]*x^2
fun(x, a) = a[1] + a[2]*x + a[3]*x^2

function ∇fun!(x, a, df) 
    df[1] = 1.0
    df[2] = x
    df[3] = x^2
end

a = gauss_newton_fit(x, y, fun, ∇fun!, [0.5, 0.5, 0.5], 1e-8, 30)

a = gaussnewtongeneric_fit(x, y, fun, ∇fun!, a0[[, eps,] maxiter])

Gauss-Newton nonlinear least squares. Given matrix x the function tries to fit parameters a to an implicit function f using least squares approximation:

$f(x₁,..., xₘ, a₁, ..., aₙ) = 0$

This function is very similar to gauss_newton_fit but more generic. It doesn't limit the number of independent variables and doesn't require a y LHS.

Arguments:

• x Matrix where each column is a variable
• `
• fun a function that is called as fun(x, a) where a is a vector of parameters.
• ∇fun! A function that calculares the derivatives with respect to parameters a
• a0 Vector with the initial guesses of the parameters
• eps Maximum residual for convergence
• maxiter Maximum number of iterations for convergence

Return value

A vector with the convrged array. If no convergence is achieved, the function throws an error.

Specification of the fitting function

The function that should be fitted shoud be specified by Julia funcion with the following signature:

fun(x::T, a::AbstractVector{T}) where {T<:Number}

The derivatives with respect to each fitting parameter a[i] should have the following signature:

∇fun!(x::T, a::AbstractVector{T}, df::AbstractVector{T}) where {T<:Number}

No return value is expected and the derivatives are returned in argument df.

Initial approximation (guess)

If the initial approximation is not good enough, divergence is possible.

Careful with parameters close to 0. The initial guess should never be 0.0 because the initial value of the parameter is used as reference value for computing resiuduals.

Convergence criteria

The argumento maxiter specifies the maximum number of iterations that should be carried out. At each iteration,

$aₖⁿ⁺¹ = aₖⁿ + δₖ$

Convergence is achieved when

$|δᵏ / aₖ⁰| < ε$

Example

x = 1.0:10.0
a = [3.0, 2.0, 1.0]
y = a[1] + a[2]*x + a[3]*x^2
fun(x, a) = a[1] + a[2]*x + a[3]*x^2

function ∇fun!(x, a, df) 
    df[1] = 1.0
    df[2] = x
    df[3] = x^2
end

a = gauss_newton_fit(x, y, fun, ∇fun!, [0.5, 0.5, 0.5], 1e-8, 30)

Uses nonlinear least squares to fit the modified King's law:

E^2 = A + B*U^n

The Original (linear) King's law is used to estimate A and B when n=1/2. This initial value is used as an initial guess for fitting the nonlinear modified King's law using the function nonlinear_fit.

Equation that computes the error of the modified King's law

Fits a straight line through a set of points, y = a₁ + a₂ * x

Original Kings law (1910) represents the relation between voltage and velocity in a hotwire anemometer. The law is given by:

E^2 = A + B*U^0.5

This function estimates A and B.

Linear Rational LeastSquares

The following curvit is done:

y = p(x) / q(x)

where p(x) and q(x) are polynomials.

The linear case is solved by doing a least square fit on

y*q(x) = p(x)

where the zero order term o q(x) is assumed to be 1.

Fits a log function through a set of points: y = a₁+ a₂*log(x)

Auxiliary function used in nonlinear least squares

Nonlinear multi-variable least squares

Uses a Newton like algorithm to compute the least squares fit of a model

Parameters:

• x an array where each colum is a variable
• fun The function that should be fitted (more later on)
• a0 An initial guess for each fitting parameter
• eps Convergence criterion
• maxiter Maximum number of iterations

Specifying the model

The model is specified in argument fun that should be callable according to the following signature:

r = fun(x, a)

where both x and a are vectors

The way the algorithm is implemented allows for the function to be implicit. This means that we are trying to fit the relationship fun(x, a) = 0.

Initial guess

The initial guess is important and should be as close as possible to the expeted values. It should, at least, give an order of magnitude of each parameter. Zero values are not recommended since in this case no order of magnitude is available and it is assumed. The initial step is assumed to be a0/10.

Convergence criterion

The initial guess a0 will probably not fit the data very well. The largest value of maxerr = maxabs(fun(x, a0)) (maximum error) is used as a reference and every time a new approximation is computed, if the change in paramaters is compared to this reference error. If it is small, maxabs(da) < eps*maxerr, the algorithm converged.

Return values

A tuple containing:

• The estimated coefficients
• A Bool stating whether the algorithm converged
• Number of iterations it took to converge.

Fits a polynomial of degree n through a set of points.

Simple algorithm, doesn't use orthogonal polynomials or any such thing and therefore unconditioned matrices are possible. Use it only for low degree polynomials.

This function returns a the coefficients of the polynomial.

Fits a power law through a set of points: y = a₁*x^a₂

Carry out a nonlinear least squares of rational polynomials

Find the polynomial coefficients that best approximate the points given by x and y.

Evaluate a rational polynomial

Create Vandermonde matrix for simple polynomial fit