Java tutorial
/** * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies * * Please see distribution for license. */ package com.opengamma.analytics.math.statistics.leastsquare; import static org.apache.commons.math.util.MathUtils.binomialCoefficient; import java.util.List; import org.apache.commons.lang.ArrayUtils; import com.google.common.collect.Lists; import com.opengamma.analytics.math.FunctionUtils; import com.opengamma.analytics.math.function.Function1D; import com.opengamma.analytics.math.linearalgebra.Decomposition; import com.opengamma.analytics.math.linearalgebra.DecompositionResult; import com.opengamma.analytics.math.linearalgebra.SVDecompositionCommons; import com.opengamma.analytics.math.matrix.ColtMatrixAlgebra; import com.opengamma.analytics.math.matrix.DoubleMatrix1D; import com.opengamma.analytics.math.matrix.DoubleMatrix2D; import com.opengamma.analytics.math.matrix.DoubleMatrixUtils; import com.opengamma.analytics.math.matrix.MatrixAlgebra; import com.opengamma.util.ArgumentChecker; /** * */ public class GeneralizedLeastSquare { private final Decomposition<?> _decomposition; private final MatrixAlgebra _algebra; public GeneralizedLeastSquare() { _decomposition = new SVDecompositionCommons(); _algebra = new ColtMatrixAlgebra(); } /** * * @param <T> The type of the independent variables (e.g. Double, double[], DoubleMatrix1D etc) * @param x independent variables * @param y dependent (scalar) variables * @param sigma (Gaussian) measurement error on dependent variables * @param basisFunctions set of basis functions - the fitting function is formed by these basis functions times a set of weights * @return the results of the least square */ public <T> GeneralizedLeastSquareResults<T> solve(final T[] x, final double[] y, final double[] sigma, final List<Function1D<T, Double>> basisFunctions) { return solve(x, y, sigma, basisFunctions, 0.0, 0); } /** * Generalised least square with penalty on (higher-order) finite differences of weights * @param <T> The type of the independent variables (e.g. Double, double[], DoubleMatrix1D etc) * @param x independent variables * @param y dependent (scalar) variables * @param sigma (Gaussian) measurement error on dependent variables * @param basisFunctions set of basis functions - the fitting function is formed by these basis functions times a set of weights * @param lambda strength of penalty function * @param differenceOrder difference order between weights used in penalty function * @return the results of the least square */ public <T> GeneralizedLeastSquareResults<T> solve(final T[] x, final double[] y, final double[] sigma, final List<Function1D<T, Double>> basisFunctions, final double lambda, final int differenceOrder) { ArgumentChecker.notNull(x, "x null"); ArgumentChecker.notNull(y, "y null"); ArgumentChecker.notNull(sigma, "sigma null"); ArgumentChecker.notEmpty(basisFunctions, "empty basisFunctions"); final int n = x.length; ArgumentChecker.isTrue(n > 0, "no data"); ArgumentChecker.isTrue(y.length == n, "y wrong length"); ArgumentChecker.isTrue(sigma.length == n, "sigma wrong length"); ArgumentChecker.isTrue(lambda >= 0.0, "negative lambda"); ArgumentChecker.isTrue(differenceOrder >= 0, "difference order"); final List<T> lx = Lists.newArrayList(x); final List<Double> ly = Lists.newArrayList(ArrayUtils.toObject(y)); final List<Double> lsigma = Lists.newArrayList(ArrayUtils.toObject(sigma)); return solveImp(lx, ly, lsigma, basisFunctions, lambda, differenceOrder); } GeneralizedLeastSquareResults<Double> solve(final double[] x, final double[] y, final double[] sigma, final List<Function1D<Double, Double>> basisFunctions, final double lambda, final int differenceOrder) { return solve(ArrayUtils.toObject(x), y, sigma, basisFunctions, lambda, differenceOrder); } /** * * @param <T> The type of the independent variables (e.g. Double, double[], DoubleMatrix1D etc) * @param x independent variables * @param y dependent (scalar) variables * @param sigma (Gaussian) measurement error on dependent variables * @param basisFunctions set of basis functions - the fitting function is formed by these basis functions times a set of weights * @return the results of the least square */ public <T> GeneralizedLeastSquareResults<T> solve(final List<T> x, final List<Double> y, final List<Double> sigma, final List<Function1D<T, Double>> basisFunctions) { return solve(x, y, sigma, basisFunctions, 0.0, 0); } /** * Generalised least square with penalty on (higher-order) finite differences of weights * @param <T> The type of the independent variables (e.g. Double, double[], DoubleMatrix1D etc) * @param x independent variables * @param y dependent (scalar) variables * @param sigma (Gaussian) measurement error on dependent variables * @param basisFunctions set of basis functions - the fitting function is formed by these basis functions times a set of weights * @param lambda strength of penalty function * @param differenceOrder difference order between weights used in penalty function * @return the results of the least square */ public <T> GeneralizedLeastSquareResults<T> solve(final List<T> x, final List<Double> y, final List<Double> sigma, final List<Function1D<T, Double>> basisFunctions, final double lambda, final int differenceOrder) { ArgumentChecker.notEmpty(x, "empty measurement points"); ArgumentChecker.notEmpty(y, "empty measurement values"); ArgumentChecker.notEmpty(sigma, "empty measurement errors"); ArgumentChecker.notEmpty(basisFunctions, "empty basisFunctions"); final int n = x.size(); ArgumentChecker.isTrue(n > 0, "no data"); ArgumentChecker.isTrue(y.size() == n, "y wrong length"); ArgumentChecker.isTrue(sigma.size() == n, "sigma wrong length"); ArgumentChecker.isTrue(lambda >= 0.0, "negative lambda"); ArgumentChecker.isTrue(differenceOrder >= 0, "difference order"); return solveImp(x, y, sigma, basisFunctions, lambda, differenceOrder); } /** * Specialist method used mainly for solving multidimensional P-spline problems where the basis functions (B-splines) span a N-dimension space, and the weights sit on an N-dimension * grid and are treated as a N-order tensor rather than a vector, so k-order differencing is done for each tensor index while varying the other indices. * @param <T> The type of the independent variables (e.g. Double, double[], DoubleMatrix1D etc) * @param x independent variables * @param y dependent (scalar) variables * @param sigma (Gaussian) measurement error on dependent variables * @param basisFunctions set of basis functions - the fitting function is formed by these basis functions times a set of weights * @param sizes The size the weights tensor in each dimension (the product of this must equal the number of basis functions) * @param lambda strength of penalty function in each dimension * @param differenceOrder difference order between weights used in penalty function for each dimension * @return the results of the least square */ public <T> GeneralizedLeastSquareResults<T> solve(final List<T> x, final List<Double> y, final List<Double> sigma, final List<Function1D<T, Double>> basisFunctions, final int[] sizes, final double[] lambda, final int[] differenceOrder) { ArgumentChecker.notEmpty(x, "empty measurement points"); ArgumentChecker.notEmpty(y, "empty measurement values"); ArgumentChecker.notEmpty(sigma, "empty measurement errors"); ArgumentChecker.notEmpty(basisFunctions, "empty basisFunctions"); final int n = x.size(); ArgumentChecker.isTrue(n > 0, "no data"); ArgumentChecker.isTrue(y.size() == n, "y wrong length"); ArgumentChecker.isTrue(sigma.size() == n, "sigma wrong length"); final int dim = sizes.length; ArgumentChecker.isTrue(dim == lambda.length, "number of penalty functions {} must be equal to number of directions {}", lambda.length, dim); ArgumentChecker.isTrue(dim == differenceOrder.length, "number of difference order {} must be equal to number of directions {}", differenceOrder.length, dim); for (int i = 0; i < dim; i++) { ArgumentChecker.isTrue(sizes[i] > 0, "sizes must be >= 1"); ArgumentChecker.isTrue(lambda[i] >= 0.0, "negative lambda"); ArgumentChecker.isTrue(differenceOrder[i] >= 0, "difference order"); } return solveImp(x, y, sigma, basisFunctions, sizes, lambda, differenceOrder); } private <T> GeneralizedLeastSquareResults<T> solveImp(final List<T> x, final List<Double> y, final List<Double> sigma, final List<Function1D<T, Double>> basisFunctions, final double lambda, final int differenceOrder) { final int n = x.size(); final int m = basisFunctions.size(); final double[] b = new double[m]; final double[] invSigmaSqr = new double[n]; final double[][] f = new double[m][n]; int i, j, k; for (i = 0; i < n; i++) { final double temp = sigma.get(i); ArgumentChecker.isTrue(temp > 0, "sigma must be greater than zero"); invSigmaSqr[i] = 1.0 / temp / temp; } for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { f[i][j] = basisFunctions.get(i).evaluate(x.get(j)); } } double sum; for (i = 0; i < m; i++) { sum = 0; for (k = 0; k < n; k++) { sum += y.get(k) * f[i][k] * invSigmaSqr[k]; } b[i] = sum; } final DoubleMatrix1D mb = new DoubleMatrix1D(b); DoubleMatrix2D ma = getAMatrix(f, invSigmaSqr); if (lambda > 0.0) { final DoubleMatrix2D d = getDiffMatrix(m, differenceOrder); ma = (DoubleMatrix2D) _algebra.add(ma, _algebra.scale(d, lambda)); } final DecompositionResult decmp = _decomposition.evaluate(ma); final DoubleMatrix1D w = decmp.solve(mb); final DoubleMatrix2D covar = decmp.solve(DoubleMatrixUtils.getIdentityMatrix2D(m)); double chiSq = 0; for (i = 0; i < n; i++) { double temp = 0; for (k = 0; k < m; k++) { temp += w.getEntry(k) * f[k][i]; } chiSq += FunctionUtils.square(y.get(i) - temp) * invSigmaSqr[i]; } return new GeneralizedLeastSquareResults<>(basisFunctions, chiSq, w, covar); } private <T> GeneralizedLeastSquareResults<T> solveImp(final List<T> x, final List<Double> y, final List<Double> sigma, final List<Function1D<T, Double>> basisFunctions, final int[] sizes, final double[] lambda, final int[] differenceOrder) { final int dim = sizes.length; final int n = x.size(); final int m = basisFunctions.size(); final double[] b = new double[m]; final double[] invSigmaSqr = new double[n]; final double[][] f = new double[m][n]; int i, j, k; for (i = 0; i < n; i++) { final double temp = sigma.get(i); ArgumentChecker.isTrue(temp > 0, "sigma must be great than zero"); invSigmaSqr[i] = 1.0 / temp / temp; } for (i = 0; i < m; i++) { for (j = 0; j < n; j++) { f[i][j] = basisFunctions.get(i).evaluate(x.get(j)); } } double sum; for (i = 0; i < m; i++) { sum = 0; for (k = 0; k < n; k++) { sum += y.get(k) * f[i][k] * invSigmaSqr[k]; } b[i] = sum; } final DoubleMatrix1D mb = new DoubleMatrix1D(b); DoubleMatrix2D ma = getAMatrix(f, invSigmaSqr); for (i = 0; i < dim; i++) { if (lambda[i] > 0.0) { final DoubleMatrix2D d = getDiffMatrix(sizes, differenceOrder[i], i); ma = (DoubleMatrix2D) _algebra.add(ma, _algebra.scale(d, lambda[i])); } } final DecompositionResult decmp = _decomposition.evaluate(ma); final DoubleMatrix1D w = decmp.solve(mb); final DoubleMatrix2D covar = decmp.solve(DoubleMatrixUtils.getIdentityMatrix2D(m)); double chiSq = 0; for (i = 0; i < n; i++) { double temp = 0; for (k = 0; k < m; k++) { temp += w.getEntry(k) * f[k][i]; } chiSq += FunctionUtils.square(y.get(i) - temp) * invSigmaSqr[i]; } return new GeneralizedLeastSquareResults<>(basisFunctions, chiSq, w, covar); } private DoubleMatrix2D getAMatrix(final double[][] funcMatrix, final double[] invSigmaSqr) { final int m = funcMatrix.length; final int n = funcMatrix[0].length; final double[][] a = new double[m][m]; for (int i = 0; i < m; i++) { double sum = 0; for (int k = 0; k < n; k++) { sum += FunctionUtils.square(funcMatrix[i][k]) * invSigmaSqr[k]; } a[i][i] = sum; for (int j = i + 1; j < m; j++) { sum = 0; for (int k = 0; k < n; k++) { sum += funcMatrix[i][k] * funcMatrix[j][k] * invSigmaSqr[k]; } a[i][j] = sum; a[j][i] = sum; } } return new DoubleMatrix2D(a); } private DoubleMatrix2D getDiffMatrix(final int m, final int k) { ArgumentChecker.isTrue(k < m, "difference order too high"); final double[][] data = new double[m][m]; if (m == 0) { for (int i = 0; i < m; i++) { data[i][i] = 1.0; } return new DoubleMatrix2D(data); } final int[] coeff = new int[k + 1]; int sign = 1; for (int i = k; i >= 0; i--) { coeff[i] = (int) (sign * binomialCoefficient(k, i)); sign *= -1; } for (int i = k; i < m; i++) { for (int j = 0; j < k + 1; j++) { data[i][j + i - k] = coeff[j]; } } final DoubleMatrix2D d = new DoubleMatrix2D(data); final DoubleMatrix2D dt = _algebra.getTranspose(d); return (DoubleMatrix2D) _algebra.multiply(dt, d); } private DoubleMatrix2D getDiffMatrix(final int[] size, final int k, final int indices) { final int dim = size.length; final DoubleMatrix2D d = getDiffMatrix(size[indices], k); int preProduct = 1; int postProduct = 1; for (int j = indices + 1; j < dim; j++) { preProduct *= size[j]; } for (int j = 0; j < indices; j++) { postProduct *= size[j]; } DoubleMatrix2D temp = d; if (preProduct != 1) { temp = (DoubleMatrix2D) _algebra.kroneckerProduct(DoubleMatrixUtils.getIdentityMatrix2D(preProduct), temp); } if (postProduct != 1) { temp = (DoubleMatrix2D) _algebra.kroneckerProduct(temp, DoubleMatrixUtils.getIdentityMatrix2D(postProduct)); } return temp; } }