Polynomial fitting

Least squares polynomial fits with one input (PolyFit) or several inputs and outputs (PolyFitMultivariate).

One input

PolyFit fits y = b0 + b1 x + ... + bd x^d by least squares, like numpy.polyfit. The columns of the Vandermonde matrix are scaled before the QR decomposition, so very large or very small x values keep their precision.

Fit a polynomialLeast squares fit of a degree 2 polynomial, like numpy.polyfit.

import 'package:scidart/numdart.dart';

void main() {
  var x = Array([0, 1, 2, 3, 4, 5]);
  var y = Array([1.1, 2.9, 9.2, 19.1, 32.8, 51.2]);

  var fit = PolyFit(x, y, 2);
  print(fit);
  print('coefficients (lowest power first): ${fit.coefficients().toList()}');
  print('R2 = ${fit.R2()}');
  print('prediction at x = 6: ${fit.predict(6)}');
}

Output

2.010714285714285 x^2 +  -0.050714285714280036 x +  1.0785714285714196   (R^2 = 0.9999313529244543)
coefficients (lowest power first): [2.010714285714285, -0.050714285714280036, 1.0785714285714196]
R2 = 0.9999313529244543
prediction at x = 6: 73.16

Note that coefficients() lists the lowest power first, the reverse of numpy.polyfit. When there are not enough distinct x values for the requested degree, or a power of x overflows, the degree is reduced: check polyDegree().

Try it: drag the degree and watch the fit and R² change.

Loading the interactive demo…

Several inputs and outputs

PolyFitMultivariate takes a matrix of samples (one row per sample, one column per feature) and a matrix of outputs. It builds every product of features up to the degree, in the same order as sklearn.preprocessing.PolynomialFeatures, and fits all outputs at once.

Polynomial regression with several featuresFit y from two inputs with interaction terms, like sklearn PolynomialFeatures plus LinearRegression.

import 'package:scidart/numdart.dart';

void main() {
  // y = 1 + 1.5 x0 + x1 + 0.5 x0^2 + 2 x0 x1
  var x = Array2d([
    Array([0, 0]),
    Array([1, 0]),
    Array([0, 1]),
    Array([1, 1]),
    Array([2, 1]),
    Array([1, 2]),
    Array([2, 2]),
  ]);
  var y = Array2d([
    for (var r in x)
      Array([1 + 1.5 * r[0] + r[1] + 0.5 * r[0] * r[0] + 2 * r[0] * r[1]])
  ]);

  var fit = PolyFitMultivariate(x, y, 2);
  for (var t = 0; t < fit.powers.length; t++) {
    var c = fit.coefficients[t][0];
    print('powers ${fit.powers[t]}: ${c.toStringAsFixed(6)}');
  }
  print('rank: ${fit.rank}, R2: ${fit.R2().toList()}');
  print('predict [3, 1]: ${fit.predict(Array([3, 1])).toList()}');
}

Output

powers [0, 0]: 1.000000
powers [1, 0]: 1.500000
powers [0, 1]: 1.000000
powers [2, 0]: 0.500000
powers [1, 1]: 2.000000
powers [0, 2]: -0.000000
rank: 6, R2: [1.0]
predict [3, 1]: [17.0]

powers[t][f] is the exponent of feature f in term t, and coefficients[t][o] is the coefficient of that term for output o. With more terms than samples it returns the minimum norm solution, like numpy.linalg.lstsq, and rank tells how many terms were really determined.