loess function

LoessResult loess(List<double> x, List<double> y, { double span = 0.75, int degree = 1, })

Computes LOESS smoothing for the given data.

span is the fraction of data to use in each local regression (0 < span <= 1). degree is the degree of the local polynomial (1 or 2).

Implementation

LoessResult loess(
  List<double> x,
  List<double> y, {
  double span = 0.75,
  int degree = 1,
}) {
  if (x.length != y.length) {
    throw ArgumentError('x and y must have the same length');
  }
  if (span <= 0 || span > 1) {
    throw ArgumentError('Span must be in (0, 1]');
  }
  if (degree != 1 && degree != 2) {
    throw ArgumentError('Degree must be 1 or 2');
  }

  final n = x.length;
  if (n == 0) return LoessResult(x: [], y: []);

  // Sort by x
  final indices = List.generate(n, (i) => i);
  indices.sort((a, b) => x[a].compareTo(x[b]));

  final sortedX = indices.map((i) => x[i]).toList();
  final sortedY = indices.map((i) => y[i]).toList();

  final k = (span * n).ceil().clamp(degree + 1, n);
  final smoothedY = List<double>.filled(n, 0);

  for (int i = 0; i < n; i++) {
    final xi = sortedX[i];

    // Find k nearest neighbors
    final distances = List.generate(n, (j) => (sortedX[j] - xi).abs());
    final neighborIndices = List.generate(n, (j) => j);
    neighborIndices.sort((a, b) => distances[a].compareTo(distances[b]));
    final neighbors = neighborIndices.take(k).toList();

    // Calculate weights using tricube function
    final maxDist = distances[neighbors.last];
    final weights = <double>[];
    for (final j in neighbors) {
      final u = maxDist > 0 ? distances[j] / maxDist : 0.0;
      weights.add(_tricube(u));
    }

    // Perform weighted polynomial regression
    if (degree == 1) {
      // Weighted linear regression
      double sumW = 0, sumWx = 0, sumWy = 0, sumWxx = 0, sumWxy = 0;
      for (int j = 0; j < k; j++) {
        final idx = neighbors[j];
        final w = weights[j];
        final xj = sortedX[idx];
        final yj = sortedY[idx];
        sumW += w;
        sumWx += w * xj;
        sumWy += w * yj;
        sumWxx += w * xj * xj;
        sumWxy += w * xj * yj;
      }

      final denom = sumW * sumWxx - sumWx * sumWx;
      if (denom.abs() < 1e-10) {
        smoothedY[i] = sumWy / sumW;
      } else {
        final slope = (sumW * sumWxy - sumWx * sumWy) / denom;
        final intercept = (sumWy - slope * sumWx) / sumW;
        smoothedY[i] = slope * xi + intercept;
      }
    } else {
      // Weighted quadratic regression
      final localX = neighbors.map((j) => sortedX[j]).toList();
      final localY = neighbors.map((j) => sortedY[j]).toList();
      final result = _weightedPolynomialRegression(localX, localY, weights, 2);
      smoothedY[i] = result[0] + result[1] * xi + result[2] * xi * xi;
    }
  }

  return LoessResult(x: sortedX, y: smoothedY);
}