umontreal.iro.lecuyer.functionfit.SmoothingCubicSpline Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation
Show all versions of ssj Show documentation
SSJ is a Java library for stochastic simulation, developed under the direction of Pierre L'Ecuyer,
in the Département d'Informatique et de Recherche Opérationnelle (DIRO), at the Université de Montréal.
It provides facilities for generating uniform and nonuniform random variates, computing different
measures related to probability distributions, performing goodness-of-fit tests, applying quasi-Monte
Carlo methods, collecting (elementary) statistics, and programming discrete-event simulations with both
events and processes.
The newest version!
/*
* Class: SmoothingCubicSpline
* Description: smoothing cubic spline algorithm of Schoenberg
* Environment: Java
* Software: SSJ
* Copyright (C) 2001 Pierre L'Ecuyer and Université de Montréal
* Organization: DIRO, Université de Montréal
* @author
* @since
* SSJ is free software: you can redistribute it and/or modify it under
* the terms of the GNU General Public License (GPL) as published by the
* Free Software Foundation, either version 3 of the License, or
* any later version.
* SSJ is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* A copy of the GNU General Public License is available at
GPL licence site.
*/
package umontreal.iro.lecuyer.functionfit;
import umontreal.iro.lecuyer.functions.*;
import umontreal.iro.lecuyer.functions.Polynomial;
import umontreal.iro.lecuyer.functions.MathFunctionWithFirstDerivative;
import umontreal.iro.lecuyer.functions.MathFunctionWithDerivative;
import umontreal.iro.lecuyer.functions.MathFunctionWithIntegral;
import umontreal.iro.lecuyer.functions.SquareMathFunction;
import umontreal.iro.lecuyer.functions.MathFunctionUtil;
/**
* Represents a cubic spline with nodes at
* (xi, yi) computed with
* the smoothing cubic spline algorithm of Schoenberg.
* A smoothing cubic spline
* is made of n + 1 cubic polynomials. The ith polynomial of such a spline,
* for
* i = 1,…, n - 1, is defined as Si(x) while the complete spline is
* defined as
*
*
*
* S(x) = Si(x), for x∈[xi-1, xi].
*
* For x < x0 and
* x > xn-1, the spline is not precisely defined, but this class performs
* extrapolation by using S0 and Sn linear polynomials.
* The algorithm which calculates the smoothing spline is a
* generalization of the algorithm for an interpolating spline.
* Si is linked to Si+1 at xi+1 and keeps continuity properties for
* first and second derivatives at this point, therefore
*
* Si(xi+1) = Si+1(xi+1),
*
* S'i(xi+1) = S'i+1(xi+1) and
* S''i(xi+1) = S''i+1(xi+1).
*
*
* The spline is computed with a smoothing parameter
* ρ∈[0, 1] which represents its accuracy with respect to the initial
* (xi, yi) nodes. The smoothing spline minimizes
*
*
*
* L = ρ∑i=0n-1wi(yi-Si(xi))2 + (1 - ρ)∫x0xn-1(S''(x))2dx
*
* In fact, by setting ρ = 1, we obtain the interpolating spline; and
* we obtain a linear function by setting ρ = 0.
* The weights wi > 0, which default to 1, can be used to change the contribution
* of each point in the error term. A large value wi will give a large weight
* to the ith point, so the spline will pass closer to it.
* Here is a small example that uses smoothing splines:
*
*
*
*
*
*
*
int n;
*
double[] X = new double[n];
*
double[] Y = new double[n];
*
// here, fill arrays X and Y with n data points (x_i, y_i)
*
// The points must be sorted with respect to x_i.
*
*
double rho = 0.1;
*
SmoothingCubicSpline fit = new SmoothingCubicSpline(X, Y, rho);
*
*
int m = 40;
*
double[] Xp = new double[m+1]; // Xp, Yp are spline points
*
double[] Yp = new double[m+1];
*
double h = (X[n-1] - X[0]) / m; // step
*
*
for (int i = 0; i <= m; i++) {
*
double z = X[0] + i * h;
*
Xp[i] = z;
*
Yp[i] = fit.evaluate(z); // evaluate spline at z
*
}
*
*
*
*
*/
public class SmoothingCubicSpline implements MathFunction,
MathFunctionWithFirstDerivative, MathFunctionWithDerivative,
MathFunctionWithIntegral {
private Polynomial[] splineVector;
private double[] x, y, weight;
private double rho;
/**
* Constructs a spline with nodes at
* (xi, yi),
* with weights wi and smoothing factor ρ = rho.
* The xi must be sorted in increasing order.
*
* @param x the xi coordinates.
*
* @param y the yi coordinates.
*
* @param w the weight for each point, must be > 0.
*
* @param rho the smoothing parameter
*
* @exception IllegalArgumentException if x, y and
* z do not have the same length, if rho has wrong value, or if
* the spline cannot be calculated.
*
*
*/
public SmoothingCubicSpline (double[] x, double[] y, double[] w,
double rho) {
if (x.length != y.length)
throw new IllegalArgumentException ("x.length != y.length");
if (w != null && x.length != w.length)
throw new IllegalArgumentException ("x.length != w.length");
if (rho < 0 || rho > 1)
throw new IllegalArgumentException ("rho not in [0, 1]");
splineVector = new Polynomial[x.length+1];
this.rho = rho;
this.x = x.clone();
this.y = y.clone();
weight = new double[x.length];
if (w == null) {
for (int i = 0; i < weight.length; i++)
weight[i] = 1.0;
} else {
for (int i = 0; i < weight.length; i++)
weight[i] = w[i];
}
resolve();
}
/**
* Constructs a spline with nodes at
* (xi, yi),
* with weights = 1 and smoothing factor ρ = rho.
* The xi must be sorted in increasing order.
*
* @param x the xi coordinates.
*
* @param y the yi coordinates.
*
* @param rho the smoothing parameter
*
* @exception IllegalArgumentException if x and y do
* not have the same length, if rho has wrong value, or if the spline
* cannot be calculated.
*
*/
public SmoothingCubicSpline (double[] x, double[] y, double rho) {
this (x, y, null, rho);
}
/**
* Evaluates and returns the value of the spline at z.
*
* @param z argument of the spline.
*
* @return value of spline.
*
*/
public double evaluate (double z) {
int i = getFitPolynomialIndex(z);
if (i == 0)
return splineVector[i].evaluate(z-x[0]);
else
return splineVector[i].evaluate(z-x[i-1]);
}
/**
* Evaluates and returns the value of the integral of the
* spline from a to b.
*
* @param a lower limit of integral.
*
* @param b upper limit of integral.
*
* @return value of integral.
*
*/
public double integral (double a, double b) {
int iA = getFitPolynomialIndex(a);
int iB = getFitPolynomialIndex(b);
double retour;
int i = 1;
// Calcule l'integrale
if (iA == iB) { // les deux valeurs sont sur le meme polynome
retour = splineVector[iB].integral(a-x[iB], b-x[iB]);
} else {
if (iA == 0)
retour = splineVector[iA].integral(a-x[iA], 0);
else
retour = splineVector[iA].integral(a-x[iA], x[iA+1]-x[iA]);
for (i = iA+1; ifirst derivative of
* the spline at z.
*
* @param z argument of the spline.
*
* @return value of first derivative.
*
*/
public double derivative (double z) {
int i = getFitPolynomialIndex(z);
if (i == 0)
return splineVector[i].derivative(z-x[0]);
else
return splineVector[i].derivative(z-x[i-1]);
}
/**
* Evaluates and returns the value of the n-th derivative of
* the spline at z.
*
* @param z argument of the spline.
*
* @param n order of the derivative.
*
* @return value of n-th derivative.
*
*/
public double derivative (double z, int n) {
int i = getFitPolynomialIndex(z);
if (i == 0)
return splineVector[i].derivative(z - x[0], n);
else
return splineVector[i].derivative(z - x[i-1], n);
}
/**
* Returns the xi coordinates for this spline.
*
* @return the xi coordinates.
*
*/
public double[] getX() {
return x.clone ();
}
/**
* Returns the yi coordinates for this spline.
*
* @return the yi coordinates.
*
*/
public double[] getY() {
return y.clone ();
}
/**
* Returns the weights of the points.
*
* @return the weights of the points.
*
*/
public double[] getWeights() {
return weight;
}
/**
* Returns the smoothing factor used to construct the spline.
*
* @return the smoothing factor.
*
*/
public double getRho() {
return rho;
}
/**
* Returns a table containing all fitting polynomials.
*
* @return Table containing the fitting polynomials.
*
*/
public Polynomial[] getSplinePolynomials() {
return splineVector.clone();
}
/**
* Returns the index of P, the {@link Polynomial} instance used to evaluate
* x, in an ArrayList table instance returned by
* getSplinePolynomials(). This index k gives also the interval in
* table X which contains the value x
* (i.e. such that
* xk < x <= xk+1).
*
* @return Index of the polynomial check with x in the Polynomial list returned by methodgetSplinePolynomials
*
*/
public int getFitPolynomialIndex (double x) {
// Algorithme de recherche binaire legerement modifie
int j = this.x.length-1;
if (x > this.x[j])
return j+1;
int tmp = 0;
int i = 0;
while (i+1 != j) {
if (x > this.x[tmp]) {
i = tmp;
tmp = i+(j-i)/2;
} else {
j = tmp;
tmp = i+(j-i)/2;
}
if (j == 0) // le point est < a x_0, on sort
i--;
}
return i+1;
}
private void resolve () {
/*
taken from D.S.G Pollock's paper, "Smoothing with Cubic Splines",
Queen Mary, University of London (1993)
http://www.qmw.ac.uk/~ugte133/PAPERS/SPLINES.PDF
*/
double[] h = new double[x.length];
double[] r = new double[x.length];
double[] u = new double[x.length];
double[] v = new double[x.length];
double[] w = new double[x.length];
double[] q = new double[x.length+1];
double[] sigma = new double[weight.length];
for (int i = 0; i < weight.length; i++) {
if (weight[i] <= 0.0)
sigma[i] = 1.0e100;
else
sigma[i] = 1.0/Math.sqrt(weight[i]);
}
double mu;
int n = x.length-1;
if (rho <= 0)
mu = 1.0e100; // arbitrary large number to avoid 1/0
else
mu = 2 * (1 - rho)/(3 * rho);
h[0] = x[1] - x[0];
r[0] = 3/h[0];
for (int i = 1; i < n; i++) {
h[i] = x[i+1] - x[i];
r[i] = 3/h[i];
q[i] = 3 * (y[i+1] - y[i])/h[i] - 3 * (y[i] - y[i-1])/h[i - 1];
}
for (int i = 1; i < n; i++) {
u[i] = r[i-1]*r[i-1] * sigma[i-1]
+ (r[i - 1] + r[i])*(r[i - 1] + r[i]) * sigma[i] +
r[i]*r[i] * sigma[i+1];
u[i] = mu * u[i] + 2 * (x[i+1] - x[i-1]);
v[i] = -(r[i - 1] + r[i]) * r[i] * sigma[i] - r[i] * (r[i] +
r[i+1]) * sigma[i+1];
v[i] = mu * v[i] + h[i];
w[i] = mu * r[i] * r[i+1] * sigma[i+1];
}
q = Quincunx(u, v, w, q);
// extrapolation a gauche
double[] params = new double[4];
double dd;
params[0] = y[0] - mu * r[0] * q[1] * sigma[0];
dd = y[1] - mu * ((-r[0] - r[1]) * q[1] + r[1] * q[2]) * sigma[1];
params[1] = (dd - params[0])/h[0] - q[1] * h[0]/3;
splineVector[0] = new Polynomial(params);
// premier polynome
params[0] = y[0] - mu * r[0] * q[1] * sigma[0];
dd = y[1] - mu * ((-r[0] - r[1]) * q[1] + r[1] * q[2]) * sigma[1];
params[3] = q[1]/(3 * h[0]);
params[2] = 0;
params[1] = (dd - params[0])/h[0] - q[1] * h[0]/3;
splineVector[1] = new Polynomial(params);
// les polynomes suivants
int j;
for (j = 1; j < n; j++) {
params[3] = (q[j + 1] - q[j])/(3 * h[j]);
params[2] = q[j];
params[1] = (q[j] + q[j - 1]) * h[j - 1] + splineVector[j].getCoefficient(1);
params[0] = r[j - 1] * q[j - 1] + (-r[j-1] - r[j]) * q[j] + r[j] * q[j + 1];
params[0] = y[j] - mu * params[0] * sigma[j];
splineVector[j+1] = new Polynomial(params);
}
// extrapolation a droite
j = n;
params[3] = 0;
params[2] = 0;
params[1] = splineVector[j].derivative(x[x.length-1]-x[x.length-2]);
params[0] = splineVector[j].evaluate(x[x.length-1]-x[x.length-2]);
splineVector[n+1] = new Polynomial(params);
}
private static double[] Quincunx (double[] u, double[] v,
double[] w, double[] q) {
u[0] = 0;
v[1] = v[1]/u[1];
w[1] = w[1]/u[1];
int j;
for (j = 2; j < u.length-1; j++) {
u[j] = u[j] - u[j - 2] * w[j-2]*w[j-2] - u[j - 1] * v[j-1]*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[1] = q[1] - v[0] * q[0];
for (j = 2; j < u.length-1; j++) {
q[j] = q[j] - v[j - 1] * q[j - 1] - w[j - 2] * q[j - 2];
}
for (j = 1; j < u.length-1; j++) {
q[j] = q[j]/u[j];
}
// back substitution
q[u.length-1] = 0;
for (j = u.length-3; j > 0; j--) {
q[j] = q[j] - v[j] * q[j + 1] - w[j] * q[j + 2];
}
return q;
}
}