PolyFit constructor

PolyFit(Array x, Array y, int degree)

Performs a polynomial reggression on the data points (y[i], x[i]). x the values of the predictor variable y the corresponding values of the response variable degree the degree of the polynomial to fit

Implementation

PolyFit(Array x, Array y, this.degree) {
  if (degree < 0) {
    throw FormatException('expected deg >= 0');
  }
  if (x.isEmpty) {
    throw FormatException('expected non-empty vector for x');
  }
  if (x.length != y.length) {
    throw FormatException('expected x and y to have same length');
  }

  for (var i = 0; i < x.length; i++) {
    if (!x[i].isFinite || !y[i].isFinite) {
      throw FormatException('expected only finite values in x and y');
    }
  }

  var n = x.length;
  // n points define at most a polynomial of degree n - 1
  if (degree > n - 1) {
    degree = n - 1;
  }
  late QR qr;
  late ({Array maxAbs, Array rootSum}) scale;

  // in case Vandermonde matrix does not have full rank, reduce degree until it does
  while (true) {
    // build Vandermonde matrix with columns scaled to unit norm, as numpy
    // does, so the rank test does not depend on the magnitude of x
    var matrixX = matrixVander(x, N: degree + 1, increasing: true);

    // a power of x that overflows cannot be used: reduce the degree
    if (isFiniteMatrix(matrixX)) {
      scale = normalizeColumns(matrixX);

      // find least squares solution
      qr = QR(matrixX);
      var tolerance = math.max(n, degree + 1) * machineEpsilon;
      if (qr.isFullRank(tolerance: tolerance)) break;
    }

    // decrease degree and try again
    degree--;

    if (degree < 0) {
      throw FormatException('Matrix is rank deficient.');
    }
  }

  // create matrix from vector
  var matrixY = Array2d.fromVector(y, n);

  // linear regression coefficients, undoing the column scaling
  beta = qr.solve(matrixY);
  for (var j = 0; j < beta.row; j++) {
    beta[j][0] = beta[j][0] / scale.maxAbs[j] / scale.rootSum[j];
  }

  // mean of y[] values
  var meany = mean(y);

  // total variation to be accounted for
  for (var i = 0; i < n; i++) {
    var dev = y[i] - meany;
    sst += dev * dev;
  }

  // variation not accounted for
  for (var i = 0; i < n; i++) {
    var residual = predict(x[i]) - y[i];
    sse += residual * residual;
  }
}