Cubic Splines are nice but I needed something more generic, i.e. smoothing splines (http://www.qmw.ac.uk/~ugte133/book/11_tsd/splines.pdf). After searching on google, I could not find anything and decided to implement it myself. I believe this would be a nice to have feature and I send you my source code if you decide to include it in your next release:
package org.apache.commons.math.analysis; import org.apache.commons.math.analysis.PolynomialFunction; import org.apache.commons.math.analysis.PolynomialSplineFunction; import org.apache.commons.math.analysis.UnivariateRealFunction; import org.apache.commons.math.analysis.UnivariateRealInterpolator; /** * Smoothing Spline Interpolator based on the algorithm described at * http://www.qmw.ac.uk/~ugte133/book/11_tsd/splines.pdf * <br> * <br> * The Smoothing Spline is used to minimize Sum(sqr((Si-yi)/sigmai)) + lambda*Sum(sqr(S''i)) * <ul> * <li>If lambda=0 (default), this is equivalent to the cubic spline interpolation * <li>If lambda=Infinity, this is equivalent to the least square fitting * </ul> * * @author freyssj */ public class SmoothingSplineInterpolator implements UnivariateRealInterpolator { private double lambda = 0; private double[] sigma = null; private double residuals; private double smoothing; /** * */ private static void quincunx(int n, double[] u, double[] v, double[] w, double[] q) { //Factorisation u[0] = 0; v[0] = 0; w[0] = 0; for (int j = 1; j <= n-1; j++) { u[j] = u[j] - (j-2>=0? u[j-2]*sqr(w[j-2]): 0) - u[j-1] *sqr(v[j-1]); v[j] = (v[j] - u[j-1] * v[j-1] * w[j-1]) / u[j]; w[j] = w[j] / u[j]; } //Forward Substitution q[0] = 0; for (int j = 1; j <= n-1; j++) { q[j] = q[j] - v[j-1] * q[j-1] - (j-2>=0? w[j-2]*q[j-2] :0); } for (int j = 1; j <= n-1; j++) { q[j] = q[j] / u[j]; } //Back Substitution q[n-1] = q[n-1]; q[n-2] = q[n-2] - v[n-2] * q[n-2+1]; for (int j = n-3; j >=1; j--) { q[j] = q[j] - v[j] * q[j+1] - w[j] * q[j+2]; } } /** * @see org.apache.commons.math.analysis.UnivariateRealInterpolator#interpolate(double [], double[]) */ public UnivariateRealFunction interpolate(double[] x, double[] y) { if (x.length != y.length) throw new IllegalArgumentException ("Dataset arrays must have same length."); if (x.length < 3) throw new IllegalArgumentException("At least 3 datapoints are required to compute a spline interpolant"); if (sigma!=null && sigma.length != x.length) throw new IllegalArgumentException("Sigma and dataset arrays must have same length"); int n = x.length; double[] h = new double[n-1]; double[] r = new double[n]; double[] f = new double[n]; double[] p = new double[n]; double[] q = new double[n]; double[] u = new double[n]; double[] v = new double[n]; double[] w = new double[n]; for (int i = 0; i < n-1; i++) { h[i] = x[i + 1] - x[i]; r[i] = 3d / h[i]; } for (int i = 1; i < n-1; i++) { f[i] = -(r[i-1] + r[i]); p[i] = 2d * (x[i + 1] - x[i - 1]); q[i] = 3d * (y[i + 1] - y[i]) / h[i] - 3d * (y[i] - y[i - 1]) / h[i - 1]; } r[n-1] = 0; f[n-1] = 0; for (int i = 1; i <= n-1; i++) { u[i] = sqr(r[i-1])*sigma(i-1) + sqr(f[i])*sigma(i) + sqr(r[i]) *sigma(i + 1); u[i] = lambda * u[i] + p[i]; if(i<n-1) v[i] = f[i] * r[i] * sigma(i) + r[i] * f[i + 1] * sigma(i + 1); if(i<n-1) v[i] = lambda * v[i] + h[i]; if(i<n-1) w[i] = lambda * r[i] * r[i + 1] * sigma(i + 1); } //Solve the system quincunx(n, u, v, w, q); //Spline Parameters residuals = smoothing = 0; double a[] = new double[n]; double b[] = new double[n]; double c[] = new double[n]; double d[] = new double[n]; d[0] = y[0] - lambda * r[0] * q[1] * sigma(0); d[1] = y[1] - lambda * (f[1] * q[1] + r[1] * q[2]) * sigma(1); a[0] = q[1] / (3d * h[0]); b[0] = 0; c[0] = (d[1] - d[0])/h[0] - q[1] * h[0]/3d; r[0] = 0; for (int j = 1; j < n-1; j++) { a[j] = (q[j + 1]-q[j]) / (3d * h[j]); b[j] = q[j]; c[j] = (q[j] + q[j-1]) * h[j-1] + c[j-1]; d[j] = r[j-1] * q[j-1] + f[j] * q[j] + r[j] * q[j+1]; d[j] = y[j] - lambda * d[j] * sigma(j); } for (int j = 0; j < n-1; j++) { if(sigma(j)>0) residuals += (d[j]-y[j]) * (d[j]-y[j]) / (sigma(j) * sigma(j)); smoothing += a[j]*a[j]*h[j]; } PolynomialFunction polynomials[] = new PolynomialFunction[n-1]; double coefficients[] = new double[4]; for (int i = 0; i < polynomials.length; i++) { coefficients[0] = d[i]; coefficients[1] = c[i]; coefficients[2] = b[i]; coefficients[3] = a[i]; polynomials[i] = new PolynomialFunction(coefficients); } return new PolynomialSplineFunction(x, polynomials); } private double getResiduals() { return residuals; } private double getSmoothing() { return smoothing; } public double getLambda() { return lambda; } public double[] getSigma() { return sigma; } public void setLambda(double d) { lambda = d; } public void setSigma(double[] ds) { sigma = ds; } private static final double sqr(double v) { return v*v; } private final double sigma(int index) { if(sigma==null) return 1; return sigma[index]; } } Joel Freyss Research Informatician Actelion Pharmaceuticals Ltd / Gewerbestrasse 16 / CH - 4123 Allschwil / Switzerland / phone +41 61 487 43 34 / [EMAIL PROTECTED] / www.actelion.com ************************************************************************* The information of this e-mail and in any file transmitted with it is strictly confidential and may be legally privileged. It is intended solely for the addressee. If you are not the intended recipient, any copying, distribution or any other use of this e-mail is prohibited and may be unlawful. In such case, you should please notify the sender immediately and destroy this e-mail. The content of this e-mail is not legally binding unless confirmed by letter. Any views expressed in this message are those of the individual sender, except where the message states otherwise and the sender is authorised to state them to be the views of the sender's company. For further information about Actelion please see our website at http://www.actelion.com. ************************************************************************** --------------------------------------------------------------------- To unsubscribe, e-mail: [EMAIL PROTECTED] For additional commands, e-mail: [EMAIL PROTECTED]
