Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

MATH-1008 

Created package "o.a.c.m.fitting.leastsquares" to contain a modified version
of the contents of (to-be-deprecated) "o.a.c.m.optim.nonlinear.vector", the
main purpose being to provide a new "fluent" API (cf. "withXxx" methods).  
Along the way, class "LevenbergMarquardtOptimizer" has been further cleaned    
up (e.g. removing protected fields and deprecated methods and using local
variables instead of instance fields).
An additional constructor was necessary in "BaseOptimizer".


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1508481 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Gilles Sadowski 2013-07-30 15:04:22 +00:00
parent 0215c2c5f5
commit acd569595e
29 changed files with 7092 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

362
src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java Normal file
View File

@ -0,0 +1,362 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MathUnsupportedOperationException;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.optim.ConvergenceChecker;
import org.apache.commons.math3.optim.BaseOptimizer;
import org.apache.commons.math3.optim.PointVectorValuePair;
import org.apache.commons.math3.optim.OptimizationData;
import org.apache.commons.math3.util.FastMath;
/**
* Base class for implementing least-squares optimizers.
* It provides methods for error estimation.
*
* @version $Id$
* @since 3.3
*/
public abstract class AbstractLeastSquaresOptimizer
extends BaseOptimizer<PointVectorValuePair> {
/** Target values for the model function at optimum. */
private final double[] target;
/** Weight matrix. */
private final RealMatrix weight;
/** Model function. */
private final MultivariateVectorFunction model;
/** Jacobian of the model function. */
private final MultivariateMatrixFunction jacobian;
/** Square-root of the weight matrix. */
private final RealMatrix weightSqrt;
/** Initial guess. */
private final double[] start;
/**
* @param target Observations.
* @param weight Weight of the observations.
* For performance, no defensive copy is performed.
* @param weightSqrt Square-root of the {@code weight} matrix.
* If {@code null}, it will be computed; otherwise it is the caller's
* responsibility that {@code weight} and {@code weightSqrt} are
* consistent.
* No defensive copy is performed.
* @param model ModelFunction.
* @param jacobian Jacobian of the model function.
* @param checker Convergence checker.
* @param start Initial guess.
* @param maxEval Maximum number of evaluations of the model
* function.
* @param maxIter Maximum number of iterations.
*/
protected AbstractLeastSquaresOptimizer(double[] target,
RealMatrix weight,
RealMatrix weightSqrt,
MultivariateVectorFunction model,
MultivariateMatrixFunction jacobian,
ConvergenceChecker<PointVectorValuePair> checker,
double[] start,
int maxEval,
int maxIter) {
super(checker, maxEval, maxIter);
this.target = target;
this.weight = weight;
this.model = model;
this.jacobian = jacobian;
this.start = start;
this.weightSqrt = weightSqrt == null ?
(weight == null ?
null : squareRoot(weight)) : weightSqrt;
}
/**
* Gets the target values.
*
* @return the target values.
*/
public double[] getTarget() {
return target == null ? null : target.clone();
}
/**
* Gets the initial guess.
*
* @return the initial guess values.
*/
public double[] getStart() {
return start == null ? null : start.clone();
}
/**
* Gets the square-root of the weight matrix.
*
* @return the square-root of the weight matrix.
*/
public RealMatrix getWeightSquareRoot() {
return weightSqrt == null ? null : weightSqrt.copy();
}
/**
* Gets the model function.
*
* @return the model function.
*/
public MultivariateVectorFunction getModel() {
return model;
}
/**
* Gets the model function's Jacobian.
*
* @return the Jacobian.
*/
public MultivariateMatrixFunction getJacobian() {
return jacobian;
}
/**
* Get the covariance matrix of the optimized parameters.
* <br/>
* Note that this operation involves the inversion of the
* <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
* Jacobian matrix.
* The {@code threshold} parameter is a way for the caller to specify
* that the result of this computation should be considered meaningless,
* and thus trigger an exception.
*
* @param params Model parameters.
* @param threshold Singularity threshold.
* @return the covariance matrix.
* @throws org.apache.commons.math3.linear.SingularMatrixException
* if the covariance matrix cannot be computed (singular problem).
*/
public double[][] computeCovariances(double[] params,
double threshold) {
// Set up the Jacobian.
final RealMatrix j = computeWeightedJacobian(params);
// Compute transpose(J)J.
final RealMatrix jTj = j.transpose().multiply(j);
// Compute the covariances matrix.
final DecompositionSolver solver
= new QRDecomposition(jTj, threshold).getSolver();
return solver.getInverse().getData();
}
/**
* Computes an estimate of the standard deviation of the parameters. The
* returned values are the square root of the diagonal coefficients of the
* covariance matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]}
* is the optimized value of the {@code i}-th parameter, and {@code C} is
* the covariance matrix.
*
* @param params Model parameters.
* @param covarianceSingularityThreshold Singularity threshold (see
* {@link #computeCovariances(double[],double) computeCovariances}).
* @return an estimate of the standard deviation of the optimized parameters
* @throws org.apache.commons.math3.linear.SingularMatrixException
* if the covariance matrix cannot be computed.
*/
public double[] computeSigma(double[] params,
double covarianceSingularityThreshold) {
final int nC = params.length;
final double[] sig = new double[nC];
final double[][] cov = computeCovariances(params, covarianceSingularityThreshold);
for (int i = 0; i < nC; ++i) {
sig[i] = FastMath.sqrt(cov[i][i]);
}
return sig;
}
/**
* Gets the weight matrix of the observations.
*
* @return the weight matrix.
*/
public RealMatrix getWeight() {
return weight.copy();
}
/**
* Computes the normalized cost.
* It is the square-root of the sum of squared of the residuals, divided
* by the number of measurements.
*
* @param params Model function parameters.
* @return the cost.
*/
public double computeRMS(double[] params) {
final double cost = computeCost(computeResiduals(getModel().value(params)));
return FastMath.sqrt(cost * cost / target.length);
}
/**
* Calling this method will raise an exception.
*
* @param optData Obsolete.
* @return nothing.
* @throws MathUnsupportedOperationException if called.
* @deprecated Do not use this method.
*/
@Deprecated
@Override
public PointVectorValuePair optimize(OptimizationData... optData)
throws MathUnsupportedOperationException {
throw new MathUnsupportedOperationException();
}
/**
* Gets a reference to the corresponding field.
* Altering it could jeopardize the consistency of this class.
*
* @return the reference.
*/
protected double[] getTargetInternal() {
return target;
}
/**
* Gets a reference to the corresponding field.
* Altering it could jeopardize the consistency of this class.
*
* @return the reference.
*/
protected RealMatrix getWeightInternal() {
return weight;
}
/**
* Gets a reference to the corresponding field.
* Altering it could jeopardize the consistency of this class.
*
* @return the reference.
*/
protected RealMatrix getWeightSquareRootInternal() {
return weightSqrt;
}
/**
* Computes the objective function value.
* This method <em>must</em> be called by subclasses to enforce the
* evaluation counter limit.
*
* @param params Point at which the objective function must be evaluated.
* @return the objective function value at the specified point.
* @throws org.apache.commons.math3.exception.TooManyEvaluationsException
* if the maximal number of evaluations (of the model vector function) is
* exceeded.
*/
protected double[] computeObjectiveValue(double[] params) {
super.incrementEvaluationCount();
return model.value(params);
}
/**
* Computes the weighted Jacobian matrix.
*
* @param params Model parameters at which to compute the Jacobian.
* @return the weighted Jacobian: W<sup>1/2</sup> J.
* @throws DimensionMismatchException if the Jacobian dimension does not
* match problem dimension.
*/
protected RealMatrix computeWeightedJacobian(double[] params) {
return weightSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian(params)));
}
/**
* Computes the Jacobian matrix.
*
* @param params Point at which the Jacobian must be evaluated.
* @return the Jacobian at the specified point.
*/
protected double[][] computeJacobian(final double[] params) {
return jacobian.value(params);
}
/**
* Computes the cost.
*
* @param residuals Residuals.
* @return the cost.
* @see #computeResiduals(double[])
*/
protected double computeCost(double[] residuals) {
final ArrayRealVector r = new ArrayRealVector(residuals);
return FastMath.sqrt(r.dotProduct(weight.operate(r)));
}
/**
* Computes the residuals.
* The residual is the difference between the observed (target)
* values and the model (objective function) value.
* There is one residual for each element of the vector-valued
* function.
*
* @param objectiveValue Value of the the objective function. This is
* the value returned from a call to
* {@link #computeObjectiveValue(double[]) computeObjectiveValue}
* (whose array argument contains the model parameters).
* @return the residuals.
* @throws DimensionMismatchException if {@code params} has a wrong
* length.
*/
protected double[] computeResiduals(double[] objectiveValue) {
if (objectiveValue.length != target.length) {
throw new DimensionMismatchException(target.length,
objectiveValue.length);
}
final double[] residuals = new double[target.length];
for (int i = 0; i < target.length; i++) {
residuals[i] = target[i] - objectiveValue[i];
}
return residuals;
}
/**
* Computes the square-root of the weight matrix.
*
* @param m Symmetric, positive-definite (weight) matrix.
* @return the square-root of the weight matrix.
*/
private RealMatrix squareRoot(RealMatrix m) {
if (m instanceof DiagonalMatrix) {
final int dim = m.getRowDimension();
final RealMatrix sqrtM = new DiagonalMatrix(dim);
for (int i = 0; i < dim; i++) {
sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
}
return sqrtM;
} else {
final EigenDecomposition dec = new EigenDecomposition(m);
return dec.getSquareRoot();
}
}
}

325
src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java Normal file
View File

@ -0,0 +1,325 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.MathInternalError;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.BlockRealMatrix;
import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.SingularMatrixException;
import org.apache.commons.math3.optim.ConvergenceChecker;
import org.apache.commons.math3.optim.PointVectorValuePair;
/**
* Gauss-Newton least-squares solver.
*
* <p>
* This class solve a least-square problem by solving the normal equations
* of the linearized problem at each iteration. Either LU decomposition or
* QR decomposition can be used to solve the normal equations. LU decomposition
* is faster but QR decomposition is more robust for difficult problems.
* </p>
*
* @version $Id$
* @since 3.3
*
*/
public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer
implements WithTarget<GaussNewtonOptimizer>,
WithWeight<GaussNewtonOptimizer>,
WithModelAndJacobian<GaussNewtonOptimizer>,
WithConvergenceChecker<GaussNewtonOptimizer>,
WithStartPoint<GaussNewtonOptimizer>,
WithMaxIterations<GaussNewtonOptimizer>,
WithMaxEvaluations<GaussNewtonOptimizer> {
/** Indicator for using LU decomposition. */
private final boolean useLU;
/**
* Constructor called by the various {@code withXxx} methods.
*
* @param target Observations.
* @param weight Weight of the observations.
* For performance, no defensive copy is performed.
* @param weightSqrt Square-root of the {@code weight} matrix.
* If {@code null}, it will be computed; otherwise it is the caller's
* responsibility that {@code weight} and {@code weightSqrt} are
* consistent.
* No defensive copy is performed.
* @param model ModelFunction.
* @param jacobian Jacobian of the model function.
* @param checker Convergence checker.
* @param start Initial guess.
* @param maxEval Maximum number of evaluations of the model
* function.
* @param maxIter Maximum number of iterations.
* @param useLU Whether to use LU decomposition.
*/
private GaussNewtonOptimizer(double[] target,
RealMatrix weight,
RealMatrix weightSqrt,
MultivariateVectorFunction model,
MultivariateMatrixFunction jacobian,
ConvergenceChecker<PointVectorValuePair> checker,
double[] start,
int maxEval,
int maxIter,
boolean useLU) {
super(target, weight, weightSqrt, model, jacobian, checker, start, maxEval, maxIter);
this.useLU = useLU;
}
/**
* Creates a bare-bones instance.
* Several calls to {@code withXxx} methods are necessary to obtain
* an object with all necessary fields set to sensible values.
* <br/>
* The default for the algorithm is to solve the normal equations
* using LU decomposition.
*
* @return an instance of this class.
*/
public static GaussNewtonOptimizer create() {
return new GaussNewtonOptimizer(null, null, null, null, null, null, null,
0, 0, true);
}
/** {@inheritDoc} */
public GaussNewtonOptimizer withTarget(double[] target) {
return new GaussNewtonOptimizer(target,
getWeightInternal(),
getWeightSquareRootInternal(),
getModel(),
getJacobian(),
getConvergenceChecker(),
getStart(),
getMaxEvaluations(),
getMaxIterations(),
useLU);
}
/** {@inheritDoc} */
public GaussNewtonOptimizer withWeight(RealMatrix weight) {
return new GaussNewtonOptimizer(getTargetInternal(),
weight,
null,
getModel(),
getJacobian(),
getConvergenceChecker(),
getStart(),
getMaxEvaluations(),
getMaxIterations(),
useLU);
}
/** {@inheritDoc} */
public GaussNewtonOptimizer withModelAndJacobian(MultivariateVectorFunction model,
MultivariateMatrixFunction jacobian) {
return new GaussNewtonOptimizer(getTargetInternal(),
getWeightInternal(),
getWeightSquareRootInternal(),
model,
jacobian,
getConvergenceChecker(),
getStart(),
getMaxEvaluations(),
getMaxIterations(),
useLU);
}
/** {@inheritDoc} */
public GaussNewtonOptimizer withConvergenceChecker(ConvergenceChecker<PointVectorValuePair> checker) {
return new GaussNewtonOptimizer(getTarget(),
getWeightInternal(),
getWeightSquareRootInternal(),
getModel(),
getJacobian(),
checker,
getStart(),
getMaxEvaluations(),
getMaxIterations(),
useLU);
}
/** {@inheritDoc} */
public GaussNewtonOptimizer withStartPoint(double[] start) {
return new GaussNewtonOptimizer(getTarget(),
getWeightInternal(),
getWeightSquareRootInternal(),
getModel(),
getJacobian(),
getConvergenceChecker(),
start,
getMaxEvaluations(),
getMaxIterations(),
useLU);
}
/** {@inheritDoc} */
public GaussNewtonOptimizer withMaxIterations(int maxIter) {
return new GaussNewtonOptimizer(getTarget(),
getWeightInternal(),
getWeightSquareRootInternal(),
getModel(),
getJacobian(),
getConvergenceChecker(),
getStart(),
getMaxEvaluations(),
maxIter,
useLU);
}
/** {@inheritDoc} */
public GaussNewtonOptimizer withMaxEvaluations(int maxEval) {
return new GaussNewtonOptimizer(getTarget(),
getWeightInternal(),
getWeightSquareRootInternal(),
getModel(),
getJacobian(),
getConvergenceChecker(),
getStart(),
maxEval,
getMaxIterations(),
useLU);
}
/**
* Creates a new instance.
*
* @param withLU Whether to use LU decomposition.
* @return a new instance with all fields identical to this instance except
* for the givens arguments.
*/
public GaussNewtonOptimizer withLU(boolean withLU) {
return new GaussNewtonOptimizer(getTarget(),
getWeightInternal(),
getWeightSquareRootInternal(),
getModel(),
getJacobian(),
getConvergenceChecker(),
getStart(),
getMaxEvaluations(),
getMaxIterations(),
withLU);
}
/** {@inheritDoc} */
@Override
public PointVectorValuePair doOptimize() {
final ConvergenceChecker<PointVectorValuePair> checker
= getConvergenceChecker();
// Computation will be useless without a checker (see "for-loop").
if (checker == null) {
throw new NullArgumentException();
}
final double[] targetValues = getTarget();
final int nR = targetValues.length; // Number of observed data.
final RealMatrix weightMatrix = getWeight();
if (weightMatrix.getRowDimension() != nR) {
throw new DimensionMismatchException(weightMatrix.getRowDimension(), nR);
}
if (weightMatrix.getColumnDimension() != nR) {
throw new DimensionMismatchException(weightMatrix.getColumnDimension(), nR);
}
// Diagonal of the weight matrix.
final double[] residualsWeights = new double[nR];
for (int i = 0; i < nR; i++) {
residualsWeights[i] = weightMatrix.getEntry(i, i);
}
final double[] currentPoint = getStart();
final int nC = currentPoint.length;
// iterate until convergence is reached
PointVectorValuePair current = null;
for (boolean converged = false; !converged;) {
incrementIterationCount();
// evaluate the objective function and its jacobian
PointVectorValuePair previous = current;
// Value of the objective function at "currentPoint".
final double[] currentObjective = computeObjectiveValue(currentPoint);
final double[] currentResiduals = computeResiduals(currentObjective);
final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint);
current = new PointVectorValuePair(currentPoint, currentObjective);
// build the linear problem
final double[] b = new double[nC];
final double[][] a = new double[nC][nC];
for (int i = 0; i < nR; ++i) {
final double[] grad = weightedJacobian.getRow(i);
final double weight = residualsWeights[i];
final double residual = currentResiduals[i];
// compute the normal equation
final double wr = weight * residual;
for (int j = 0; j < nC; ++j) {
b[j] += wr * grad[j];
}
// build the contribution matrix for measurement i
for (int k = 0; k < nC; ++k) {
double[] ak = a[k];
double wgk = weight * grad[k];
for (int l = 0; l < nC; ++l) {
ak[l] += wgk * grad[l];
}
}
}
// Check convergence.
if (previous != null) {
converged = checker.converged(getIterations(), previous, current);
if (converged) {
return current;
}
}
try {
// solve the linearized least squares problem
RealMatrix mA = new BlockRealMatrix(a);
DecompositionSolver solver = useLU ?
new LUDecomposition(mA).getSolver() :
new QRDecomposition(mA).getSolver();
final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
// update the estimated parameters
for (int i = 0; i < nC; ++i) {
currentPoint[i] += dX[i];
}
} catch (SingularMatrixException e) {
throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
}
}
// Must never happen.
throw new MathInternalError();
}
}

1172
src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java Normal file
View File

File diff suppressed because it is too large Load Diff

39
src/main/java/org/apache/commons/math3/fitting/leastsquares/WithConvergenceChecker.java Normal file
View File

@ -0,0 +1,39 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import org.apache.commons.math3.optim.ConvergenceChecker;
import org.apache.commons.math3.optim.PointVectorValuePair;
/**
* Interface for "fluent-API" that advertizes a capability of the optimizer.
*
* @param <T> Concrete optimizer implementation.
*
* @version $Id$
* @since 3.3
*/
public interface WithConvergenceChecker<T> {
/**
* Creates a new instance with the specified parameter.
*
* @param checker Convergence checker.
* @return a new optimizer instance with all fields identical to this
* instance except for the given argument.
*/
T withConvergenceChecker(ConvergenceChecker<PointVectorValuePair> checker);
}

36
src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxEvaluations.java Normal file
View File

@ -0,0 +1,36 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
/**
* Interface for "fluent-API" that advertizes a capability of the optimizer.
*
* @param <T> Concrete optimizer implementation.
*
* @version $Id$
* @since 3.3
*/
public interface WithMaxEvaluations<T> {
/**
* Creates a new instance with the specified parameter.
*
* @param maxEval Maximum number of evaluations of the model function.
* @return a new optimizer instance with all fields identical to this
* instance except for the given argument.
*/
T withMaxEvaluations(int maxEval);
}

36
src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxIterations.java Normal file
View File

@ -0,0 +1,36 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
/**
* Interface for "fluent-API" that advertizes a capability of the optimizer.
*
* @param <T> Concrete optimizer implementation.
*
* @version $Id$
* @since 3.3
*/
public interface WithMaxIterations<T> {
/**
* Creates a new instance with the specified parameter.
*
* @param maxIter Maximum number of iterations.
* @return a new optimizer instance with all fields identical to this
* instance except for the given argument.
*/
T withMaxIterations(int maxIter);
}

41
src/main/java/org/apache/commons/math3/fitting/leastsquares/WithModelAndJacobian.java Normal file
View File

@ -0,0 +1,41 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
/**
* Interface for "fluent-API" that advertizes a capability of the optimizer.
*
* @param <T> Concrete optimizer implementation.
*
* @version $Id$
* @since 3.3
*/
public interface WithModelAndJacobian<T> {
/**
* Creates a new instance with the specified parameters.
*
* @param model ModelFunction.
* @param jacobian Jacobian of the model function.
* @return a new optimizer instance with all fields identical to this
* instance except for the given arguments.
*/
T withModelAndJacobian(MultivariateVectorFunction model,
MultivariateMatrixFunction jacobian);
}

36
src/main/java/org/apache/commons/math3/fitting/leastsquares/WithStartPoint.java Normal file
View File

@ -0,0 +1,36 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
/**
* Interface for "fluent-API" that advertizes a capability of the optimizer.
*
* @param <T> Concrete optimizer implementation.
*
* @version $Id$
* @since 3.3
*/
public interface WithStartPoint<T> {
/**
* Creates a new instance with the specified parameter.
*
* @param start Initial guess for the parameters of the model function.
* @return a new optimizer instance with all fields identical to this
* instance except for the given argument.
*/
T withStartPoint(double[] start);
}

36
src/main/java/org/apache/commons/math3/fitting/leastsquares/WithTarget.java Normal file
View File

@ -0,0 +1,36 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
/**
* Interface for "fluent-API" that advertizes a capability of the optimizer.
*
* @param <T> Concrete optimizer implementation.
*
* @version $Id$
* @since 3.3
*/
public interface WithTarget<T> {
/**
* Creates a new instance with the specified parameter.
*
* @param target Objective points of the model function.
* @return a new optimizer instance with all fields identical to this
* instance except for the given argument.
*/
T withTarget(double[] target);
}

38
src/main/java/org/apache/commons/math3/fitting/leastsquares/WithWeight.java Normal file
View File

@ -0,0 +1,38 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import org.apache.commons.math3.linear.RealMatrix;
/**
* Interface for "fluent-API" that advertizes a capability of the optimizer.
*
* @param <T> Concrete optimizer implementation.
*
* @version $Id$
* @since 3.3
*/
public interface WithWeight<T> {
/**
* Creates a new instance with the specified parameter.
*
* @param weight Weight matrix of the observations.
* @return a new optimizer instance with all fields identical to this
* instance except for the given argument.
*/
T withWeight(RealMatrix weight);
}

37
src/main/java/org/apache/commons/math3/fitting/leastsquares/package-info.java Normal file
View File

@ -0,0 +1,37 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
/**
* This package provides algorithms that minimize the residuals
* between observations and model values.
* The {@link org.apache.commons.math3.fitting.leastsquares.AbstractLeastSquaresOptimizer
* non-linear least-squares optimizers} minimize the distance (called <em>cost</em> or
* <em>&chi;<sup>2</sup></em>) between model and observations.
*
* <br/>
* Algorithms in this category need access to a <em>model function</em>
* (represented by a {@link org.apache.commons.math3.analysis.MultivariateVectorFunction
* MultivariateVectorFunction}).
* Such a model predicts a set of values which the algorithm tries to match with a set
* of given set of {@link WithTarget observed values}.
* <br/>
* The algorithms implemented in this package also require that the user specifies the
* Jacobian matrix of the model (represented by a
* {@link org.apache.commons.math3.analysis.MultivariateMatrixFunction
* MultivariateMatrixFunction}).
*/
package org.apache.commons.math3.fitting.leastsquares;

34
src/main/java/org/apache/commons/math3/optim/BaseOptimizer.java
View File

@ -45,10 +45,21 @@ public abstract class BaseOptimizer<PAIR> {
* @param checker Convergence checker.
*/
protected BaseOptimizer(ConvergenceChecker<PAIR> checker) {
this(checker, 0, Integer.MAX_VALUE);
}
/**
* @param checker Convergence checker.
* @param maxEval Maximum number of objective function evaluations.
* @param maxIter Maximum number of algorithm iterations.
*/
protected BaseOptimizer(ConvergenceChecker<PAIR> checker,
int maxEval,
int maxIter) {
this.checker = checker;
evaluations = new Incrementor(0, new MaxEvalCallback());
iterations = new Incrementor(Integer.MAX_VALUE, new MaxIterCallback());
evaluations = new Incrementor(maxEval, new MaxEvalCallback());
iterations = new Incrementor(maxIter, new MaxIterCallback());
}
/**
@ -143,6 +154,25 @@ public abstract class BaseOptimizer<PAIR> {
return doOptimize();
}
/**
* Performs the optimization.
*
* @return a point/value pair that satifies the convergence criteria.
* @throws TooManyEvaluationsException if the maximal number of
* evaluations is exceeded.
* @throws TooManyIterationsException if the maximal number of
* iterations is exceeded.
*/
public PAIR optimize()
throws TooManyEvaluationsException,
TooManyIterationsException {
// Reset counters.
evaluations.resetCount();
iterations.resetCount();
// Perform optimization.
return doOptimize();
}
/**
* Performs the bulk of the optimization algorithm.
*

643
src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerAbstractTest.java Normal file
View File

@ -0,0 +1,643 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.io.IOException;
import java.io.Serializable;
import java.util.Arrays;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
import org.apache.commons.math3.linear.BlockRealMatrix;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.optim.PointVectorValuePair;
import org.apache.commons.math3.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
/**
* Some of the unit tests are re-implementations of the MINPACK
* <a href="http://www.netlib.org/minpack/ex/file17">file17</a> and
* <a href="http://www.netlib.org/minpack/ex/file22">file22</a> test files.
* The redistribution policy for MINPACK is available
* <a href="http://www.netlib.org/minpack/disclaimer">here</a>.
*
* <T> Concrete implementation of an optimizer.
*
* @version $Id$
*/
public abstract class AbstractLeastSquaresOptimizerAbstractTest<T extends AbstractLeastSquaresOptimizer &
WithTarget<T> &
WithWeight<T> &
WithModelAndJacobian<T> &
WithConvergenceChecker<T> &
WithStartPoint<T> &
WithMaxIterations<T> &
WithMaxEvaluations<T>> {
/**
* @return a concrete optimizer.
*/
public abstract T createOptimizer();
/**
* @return the default number of allowed iterations (which will be
* used when not specified otherwise).
*/
public abstract int getMaxIterations();
@Test
public void testGetIterations() {
T optim = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withTarget(new double[] { 1 })
.withWeight(new DiagonalMatrix(new double[] { 1 }))
.withStartPoint(new double[] { 3 })
.withModelAndJacobian(new MultivariateVectorFunction() {
public double[] value(double[] point) {
return new double[] {
FastMath.pow(point[0], 4)
};
}},
new MultivariateMatrixFunction() {
public double[][] value(double[] point) {
return new double[][] {
{ 0.25 * FastMath.pow(point[0], 3) }
};
}
});
optim.optimize();
Assert.assertTrue(optim.getIterations() > 0);
}
@Test
public void testTrivial() {
LinearProblem problem
= new LinearProblem(new double[][] { { 2 } },
new double[] { 3 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1 }))
.withStartPoint(new double[] { 0 });
PointVectorValuePair optimum = optimizer.optimize();
Assert.assertEquals(0, optimizer.computeRMS(optimum.getPoint()), 1e-10);
Assert.assertEquals(1.5, optimum.getPoint()[0], 1e-10);
Assert.assertEquals(3.0, optimum.getValue()[0], 1e-10);
}
@Test
public void testQRColumnsPermutation() {
LinearProblem problem
= new LinearProblem(new double[][] { { 1, -1 }, { 0, 2 }, { 1, -2 } },
new double[] { 4, 6, 1 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 }))
.withStartPoint(new double[] { 0, 0 });
PointVectorValuePair optimum = optimizer.optimize();
Assert.assertEquals(0, optimizer.computeRMS(optimum.getPoint()), 1e-10);
Assert.assertEquals(7, optimum.getPoint()[0], 1e-10);
Assert.assertEquals(3, optimum.getPoint()[1], 1e-10);
Assert.assertEquals(4, optimum.getValue()[0], 1e-10);
Assert.assertEquals(6, optimum.getValue()[1], 1e-10);
Assert.assertEquals(1, optimum.getValue()[2], 1e-10);
}
@Test
public void testNoDependency() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 2, 0, 0, 0, 0, 0 },
{ 0, 2, 0, 0, 0, 0 },
{ 0, 0, 2, 0, 0, 0 },
{ 0, 0, 0, 2, 0, 0 },
{ 0, 0, 0, 0, 2, 0 },
{ 0, 0, 0, 0, 0, 2 }
}, new double[] { 0, 1.1, 2.2, 3.3, 4.4, 5.5 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1, 1 }))
.withStartPoint(new double[] { 0, 0, 0, 0, 0, 0 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10);
for (int i = 0; i < problem.target.length; ++i) {
Assert.assertEquals(0.55 * i, optimum[i], 1e-10);
}
}
@Test
public void testOneSet() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 0, 0 },
{ -1, 1, 0 },
{ 0, -1, 1 }
}, new double[] { 1, 1, 1});
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 }))
.withStartPoint(new double[] { 0, 0, 0 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10);
Assert.assertEquals(1, optimum[0], 1e-10);
Assert.assertEquals(2, optimum[1], 1e-10);
Assert.assertEquals(3, optimum[2], 1e-10);
}
@Test
public void testTwoSets() {
double epsilon = 1e-7;
LinearProblem problem = new LinearProblem(new double[][] {
{ 2, 1, 0, 4, 0, 0 },
{ -4, -2, 3, -7, 0, 0 },
{ 4, 1, -2, 8, 0, 0 },
{ 0, -3, -12, -1, 0, 0 },
{ 0, 0, 0, 0, epsilon, 1 },
{ 0, 0, 0, 0, 1, 1 }
}, new double[] { 2, -9, 2, 2, 1 + epsilon * epsilon, 2});
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1, 1 }))
.withStartPoint(new double[] { 0, 0, 0, 0, 0, 0 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10);
Assert.assertEquals(3, optimum[0], 1e-10);
Assert.assertEquals(4, optimum[1], 1e-10);
Assert.assertEquals(-1, optimum[2], 1e-10);
Assert.assertEquals(-2, optimum[3], 1e-10);
Assert.assertEquals(1 + epsilon, optimum[4], 1e-10);
Assert.assertEquals(1 - epsilon, optimum[5], 1e-10);
}
@Test(expected=ConvergenceException.class)
public void testNonInvertible() throws Exception {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 2, -3 },
{ 2, 1, 3 },
{ -3, 0, -9 }
}, new double[] { 1, 1, 1 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 }))
.withStartPoint(new double[] { 0, 0, 0 });
optimizer.optimize();
}
@Test
public void testIllConditioned() {
LinearProblem problem1 = new LinearProblem(new double[][] {
{ 10, 7, 8, 7 },
{ 7, 5, 6, 5 },
{ 8, 6, 10, 9 },
{ 7, 5, 9, 10 }
}, new double[] { 32, 23, 33, 31 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem1.getModelFunction(),
problem1.getModelFunctionJacobian())
.withTarget(problem1.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1 }))
.withStartPoint(new double[] { 0, 1, 2, 3 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10);
Assert.assertEquals(1, optimum[0], 1e-10);
Assert.assertEquals(1, optimum[1], 1e-10);
Assert.assertEquals(1, optimum[2], 1e-10);
Assert.assertEquals(1, optimum[3], 1e-10);
LinearProblem problem2 = new LinearProblem(new double[][] {
{ 10.00, 7.00, 8.10, 7.20 },
{ 7.08, 5.04, 6.00, 5.00 },
{ 8.00, 5.98, 9.89, 9.00 },
{ 6.99, 4.99, 9.00, 9.98 }
}, new double[] { 32, 23, 33, 31 });
optimizer = optimizer
.withModelAndJacobian(problem2.getModelFunction(),
problem2.getModelFunctionJacobian())
.withTarget(problem2.getTarget());
optimum = optimizer.optimize().getPoint();
Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10);
Assert.assertEquals(-81, optimum[0], 1e-8);
Assert.assertEquals(137, optimum[1], 1e-8);
Assert.assertEquals(-34, optimum[2], 1e-8);
Assert.assertEquals( 22, optimum[3], 1e-8);
}
@Test
public void testMoreEstimatedParametersSimple() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 3, 2, 0, 0 },
{ 0, 1, -1, 1 },
{ 2, 0, 1, 0 }
}, new double[] { 7, 3, 5 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 }))
.withStartPoint(new double[] { 7, 6, 5, 4 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10);
}
@Test
public void testMoreEstimatedParametersUnsorted() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 1, 0, 0, 0, 0 },
{ 0, 0, 1, 1, 1, 0 },
{ 0, 0, 0, 0, 1, -1 },
{ 0, 0, -1, 1, 0, 1 },
{ 0, 0, 0, -1, 1, 0 }
}, new double[] { 3, 12, -1, 7, 1 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1 }))
.withStartPoint(new double[] { 2, 2, 2, 2, 2, 2 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10);
Assert.assertEquals(3, optimum[2], 1e-10);
Assert.assertEquals(4, optimum[3], 1e-10);
Assert.assertEquals(5, optimum[4], 1e-10);
Assert.assertEquals(6, optimum[5], 1e-10);
}
@Test
public void testRedundantEquations() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 1 },
{ 1, -1 },
{ 1, 3 }
}, new double[] { 3, 1, 5 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 }))
.withStartPoint(new double[] { 1, 1 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10);
Assert.assertEquals(2, optimum[0], 1e-10);
Assert.assertEquals(1, optimum[1], 1e-10);
}
@Test
public void testInconsistentEquations() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 1 },
{ 1, -1 },
{ 1, 3 }
}, new double[] { 3, 1, 4 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 }))
.withStartPoint(new double[] { 1, 1 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertTrue(optimizer.computeRMS(optimum) > 0.1);
}
@Test(expected=DimensionMismatchException.class)
public void testInconsistentSizes1() {
LinearProblem problem
= new LinearProblem(new double[][] { { 1, 0 },
{ 0, 1 } },
new double[] { -1, 1 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1 }))
.withStartPoint(new double[] { 0, 0 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10);
Assert.assertEquals(-1, optimum[0], 1e-10);
Assert.assertEquals(1, optimum[1], 1e-10);
optimizer.withWeight(new DiagonalMatrix(new double[] { 1 })).optimize();
}
@Test(expected=DimensionMismatchException.class)
public void testInconsistentSizes2() {
LinearProblem problem
= new LinearProblem(new double[][] { { 1, 0 }, { 0, 1 } },
new double[] { -1, 1 });
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1 }))
.withStartPoint(new double[] { 0, 0 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10);
Assert.assertEquals(-1, optimum[0], 1e-10);
Assert.assertEquals(1, optimum[1], 1e-10);
optimizer
.withTarget(new double[] { 1 })
.withWeight(new DiagonalMatrix(new double[] { 1 }))
.optimize();
}
@Test
public void testCircleFitting() {
CircleVectorial circle = new CircleVectorial();
circle.addPoint( 30, 68);
circle.addPoint( 50, -6);
circle.addPoint(110, -20);
circle.addPoint( 35, 15);
circle.addPoint( 45, 97);
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(circle.getModelFunction(),
circle.getModelFunctionJacobian())
.withTarget(new double[] { 0, 0, 0, 0, 0 })
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1 }))
.withStartPoint(new double[] { 98.680, 47.345 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertTrue(optimizer.getEvaluations() < 10);
double rms = optimizer.computeRMS(optimum);
Assert.assertEquals(1.768262623567235, FastMath.sqrt(circle.getN()) * rms, 1e-10);
Vector2D center = new Vector2D(optimum[0], optimum[1]);
Assert.assertEquals(69.96016176931406, circle.getRadius(center), 1e-6);
Assert.assertEquals(96.07590211815305, center.getX(), 1e-6);
Assert.assertEquals(48.13516790438953, center.getY(), 1e-6);
double[][] cov = optimizer.computeCovariances(optimum, 1e-14);
Assert.assertEquals(1.839, cov[0][0], 0.001);
Assert.assertEquals(0.731, cov[0][1], 0.001);
Assert.assertEquals(cov[0][1], cov[1][0], 1e-14);
Assert.assertEquals(0.786, cov[1][1], 0.001);
// add perfect measurements and check errors are reduced
double r = circle.getRadius(center);
for (double d= 0; d < 2 * FastMath.PI; d += 0.01) {
circle.addPoint(center.getX() + r * FastMath.cos(d), center.getY() + r * FastMath.sin(d));
}
double[] target = new double[circle.getN()];
Arrays.fill(target, 0);
double[] weights = new double[circle.getN()];
Arrays.fill(weights, 2);
optimizer = optimizer.withTarget(target).withWeight(new DiagonalMatrix(weights));
optimum = optimizer.optimize().getPoint();
cov = optimizer.computeCovariances(optimum, 1e-14);
Assert.assertEquals(0.0016, cov[0][0], 0.001);
Assert.assertEquals(3.2e-7, cov[0][1], 1e-9);
Assert.assertEquals(cov[0][1], cov[1][0], 1e-14);
Assert.assertEquals(0.0016, cov[1][1], 0.001);
}
@Test
public void testCircleFittingBadInit() {
CircleVectorial circle = new CircleVectorial();
double[][] points = circlePoints;
double[] target = new double[points.length];
Arrays.fill(target, 0);
double[] weights = new double[points.length];
Arrays.fill(weights, 2);
for (int i = 0; i < points.length; ++i) {
circle.addPoint(points[i][0], points[i][1]);
}
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(circle.getModelFunction(),
circle.getModelFunctionJacobian())
.withTarget(target)
.withWeight(new DiagonalMatrix(weights))
.withStartPoint(new double[] { -12, -12 });
double[] optimum = optimizer.optimize().getPoint();
Vector2D center = new Vector2D(optimum[0], optimum[1]);
Assert.assertTrue(optimizer.getEvaluations() < 25);
Assert.assertEquals( 0.043, optimizer.computeRMS(optimum), 1e-3);
Assert.assertEquals( 0.292235, circle.getRadius(center), 1e-6);
Assert.assertEquals(-0.151738, center.getX(), 1e-6);
Assert.assertEquals( 0.2075001, center.getY(), 1e-6);
}
@Test
public void testCircleFittingGoodInit() {
CircleVectorial circle = new CircleVectorial();
double[][] points = circlePoints;
double[] target = new double[points.length];
Arrays.fill(target, 0);
double[] weights = new double[points.length];
Arrays.fill(weights, 2);
for (int i = 0; i < points.length; ++i) {
circle.addPoint(points[i][0], points[i][1]);
}
T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(circle.getModelFunction(),
circle.getModelFunctionJacobian())
.withTarget(target)
.withWeight(new DiagonalMatrix(weights))
.withStartPoint(new double[] { 0, 0 });
double[] optimum = optimizer.optimize().getPoint();
Assert.assertEquals(-0.1517383071957963, optimum[0], 1e-6);
Assert.assertEquals(0.2074999736353867, optimum[1], 1e-6);
Assert.assertEquals(0.04268731682389561, optimizer.computeRMS(optimum), 1e-8);
}
private final double[][] circlePoints = new double[][] {
{-0.312967, 0.072366}, {-0.339248, 0.132965}, {-0.379780, 0.202724},
{-0.390426, 0.260487}, {-0.361212, 0.328325}, {-0.346039, 0.392619},
{-0.280579, 0.444306}, {-0.216035, 0.470009}, {-0.149127, 0.493832},
{-0.075133, 0.483271}, {-0.007759, 0.452680}, { 0.060071, 0.410235},
{ 0.103037, 0.341076}, { 0.118438, 0.273884}, { 0.131293, 0.192201},
{ 0.115869, 0.129797}, { 0.072223, 0.058396}, { 0.022884, 0.000718},
{-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862},
{-0.278592, -0.005008}, {-0.337655, 0.056658}, {-0.385899, 0.112526},
{-0.405517, 0.186957}, {-0.415374, 0.262071}, {-0.387482, 0.343398},
{-0.347322, 0.397943}, {-0.287623, 0.458425}, {-0.223502, 0.475513},
{-0.135352, 0.478186}, {-0.061221, 0.483371}, { 0.003711, 0.422737},
{ 0.065054, 0.375830}, { 0.108108, 0.297099}, { 0.123882, 0.222850},
{ 0.117729, 0.134382}, { 0.085195, 0.056820}, { 0.029800, -0.019138},
{-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578},
{-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561, 0.014926},
{-0.471036, 0.074716}, {-0.488638, 0.182508}, {-0.485990, 0.254068},
{-0.463943, 0.338438}, {-0.406453, 0.404704}, {-0.334287, 0.466119},
{-0.254244, 0.503188}, {-0.161548, 0.495769}, {-0.075733, 0.495560},
{ 0.001375, 0.434937}, { 0.082787, 0.385806}, { 0.115490, 0.323807},
{ 0.141089, 0.223450}, { 0.138693, 0.131703}, { 0.126415, 0.049174},
{ 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635},
{-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251},
{-0.405195, -0.000895}, {-0.444937, 0.085456}, {-0.484357, 0.175597},
{-0.472453, 0.248681}, {-0.438580, 0.347463}, {-0.402304, 0.422428},
{-0.326777, 0.479438}, {-0.247797, 0.505581}, {-0.152676, 0.519380},
{-0.071754, 0.516264}, { 0.015942, 0.472802}, { 0.076608, 0.419077},
{ 0.127673, 0.330264}, { 0.159951, 0.262150}, { 0.153530, 0.172681},
{ 0.140653, 0.089229}, { 0.078666, 0.024981}, { 0.023807, -0.037022},
{-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526}
};
public void doTestStRD(final StatisticalReferenceDataset dataset,
final double errParams,
final double errParamsSd) {
final double[] w = new double[dataset.getNumObservations()];
Arrays.fill(w, 1);
final double[][] data = dataset.getData();
final double[] initial = dataset.getStartingPoint(0);
final StatisticalReferenceDataset.LeastSquaresProblem problem = dataset.getLeastSquaresProblem();
final T optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(getMaxIterations())
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(data[1])
.withWeight(new DiagonalMatrix(w))
.withStartPoint(initial);
final double[] actual = optimizer.optimize().getPoint();
for (int i = 0; i < actual.length; i++) {
double expected = dataset.getParameter(i);
double delta = FastMath.abs(errParams * expected);
Assert.assertEquals(dataset.getName() + ", param #" + i,
expected, actual[i], delta);
}
}
@Test
public void testKirby2() throws IOException {
doTestStRD(StatisticalReferenceDatasetFactory.createKirby2(), 1E-7, 1E-7);
}
@Test
public void testHahn1() throws IOException {
doTestStRD(StatisticalReferenceDatasetFactory.createHahn1(), 1E-7, 1E-4);
}
static class LinearProblem {
private final RealMatrix factors;
private final double[] target;
public LinearProblem(double[][] factors, double[] target) {
this.factors = new BlockRealMatrix(factors);
this.target = target;
}
public double[] getTarget() {
return target;
}
public MultivariateVectorFunction getModelFunction() {
return new MultivariateVectorFunction() {
public double[] value(double[] params) {
return factors.operate(params);
}
};
}
public MultivariateMatrixFunction getModelFunctionJacobian() {
return new MultivariateMatrixFunction() {
public double[][] value(double[] params) {
return factors.getData();
}
};
}
}
}

109
src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTest.java Normal file
View File

@ -0,0 +1,109 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with this
* work for additional information regarding copyright ownership. The ASF
* licenses this file to You under the Apache License, Version 2.0 (the
* "License"); you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
* http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law
* or agreed to in writing, software distributed under the License is
* distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
* KIND, either express or implied. See the License for the specific language
* governing permissions and limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.io.IOException;
import java.util.Arrays;
import org.apache.commons.math3.optim.PointVectorValuePair;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.util.FastMath;
import org.junit.Test;
import org.junit.Assert;
/**
* The only features tested here are utility methods defined
* in {@link AbstractLeastSquaresOptimizer} that compute the
* chi-square and parameters standard-deviations.
*/
public class AbstractLeastSquaresOptimizerTest {
@Test
public void testComputeCost() throws IOException {
final StatisticalReferenceDataset dataset
= StatisticalReferenceDatasetFactory.createKirby2();
final double[] a = dataset.getParameters();
final double[] y = dataset.getData()[1];
final double[] w = new double[y.length];
Arrays.fill(w, 1d);
StatisticalReferenceDataset.LeastSquaresProblem problem
= dataset.getLeastSquaresProblem();
final LevenbergMarquardtOptimizer optim = LevenbergMarquardtOptimizer.create()
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(y)
.withWeight(new DiagonalMatrix(w))
.withStartPoint(a);
final double expected = dataset.getResidualSumOfSquares();
final double cost = optim.computeCost(optim.computeResiduals(optim.getModel().value(optim.getStart())));
final double actual = cost * cost;
Assert.assertEquals(dataset.getName(), expected, actual, 1e-11 * expected);
}
@Test
public void testComputeRMS() throws IOException {
final StatisticalReferenceDataset dataset
= StatisticalReferenceDatasetFactory.createKirby2();
final double[] a = dataset.getParameters();
final double[] y = dataset.getData()[1];
final double[] w = new double[y.length];
Arrays.fill(w, 1d);
StatisticalReferenceDataset.LeastSquaresProblem problem
= dataset.getLeastSquaresProblem();
final LevenbergMarquardtOptimizer optim = LevenbergMarquardtOptimizer.create()
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(y)
.withWeight(new DiagonalMatrix(w))
.withStartPoint(a);
final double expected = FastMath.sqrt(dataset.getResidualSumOfSquares() /
dataset.getNumObservations());
final double actual = optim.computeRMS(optim.getStart());
Assert.assertEquals(dataset.getName(), expected, actual, 1e-11 * expected);
}
@Test
public void testComputeSigma() throws IOException {
final StatisticalReferenceDataset dataset
= StatisticalReferenceDatasetFactory.createKirby2();
final double[] a = dataset.getParameters();
final double[] y = dataset.getData()[1];
final double[] w = new double[y.length];
Arrays.fill(w, 1d);
StatisticalReferenceDataset.LeastSquaresProblem problem
= dataset.getLeastSquaresProblem();
final LevenbergMarquardtOptimizer optim = LevenbergMarquardtOptimizer.create()
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(y)
.withWeight(new DiagonalMatrix(w))
.withStartPoint(a);
final double[] expected = dataset.getParametersStandardDeviations();
final double cost = optim.computeCost(optim.computeResiduals(optim.getModel().value(optim.getStart())));
final double[] sig = optim.computeSigma(optim.getStart(), 1e-14);
final int dof = y.length - a.length;
for (int i = 0; i < sig.length; i++) {
final double actual = FastMath.sqrt(cost * cost / dof) * sig[i];
Assert.assertEquals(dataset.getName() + ", parameter #" + i,
expected[i], actual, 1e-6 * expected[i]);
}
}
}

308
src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTestValidation.java Normal file
View File

@ -0,0 +1,308 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with this
* work for additional information regarding copyright ownership. The ASF
* licenses this file to You under the Apache License, Version 2.0 (the
* "License"); you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
* http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law
* or agreed to in writing, software distributed under the License is
* distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
* KIND, either express or implied. See the License for the specific language
* governing permissions and limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.util.Arrays;
import java.util.List;
import java.util.ArrayList;
import java.awt.geom.Point2D;
import org.apache.commons.math3.optim.PointVectorValuePair;
import org.apache.commons.math3.stat.descriptive.SummaryStatistics;
import org.apache.commons.math3.stat.descriptive.StatisticalSummary;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.util.FastMath;
import org.junit.Test;
import org.junit.Assert;
/**
* This class demonstrates the main functionality of the
* {@link AbstractLeastSquaresOptimizer}, common to the
* optimizer implementations in package
* {@link org.apache.commons.math3.fitting.leastsquares}.
* <br/>
* Not enabled by default, as the class name does not end with "Test".
* <br/>
* Invoke by running
* <pre><code>
* mvn test -Dtest=AbstractLeastSquaresOptimizerTestValidation
* </code></pre>
* or by running
* <pre><code>
* mvn test -Dtest=AbstractLeastSquaresOptimizerTestValidation -DargLine="-DmcRuns=1234 -server"
* </code></pre>
*/
public class AbstractLeastSquaresOptimizerTestValidation {
/** Number of runs. */
private static final int MONTE_CARLO_RUNS = Integer.parseInt(System.getProperty("mcRuns",
"100"));
/**
* Using a Monte-Carlo procedure, this test checks the error estimations
* as provided by the square-root of the diagonal elements of the
* covariance matrix.
* <br/>
* The test generates sets of observations, each sampled from
* a Gaussian distribution.
* <br/>
* The optimization problem solved is defined in class
* {@link StraightLineProblem}.
* <br/>
* The output (on stdout) will be a table summarizing the distribution
* of parameters generated by the Monte-Carlo process and by the direct
* estimation provided by the diagonal elements of the covariance matrix.
*/
@Test
public void testParametersErrorMonteCarloObservations() {
// Error on the observations.
final double yError = 15;
// True values of the parameters.
final double slope = 123.456;
final double offset = -98.765;
// Samples generator.
final RandomStraightLinePointGenerator lineGenerator
= new RandomStraightLinePointGenerator(slope, offset,
yError,
-1e3, 1e4,
138577L);
// Number of observations.
final int numObs = 100; // XXX Should be a command-line option.
// number of parameters.
final int numParams = 2;
// Parameters found for each of Monte-Carlo run.
final SummaryStatistics[] paramsFoundByDirectSolution = new SummaryStatistics[numParams];
// Sigma estimations (square-root of the diagonal elements of the
// covariance matrix), for each Monte-Carlo run.
final SummaryStatistics[] sigmaEstimate = new SummaryStatistics[numParams];
// Initialize statistics accumulators.
for (int i = 0; i < numParams; i++) {
paramsFoundByDirectSolution[i] = new SummaryStatistics();
sigmaEstimate[i] = new SummaryStatistics();
}
final double[] init = { slope, offset };
// Monte-Carlo (generates many sets of observations).
final int mcRepeat = MONTE_CARLO_RUNS;
int mcCount = 0;
while (mcCount < mcRepeat) {
// Observations.
final Point2D.Double[] obs = lineGenerator.generate(numObs);
final StraightLineProblem problem = new StraightLineProblem(yError);
for (int i = 0; i < numObs; i++) {
final Point2D.Double p = obs[i];
problem.addPoint(p.x, p.y);
}
// Direct solution (using simple regression).
final double[] regress = problem.solve();
// Estimation of the standard deviation (diagonal elements of the
// covariance matrix).
// Dummy optimizer (to compute the covariance matrix).
final AbstractLeastSquaresOptimizer optim = LevenbergMarquardtOptimizer.create()
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.target())
.withWeight(new DiagonalMatrix(problem.weight()));
final double[] sigma = optim.computeSigma(init, 1e-14);
// Accumulate statistics.
for (int i = 0; i < numParams; i++) {
paramsFoundByDirectSolution[i].addValue(regress[i]);
sigmaEstimate[i].addValue(sigma[i]);
}
// Next Monte-Carlo.
++mcCount;
}
// Print statistics.
final String line = "--------------------------------------------------------------";
System.out.println(" True value Mean Std deviation");
for (int i = 0; i < numParams; i++) {
System.out.println(line);
System.out.println("Parameter #" + i);
StatisticalSummary s = paramsFoundByDirectSolution[i].getSummary();
System.out.printf(" %+.6e %+.6e %+.6e\n",
init[i],
s.getMean(),
s.getStandardDeviation());
s = sigmaEstimate[i].getSummary();
System.out.printf("sigma: %+.6e (%+.6e)\n",
s.getMean(),
s.getStandardDeviation());
}
System.out.println(line);
// Check the error estimation.
for (int i = 0; i < numParams; i++) {
Assert.assertEquals(paramsFoundByDirectSolution[i].getSummary().getStandardDeviation(),
sigmaEstimate[i].getSummary().getMean(),
8e-2);
}
}
/**
* In this test, the set of observations is fixed.
* Using a Monte-Carlo procedure, it generates sets of parameters,
* and determine the parameter change that will result in the
* normalized chi-square becoming larger by one than the value from
* the best fit solution.
* <br/>
* The optimization problem solved is defined in class
* {@link StraightLineProblem}.
* <br/>
* The output (on stdout) will be a list of lines containing:
* <ul>
* <li>slope of the straight line,</li>
* <li>intercept of the straight line,</li>
* <li>chi-square of the solution defined by the above two values.</li>
* </ul>
* The output is separated into two blocks (with a blank line between
* them); the first block will contain all parameter sets for which
* {@code chi2 < chi2_b + 1}
* and the second block, all sets for which
* {@code chi2 >= chi2_b + 1}
* where {@code chi2_b} is the lowest chi-square (corresponding to the
* best solution).
*/
@Test
public void testParametersErrorMonteCarloParameters() {
// Error on the observations.
final double yError = 15;
// True values of the parameters.
final double slope = 123.456;
final double offset = -98.765;
// Samples generator.
final RandomStraightLinePointGenerator lineGenerator
= new RandomStraightLinePointGenerator(slope, offset,
yError,
-1e3, 1e4,
13839013L);
// Number of observations.
final int numObs = 10;
// number of parameters.
final int numParams = 2;
// Create a single set of observations.
final Point2D.Double[] obs = lineGenerator.generate(numObs);
final StraightLineProblem problem = new StraightLineProblem(yError);
for (int i = 0; i < numObs; i++) {
final Point2D.Double p = obs[i];
problem.addPoint(p.x, p.y);
}
// Direct solution (using simple regression).
final double[] regress = problem.solve();
// Dummy optimizer (to compute the chi-square).
final AbstractLeastSquaresOptimizer optim = LevenbergMarquardtOptimizer.create()
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.target())
.withWeight(new DiagonalMatrix(problem.weight()));
final double[] init = { slope, offset };
// Get chi-square of the best parameters set for the given set of
// observations.
final double bestChi2N = getChi2N(optim, problem, regress);
final double[] sigma = optim.computeSigma(regress, 1e-14);
// Monte-Carlo (generates a grid of parameters).
final int mcRepeat = MONTE_CARLO_RUNS;
final int gridSize = (int) FastMath.sqrt(mcRepeat);
// Parameters found for each of Monte-Carlo run.
// Index 0 = slope
// Index 1 = offset
// Index 2 = normalized chi2
final List<double[]> paramsAndChi2 = new ArrayList<double[]>(gridSize * gridSize);
final double slopeRange = 10 * sigma[0];
final double offsetRange = 10 * sigma[1];
final double minSlope = slope - 0.5 * slopeRange;
final double minOffset = offset - 0.5 * offsetRange;
final double deltaSlope = slopeRange/ gridSize;
final double deltaOffset = offsetRange / gridSize;
for (int i = 0; i < gridSize; i++) {
final double s = minSlope + i * deltaSlope;
for (int j = 0; j < gridSize; j++) {
final double o = minOffset + j * deltaOffset;
final double chi2N = getChi2N(optim, problem, new double[] {s, o});
paramsAndChi2.add(new double[] {s, o, chi2N});
}
}
// Output (for use with "gnuplot").
// Some info.
// For plotting separately sets of parameters that have a large chi2.
final double chi2NPlusOne = bestChi2N + 1;
int numLarger = 0;
final String lineFmt = "%+.10e %+.10e %.8e\n";
// Point with smallest chi-square.
System.out.printf(lineFmt, regress[0], regress[1], bestChi2N);
System.out.println(); // Empty line.
// Points within the confidence interval.
for (double[] d : paramsAndChi2) {
if (d[2] <= chi2NPlusOne) {
System.out.printf(lineFmt, d[0], d[1], d[2]);
}
}
System.out.println(); // Empty line.
// Points outside the confidence interval.
for (double[] d : paramsAndChi2) {
if (d[2] > chi2NPlusOne) {
++numLarger;
System.out.printf(lineFmt, d[0], d[1], d[2]);
}
}
System.out.println(); // Empty line.
System.out.println("# sigma=" + Arrays.toString(sigma));
System.out.println("# " + numLarger + " sets filtered out");
}
/**
* @return the normalized chi-square.
*/
private double getChi2N(AbstractLeastSquaresOptimizer optim,
StraightLineProblem problem,
double[] params) {
final double[] t = problem.target();
final double[] w = problem.weight();
final double cost = optim.computeCost(optim.computeResiduals(optim.getModel().value(params)));
return cost * cost / (t.length - params.length);
}
}

175
src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleProblem.java Normal file
View File

@ -0,0 +1,175 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.util.ArrayList;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.util.MathUtils;
import org.apache.commons.math3.util.FastMath;
/**
* Class that models a circle.
* The parameters of problem are:
* <ul>
* <li>the x-coordinate of the circle center,</li>
* <li>the y-coordinate of the circle center,</li>
* <li>the radius of the circle.</li>
* </ul>
* The model functions are:
* <ul>
* <li>for each triplet (cx, cy, r), the (x, y) coordinates of a point on the
* corresponding circle.</li>
* </ul>
*/
class CircleProblem {
/** Cloud of points assumed to be fitted by a circle. */
private final ArrayList<double[]> points;
/** Error on the x-coordinate of the points. */
private final double xSigma;
/** Error on the y-coordinate of the points. */
private final double ySigma;
/** Number of points on the circumference (when searching which
model point is closest to a given "observation". */
private final int resolution;
/**
* @param xError Assumed error for the x-coordinate of the circle points.
* @param yError Assumed error for the y-coordinate of the circle points.
* @param searchResolution Number of points to try when searching the one
* that is closest to a given "observed" point.
*/
public CircleProblem(double xError,
double yError,
int searchResolution) {
points = new ArrayList<double[]>();
xSigma = xError;
ySigma = yError;
resolution = searchResolution;
}
/**
* @param xError Assumed error for the x-coordinate of the circle points.
* @param yError Assumed error for the y-coordinate of the circle points.
*/
public CircleProblem(double xError,
double yError) {
this(xError, yError, 500);
}
public void addPoint(double px, double py) {
points.add(new double[] { px, py });
}
public double[] target() {
final double[] t = new double[points.size() * 2];
for (int i = 0; i < points.size(); i++) {
final double[] p = points.get(i);
final int index = i * 2;
t[index] = p[0];
t[index + 1] = p[1];
}
return t;
}
public double[] weight() {
final double wX = 1 / (xSigma * xSigma);
final double wY = 1 / (ySigma * ySigma);
final double[] w = new double[points.size() * 2];
for (int i = 0; i < points.size(); i++) {
final int index = i * 2;
w[index] = wX;
w[index + 1] = wY;
}
return w;
}
public MultivariateVectorFunction getModelFunction() {
return new MultivariateVectorFunction() {
public double[] value(double[] params) {
final double cx = params[0];
final double cy = params[1];
final double r = params[2];
final double[] model = new double[points.size() * 2];
final double deltaTheta = MathUtils.TWO_PI / resolution;
for (int i = 0; i < points.size(); i++) {
final double[] p = points.get(i);
final double px = p[0];
final double py = p[1];
double bestX = 0;
double bestY = 0;
double dMin = Double.POSITIVE_INFINITY;
// Find the angle for which the circle passes closest to the
// current point (using a resolution of 100 points along the
// circumference).
for (double theta = 0; theta <= MathUtils.TWO_PI; theta += deltaTheta) {
final double currentX = cx + r * FastMath.cos(theta);
final double currentY = cy + r * FastMath.sin(theta);
final double dX = currentX - px;
final double dY = currentY - py;
final double d = dX * dX + dY * dY;
if (d < dMin) {
dMin = d;
bestX = currentX;
bestY = currentY;
}
}
final int index = i * 2;
model[index] = bestX;
model[index + 1] = bestY;
}
return model;
}
};
}
public MultivariateMatrixFunction getModelFunctionJacobian() {
return new MultivariateMatrixFunction() {
public double[][] value(double[] point) {
return jacobian(point);
}
};
}
private double[][] jacobian(double[] params) {
final double[][] jacobian = new double[points.size() * 2][3];
for (int i = 0; i < points.size(); i++) {
final int index = i * 2;
// Partial derivative wrt x-coordinate of center.
jacobian[index][0] = 1;
jacobian[index + 1][0] = 0;
// Partial derivative wrt y-coordinate of center.
jacobian[index][1] = 0;
jacobian[index + 1][1] = 1;
// Partial derivative wrt radius.
final double[] p = points.get(i);
jacobian[index][2] = (p[0] - params[0]) / params[2];
jacobian[index + 1][2] = (p[1] - params[1]) / params[2];
}
return jacobian;
}
}

94
src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleVectorial.java Normal file
View File

@ -0,0 +1,94 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.util.ArrayList;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
/**
* Class used in the tests.
*/
class CircleVectorial {
private ArrayList<Vector2D> points;
public CircleVectorial() {
points = new ArrayList<Vector2D>();
}
public void addPoint(double px, double py) {
points.add(new Vector2D(px, py));
}
public int getN() {
return points.size();
}
public double getRadius(Vector2D center) {
double r = 0;
for (Vector2D point : points) {
r += point.distance(center);
}
return r / points.size();
}
public MultivariateVectorFunction getModelFunction() {
return new MultivariateVectorFunction() {
public double[] value(double[] params) {
Vector2D center = new Vector2D(params[0], params[1]);
double radius = getRadius(center);
double[] residuals = new double[points.size()];
for (int i = 0; i < residuals.length; i++) {
residuals[i] = points.get(i).distance(center) - radius;
}
return residuals;
}
};
}
public MultivariateMatrixFunction getModelFunctionJacobian() {
return new MultivariateMatrixFunction() {
public double[][] value(double[] params) {
final int n = points.size();
final Vector2D center = new Vector2D(params[0], params[1]);
double dRdX = 0;
double dRdY = 0;
for (Vector2D pk : points) {
double dk = pk.distance(center);
dRdX += (center.getX() - pk.getX()) / dk;
dRdY += (center.getY() - pk.getY()) / dk;
}
dRdX /= n;
dRdY /= n;
// Jacobian of the radius residuals.
double[][] jacobian = new double[n][2];
for (int i = 0; i < n; i++) {
final Vector2D pi = points.get(i);
final double di = pi.distance(center);
jacobian[i][0] = (center.getX() - pi.getX()) / di - dRdX;
jacobian[i][1] = (center.getY() - pi.getY()) / di - dRdY;
}
return jacobian;
}
};
}
}

109
src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerTest.java Normal file
View File

@ -0,0 +1,109 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.io.IOException;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.exception.TooManyEvaluationsException;
import org.apache.commons.math3.exception.MathUnsupportedOperationException;
import org.apache.commons.math3.optim.SimpleVectorValueChecker;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.junit.Test;
/**
* <p>Some of the unit tests are re-implementations of the MINPACK <a
* href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
* href="http://www.netlib.org/minpack/ex/file22">file22</a> test files.
* The redistribution policy for MINPACK is available <a
* href="http://www.netlib.org/minpack/disclaimer">here</a>/
*
* @version $Id$
*/
public class GaussNewtonOptimizerTest
extends AbstractLeastSquaresOptimizerAbstractTest<GaussNewtonOptimizer> {
@Override
public GaussNewtonOptimizer createOptimizer() {
return GaussNewtonOptimizer.create()
.withConvergenceChecker(new SimpleVectorValueChecker(1e-6, 1e-6));
}
@Override
public int getMaxIterations() {
return 1000;
}
@Override
@Test(expected=ConvergenceException.class)
public void testMoreEstimatedParametersSimple() {
/*
* Exception is expected with this optimizer
*/
super.testMoreEstimatedParametersSimple();
}
@Override
@Test(expected=ConvergenceException.class)
public void testMoreEstimatedParametersUnsorted() {
/*
* Exception is expected with this optimizer
*/
super.testMoreEstimatedParametersUnsorted();
}
@Test(expected=TooManyEvaluationsException.class)
public void testMaxEvaluations() throws Exception {
CircleVectorial circle = new CircleVectorial();
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
circle.addPoint( 35.0, 15.0);
circle.addPoint( 45.0, 97.0);
GaussNewtonOptimizer optimizer = createOptimizer()
.withConvergenceChecker(new SimpleVectorValueChecker(1e-30, 1e-30))
.withMaxIterations(Integer.MAX_VALUE)
.withMaxEvaluations(100)
.withModelAndJacobian(circle.getModelFunction(),
circle.getModelFunctionJacobian())
.withTarget(new double[] { 0, 0, 0, 0, 0 })
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1 }))
.withStartPoint(new double[] { 98.680, 47.345 });
optimizer.optimize();
}
@Override
@Test(expected=ConvergenceException.class)
public void testCircleFittingBadInit() {
/*
* This test does not converge with this optimizer.
*/
super.testCircleFittingBadInit();
}
@Override
@Test(expected=ConvergenceException.class)
public void testHahn1()
throws IOException {
/*
* TODO This test leads to a singular problem with the Gauss-Newton
* optimizer. This should be inquired.
*/
super.testHahn1();
}
}

356
src/test/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizerTest.java Normal file
View File

@ -0,0 +1,356 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.List;
import org.apache.commons.math3.optim.PointVectorValuePair;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.TooManyEvaluationsException;
import org.apache.commons.math3.exception.MathUnsupportedOperationException;
import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
import org.apache.commons.math3.linear.SingularMatrixException;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.Precision;
import org.junit.Assert;
import org.junit.Test;
import org.junit.Ignore;
/**
* <p>Some of the unit tests are re-implementations of the MINPACK <a
* href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
* href="http://www.netlib.org/minpack/ex/file22">file22</a> test files.
* The redistribution policy for MINPACK is available <a
* href="http://www.netlib.org/minpack/disclaimer">here</a>.
*
* @version $Id$
*/
public class LevenbergMarquardtOptimizerTest
extends AbstractLeastSquaresOptimizerAbstractTest<LevenbergMarquardtOptimizer> {
@Override
public LevenbergMarquardtOptimizer createOptimizer() {
return LevenbergMarquardtOptimizer.create();
}
@Override
public int getMaxIterations() {
return 25;
}
@Override
@Test(expected=SingularMatrixException.class)
public void testNonInvertible() {
/*
* Overrides the method from parent class, since the default singularity
* threshold (1e-14) does not trigger the expected exception.
*/
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 2, -3 },
{ 2, 1, 3 },
{ -3, 0, -9 }
}, new double[] { 1, 1, 1 });
final LevenbergMarquardtOptimizer optimizer = createOptimizer()
.withMaxEvaluations(100)
.withMaxIterations(20)
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(problem.getTarget())
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 }))
.withStartPoint(new double[] { 0, 0, 0 });
final double[] optimum = optimizer.optimize().getPoint();
Assert.assertTrue(FastMath.sqrt(optimizer.getTarget().length) * optimizer.computeRMS(optimum) > 0.6);
optimizer.computeCovariances(optimum, 1.5e-14);
}
@Test
public void testControlParameters() {
CircleVectorial circle = new CircleVectorial();
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
circle.addPoint( 35.0, 15.0);
circle.addPoint( 45.0, 97.0);
checkEstimate(circle.getModelFunction(),
circle.getModelFunctionJacobian(),
0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false);
checkEstimate(circle.getModelFunction(),
circle.getModelFunctionJacobian(),
0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true);
checkEstimate(circle.getModelFunction(),
circle.getModelFunctionJacobian(),
0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true);
circle.addPoint(300, -300);
checkEstimate(circle.getModelFunction(),
circle.getModelFunctionJacobian(),
0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true);
}
private void checkEstimate(MultivariateVectorFunction problem,
MultivariateMatrixFunction problemJacobian,
double initialStepBoundFactor, int maxCostEval,
double costRelativeTolerance, double parRelativeTolerance,
double orthoTolerance, boolean shouldFail) {
try {
final LevenbergMarquardtOptimizer optimizer = LevenbergMarquardtOptimizer.create()
.withTuningParameters(initialStepBoundFactor,
costRelativeTolerance,
parRelativeTolerance,
orthoTolerance,
Precision.SAFE_MIN)
.withMaxEvaluations(maxCostEval)
.withMaxIterations(100)
.withModelAndJacobian(problem, problemJacobian)
.withTarget(new double[] { 0, 0, 0, 0, 0 })
.withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1 }))
.withStartPoint(new double[] { 98.680, 47.345 });
optimizer.optimize();
Assert.assertTrue(!shouldFail);
} catch (DimensionMismatchException ee) {
Assert.assertTrue(shouldFail);
} catch (TooManyEvaluationsException ee) {
Assert.assertTrue(shouldFail);
}
}
/**
* Non-linear test case: fitting of decay curve (from Chapter 8 of
* Bevington's textbook, "Data reduction and analysis for the physical sciences").
* XXX The expected ("reference") values may not be accurate and the tolerance too
* relaxed for this test to be currently really useful (the issue is under
* investigation).
*/
@Test
public void testBevington() {
final double[][] dataPoints = {
// column 1 = times
{ 15, 30, 45, 60, 75, 90, 105, 120, 135, 150,
165, 180, 195, 210, 225, 240, 255, 270, 285, 300,
315, 330, 345, 360, 375, 390, 405, 420, 435, 450,
465, 480, 495, 510, 525, 540, 555, 570, 585, 600,
615, 630, 645, 660, 675, 690, 705, 720, 735, 750,
765, 780, 795, 810, 825, 840, 855, 870, 885, },
// column 2 = measured counts
{ 775, 479, 380, 302, 185, 157, 137, 119, 110, 89,
74, 61, 66, 68, 48, 54, 51, 46, 55, 29,
28, 37, 49, 26, 35, 29, 31, 24, 25, 35,
24, 30, 26, 28, 21, 18, 20, 27, 17, 17,
14, 17, 24, 11, 22, 17, 12, 10, 13, 16,
9, 9, 14, 21, 17, 13, 12, 18, 10, },
};
final BevingtonProblem problem = new BevingtonProblem();
final int len = dataPoints[0].length;
final double[] weights = new double[len];
for (int i = 0; i < len; i++) {
problem.addPoint(dataPoints[0][i],
dataPoints[1][i]);
weights[i] = 1 / dataPoints[1][i];
}
final LevenbergMarquardtOptimizer optimizer = LevenbergMarquardtOptimizer.create()
.withMaxEvaluations(100)
.withMaxIterations(20)
.withModelAndJacobian(problem.getModelFunction(),
problem.getModelFunctionJacobian())
.withTarget(dataPoints[1])
.withWeight(new DiagonalMatrix(weights))
.withStartPoint(new double[] { 10, 900, 80, 27, 225 });
final PointVectorValuePair optimum = optimizer.optimize();
final double[] solution = optimum.getPoint();
final double[] expectedSolution = { 10.4, 958.3, 131.4, 33.9, 205.0 };
final double[][] covarMatrix = optimizer.computeCovariances(solution, 1e-14);
final double[][] expectedCovarMatrix = {
{ 3.38, -3.69, 27.98, -2.34, -49.24 },
{ -3.69, 2492.26, 81.89, -69.21, -8.9 },
{ 27.98, 81.89, 468.99, -44.22, -615.44 },
{ -2.34, -69.21, -44.22, 6.39, 53.80 },
{ -49.24, -8.9, -615.44, 53.8, 929.45 }
};
final int numParams = expectedSolution.length;
// Check that the computed solution is within the reference error range.
for (int i = 0; i < numParams; i++) {
final double error = FastMath.sqrt(expectedCovarMatrix[i][i]);
Assert.assertEquals("Parameter " + i, expectedSolution[i], solution[i], error);
}
// Check that each entry of the computed covariance matrix is within 10%
// of the reference matrix entry.
for (int i = 0; i < numParams; i++) {
for (int j = 0; j < numParams; j++) {
Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]",
expectedCovarMatrix[i][j],
covarMatrix[i][j],
FastMath.abs(0.1 * expectedCovarMatrix[i][j]));
}
}
}
@Test
public void testCircleFitting2() {
final double xCenter = 123.456;
final double yCenter = 654.321;
final double xSigma = 10;
final double ySigma = 15;
final double radius = 111.111;
// The test is extremely sensitive to the seed.
final long seed = 59421061L;
final RandomCirclePointGenerator factory
= new RandomCirclePointGenerator(xCenter, yCenter, radius,
xSigma, ySigma,
seed);
final CircleProblem circle = new CircleProblem(xSigma, ySigma);
final int numPoints = 10;
for (Vector2D p : factory.generate(numPoints)) {
circle.addPoint(p.getX(), p.getY());
}
// First guess for the center's coordinates and radius.
final double[] init = { 90, 659, 115 };
final LevenbergMarquardtOptimizer optimizer = LevenbergMarquardtOptimizer.create()
.withMaxEvaluations(100)
.withMaxIterations(50)
.withModelAndJacobian(circle.getModelFunction(),
circle.getModelFunctionJacobian())
.withTarget(circle.target())
.withWeight(new DiagonalMatrix(circle.weight()))
.withStartPoint(init);
final PointVectorValuePair optimum = optimizer.optimize();
final double[] paramFound = optimum.getPoint();
// Retrieve errors estimation.
final double[] asymptoticStandardErrorFound = optimizer.computeSigma(paramFound, 1e-14);
// Check that the parameters are found within the assumed error bars.
Assert.assertEquals(xCenter, paramFound[0], asymptoticStandardErrorFound[0]);
Assert.assertEquals(yCenter, paramFound[1], asymptoticStandardErrorFound[1]);
Assert.assertEquals(radius, paramFound[2], asymptoticStandardErrorFound[2]);
}
private static class QuadraticProblem {
private List<Double> x;
private List<Double> y;
public QuadraticProblem() {
x = new ArrayList<Double>();
y = new ArrayList<Double>();
}
public void addPoint(double x, double y) {
this.x.add(x);
this.y.add(y);
}
public MultivariateVectorFunction getModelFunction() {
return new MultivariateVectorFunction() {
public double[] value(double[] variables) {
double[] values = new double[x.size()];
for (int i = 0; i < values.length; ++i) {
values[i] = (variables[0] * x.get(i) + variables[1]) * x.get(i) + variables[2];
}
return values;
}
};
}
public MultivariateMatrixFunction getModelFunctionJacobian() {
return new MultivariateMatrixFunction() {
public double[][] value(double[] params) {
double[][] jacobian = new double[x.size()][3];
for (int i = 0; i < jacobian.length; ++i) {
jacobian[i][0] = x.get(i) * x.get(i);
jacobian[i][1] = x.get(i);
jacobian[i][2] = 1.0;
}
return jacobian;
}
};
}
}
private static class BevingtonProblem {
private List<Double> time;
private List<Double> count;
public BevingtonProblem() {
time = new ArrayList<Double>();
count = new ArrayList<Double>();
}
public void addPoint(double t, double c) {
time.add(t);
count.add(c);
}
public MultivariateVectorFunction getModelFunction() {
return new MultivariateVectorFunction() {
public double[] value(double[] params) {
double[] values = new double[time.size()];
for (int i = 0; i < values.length; ++i) {
final double t = time.get(i);
values[i] = params[0] +
params[1] * Math.exp(-t / params[3]) +
params[2] * Math.exp(-t / params[4]);
}
return values;
}
};
}
public MultivariateMatrixFunction getModelFunctionJacobian() {
return new MultivariateMatrixFunction() {
public double[][] value(double[] params) {
double[][] jacobian = new double[time.size()][5];
for (int i = 0; i < jacobian.length; ++i) {
final double t = time.get(i);
jacobian[i][0] = 1;
final double p3 = params[3];
final double p4 = params[4];
final double tOp3 = t / p3;
final double tOp4 = t / p4;
jacobian[i][1] = Math.exp(-tOp3);
jacobian[i][2] = Math.exp(-tOp4);
jacobian[i][3] = params[1] * Math.exp(-tOp3) * tOp3 / p3;
jacobian[i][4] = params[2] * Math.exp(-tOp4) * tOp4 / p4;
}
return jacobian;
}
};
}
}
}

1499
src/test/java/org/apache/commons/math3/fitting/leastsquares/MinpackTest.java Normal file
View File

File diff suppressed because it is too large Load Diff

91
src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomCirclePointGenerator.java Normal file
View File

@ -0,0 +1,91 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well44497b;
import org.apache.commons.math3.util.MathUtils;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.distribution.RealDistribution;
import org.apache.commons.math3.distribution.UniformRealDistribution;
import org.apache.commons.math3.distribution.NormalDistribution;
import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
/**
* Factory for generating a cloud of points that approximate a circle.
*/
public class RandomCirclePointGenerator {
/** RNG for the x-coordinate of the center. */
private final RealDistribution cX;
/** RNG for the y-coordinate of the center. */
private final RealDistribution cY;
/** RNG for the parametric position of the point. */
private final RealDistribution tP;
/** Radius of the circle. */
private final double radius;
/**
* @param x Abscissa of the circle center.
* @param y Ordinate of the circle center.
* @param radius Radius of the circle.
* @param xSigma Error on the x-coordinate of the circumference points.
* @param ySigma Error on the y-coordinate of the circumference points.
* @param seed RNG seed.
*/
public RandomCirclePointGenerator(double x,
double y,
double radius,
double xSigma,
double ySigma,
long seed) {
final RandomGenerator rng = new Well44497b(seed);
this.radius = radius;
cX = new NormalDistribution(rng, x, xSigma,
NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
cY = new NormalDistribution(rng, y, ySigma,
NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
tP = new UniformRealDistribution(rng, 0, MathUtils.TWO_PI,
UniformRealDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
}
/**
* Point generator.
*
* @param n Number of points to create.
* @return the cloud of {@code n} points.
*/
public Vector2D[] generate(int n) {
final Vector2D[] cloud = new Vector2D[n];
for (int i = 0; i < n; i++) {
cloud[i] = create();
}
return cloud;
}
/**
* Create one point.
*
* @return a point.
*/
private Vector2D create() {
final double t = tP.sample();
final double pX = cX.sample() + radius * FastMath.cos(t);
final double pY = cY.sample() + radius * FastMath.sin(t);
return new Vector2D(pX, pY);
}
}

98
src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomStraightLinePointGenerator.java Normal file
View File

@ -0,0 +1,98 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.awt.geom.Point2D;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well44497b;
import org.apache.commons.math3.distribution.RealDistribution;
import org.apache.commons.math3.distribution.UniformRealDistribution;
import org.apache.commons.math3.distribution.NormalDistribution;
/**
* Factory for generating a cloud of points that approximate a straight line.
*/
public class RandomStraightLinePointGenerator {
/** Slope. */
private final double slope;
/** Intercept. */
private final double intercept;
/** RNG for the x-coordinate. */
private final RealDistribution x;
/** RNG for the error on the y-coordinate. */
private final RealDistribution error;
/**
* The generator will create a cloud of points whose x-coordinates
* will be randomly sampled between {@code xLo} and {@code xHi}, and
* the corresponding y-coordinates will be computed as
* <pre><code>
* y = a x + b + N(0, error)
* </code></pre>
* where {@code N(mean, sigma)} is a Gaussian distribution with the
* given mean and standard deviation.
*
* @param a Slope.
* @param b Intercept.
* @param sigma Standard deviation on the y-coordinate of the point.
* @param lo Lowest value of the x-coordinate.
* @param hi Highest value of the x-coordinate.
* @param seed RNG seed.
*/
public RandomStraightLinePointGenerator(double a,
double b,
double sigma,
double lo,
double hi,
long seed) {
final RandomGenerator rng = new Well44497b(seed);
slope = a;
intercept = b;
error = new NormalDistribution(rng, 0, sigma,
NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
x = new UniformRealDistribution(rng, lo, hi,
UniformRealDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
}
/**
* Point generator.
*
* @param n Number of points to create.
* @return the cloud of {@code n} points.
*/
public Point2D.Double[] generate(int n) {
final Point2D.Double[] cloud = new Point2D.Double[n];
for (int i = 0; i < n; i++) {
cloud[i] = create();
}
return cloud;
}
/**
* Create one point.
*
* @return a point.
*/
private Point2D.Double create() {
final double abscissa = x.sample();
final double yModel = slope * abscissa + intercept;
final double ordinate = yModel + error.sample();
return new Point2D.Double(abscissa, ordinate);
}
}

370
src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDataset.java Normal file
View File

@ -0,0 +1,370 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.io.BufferedReader;
import java.io.IOException;
import java.util.ArrayList;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.util.MathArrays;
/**
* This class gives access to the statistical reference datasets provided by the
* NIST (available
* <a href="http://www.itl.nist.gov/div898/strd/general/dataarchive.html">here</a>).
* Instances of this class can be created by invocation of the
* {@link StatisticalReferenceDatasetFactory}.
*/
public abstract class StatisticalReferenceDataset {
/** The name of this dataset. */
private final String name;
/** The total number of observations (data points). */
private final int numObservations;
/** The total number of parameters. */
private final int numParameters;
/** The total number of starting points for the optimizations. */
private final int numStartingPoints;
/** The values of the predictor. */
private final double[] x;
/** The values of the response. */
private final double[] y;
/**
* The starting values. {@code startingValues[j][i]} is the value of the
* {@code i}-th parameter in the {@code j}-th set of starting values.
*/
private final double[][] startingValues;
/** The certified values of the parameters. */
private final double[] a;
/** The certified values of the standard deviation of the parameters. */
private final double[] sigA;
/** The certified value of the residual sum of squares. */
private double residualSumOfSquares;
/** The least-squares problem. */
private final LeastSquaresProblem problem;
/**
* Creates a new instance of this class from the specified data file. The
* file must follow the StRD format.
*
* @param in the data file
* @throws IOException if an I/O error occurs
*/
public StatisticalReferenceDataset(final BufferedReader in)
throws IOException {
final ArrayList<String> lines = new ArrayList<String>();
for (String line = in.readLine(); line != null; line = in.readLine()) {
lines.add(line);
}
int[] index = findLineNumbers("Data", lines);
if (index == null) {
throw new AssertionError("could not find line indices for data");
}
this.numObservations = index[1] - index[0] + 1;
this.x = new double[this.numObservations];
this.y = new double[this.numObservations];
for (int i = 0; i < this.numObservations; i++) {
final String line = lines.get(index[0] + i - 1);
final String[] tokens = line.trim().split(" ++");
// Data columns are in reverse order!!!
this.y[i] = Double.parseDouble(tokens[0]);
this.x[i] = Double.parseDouble(tokens[1]);
}
index = findLineNumbers("Starting Values", lines);
if (index == null) {
throw new AssertionError(
"could not find line indices for starting values");
}
this.numParameters = index[1] - index[0] + 1;
double[][] start = null;
this.a = new double[numParameters];
this.sigA = new double[numParameters];
for (int i = 0; i < numParameters; i++) {
final String line = lines.get(index[0] + i - 1);
final String[] tokens = line.trim().split(" ++");
if (start == null) {
start = new double[tokens.length - 4][numParameters];
}
for (int j = 2; j < tokens.length - 2; j++) {
start[j - 2][i] = Double.parseDouble(tokens[j]);
}
this.a[i] = Double.parseDouble(tokens[tokens.length - 2]);
this.sigA[i] = Double.parseDouble(tokens[tokens.length - 1]);
}
if (start == null) {
throw new IOException("could not find starting values");
}
this.numStartingPoints = start.length;
this.startingValues = start;
double dummyDouble = Double.NaN;
String dummyString = null;
for (String line : lines) {
if (line.contains("Dataset Name:")) {
dummyString = line
.substring(line.indexOf("Dataset Name:") + 13,
line.indexOf("(")).trim();
}
if (line.contains("Residual Sum of Squares")) {
final String[] tokens = line.split(" ++");
dummyDouble = Double.parseDouble(tokens[4].trim());
}
}
if (Double.isNaN(dummyDouble)) {
throw new IOException(
"could not find certified value of residual sum of squares");
}
this.residualSumOfSquares = dummyDouble;
if (dummyString == null) {
throw new IOException("could not find dataset name");
}
this.name = dummyString;
this.problem = new LeastSquaresProblem();
}
class LeastSquaresProblem {
public MultivariateVectorFunction getModelFunction() {
return new MultivariateVectorFunction() {
public double[] value(final double[] a) {
final int n = getNumObservations();
final double[] yhat = new double[n];
for (int i = 0; i < n; i++) {
yhat[i] = getModelValue(getX(i), a);
}
return yhat;
}
};
}
public MultivariateMatrixFunction getModelFunctionJacobian() {
return new MultivariateMatrixFunction() {
public double[][] value(final double[] a)
throws IllegalArgumentException {
final int n = getNumObservations();
final double[][] j = new double[n][];
for (int i = 0; i < n; i++) {
j[i] = getModelDerivatives(getX(i), a);
}
return j;
}
};
}
}
/**
* Returns the name of this dataset.
*
* @return the name of the dataset
*/
public String getName() {
return name;
}
/**
* Returns the total number of observations (data points).
*
* @return the number of observations
*/
public int getNumObservations() {
return numObservations;
}
/**
* Returns a copy of the data arrays. The data is laid out as follows <li>
* {@code data[0][i] = x[i]},</li> <li>{@code data[1][i] = y[i]},</li>
*
* @return the array of data points.
*/
public double[][] getData() {
return new double[][] {
MathArrays.copyOf(x), MathArrays.copyOf(y)
};
}
/**
* Returns the x-value of the {@code i}-th data point.
*
* @param i the index of the data point
* @return the x-value
*/
public double getX(final int i) {
return x[i];
}
/**
* Returns the y-value of the {@code i}-th data point.
*
* @param i the index of the data point
* @return the y-value
*/
public double getY(final int i) {
return y[i];
}
/**
* Returns the total number of parameters.
*
* @return the number of parameters
*/
public int getNumParameters() {
return numParameters;
}
/**
* Returns the certified values of the paramters.
*
* @return the values of the parameters
*/
public double[] getParameters() {
return MathArrays.copyOf(a);
}
/**
* Returns the certified value of the {@code i}-th parameter.
*
* @param i the index of the parameter
* @return the value of the parameter
*/
public double getParameter(final int i) {
return a[i];
}
/**
* Reurns the certified values of the standard deviations of the parameters.
*
* @return the standard deviations of the parameters
*/
public double[] getParametersStandardDeviations() {
return MathArrays.copyOf(sigA);
}
/**
* Returns the certified value of the standard deviation of the {@code i}-th
* parameter.
*
* @param i the index of the parameter
* @return the standard deviation of the parameter
*/
public double getParameterStandardDeviation(final int i) {
return sigA[i];
}
/**
* Returns the certified value of the residual sum of squares.
*
* @return the residual sum of squares
*/
public double getResidualSumOfSquares() {
return residualSumOfSquares;
}
/**
* Returns the total number of starting points (initial guesses for the
* optimization process).
*
* @return the number of starting points
*/
public int getNumStartingPoints() {
return numStartingPoints;
}
/**
* Returns the {@code i}-th set of initial values of the parameters.
*
* @param i the index of the starting point
* @return the starting point
*/
public double[] getStartingPoint(final int i) {
return MathArrays.copyOf(startingValues[i]);
}
/**
* Returns the least-squares problem corresponding to fitting the model to
* the specified data.
*
* @return the least-squares problem
*/
public LeastSquaresProblem getLeastSquaresProblem() {
return problem;
}
/**
* Returns the value of the model for the specified values of the predictor
* variable and the parameters.
*
* @param x the predictor variable
* @param a the parameters
* @return the value of the model
*/
public abstract double getModelValue(final double x, final double[] a);
/**
* Returns the values of the partial derivatives of the model with respect
* to the parameters.
*
* @param x the predictor variable
* @param a the parameters
* @return the partial derivatives
*/
public abstract double[] getModelDerivatives(final double x,
final double[] a);
/**
* <p>
* Parses the specified text lines, and extracts the indices of the first
* and last lines of the data defined by the specified {@code key}. This key
* must be one of
* </p>
* <ul>
* <li>{@code "Starting Values"},</li>
* <li>{@code "Certified Values"},</li>
* <li>{@code "Data"}.</li>
* </ul>
* <p>
* In the NIST data files, the line indices are separated by the keywords
* {@code "lines"} and {@code "to"}.
* </p>
*
* @param lines the line of text to be parsed
* @return an array of two {@code int}s. First value is the index of the
* first line, second value is the index of the last line.
* {@code null} if the line could not be parsed.
*/
private static int[] findLineNumbers(final String key,
final Iterable<String> lines) {
for (String text : lines) {
boolean flag = text.contains(key) && text.contains("lines") &&
text.contains("to") && text.contains(")");
if (flag) {
final int[] numbers = new int[2];
final String from = text.substring(text.indexOf("lines") + 5,
text.indexOf("to"));
numbers[0] = Integer.parseInt(from.trim());
final String to = text.substring(text.indexOf("to") + 2,
text.indexOf(")"));
numbers[1] = Integer.parseInt(to.trim());
return numbers;
}
}
return null;
}
}

201
src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDatasetFactory.java Normal file
View File

@ -0,0 +1,201 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import org.apache.commons.math3.util.FastMath;
/**
* A factory to create instances of {@link StatisticalReferenceDataset} from
* available resources.
*/
public class StatisticalReferenceDatasetFactory {
private StatisticalReferenceDatasetFactory() {
// Do nothing
}
/**
* Creates a new buffered reader from the specified resource name.
*
* @param name the name of the resource
* @return a buffered reader
* @throws IOException if an I/O error occured
*/
public static BufferedReader createBufferedReaderFromResource(final String name)
throws IOException {
final InputStream resourceAsStream;
resourceAsStream = StatisticalReferenceDatasetFactory.class
.getResourceAsStream(name);
if (resourceAsStream == null) {
throw new IOException("could not find resource " + name);
}
return new BufferedReader(new InputStreamReader(resourceAsStream));
}
public static StatisticalReferenceDataset createKirby2()
throws IOException {
final BufferedReader in = createBufferedReaderFromResource("Kirby2.dat");
StatisticalReferenceDataset dataset = null;
try {
dataset = new StatisticalReferenceDataset(in) {
@Override
public double getModelValue(final double x, final double[] a) {
final double p = a[0] + x * (a[1] + x * a[2]);
final double q = 1.0 + x * (a[3] + x * a[4]);
return p / q;
}
@Override
public double[] getModelDerivatives(final double x,
final double[] a) {
final double[] dy = new double[5];
final double p = a[0] + x * (a[1] + x * a[2]);
final double q = 1.0 + x * (a[3] + x * a[4]);
dy[0] = 1.0 / q;
dy[1] = x / q;
dy[2] = x * dy[1];
dy[3] = -x * p / (q * q);
dy[4] = x * dy[3];
return dy;
}
};
} finally {
in.close();
}
return dataset;
}
public static StatisticalReferenceDataset createHahn1()
throws IOException {
final BufferedReader in = createBufferedReaderFromResource("Hahn1.dat");
StatisticalReferenceDataset dataset = null;
try {
dataset = new StatisticalReferenceDataset(in) {
@Override
public double getModelValue(final double x, final double[] a) {
final double p = a[0] + x * (a[1] + x * (a[2] + x * a[3]));
final double q = 1.0 + x * (a[4] + x * (a[5] + x * a[6]));
return p / q;
}
@Override
public double[] getModelDerivatives(final double x,
final double[] a) {
final double[] dy = new double[7];
final double p = a[0] + x * (a[1] + x * (a[2] + x * a[3]));
final double q = 1.0 + x * (a[4] + x * (a[5] + x * a[6]));
dy[0] = 1.0 / q;
dy[1] = x * dy[0];
dy[2] = x * dy[1];
dy[3] = x * dy[2];
dy[4] = -x * p / (q * q);
dy[5] = x * dy[4];
dy[6] = x * dy[5];
return dy;
}
};
} finally {
in.close();
}
return dataset;
}
public static StatisticalReferenceDataset createMGH17()
throws IOException {
final BufferedReader in = createBufferedReaderFromResource("MGH17.dat");
StatisticalReferenceDataset dataset = null;
try {
dataset = new StatisticalReferenceDataset(in) {
@Override
public double getModelValue(final double x, final double[] a) {
return a[0] + a[1] * FastMath.exp(-a[3] * x) + a[2] *
FastMath.exp(-a[4] * x);
}
@Override
public double[] getModelDerivatives(final double x,
final double[] a) {
final double[] dy = new double[5];
dy[0] = 1.0;
dy[1] = FastMath.exp(-x * a[3]);
dy[2] = FastMath.exp(-x * a[4]);
dy[3] = -x * a[1] * dy[1];
dy[4] = -x * a[2] * dy[2];
return dy;
}
};
} finally {
in.close();
}
return dataset;
}
public static StatisticalReferenceDataset createLanczos1()
throws IOException {
final BufferedReader in =
createBufferedReaderFromResource("Lanczos1.dat");
StatisticalReferenceDataset dataset = null;
try {
dataset = new StatisticalReferenceDataset(in) {
@Override
public double getModelValue(final double x, final double[] a) {
System.out.println(a[0]+", "+a[1]+", "+a[2]+", "+a[3]+", "+a[4]+", "+a[5]);
return a[0] * FastMath.exp(-a[3] * x) +
a[1] * FastMath.exp(-a[4] * x) +
a[2] * FastMath.exp(-a[5] * x);
}
@Override
public double[] getModelDerivatives(final double x,
final double[] a) {
final double[] dy = new double[6];
dy[0] = FastMath.exp(-x * a[3]);
dy[1] = FastMath.exp(-x * a[4]);
dy[2] = FastMath.exp(-x * a[5]);
dy[3] = -x * a[0] * dy[0];
dy[4] = -x * a[1] * dy[1];
dy[5] = -x * a[2] * dy[2];
return dy;
}
};
} finally {
in.close();
}
return dataset;
}
/**
* Returns an array with all available reference datasets.
*
* @return the array of datasets
* @throws IOException if an I/O error occurs
*/
public StatisticalReferenceDataset[] createAll()
throws IOException {
return new StatisticalReferenceDataset[] {
createKirby2(), createMGH17()
};
}
}

165
src/test/java/org/apache/commons/math3/fitting/leastsquares/StraightLineProblem.java Normal file
View File

@ -0,0 +1,165 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.fitting.leastsquares;
import java.util.ArrayList;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.stat.regression.SimpleRegression;
/**
* Class that models a straight line defined as {@code y = a x + b}.
* The parameters of problem are:
* <ul>
* <li>{@code a}</li>
* <li>{@code b}</li>
* </ul>
* The model functions are:
* <ul>
* <li>for each pair (a, b), the y-coordinate of the line.</li>
* </ul>
*/
class StraightLineProblem {
/** Cloud of points assumed to be fitted by a straight line. */
private final ArrayList<double[]> points;
/** Error (on the y-coordinate of the points). */
private final double sigma;
/**
* @param error Assumed error for the y-coordinate.
*/
public StraightLineProblem(double error) {
points = new ArrayList<double[]>();
sigma = error;
}
public void addPoint(double px, double py) {
points.add(new double[] { px, py });
}
/**
* @return the list of x-coordinates.
*/
public double[] x() {
final double[] v = new double[points.size()];
for (int i = 0; i < points.size(); i++) {
final double[] p = points.get(i);
v[i] = p[0]; // x-coordinate.
}
return v;
}
/**
* @return the list of y-coordinates.
*/
public double[] y() {
final double[] v = new double[points.size()];
for (int i = 0; i < points.size(); i++) {
final double[] p = points.get(i);
v[i] = p[1]; // y-coordinate.
}
return v;
}
public double[] target() {
return y();
}
public double[] weight() {
final double weight = 1 / (sigma * sigma);
final double[] w = new double[points.size()];
for (int i = 0; i < points.size(); i++) {
w[i] = weight;
}
return w;
}
public MultivariateVectorFunction getModelFunction() {
return new MultivariateVectorFunction() {
public double[] value(double[] params) {
final Model line = new Model(params[0], params[1]);
final double[] model = new double[points.size()];
for (int i = 0; i < points.size(); i++) {
final double[] p = points.get(i);
model[i] = line.value(p[0]);
}
return model;
}
};
}
public MultivariateMatrixFunction getModelFunctionJacobian() {
return new MultivariateMatrixFunction() {
public double[][] value(double[] point) {
return jacobian(point);
}
};
}
/**
* Directly solve the linear problem, using the {@link SimpleRegression}
* class.
*/
public double[] solve() {
final SimpleRegression regress = new SimpleRegression(true);
for (double[] d : points) {
regress.addData(d[0], d[1]);
}
final double[] result = { regress.getSlope(), regress.getIntercept() };
return result;
}
private double[][] jacobian(double[] params) {
final double[][] jacobian = new double[points.size()][2];
for (int i = 0; i < points.size(); i++) {
final double[] p = points.get(i);
// Partial derivative wrt "a".
jacobian[i][0] = p[0];
// Partial derivative wrt "b".
jacobian[i][1] = 1;
}
return jacobian;
}
/**
* Linear function.
*/
public static class Model implements UnivariateFunction {
final double a;
final double b;
public Model(double a,
double b) {
this.a = a;
this.b = b;
}
public double value(double x) {
return a * x + b;
}
}
}

296
src/test/resources/org/apache/commons/math3/fitting/leastsquares/Hahn1.dat Normal file
View File

@ -0,0 +1,296 @@
NIST/ITL StRD
Dataset Name: Hahn1 (Hahn1.dat)
File Format: ASCII
Starting Values (lines 41 to 47)
Certified Values (lines 41 to 52)
Data (lines 61 to 296)
Procedure: Nonlinear Least Squares Regression
Description: These data are the result of a NIST study involving
the thermal expansion of copper. The response
variable is the coefficient of thermal expansion, and
the predictor variable is temperature in degrees
kelvin.
Reference: Hahn, T., NIST (197?).
Copper Thermal Expansion Study.
Data: 1 Response (y = coefficient of thermal expansion)
1 Predictor (x = temperature, degrees kelvin)
236 Observations
Average Level of Difficulty
Observed Data
Model: Rational Class (cubic/cubic)
7 Parameters (b1 to b7)
y = (b1+b2*x+b3*x**2+b4*x**3) /
(1+b5*x+b6*x**2+b7*x**3) + e
Starting values Certified Values
Start 1 Start 2 Parameter Standard Deviation
b1 = 10 1 1.0776351733E+00 1.7070154742E-01
b2 = -1 -0.1 -1.2269296921E-01 1.2000289189E-02
b3 = 0.05 0.005 4.0863750610E-03 2.2508314937E-04
b4 = -0.00001 -0.000001 -1.4262662514E-06 2.7578037666E-07
b5 = -0.05 -0.005 -5.7609940901E-03 2.4712888219E-04
b6 = 0.001 0.0001 2.4053735503E-04 1.0449373768E-05
b7 = -0.000001 -0.0000001 -1.2314450199E-07 1.3027335327E-08
Residual Sum of Squares: 1.5324382854E+00
Residual Standard Deviation: 8.1803852243E-02
Degrees of Freedom: 229
Number of Observations: 236
Data: y x
.591E0 24.41E0
1.547E0 34.82E0
2.902E0 44.09E0
2.894E0 45.07E0
4.703E0 54.98E0
6.307E0 65.51E0
7.03E0 70.53E0
7.898E0 75.70E0
9.470E0 89.57E0
9.484E0 91.14E0
10.072E0 96.40E0
10.163E0 97.19E0
11.615E0 114.26E0
12.005E0 120.25E0
12.478E0 127.08E0
12.982E0 133.55E0
12.970E0 133.61E0
13.926E0 158.67E0
14.452E0 172.74E0
14.404E0 171.31E0
15.190E0 202.14E0
15.550E0 220.55E0
15.528E0 221.05E0
15.499E0 221.39E0
16.131E0 250.99E0
16.438E0 268.99E0
16.387E0 271.80E0
16.549E0 271.97E0
16.872E0 321.31E0
16.830E0 321.69E0
16.926E0 330.14E0
16.907E0 333.03E0
16.966E0 333.47E0
17.060E0 340.77E0
17.122E0 345.65E0
17.311E0 373.11E0
17.355E0 373.79E0
17.668E0 411.82E0
17.767E0 419.51E0
17.803E0 421.59E0
17.765E0 422.02E0
17.768E0 422.47E0
17.736E0 422.61E0
17.858E0 441.75E0
17.877E0 447.41E0
17.912E0 448.7E0
18.046E0 472.89E0
18.085E0 476.69E0
18.291E0 522.47E0
18.357E0 522.62E0
18.426E0 524.43E0
18.584E0 546.75E0
18.610E0 549.53E0
18.870E0 575.29E0
18.795E0 576.00E0
19.111E0 625.55E0
.367E0 20.15E0
.796E0 28.78E0
0.892E0 29.57E0
1.903E0 37.41E0
2.150E0 39.12E0
3.697E0 50.24E0
5.870E0 61.38E0
6.421E0 66.25E0
7.422E0 73.42E0
9.944E0 95.52E0
11.023E0 107.32E0
11.87E0 122.04E0
12.786E0 134.03E0
14.067E0 163.19E0
13.974E0 163.48E0
14.462E0 175.70E0
14.464E0 179.86E0
15.381E0 211.27E0
15.483E0 217.78E0
15.59E0 219.14E0
16.075E0 262.52E0
16.347E0 268.01E0
16.181E0 268.62E0
16.915E0 336.25E0
17.003E0 337.23E0
16.978E0 339.33E0
17.756E0 427.38E0
17.808E0 428.58E0
17.868E0 432.68E0
18.481E0 528.99E0
18.486E0 531.08E0
19.090E0 628.34E0
16.062E0 253.24E0
16.337E0 273.13E0
16.345E0 273.66E0
16.388E0 282.10E0
17.159E0 346.62E0
17.116E0 347.19E0
17.164E0 348.78E0
17.123E0 351.18E0
17.979E0 450.10E0
17.974E0 450.35E0
18.007E0 451.92E0
17.993E0 455.56E0
18.523E0 552.22E0
18.669E0 553.56E0
18.617E0 555.74E0
19.371E0 652.59E0
19.330E0 656.20E0
0.080E0 14.13E0
0.248E0 20.41E0
1.089E0 31.30E0
1.418E0 33.84E0
2.278E0 39.70E0
3.624E0 48.83E0
4.574E0 54.50E0
5.556E0 60.41E0
7.267E0 72.77E0
7.695E0 75.25E0
9.136E0 86.84E0
9.959E0 94.88E0
9.957E0 96.40E0
11.600E0 117.37E0
13.138E0 139.08E0
13.564E0 147.73E0
13.871E0 158.63E0
13.994E0 161.84E0
14.947E0 192.11E0
15.473E0 206.76E0
15.379E0 209.07E0
15.455E0 213.32E0
15.908E0 226.44E0
16.114E0 237.12E0
17.071E0 330.90E0
17.135E0 358.72E0
17.282E0 370.77E0
17.368E0 372.72E0
17.483E0 396.24E0
17.764E0 416.59E0
18.185E0 484.02E0
18.271E0 495.47E0
18.236E0 514.78E0
18.237E0 515.65E0
18.523E0 519.47E0
18.627E0 544.47E0
18.665E0 560.11E0
19.086E0 620.77E0
0.214E0 18.97E0
0.943E0 28.93E0
1.429E0 33.91E0
2.241E0 40.03E0
2.951E0 44.66E0
3.782E0 49.87E0
4.757E0 55.16E0
5.602E0 60.90E0
7.169E0 72.08E0
8.920E0 85.15E0
10.055E0 97.06E0
12.035E0 119.63E0
12.861E0 133.27E0
13.436E0 143.84E0
14.167E0 161.91E0
14.755E0 180.67E0
15.168E0 198.44E0
15.651E0 226.86E0
15.746E0 229.65E0
16.216E0 258.27E0
16.445E0 273.77E0
16.965E0 339.15E0
17.121E0 350.13E0
17.206E0 362.75E0
17.250E0 371.03E0
17.339E0 393.32E0
17.793E0 448.53E0
18.123E0 473.78E0
18.49E0 511.12E0
18.566E0 524.70E0
18.645E0 548.75E0
18.706E0 551.64E0
18.924E0 574.02E0
19.1E0 623.86E0
0.375E0 21.46E0
0.471E0 24.33E0
1.504E0 33.43E0
2.204E0 39.22E0
2.813E0 44.18E0
4.765E0 55.02E0
9.835E0 94.33E0
10.040E0 96.44E0
11.946E0 118.82E0
12.596E0 128.48E0
13.303E0 141.94E0
13.922E0 156.92E0
14.440E0 171.65E0
14.951E0 190.00E0
15.627E0 223.26E0
15.639E0 223.88E0
15.814E0 231.50E0
16.315E0 265.05E0
16.334E0 269.44E0
16.430E0 271.78E0
16.423E0 273.46E0
17.024E0 334.61E0
17.009E0 339.79E0
17.165E0 349.52E0
17.134E0 358.18E0
17.349E0 377.98E0
17.576E0 394.77E0
17.848E0 429.66E0
18.090E0 468.22E0
18.276E0 487.27E0
18.404E0 519.54E0
18.519E0 523.03E0
19.133E0 612.99E0
19.074E0 638.59E0
19.239E0 641.36E0
19.280E0 622.05E0
19.101E0 631.50E0
19.398E0 663.97E0
19.252E0 646.9E0
19.89E0 748.29E0
20.007E0 749.21E0
19.929E0 750.14E0
19.268E0 647.04E0
19.324E0 646.89E0
20.049E0 746.9E0
20.107E0 748.43E0
20.062E0 747.35E0
20.065E0 749.27E0
19.286E0 647.61E0
19.972E0 747.78E0
20.088E0 750.51E0
20.743E0 851.37E0
20.83E0 845.97E0
20.935E0 847.54E0
21.035E0 849.93E0
20.93E0 851.61E0
21.074E0 849.75E0
21.085E0 850.98E0
20.935E0 848.23E0

211
src/test/resources/org/apache/commons/math3/fitting/leastsquares/Kirby2.dat Normal file
View File

@ -0,0 +1,211 @@
NIST/ITL StRD
Dataset Name: Kirby2 (Kirby2.dat)
File Format: ASCII
Starting Values (lines 41 to 45)
Certified Values (lines 41 to 50)
Data (lines 61 to 211)
Procedure: Nonlinear Least Squares Regression
Description: These data are the result of a NIST study involving
scanning electron microscope line with standards.
Reference: Kirby, R., NIST (197?).
Scanning electron microscope line width standards.
Data: 1 Response (y)
1 Predictor (x)
151 Observations
Average Level of Difficulty
Observed Data
Model: Rational Class (quadratic/quadratic)
5 Parameters (b1 to b5)
y = (b1 + b2*x + b3*x**2) /
(1 + b4*x + b5*x**2) + e
Starting values Certified Values
Start 1 Start 2 Parameter Standard Deviation
b1 = 2 1.5 1.6745063063E+00 8.7989634338E-02
b2 = -0.1 -0.15 -1.3927397867E-01 4.1182041386E-03
b3 = 0.003 0.0025 2.5961181191E-03 4.1856520458E-05
b4 = -0.001 -0.0015 -1.7241811870E-03 5.8931897355E-05
b5 = 0.00001 0.00002 2.1664802578E-05 2.0129761919E-07
Residual Sum of Squares: 3.9050739624E+00
Residual Standard Deviation: 1.6354535131E-01
Degrees of Freedom: 146
Number of Observations: 151
Data: y x
0.0082E0 9.65E0
0.0112E0 10.74E0
0.0149E0 11.81E0
0.0198E0 12.88E0
0.0248E0 14.06E0
0.0324E0 15.28E0
0.0420E0 16.63E0
0.0549E0 18.19E0
0.0719E0 19.88E0
0.0963E0 21.84E0
0.1291E0 24.00E0
0.1710E0 26.25E0
0.2314E0 28.86E0
0.3227E0 31.85E0
0.4809E0 35.79E0
0.7084E0 40.18E0
1.0220E0 44.74E0
1.4580E0 49.53E0
1.9520E0 53.94E0
2.5410E0 58.29E0
3.2230E0 62.63E0
3.9990E0 67.03E0
4.8520E0 71.25E0
5.7320E0 75.22E0
6.7270E0 79.33E0
7.8350E0 83.56E0
9.0250E0 87.75E0
10.2670E0 91.93E0
11.5780E0 96.10E0
12.9440E0 100.28E0
14.3770E0 104.46E0
15.8560E0 108.66E0
17.3310E0 112.71E0
18.8850E0 116.88E0
20.5750E0 121.33E0
22.3200E0 125.79E0
22.3030E0 125.79E0
23.4600E0 128.74E0
24.0600E0 130.27E0
25.2720E0 133.33E0
25.8530E0 134.79E0
27.1100E0 137.93E0
27.6580E0 139.33E0
28.9240E0 142.46E0
29.5110E0 143.90E0
30.7100E0 146.91E0
31.3500E0 148.51E0
32.5200E0 151.41E0
33.2300E0 153.17E0
34.3300E0 155.97E0
35.0600E0 157.76E0
36.1700E0 160.56E0
36.8400E0 162.30E0
38.0100E0 165.21E0
38.6700E0 166.90E0
39.8700E0 169.92E0
40.0300E0 170.32E0
40.5000E0 171.54E0
41.3700E0 173.79E0
41.6700E0 174.57E0
42.3100E0 176.25E0
42.7300E0 177.34E0
43.4600E0 179.19E0
44.1400E0 181.02E0
44.5500E0 182.08E0
45.2200E0 183.88E0
45.9200E0 185.75E0
46.3000E0 186.80E0
47.0000E0 188.63E0
47.6800E0 190.45E0
48.0600E0 191.48E0
48.7400E0 193.35E0
49.4100E0 195.22E0
49.7600E0 196.23E0
50.4300E0 198.05E0
51.1100E0 199.97E0
51.5000E0 201.06E0
52.1200E0 202.83E0
52.7600E0 204.69E0
53.1800E0 205.86E0
53.7800E0 207.58E0
54.4600E0 209.50E0
54.8300E0 210.65E0
55.4000E0 212.33E0
56.4300E0 215.43E0
57.0300E0 217.16E0
58.0000E0 220.21E0
58.6100E0 221.98E0
59.5800E0 225.06E0
60.1100E0 226.79E0
61.1000E0 229.92E0
61.6500E0 231.69E0
62.5900E0 234.77E0
63.1200E0 236.60E0
64.0300E0 239.63E0
64.6200E0 241.50E0
65.4900E0 244.48E0
66.0300E0 246.40E0
66.8900E0 249.35E0
67.4200E0 251.32E0
68.2300E0 254.22E0
68.7700E0 256.24E0
69.5900E0 259.11E0
70.1100E0 261.18E0
70.8600E0 264.02E0
71.4300E0 266.13E0
72.1600E0 268.94E0
72.7000E0 271.09E0
73.4000E0 273.87E0
73.9300E0 276.08E0
74.6000E0 278.83E0
75.1600E0 281.08E0
75.8200E0 283.81E0
76.3400E0 286.11E0
76.9800E0 288.81E0
77.4800E0 291.08E0
78.0800E0 293.75E0
78.6000E0 295.99E0
79.1700E0 298.64E0
79.6200E0 300.84E0
79.8800E0 302.02E0
80.1900E0 303.48E0
80.6600E0 305.65E0
81.2200E0 308.27E0
81.6600E0 310.41E0
82.1600E0 313.01E0
82.5900E0 315.12E0
83.1400E0 317.71E0
83.5000E0 319.79E0
84.0000E0 322.36E0
84.4000E0 324.42E0
84.8900E0 326.98E0
85.2600E0 329.01E0
85.7400E0 331.56E0
86.0700E0 333.56E0
86.5400E0 336.10E0
86.8900E0 338.08E0
87.3200E0 340.60E0
87.6500E0 342.57E0
88.1000E0 345.08E0
88.4300E0 347.02E0
88.8300E0 349.52E0
89.1200E0 351.44E0
89.5400E0 353.93E0
89.8500E0 355.83E0
90.2500E0 358.32E0
90.5500E0 360.20E0
90.9300E0 362.67E0
91.2000E0 364.53E0
91.5500E0 367.00E0
92.2000E0 371.30E0

84
src/test/resources/org/apache/commons/math3/fitting/leastsquares/Lanczos1.dat Normal file
View File

@ -0,0 +1,84 @@
NIST/ITL StRD
Dataset Name: Lanczos1 (Lanczos1.dat)
File Format: ASCII
Starting Values (lines 41 to 46)
Certified Values (lines 41 to 51)
Data (lines 61 to 84)
Procedure: Nonlinear Least Squares Regression
Description: These data are taken from an example discussed in
Lanczos (1956). The data were generated to 14-digits
of accuracy using
f(x) = 0.0951*exp(-x) + 0.8607*exp(-3*x)
+ 1.5576*exp(-5*x).
Reference: Lanczos, C. (1956).
Applied Analysis.
Englewood Cliffs, NJ: Prentice Hall, pp. 272-280.
Data: 1 Response (y)
1 Predictor (x)
24 Observations
Average Level of Difficulty
Generated Data
Model: Exponential Class
6 Parameters (b1 to b6)
y = b1*exp(-b2*x) + b3*exp(-b4*x) + b5*exp(-b6*x) + e
Starting values Certified Values
Start 1 Start 2 Parameter Standard Deviation
b1 = 1.2 0.5 9.5100000027E-02 5.3347304234E-11
b2 = 0.3 0.7 1.0000000001E+00 2.7473038179E-10
b3 = 5.6 3.6 8.6070000013E-01 1.3576062225E-10
b4 = 5.5 4.2 3.0000000002E+00 3.3308253069E-10
b5 = 6.5 4 1.5575999998E+00 1.8815731448E-10
b6 = 7.6 6.3 5.0000000001E+00 1.1057500538E-10
Residual Sum of Squares: 1.4307867721E-25
Residual Standard Deviation: 8.9156129349E-14
Degrees of Freedom: 18
Number of Observations: 24
Data: y x
2.513400000000E+00 0.000000000000E+00
2.044333373291E+00 5.000000000000E-02
1.668404436564E+00 1.000000000000E-01
1.366418021208E+00 1.500000000000E-01
1.123232487372E+00 2.000000000000E-01
9.268897180037E-01 2.500000000000E-01
7.679338563728E-01 3.000000000000E-01
6.388775523106E-01 3.500000000000E-01
5.337835317402E-01 4.000000000000E-01
4.479363617347E-01 4.500000000000E-01
3.775847884350E-01 5.000000000000E-01
3.197393199326E-01 5.500000000000E-01
2.720130773746E-01 6.000000000000E-01
2.324965529032E-01 6.500000000000E-01
1.996589546065E-01 7.000000000000E-01
1.722704126914E-01 7.500000000000E-01
1.493405660168E-01 8.000000000000E-01
1.300700206922E-01 8.500000000000E-01
1.138119324644E-01 9.000000000000E-01
1.000415587559E-01 9.500000000000E-01
8.833209084540E-02 1.000000000000E+00
7.833544019350E-02 1.050000000000E+00
6.976693743449E-02 1.100000000000E+00
6.239312536719E-02 1.150000000000E+00

93
src/test/resources/org/apache/commons/math3/fitting/leastsquares/MGH17.dat Normal file
View File

@ -0,0 +1,93 @@
NIST/ITL StRD
Dataset Name: MGH17 (MGH17.dat)
File Format: ASCII
Starting Values (lines 41 to 45)
Certified Values (lines 41 to 50)
Data (lines 61 to 93)
Procedure: Nonlinear Least Squares Regression
Description: This problem was found to be difficult for some very
good algorithms.
See More, J. J., Garbow, B. S., and Hillstrom, K. E.
(1981). Testing unconstrained optimization software.
ACM Transactions on Mathematical Software. 7(1):
pp. 17-41.
Reference: Osborne, M. R. (1972).
Some aspects of nonlinear least squares
calculations. In Numerical Methods for Nonlinear
Optimization, Lootsma (Ed).
New York, NY: Academic Press, pp. 171-189.
Data: 1 Response (y)
1 Predictor (x)
33 Observations
Average Level of Difficulty
Generated Data
Model: Exponential Class
5 Parameters (b1 to b5)
y = b1 + b2*exp[-x*b4] + b3*exp[-x*b5] + e
Starting values Certified Values
Start 1 Start 2 Parameter Standard Deviation
b1 = 50 0.5 3.7541005211E-01 2.0723153551E-03
b2 = 150 1.5 1.9358469127E+00 2.2031669222E-01
b3 = -100 -1 -1.4646871366E+00 2.2175707739E-01
b4 = 1 0.01 1.2867534640E-02 4.4861358114E-04
b5 = 2 0.02 2.2122699662E-02 8.9471996575E-04
Residual Sum of Squares: 5.4648946975E-05
Residual Standard Deviation: 1.3970497866E-03
Degrees of Freedom: 28
Number of Observations: 33
Data: y x
8.440000E-01 0.000000E+00
9.080000E-01 1.000000E+01
9.320000E-01 2.000000E+01
9.360000E-01 3.000000E+01
9.250000E-01 4.000000E+01
9.080000E-01 5.000000E+01
8.810000E-01 6.000000E+01
8.500000E-01 7.000000E+01
8.180000E-01 8.000000E+01
7.840000E-01 9.000000E+01
7.510000E-01 1.000000E+02
7.180000E-01 1.100000E+02
6.850000E-01 1.200000E+02
6.580000E-01 1.300000E+02
6.280000E-01 1.400000E+02
6.030000E-01 1.500000E+02
5.800000E-01 1.600000E+02
5.580000E-01 1.700000E+02
5.380000E-01 1.800000E+02
5.220000E-01 1.900000E+02
5.060000E-01 2.000000E+02
4.900000E-01 2.100000E+02
4.780000E-01 2.200000E+02
4.670000E-01 2.300000E+02
4.570000E-01 2.400000E+02
4.480000E-01 2.500000E+02
4.380000E-01 2.600000E+02
4.310000E-01 2.700000E+02
4.240000E-01 2.800000E+02
4.200000E-01 2.900000E+02
4.140000E-01 3.000000E+02
4.110000E-01 3.100000E+02
4.060000E-01 3.200000E+02