jdistlib.math.spline

Class SmoothSpline

java.lang.Object

jdistlib.math.spline.SmoothSpline

public class SmoothSpline
extends java.lang.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.

Author:

Roby Joehanes

Field Summary

Fields

Modifier and Type	Field and Description
static double	kDefaultEpsilon
static int	kDefaultMaxNumIterations
static double	kDefaultSmoothingParamLowerBound
static double	kDefaultSmoothingParamUpperBound
static double	kDefaultTolerance

Constructor Summary

Constructors

Constructor and Description
SmoothSpline()

Method Summary

Modifier and Type	Method and Description
static SmoothSplineResult	fit(double[] x, double[] y)
static SmoothSplineResult	fit(double[] x, double[] y, double[] weights)
static SmoothSplineResult	fit(double[] x, double[] y, double[] weights, SmoothSplineCriterion criterion, double penalty, double df)
static SmoothSplineResult	fit(double[] x, double[] y, double[] weights, SmoothSplineCriterion criterion, double penalty, double df, double smoothingParam) Fit a smooth spline
static 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
static SmoothSplineResult	fitDFMatch(double[] x, double[] y, double df)
static SmoothSplineResult	fitDFMatch(double[] x, double[] y, double[] weights, double df)
static double	predict(double[] knots, double[] coefs, double xmin, double xmax, double val, int deriv)
static double	predict(SmoothSplineResult result, double val, int deriv)

Methods inherited from class java.lang.Object

equals, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Field Detail

kDefaultTolerance

public static final double kDefaultTolerance

See Also:

Constant Field Values

kDefaultEpsilon

public static final double kDefaultEpsilon

See Also:

Constant Field Values

kDefaultSmoothingParamLowerBound

public static final double kDefaultSmoothingParamLowerBound

See Also:

Constant Field Values

kDefaultSmoothingParamUpperBound

public static final double kDefaultSmoothingParamUpperBound

See Also:

Constant Field Values

kDefaultMaxNumIterations

public static final int kDefaultMaxNumIterations

See Also:

Constant Field Values

Constructor Detail

SmoothSpline

public SmoothSpline()

Method Detail

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.

Returns: