Nov 21, 2013 at 9:52 AM
Edited Nov 21, 2013 at 9:56 AM

Using Fit.Polynomial() works quite well for nonweighted 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 6:20 PM
Edited Nov 27, 2013 at 6: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),
Matrix<double>.Build.DenseOfDiagonalArray(w)).ToArray();
}
Example:
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 LevenbergMarquardt internally which we do not provide yet either, although we're currently working on integration nonlinear fitting and minimization routines.
Thanks,
Christoph



Update: I've actually pushed PolynomialWeighted to master, so it will be available in the next alpha build.



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 9:28 PM
Edited Nov 27, 2013 at 10:29 PM

Good catch, yes, we should rewrite it that way!
Thanks,
Christoph



Optimizations now applied in master.
Thanks,
Christoph

