PolyFit

This package extends Polynomials.jl and provides methods to fit a polynomial to data with error bars. The solution is obtained through an L2 minimization, where the covariance matrix is used to define the misfit that is minimized.

We solve the system

We can express this as

Solving the equation system is equivalent to minimizing the misfit

where C is the covariance matrix. The diagonal elements of the covariance matrix are the variances in y. Fitting data without error bars is equivalent to setting C = σ^2 I for a constant σ.

Installation

Install the package using

pkg> add https://github.com/jishnub/PolyFit.jl

julia> using PolyFit

Usage

This package adds methods to polyfit, and the usage is similar to Polynomials.jl

julia> x=1:50;y=x.^2;

julia> polyfit(x,y,2)
Poly(3.1065291710470317e-14 - 2.6679598459471846e-15*x + 1.0000000000000002*x^2)

The extra functionality this adds is the ability to pass error bars on y to the fit. This can be done either by passing error bars σ(y) or a covariance matrix C(yi,yj).

Error bars

julia> σ=10rand(length(y));

julia> polyfit(x,y,2,σ)
Poly(2.643961948331425e-14 - 2.7668654490519844e-15*x + 1.0000000000000002*x^2)

Covariance Matrix

A diagonal covariance matrix is equiavlent to passing error bars. The diagonal elements are the variances.

julia> polyfit(x,y,2,Diagonal(σ.^2))
Poly(2.6426223525629683e-14 - 2.7636271626091054e-15*x + 1.0000000000000002*x^2)

A non-diagonal covariance matrix has to be positive-definite

julia> C = Symmetric(Diagonal(σ.^2) + 0.1*rand(length(y),length(y)));

julia> isposdef(C)
true

julia> polyfit(x,y,2,C)
Poly(5.517934888038756e-14 + 5.126632798188113e-15*x + 0.9999999999999998*x^2)

SpecialMatrices

This package uses SpecialMatrices.jl to solve Van der Monde matrices, so the special case of fitting a degree-n polynomial to n+1 points without error bars is faster than Polynomials.jl. This method can be called using the special tag SquareSystem().

julia> @btime polyfit($x,$y,length($x)-1);
  168.393 μs (59 allocations: 65.02 KiB)

julia> @btime polyfit(SquareSystem(),$x,$y,length($x)-1);
  2.071 μs (4 allocations: 1.11 KiB)

julia> polyfit(x,y,length(x)-1)≈polyfit(SquareSystem(),x,y,length(x)-1)
true