This project has moved. For the latest updates, please go here.

How to do a weighted polynomial fit?

Nov 21, 2013 at 8:52 AM
Edited Nov 21, 2013 at 8:56 AM
Using Fit.Polynomial() works quite well for non-weighted fitting.

I'm translating some python/scipy code that uses optimize.curve_fit() with relative weights (sigma).
I saw that there is WeightedRegression class, but I don't know how to use it.
Nov 27, 2013 at 5:20 PM
Edited Nov 27, 2013 at 5:26 PM
Hi Thomas,

A weighted version of Fit.Polynomial using WeightedRegression (see also) could look like this:
double[] PolynomialWeighted(double[] x, double[] y, double[] w, int order)
    var design = Matrix<double>.Build.Dense(
        x.Length, order + 1, (i, j) => Math.Pow(x[i], j));

    return WeightedRegression.Weighted(
        design, Vector<double>.Build.Dense(y),
var x = Enumerable.Range(1, 6).Select(Convert.ToDouble).ToArray();
var y = new[] { 4.986, 2.347, 2.061, -2.995, -2.352, -5.782 };
var w = new[] { 0.5, 1.0, 1.0, 1.0, 1.0, 0.5 };
var order = 2;

var resp = Fit.Polynomial(x, y, order).Dump();
var resp2 = Fit.PolynomialWeighted(x, y, w, order).Dump();
However, note that WeightedRegression actually uses normal equations (with Cholesky decomposition) internally, so it is not very robust and not a good choice for polynomials of high order. Fit.Polynomial uses a QR decomposition which is a bit better, but still not very good for polynomials. At some point we'll want to replace the polynomial case with a specialized algorithm.

To my understanding, scypi optimize.curve_fit uses Levenberg-Marquardt internally which we do not provide yet either, although we're currently working on integration non-linear fitting and minimization routines.

Marked as answer by tibel on 11/27/2013 at 12:16 PM
Nov 27, 2013 at 6:32 PM
Update: I've actually pushed PolynomialWeighted to master, so it will be available in the next alpha build.
Nov 27, 2013 at 7:06 PM
Thanks Christoph,
that is exactly what I needed ;-)

Optimization hint (not sure):
In WeightedRegression the lines like
return x.TransposeThisAndMultiply(w*x).Cholesky().Solve(x.Transpose()*(w*y));
can be rewritten to
return x.TransposeThisAndMultiply(w*x).Cholesky().Solve(x.TransposeThisAndMultiply(w*y));
Also DiagonalMatrixStorage shouldn't be used at the moment, as they are not optimized and therefor slow.
Nov 27, 2013 at 8:28 PM
Edited Nov 27, 2013 at 9:29 PM
Good catch, yes, we should rewrite it that way!

Nov 28, 2013 at 7:36 AM
Optimizations now applied in master.