Package jdistlib.math.spline

Class SmoothSpline

java.lang.Object
jdistlib.math.spline.SmoothSpline
public class SmoothSpline extends Object

This class deals with smoothing of cubic B-Splines. The algorithm is explained in the following paper:

Silverman, B. W. (1985) Some aspects of the spline smoothing approach to non-parametric regression curve-fitting. Journal of the Royal Statistical Society. Series B (Methodological), 47(1): 1--52.
JSTOR link: http://www.jstor.org/stable/2345542

And the comments of Finbarr O'Sullivan at page 39--40 of the paper above.

The code is in FORTRAN 77 and it is readily available in Netlib:
sbart.f in GCV package.

This code uses the idea explained in the following paper:

de Boor, C. (1977) Package for calculating with B-Splines, SIAM Journal on Numerical Analysis, 14(3): 441--472.
JSTOR link: http://www.jstor.org/stable/2156696

De Boor's code is also in FORTRAN 77 and it is also readily available in Netlib as PPPack package.

De Boor uses LINPACK for Cholesky decomposition and solver for band matrices (dpbfa and dpbsl). LINPACK is also available in Netlib here:
http://www.netlib.org/linpack/

All routines above are available as public domain. They are also included in R, a popular statistical framework. In R, some routines are translated into C and hand-polished. The codes included in R is in GPL. My translation is also in GPL.

Relevant methods:

  • fit = Essentially a call akin to R's smooth.spline. Returns a SmoothSplineResult object that has all outputs you need for prediction.
  • predict = Essentially like R's predict.smooth.spline. Returns the predicted value of the B-Spline

Here's what's different here from the original code and/or R:

  • In the original codes, some of the routines are in single precision (real). In R, they are already converted to double precision. The ones in here are all in double-precision.
  • The original codes and R's code contains a questionable constant value, namely 0.3330 (in sgram). The original paper states that it should be 1/3. So, I used that instead. (See kAThird) Inevitably, the output of the original code and R loses some precision. Note that the output starts to differ only at the fourth or fifth decimal place.
  • Removed superfluous parameters. FORTRAN 77 doesn't allow dynamic arrays. So, all intermediary / buffer arrays have to be passed from the caller. This leads to code ugliness. I created arrays just where needed without seriously affecting the performance. Also, array lengths / dimensions are no longer passed here since Java allows .length for all arrays.
  • Some matrices are better put as double-dimension arrays (most notably sg0, sg1, sg2, sg3, which presents Sigma matrix, and hs0, hs1, hs2, hs3, which represents X'WX matrix)
  • Some loop unrolling is reversed (most notably in sgram). Some of the loops are unrolled to avoid costly if() statements (e.g. in sinerp). Some are left intact (e.g. in stxwx)
  • I changed the derivative array to ROW major. FORTRAN uses column major array convention and it's a pain to pass the array around.
  • Variable / parameter / method renaming. FORTRAN 77's harsh 6-character names no longer applies in Java. So, I put better names around. But I keep the original comments intact, unless it's either superfluous or no longer applies.
  • I uses R's default for epsilon, tolerance, the lower- and upper-bounds for smoothing parameter search, and the maximum number of iterations.

Here's the gotchas:

  • I have tested this routine quite a lot, but far from extensively. I think it should be relatively free of bugs / errors. But, let me know if you do find some discrepancies.
  • Also, despite the testing, some branches of the code have never gotten tested. I put some markers at the code. I've never encountered errors so far.
  • Note that the original comment is still intact. The call parameters are different now since some have been deleted.

  • Field Details

    • kDefaultTolerance

      public static final double kDefaultTolerance
      See Also:

    • kDefaultEpsilon

      public static final double kDefaultEpsilon
      See Also:

    • kDefaultSmoothingParamLowerBound

      public static final double kDefaultSmoothingParamLowerBound
      See Also:

    • kDefaultSmoothingParamUpperBound

      public static final double kDefaultSmoothingParamUpperBound
      See Also:

    • kDefaultMaxNumIterations

      public static final int kDefaultMaxNumIterations
      See Also:

  • Constructor Details

    • SmoothSpline

      public SmoothSpline()

  • Method Details

    • predict

      public static final double predict(SmoothSplineResult result, double val, int deriv)

    • predict

      public static final double predict(double[] knots, double[] coefs, double xmin, double xmax, double val, int deriv)

    • fitDFMatch

      public static final SmoothSplineResult fitDFMatch(double[] x, double[] y, double df)

    • fitDFMatch

      public static final SmoothSplineResult fitDFMatch(double[] x, double[] y, double[] weights, double df)

    • fit

      public static final SmoothSplineResult fit(double[] x, double[] y)

    • fit

      public static final SmoothSplineResult fit(double[] x, double[] y, double[] weights)

    • fit

      public static final SmoothSplineResult fit(double[] x, double[] y, double[] weights, SmoothSplineCriterion criterion, double penalty, double df)

    • fit

      public static final SmoothSplineResult fit(double[] x, double[] y, double[] weights, SmoothSplineCriterion criterion, double penalty, double df, double smoothingParam)
      Fit a smooth spline
      Parameters:
      x - array of n observations
      y - array of n observations
      weights - array of n observations. Null for default weights.
      criterion - one of NO_CRITERION, GCV, CV, and DF_MATCHING
      penalty - Must be 0 < penalty <= 1. Default to 1
      df - Supply if you want DF matching. Else, it will be used for DF offset.
      smoothingParam - Put Double.NaN to estimate it.

    • fit

      public static final SmoothSplineResult fit(double[] x, double[] y, double[] weights, SmoothSplineCriterion criterion, double penalty, double df, double smoothingParam, double smoothingParamLBound, double smoothingParamUBound, double tolerance, int maxNumIterations)
      Fit a smooth spline
      Parameters:
      x - array of n observations
      y - array of n observations
      weights - array of n observations. Null for default weights.
      criterion - one of NO_CRITERION, GCV, CV, and DF_MATCHING
      penalty - Must be 0 < penalty <= 1. Default to 1
      df - Supply if you want DF matching. Else, it will be used for DF offset.
      smoothingParam - Put Double.NaN to estimate it.
      smoothingParamLBound - Will be ignored if smoothingParam is not NaN. Default = -1.5
      smoothingParamUBound - Will be ignored if smoothingParam is not NaN. Default = +1.5
      tolerance - Tolerance threshold. Default is 0.0001.
      maxNumIterations - maximum number of iterations. Default is 500.