public final double[] smooth(final double[] xval,final double[] yval,final double[] weights) throws MathIllegalArgumentException {
  if (xval.length != yval.length) {
    throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,xval.length,yval.length);
  }
  final int n=xval.length;
  if (n == 0) {
    throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
  }
  checkAllFiniteReal(xval);
  checkAllFiniteReal(yval);
  checkAllFiniteReal(weights);
  MathArrays.checkOrder(xval);
  if (n == 1) {
    return new double[]{yval[0]};
  }
  if (n == 2) {
    return new double[]{yval[0],yval[1]};
  }
  int bandwidthInPoints=(int)(bandwidth * n);
  if (bandwidthInPoints < 2) {
    throw new MathIllegalArgumentException(LocalizedCoreFormats.BANDWIDTH,bandwidthInPoints,2,true);
  }
  final double[] res=new double[n];
  final double[] residuals=new double[n];
  final double[] sortedResiduals=new double[n];
  final double[] robustnessWeights=new double[n];
  Arrays.fill(robustnessWeights,1);
  for (int iter=0; iter <= robustnessIters; ++iter) {
    final int[] bandwidthInterval={0,bandwidthInPoints - 1};
    for (int i=0; i < n; ++i) {
      final double x=xval[i];
      if (i > 0) {
        updateBandwidthInterval(xval,weights,i,bandwidthInterval);
      }
      final int ileft=bandwidthInterval[0];
      final int iright=bandwidthInterval[1];
      final int edge;
      if (xval[i] - xval[ileft] > xval[iright] - xval[i]) {
        edge=ileft;
      }
 else {
        edge=iright;
      }
      double sumWeights=0;
      double sumX=0;
      double sumXSquared=0;
      double sumY=0;
      double sumXY=0;
      double denom=FastMath.abs(1.0 / (xval[edge] - x));
      for (int k=ileft; k <= iright; ++k) {
        final double xk=xval[k];
        final double yk=yval[k];
        final double dist=(k < i) ? x - xk : xk - x;
        final double w=tricube(dist * denom) * robustnessWeights[k] * weights[k];
        final double xkw=xk * w;
        sumWeights+=w;
        sumX+=xkw;
        sumXSquared+=xk * xkw;
        sumY+=yk * w;
        sumXY+=yk * xkw;
      }
      final double meanX=sumX / sumWeights;
      final double meanY=sumY / sumWeights;
      final double meanXY=sumXY / sumWeights;
      final double meanXSquared=sumXSquared / sumWeights;
      final double beta;
      if (FastMath.sqrt(FastMath.abs(meanXSquared - meanX * meanX)) < accuracy) {
        beta=0;
      }
 else {
        beta=(meanXY - meanX * meanY) / (meanXSquared - meanX * meanX);
      }
      final double alpha=meanY - beta * meanX;
      res[i]=beta * x + alpha;
      residuals[i]=FastMath.abs(yval[i] - res[i]);
    }
    if (iter == robustnessIters) {
      break;
    }
    System.arraycopy(residuals,0,sortedResiduals,0,n);
    Arrays.sort(sortedResiduals);
    final double medianResidual=sortedResiduals[n / 2];
    if (FastMath.abs(medianResidual) < accuracy) {
      break;
    }
    for (int i=0; i < n; ++i) {
      final double arg=residuals[i] / (6 * medianResidual);
      if (arg >= 1) {
        robustnessWeights[i]=0;
      }
 else {
        final double w=1 - arg * arg;
        robustnessWeights[i]=w * w;
      }
    }
  }
  return res;
}

public final double [ ] smooth ( final double [ ] xval , final double [ ] yval , final double [ ] weights ) throws MathIllegalArgumentException { if ( xval . length != yval . length ) { throw new MathIllegalArgumentException ( LocalizedCoreFormats . DIMENSIONS_MISMATCH , xval . length , yval . length ) ; } final int n = xval . length ; if ( n == 0 ) { throw new MathIllegalArgumentException ( LocalizedCoreFormats . NO_DATA ) ; } checkAllFiniteReal ( xval ) ; checkAllFiniteReal ( yval ) ; checkAllFiniteReal ( weights ) ; MathArrays . checkOrder ( xval ) ; if ( n == 1 ) { return new double [ ] { yval [ 0 ] } ; } if ( n == 2 ) { return new double [ ] { yval [ 0 ] , yval [ 1 ] } ; } int bandwidthInPoints = ( int ) ( bandwidth * n ) ; if ( bandwidthInPoints < 2 ) { throw new MathIllegalArgumentException ( LocalizedCoreFormats . BANDWIDTH , bandwidthInPoints , 2 , true ) ; } final double [ ] res = new double [ n ] ; final double [ ] residuals = new double [ n ] ; final double [ ] sortedResiduals = new double [ n ] ; final double [ ] robustnessWeights = new double [ n ] ; Arrays . fill ( robustnessWeights , 1 ) ; for ( int iter = 0 ; iter <= robustnessIters ; ++ iter ) { final int [ ] bandwidthInterval = { 0 , bandwidthInPoints - 1 } ; for ( int i = 0 ; i < n ; ++ i ) { final double x = xval [ i ] ; if ( i > 0 ) { updateBandwidthInterval ( xval , weights , i , bandwidthInterval ) ; } final int ileft = bandwidthInterval [ 0 ] ; final int iright = bandwidthInterval [ 1 ] ; final int edge ; if ( xval [ i ] - xval [ ileft ] > xval [ iright ] - xval [ i ] ) { edge = ileft ; } else { edge = iright ; } double sumWeights = 0 ; double sumX = 0 ; double sumXSquared = 0 ; double sumY = 0 ; double sumXY = 0 ; double denom = FastMath . abs ( 1.0 / ( xval [ edge ] - x ) ) ; for ( int k = ileft ; k <= iright ; ++ k ) { final double xk = xval [ k ] ; final double yk = yval [ k ] ; final double dist = ( k < i ) ? x - xk : xk - x ; final double w = tricube ( dist * denom ) * robustnessWeights [ k ] * weights [ k ] ; final double xkw = xk * w ; sumWeights += w ; sumX += xkw ; sumXSquared += xk * xkw ; sumY += yk * w ; sumXY += yk * xkw ; } final double meanX = sumX / sumWeights ; final double meanY = sumY / sumWeights ; final double meanXY = sumXY / sumWeights ; final double meanXSquared = sumXSquared / sumWeights ; final double beta ; if ( FastMath . sqrt ( FastMath . abs ( meanXSquared - meanX * meanX ) ) < accuracy ) { beta = 0 ; } else { beta = ( meanXY - meanX * meanY ) / ( meanXSquared - meanX * meanX ) ; } final double alpha = meanY - beta * meanX ; res [ i ] = beta * x + alpha ; residuals [ i ] = FastMath . abs ( yval [ i ] - res [ i ] ) ; } if ( iter == robustnessIters ) { break ; } System . arraycopy ( residuals , 0 , sortedResiduals , 0 , n ) ; Arrays . sort ( sortedResiduals ) ; final double medianResidual = sortedResiduals [ n / 2 ] ; if ( FastMath . abs ( medianResidual ) < accuracy ) { break ; } for ( int i = 0 ; i < n ; ++ i ) { final double arg = residuals [ i ] / ( 6 * medianResidual ) ; if ( arg >= 1 ) { robustnessWeights [ i ] = 0 ; } else { final double w = 1 - arg * arg ; robustnessWeights [ i ] = w * w ; } } } return res ; }
compute a weighted loess fit on the data at the original abscissae .
