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

resurrected the estimation package but with all interfaces and classes deprecated 

this will help users switch smoothly to the new optimization.general package

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@754732 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-03-15 19:30:44 +00:00
parent c37f06ed3a
commit 61ca8c04df
15 changed files with 5356 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

314
src/java/org/apache/commons/math/estimation/AbstractEstimator.java Normal file
View File

@ -0,0 +1,314 @@
/*
* 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.math.estimation;
import java.util.Arrays;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.decomposition.LUDecompositionImpl;
/**
* Base class for implementing estimators.
* <p>This base class handles the boilerplates methods associated to thresholds
* settings, jacobian and error estimation.</p>
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public abstract class AbstractEstimator implements Estimator {
/** Default maximal number of cost evaluations allowed. */
public static final int DEFAULT_MAX_COST_EVALUATIONS = 100;
/**
* Build an abstract estimator for least squares problems.
* <p>The maximal number of cost evaluations allowed is set
* to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.</p>
*/
protected AbstractEstimator() {
setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS);
}
/**
* Set the maximal number of cost evaluations allowed.
*
* @param maxCostEval maximal number of cost evaluations allowed
* @see #estimate
*/
public final void setMaxCostEval(int maxCostEval) {
this.maxCostEval = maxCostEval;
}
/**
* Get the number of cost evaluations.
*
* @return number of cost evaluations
* */
public final int getCostEvaluations() {
return costEvaluations;
}
/**
* Get the number of jacobian evaluations.
*
* @return number of jacobian evaluations
* */
public final int getJacobianEvaluations() {
return jacobianEvaluations;
}
/**
* Update the jacobian matrix.
*/
protected void updateJacobian() {
incrementJacobianEvaluationsCounter();
Arrays.fill(jacobian, 0);
for (int i = 0, index = 0; i < rows; i++) {
WeightedMeasurement wm = measurements[i];
double factor = -Math.sqrt(wm.getWeight());
for (int j = 0; j < cols; ++j) {
jacobian[index++] = factor * wm.getPartial(parameters[j]);
}
}
}
/**
* Increment the jacobian evaluations counter.
*/
protected final void incrementJacobianEvaluationsCounter() {
++jacobianEvaluations;
}
/**
* Update the residuals array and cost function value.
* @exception EstimationException if the number of cost evaluations
* exceeds the maximum allowed
*/
protected void updateResidualsAndCost()
throws EstimationException {
if (++costEvaluations > maxCostEval) {
throw new EstimationException("maximal number of evaluations exceeded ({0})",
maxCostEval);
}
cost = 0;
for (int i = 0, index = 0; i < rows; i++, index += cols) {
WeightedMeasurement wm = measurements[i];
double residual = wm.getResidual();
residuals[i] = Math.sqrt(wm.getWeight()) * residual;
cost += wm.getWeight() * residual * residual;
}
cost = Math.sqrt(cost);
}
/**
* Get the Root Mean Square value.
* Get the Root Mean Square value, i.e. the root of the arithmetic
* mean of the square of all weighted residuals. This is related to the
* criterion that is minimized by the estimator as follows: if
* <em>c</em> if the criterion, and <em>n</em> is the number of
* measurements, then the RMS is <em>sqrt (c/n)</em>.
*
* @param problem estimation problem
* @return RMS value
*/
public double getRMS(EstimationProblem problem) {
WeightedMeasurement[] wm = problem.getMeasurements();
double criterion = 0;
for (int i = 0; i < wm.length; ++i) {
double residual = wm[i].getResidual();
criterion += wm[i].getWeight() * residual * residual;
}
return Math.sqrt(criterion / wm.length);
}
/**
* Get the Chi-Square value.
* @param problem estimation problem
* @return chi-square value
*/
public double getChiSquare(EstimationProblem problem) {
WeightedMeasurement[] wm = problem.getMeasurements();
double chiSquare = 0;
for (int i = 0; i < wm.length; ++i) {
double residual = wm[i].getResidual();
chiSquare += residual * residual / wm[i].getWeight();
}
return chiSquare;
}
/**
* Get the covariance matrix of unbound estimated parameters.
* @param problem estimation problem
* @return covariance matrix
* @exception EstimationException if the covariance matrix
* cannot be computed (singular problem)
*/
public double[][] getCovariances(EstimationProblem problem)
throws EstimationException {
// set up the jacobian
updateJacobian();
// compute transpose(J).J, avoiding building big intermediate matrices
final int rows = problem.getMeasurements().length;
final int cols = problem.getUnboundParameters().length;
final int max = cols * rows;
double[][] jTj = new double[cols][cols];
for (int i = 0; i < cols; ++i) {
for (int j = i; j < cols; ++j) {
double sum = 0;
for (int k = 0; k < max; k += cols) {
sum += jacobian[k + i] * jacobian[k + j];
}
jTj[i][j] = sum;
jTj[j][i] = sum;
}
}
try {
// compute the covariances matrix
RealMatrix inverse =
new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
return inverse.getData();
} catch (InvalidMatrixException ime) {
throw new EstimationException("unable to compute covariances: singular problem");
}
}
/**
* Guess the errors in unbound estimated parameters.
* <p>Guessing is covariance-based, it only gives rough order of magnitude.</p>
* @param problem estimation problem
* @return errors in estimated parameters
* @exception EstimationException if the covariances matrix cannot be computed
* or the number of degrees of freedom is not positive (number of measurements
* lesser or equal to number of parameters)
*/
public double[] guessParametersErrors(EstimationProblem problem)
throws EstimationException {
int m = problem.getMeasurements().length;
int p = problem.getUnboundParameters().length;
if (m <= p) {
throw new EstimationException(
"no degrees of freedom ({0} measurements, {1} parameters)",
m, p);
}
double[] errors = new double[problem.getUnboundParameters().length];
final double c = Math.sqrt(getChiSquare(problem) / (m - p));
double[][] covar = getCovariances(problem);
for (int i = 0; i < errors.length; ++i) {
errors[i] = Math.sqrt(covar[i][i]) * c;
}
return errors;
}
/**
* Initialization of the common parts of the estimation.
* <p>This method <em>must</em> be called at the start
* of the {@link #estimate(EstimationProblem) estimate}
* method.</p>
* @param problem estimation problem to solve
*/
protected void initializeEstimate(EstimationProblem problem) {
// reset counters
costEvaluations = 0;
jacobianEvaluations = 0;
// retrieve the equations and the parameters
measurements = problem.getMeasurements();
parameters = problem.getUnboundParameters();
// arrays shared with the other private methods
rows = measurements.length;
cols = parameters.length;
jacobian = new double[rows * cols];
residuals = new double[rows];
cost = Double.POSITIVE_INFINITY;
}
/**
* Solve an estimation problem.
*
* <p>The method should set the parameters of the problem to several
* trial values until it reaches convergence. If this method returns
* normally (i.e. without throwing an exception), then the best
* estimate of the parameters can be retrieved from the problem
* itself, through the {@link EstimationProblem#getAllParameters
* EstimationProblem.getAllParameters} method.</p>
*
* @param problem estimation problem to solve
* @exception EstimationException if the problem cannot be solved
*
*/
public abstract void estimate(EstimationProblem problem)
throws EstimationException;
/** Array of measurements. */
protected WeightedMeasurement[] measurements;
/** Array of parameters. */
protected EstimatedParameter[] parameters;
/**
* Jacobian matrix.
* <p>This matrix is in canonical form just after the calls to
* {@link #updateJacobian()}, but may be modified by the solver
* in the derived class (the {@link LevenbergMarquardtEstimator
* Levenberg-Marquardt estimator} does this).</p>
*/
protected double[] jacobian;
/** Number of columns of the jacobian matrix. */
protected int cols;
/** Number of rows of the jacobian matrix. */
protected int rows;
/** Residuals array.
* <p>This array is in canonical form just after the calls to
* {@link #updateJacobian()}, but may be modified by the solver
* in the derived class (the {@link LevenbergMarquardtEstimator
* Levenberg-Marquardt estimator} does this).</p>
*/
protected double[] residuals;
/** Cost value (square root of the sum of the residuals). */
protected double cost;
/** Maximal allowed number of cost evaluations. */
private int maxCostEval;
/** Number of cost evaluations. */
private int costEvaluations;
/** Number of jacobian evaluations. */
private int jacobianEvaluations;
}

126
src/java/org/apache/commons/math/estimation/EstimatedParameter.java Normal file
View File

@ -0,0 +1,126 @@
/*
* 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.math.estimation;
import java.io.Serializable;
/** This class represents the estimated parameters of an estimation problem.
*
* <p>The parameters of an estimation problem have a name, a value and
* a bound flag. The value of bound parameters is considered trusted
* and the solvers should not adjust them. On the other hand, the
* solvers should adjust the value of unbounds parameters until they
* satisfy convergence criterions specific to each solver.</p>
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public class EstimatedParameter
implements Serializable {
/** Simple constructor.
* Build an instance from a first estimate of the parameter,
* initially considered unbound.
* @param name name of the parameter
* @param firstEstimate first estimate of the parameter
*/
public EstimatedParameter(String name, double firstEstimate) {
this.name = name;
estimate = firstEstimate;
bound = false;
}
/** Simple constructor.
* Build an instance from a first estimate of the parameter and a
* bound flag
* @param name name of the parameter
* @param firstEstimate first estimate of the parameter
* @param bound flag, should be true if the parameter is bound
*/
public EstimatedParameter(String name,
double firstEstimate,
boolean bound) {
this.name = name;
estimate = firstEstimate;
this.bound = bound;
}
/** Copy constructor.
* Build a copy of a parameter
* @param parameter instance to copy
*/
public EstimatedParameter(EstimatedParameter parameter) {
name = parameter.name;
estimate = parameter.estimate;
bound = parameter.bound;
}
/** Set a new estimated value for the parameter.
* @param estimate new estimate for the parameter
*/
public void setEstimate(double estimate) {
this.estimate = estimate;
}
/** Get the current estimate of the parameter
* @return current estimate
*/
public double getEstimate() {
return estimate;
}
/** get the name of the parameter
* @return parameter name
*/
public String getName() {
return name;
}
/** Set the bound flag of the parameter
* @param bound this flag should be set to true if the parameter is
* bound (i.e. if it should not be adjusted by the solver).
*/
public void setBound(boolean bound) {
this.bound = bound;
}
/** Check if the parameter is bound
* @return true if the parameter is bound */
public boolean isBound() {
return bound;
}
/** Name of the parameter */
private String name;
/** Current value of the parameter */
protected double estimate;
/** Indicator for bound parameters
* (ie parameters that should not be estimated)
*/
private boolean bound;
/** Serializable version identifier */
private static final long serialVersionUID = -555440800213416949L;
}

48
src/java/org/apache/commons/math/estimation/EstimationException.java Normal file
View File

@ -0,0 +1,48 @@
/*
* 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.math.estimation;
import org.apache.commons.math.MathException;
/**
* This class represents exceptions thrown by the estimation solvers.
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public class EstimationException
extends MathException {
/** Serializable version identifier. */
private static final long serialVersionUID = -573038581493881337L;
/**
* Simple constructor.
* Build an exception by translating and formating a message
* @param specifier format specifier (to be translated)
* @param parts to insert in the format (no translation)
*/
public EstimationException(String specifier, Object ... parts) {
super(specifier, parts);
}
}

67
src/java/org/apache/commons/math/estimation/EstimationProblem.java Normal file
View File

@ -0,0 +1,67 @@
/*
* 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.math.estimation;
/**
* This interface represents an estimation problem.
*
* <p>This interface should be implemented by all real estimation
* problems before they can be handled by the estimators through the
* {@link Estimator#estimate Estimator.estimate} method.</p>
*
* <p>An estimation problem, as seen by a solver is a set of
* parameters and a set of measurements. The parameters are adjusted
* during the estimation through the {@link #getUnboundParameters
* getUnboundParameters} and {@link EstimatedParameter#setEstimate
* EstimatedParameter.setEstimate} methods. The measurements both have
* a measured value which is generally fixed at construction and a
* theoretical value which depends on the model and hence varies as
* the parameters are adjusted. The purpose of the solver is to reduce
* the residual between these values, it can retrieve the measurements
* through the {@link #getMeasurements getMeasurements} method.</p>
*
* @see Estimator
* @see WeightedMeasurement
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public interface EstimationProblem {
/**
* Get the measurements of an estimation problem.
* @return measurements
*/
public WeightedMeasurement[] getMeasurements();
/**
* Get the unbound parameters of the problem.
* @return unbound parameters
*/
public EstimatedParameter[] getUnboundParameters();
/**
* Get all the parameters of the problem.
* @return parameters
*/
public EstimatedParameter[] getAllParameters();
}

93
src/java/org/apache/commons/math/estimation/Estimator.java Normal file
View File

@ -0,0 +1,93 @@
/*
* 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.math.estimation;
/**
* This interface represents solvers for estimation problems.
*
* <p>The classes which are devoted to solve estimation problems
* should implement this interface. The problems which can be handled
* should implement the {@link EstimationProblem} interface which
* gather all the information needed by the solver.</p>
*
* <p>The interface is composed only of the {@link #estimate estimate}
* method.</p>
*
* @see EstimationProblem
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public interface Estimator {
/**
* Solve an estimation problem.
*
* <p>The method should set the parameters of the problem to several
* trial values until it reaches convergence. If this method returns
* normally (i.e. without throwing an exception), then the best
* estimate of the parameters can be retrieved from the problem
* itself, through the {@link EstimationProblem#getAllParameters
* EstimationProblem.getAllParameters} method.</p>
*
* @param problem estimation problem to solve
* @exception EstimationException if the problem cannot be solved
*
*/
public void estimate(EstimationProblem problem)
throws EstimationException;
/**
* Get the Root Mean Square value.
* Get the Root Mean Square value, i.e. the root of the arithmetic
* mean of the square of all weighted residuals. This is related to the
* criterion that is minimized by the estimator as follows: if
* <em>c</em> is the criterion, and <em>n</em> is the number of
* measurements, then the RMS is <em>sqrt (c/n)</em>.
* @see #guessParametersErrors(EstimationProblem)
*
* @param problem estimation problem
* @return RMS value
*/
public double getRMS(EstimationProblem problem);
/**
* Get the covariance matrix of estimated parameters.
* @param problem estimation problem
* @return covariance matrix
* @exception EstimationException if the covariance matrix
* cannot be computed (singular problem)
*/
public double[][] getCovariances(EstimationProblem problem)
throws EstimationException;
/**
* Guess the errors in estimated parameters.
* @see #getRMS(EstimationProblem)
* @param problem estimation problem
* @return errors in estimated parameters
* @exception EstimationException if the error cannot be guessed
*/
public double[] guessParametersErrors(EstimationProblem problem)
throws EstimationException;
}

228
src/java/org/apache/commons/math/estimation/GaussNewtonEstimator.java Normal file
View File

@ -0,0 +1,228 @@
/*
* 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.math.estimation;
import java.io.Serializable;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealVector;
import org.apache.commons.math.linear.RealVectorImpl;
import org.apache.commons.math.linear.decomposition.LUDecompositionImpl;
/**
* This class implements a solver for estimation problems.
*
* <p>This class solves estimation problems using a weighted least
* squares criterion on the measurement residuals. It uses a
* Gauss-Newton algorithm.</p>
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public class GaussNewtonEstimator extends AbstractEstimator implements Serializable {
/** Serializable version identifier */
private static final long serialVersionUID = 5485001826076289109L;
/** Default threshold for cost steady state detection. */
private static final double DEFAULT_STEADY_STATE_THRESHOLD = 1.0e-6;
/** Default threshold for cost convergence. */
private static final double DEFAULT_CONVERGENCE = 1.0e-6;
/** Threshold for cost steady state detection. */
private double steadyStateThreshold;
/** Threshold for cost convergence. */
private double convergence;
/** Simple constructor with default settings.
* <p>
* The estimator is built with default values for all settings.
* </p>
* @see #DEFAULT_STEADY_STATE_THRESHOLD
* @see #DEFAULT_CONVERGENCE
* @see AbstractEstimator#DEFAULT_MAX_COST_EVALUATIONS
*/
public GaussNewtonEstimator() {
this.steadyStateThreshold = DEFAULT_STEADY_STATE_THRESHOLD;
this.convergence = DEFAULT_CONVERGENCE;
}
/**
* Simple constructor.
*
* <p>This constructor builds an estimator and stores its convergence
* characteristics.</p>
*
* <p>An estimator is considered to have converged whenever either
* the criterion goes below a physical threshold under which
* improvements are considered useless or when the algorithm is
* unable to improve it (even if it is still high). The first
* condition that is met stops the iterations.</p>
*
* <p>The fact an estimator has converged does not mean that the
* model accurately fits the measurements. It only means no better
* solution can be found, it does not mean this one is good. Such an
* analysis is left to the caller.</p>
*
* <p>If neither conditions are fulfilled before a given number of
* iterations, the algorithm is considered to have failed and an
* {@link EstimationException} is thrown.</p>
*
* @param maxCostEval maximal number of cost evaluations allowed
* @param convergence criterion threshold below which we do not need
* to improve the criterion anymore
* @param steadyStateThreshold steady state detection threshold, the
* problem has converged has reached a steady state if
* <code>Math.abs(J<sub>n</sub> - J<sub>n-1</sub>) &lt;
* J<sub>n</sub> &times convergence</code>, where <code>J<sub>n</sub></code>
* and <code>J<sub>n-1</sub></code> are the current and preceding criterion
* values (square sum of the weighted residuals of considered measurements).
*/
public GaussNewtonEstimator(final int maxCostEval, final double convergence,
final double steadyStateThreshold) {
setMaxCostEval(maxCostEval);
this.steadyStateThreshold = steadyStateThreshold;
this.convergence = convergence;
}
/**
* Set the convergence criterion threshold.
* @param convergence criterion threshold below which we do not need
* to improve the criterion anymore
*/
public void setConvergence(final double convergence) {
this.convergence = convergence;
}
/**
* Set the steady state detection threshold.
* <p>
* The problem has converged has reached a steady state if
* <code>Math.abs(J<sub>n</sub> - J<sub>n-1</sub>) &lt;
* J<sub>n</sub> &times convergence</code>, where <code>J<sub>n</sub></code>
* and <code>J<sub>n-1</sub></code> are the current and preceding criterion
* values (square sum of the weighted residuals of considered measurements).
* </p>
* @param steadyStateThreshold steady state detection threshold
*/
public void setSteadyStateThreshold(final double steadyStateThreshold) {
this.steadyStateThreshold = steadyStateThreshold;
}
/**
* Solve an estimation problem using a least squares criterion.
*
* <p>This method set the unbound parameters of the given problem
* starting from their current values through several iterations. At
* each step, the unbound parameters are changed in order to
* minimize a weighted least square criterion based on the
* measurements of the problem.</p>
*
* <p>The iterations are stopped either when the criterion goes
* below a physical threshold under which improvement are considered
* useless or when the algorithm is unable to improve it (even if it
* is still high). The first condition that is met stops the
* iterations. If the convergence it not reached before the maximum
* number of iterations, an {@link EstimationException} is
* thrown.</p>
*
* @param problem estimation problem to solve
* @exception EstimationException if the problem cannot be solved
*
* @see EstimationProblem
*
*/
public void estimate(EstimationProblem problem)
throws EstimationException {
initializeEstimate(problem);
// work matrices
double[] grad = new double[parameters.length];
RealVectorImpl bDecrement = new RealVectorImpl(parameters.length);
double[] bDecrementData = bDecrement.getDataRef();
RealMatrix wGradGradT = MatrixUtils.createRealMatrix(parameters.length, parameters.length);
// iterate until convergence is reached
double previous = Double.POSITIVE_INFINITY;
do {
// build the linear problem
incrementJacobianEvaluationsCounter();
RealVector b = new RealVectorImpl(parameters.length);
RealMatrix a = MatrixUtils.createRealMatrix(parameters.length, parameters.length);
for (int i = 0; i < measurements.length; ++i) {
if (! measurements [i].isIgnored()) {
double weight = measurements[i].getWeight();
double residual = measurements[i].getResidual();
// compute the normal equation
for (int j = 0; j < parameters.length; ++j) {
grad[j] = measurements[i].getPartial(parameters[j]);
bDecrementData[j] = weight * residual * grad[j];
}
// build the contribution matrix for measurement i
for (int k = 0; k < parameters.length; ++k) {
double gk = grad[k];
for (int l = 0; l < parameters.length; ++l) {
wGradGradT.setEntry(k, l, weight * gk * grad[l]);
}
}
// update the matrices
a = a.add(wGradGradT);
b = b.add(bDecrement);
}
}
try {
// solve the linearized least squares problem
RealVector dX = new LUDecompositionImpl(a).getSolver().solve(b);
// update the estimated parameters
for (int i = 0; i < parameters.length; ++i) {
parameters[i].setEstimate(parameters[i].getEstimate() + dX.getEntry(i));
}
} catch(InvalidMatrixException e) {
throw new EstimationException("unable to solve: singular problem");
}
previous = cost;
updateResidualsAndCost();
} while ((getCostEvaluations() < 2) ||
(Math.abs(previous - cost) > (cost * steadyStateThreshold) &&
(Math.abs(cost) > convergence)));
}
}

874
src/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java Normal file
View File

@ -0,0 +1,874 @@
/*
* 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.math.estimation;
import java.io.Serializable;
import java.util.Arrays;
/**
* This class solves a least squares problem.
*
* <p>This implementation <em>should</em> work even for over-determined systems
* (i.e. systems having more variables than equations). Over-determined systems
* are solved by ignoring the variables which have the smallest impact according
* to their jacobian column norm. Only the rank of the matrix and some loop bounds
* are changed to implement this.</p>
*
* <p>The resolution engine is a simple translation of the MINPACK <a
* href="http://www.netlib.org/minpack/lmder.f">lmder</a> routine with minor
* changes. The changes include the over-determined resolution and the Q.R.
* decomposition which has been rewritten following the algorithm described in the
* P. Lascaux and R. Theodor book <i>Analyse num&eacute;rique matricielle
* appliqu&eacute;e &agrave; l'art de l'ing&eacute;nieur</i>, Masson 1986. The
* redistribution policy for MINPACK is available <a
* href="http://www.netlib.org/minpack/disclaimer">here</a>, for convenience, it
* is reproduced below.</p>
*
* <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
* <tr><td>
* Minpack Copyright Notice (1999) University of Chicago.
* All rights reserved
* </td></tr>
* <tr><td>
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* <ol>
* <li>Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.</li>
* <li>Redistributions in binary form must reproduce the above
* copyright notice, this list of conditions and the following
* disclaimer in the documentation and/or other materials provided
* with the distribution.</li>
* <li>The end-user documentation included with the redistribution, if any,
* must include the following acknowledgment:
* <code>This product includes software developed by the University of
* Chicago, as Operator of Argonne National Laboratory.</code>
* Alternately, this acknowledgment may appear in the software itself,
* if and wherever such third-party acknowledgments normally appear.</li>
* <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
* WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
* UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
* THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
* OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
* OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
* USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
* THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
* DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
* UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
* BE CORRECTED.</strong></li>
* <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
* HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
* ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
* INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
* ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
* PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
* SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
* (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
* EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
* POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
* <ol></td></tr>
* </table>
* @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran)
* @author Burton S. Garbow (original fortran)
* @author Kenneth E. Hillstrom (original fortran)
* @author Jorge J. More (original fortran)
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public class LevenbergMarquardtEstimator extends AbstractEstimator implements Serializable {
/**
* Build an estimator for least squares problems.
* <p>The default values for the algorithm settings are:
* <ul>
* <li>{@link #setInitialStepBoundFactor initial step bound factor}: 100.0</li>
* <li>{@link #setMaxCostEval maximal cost evaluations}: 1000</li>
* <li>{@link #setCostRelativeTolerance cost relative tolerance}: 1.0e-10</li>
* <li>{@link #setParRelativeTolerance parameters relative tolerance}: 1.0e-10</li>
* <li>{@link #setOrthoTolerance orthogonality tolerance}: 1.0e-10</li>
* </ul>
* </p>
*/
public LevenbergMarquardtEstimator() {
// set up the superclass with a default max cost evaluations setting
setMaxCostEval(1000);
// default values for the tuning parameters
setInitialStepBoundFactor(100.0);
setCostRelativeTolerance(1.0e-10);
setParRelativeTolerance(1.0e-10);
setOrthoTolerance(1.0e-10);
}
/**
* Set the positive input variable used in determining the initial step bound.
* This bound is set to the product of initialStepBoundFactor and the euclidean norm of diag*x if nonzero,
* or else to initialStepBoundFactor itself. In most cases factor should lie
* in the interval (0.1, 100.0). 100.0 is a generally recommended value
*
* @param initialStepBoundFactor initial step bound factor
* @see #estimate
*/
public void setInitialStepBoundFactor(double initialStepBoundFactor) {
this.initialStepBoundFactor = initialStepBoundFactor;
}
/**
* Set the desired relative error in the sum of squares.
*
* @param costRelativeTolerance desired relative error in the sum of squares
* @see #estimate
*/
public void setCostRelativeTolerance(double costRelativeTolerance) {
this.costRelativeTolerance = costRelativeTolerance;
}
/**
* Set the desired relative error in the approximate solution parameters.
*
* @param parRelativeTolerance desired relative error
* in the approximate solution parameters
* @see #estimate
*/
public void setParRelativeTolerance(double parRelativeTolerance) {
this.parRelativeTolerance = parRelativeTolerance;
}
/**
* Set the desired max cosine on the orthogonality.
*
* @param orthoTolerance desired max cosine on the orthogonality
* between the function vector and the columns of the jacobian
* @see #estimate
*/
public void setOrthoTolerance(double orthoTolerance) {
this.orthoTolerance = orthoTolerance;
}
/**
* Solve an estimation problem using the Levenberg-Marquardt algorithm.
* <p>The algorithm used is a modified Levenberg-Marquardt one, based
* on the MINPACK <a href="http://www.netlib.org/minpack/lmder.f">lmder</a>
* routine. The algorithm settings must have been set up before this method
* is called with the {@link #setInitialStepBoundFactor},
* {@link #setMaxCostEval}, {@link #setCostRelativeTolerance},
* {@link #setParRelativeTolerance} and {@link #setOrthoTolerance} methods.
* If these methods have not been called, the default values set up by the
* {@link #LevenbergMarquardtEstimator() constructor} will be used.</p>
* <p>The authors of the original fortran function are:</p>
* <ul>
* <li>Argonne National Laboratory. MINPACK project. March 1980</li>
* <li>Burton S. Garbow</li>
* <li>Kenneth E. Hillstrom</li>
* <li>Jorge J. More</li>
* </ul>
* <p>Luc Maisonobe did the Java translation.</p>
*
* @param problem estimation problem to solve
* @exception EstimationException if convergence cannot be
* reached with the specified algorithm settings or if there are more variables
* than equations
* @see #setInitialStepBoundFactor
* @see #setCostRelativeTolerance
* @see #setParRelativeTolerance
* @see #setOrthoTolerance
*/
public void estimate(EstimationProblem problem)
throws EstimationException {
initializeEstimate(problem);
// arrays shared with the other private methods
solvedCols = Math.min(rows, cols);
diagR = new double[cols];
jacNorm = new double[cols];
beta = new double[cols];
permutation = new int[cols];
lmDir = new double[cols];
// local variables
double delta = 0, xNorm = 0;
double[] diag = new double[cols];
double[] oldX = new double[cols];
double[] oldRes = new double[rows];
double[] work1 = new double[cols];
double[] work2 = new double[cols];
double[] work3 = new double[cols];
// evaluate the function at the starting point and calculate its norm
updateResidualsAndCost();
// outer loop
lmPar = 0;
boolean firstIteration = true;
while (true) {
// compute the Q.R. decomposition of the jacobian matrix
updateJacobian();
qrDecomposition();
// compute Qt.res
qTy(residuals);
// now we don't need Q anymore,
// so let jacobian contain the R matrix with its diagonal elements
for (int k = 0; k < solvedCols; ++k) {
int pk = permutation[k];
jacobian[k * cols + pk] = diagR[pk];
}
if (firstIteration) {
// scale the variables according to the norms of the columns
// of the initial jacobian
xNorm = 0;
for (int k = 0; k < cols; ++k) {
double dk = jacNorm[k];
if (dk == 0) {
dk = 1.0;
}
double xk = dk * parameters[k].getEstimate();
xNorm += xk * xk;
diag[k] = dk;
}
xNorm = Math.sqrt(xNorm);
// initialize the step bound delta
delta = (xNorm == 0) ? initialStepBoundFactor : (initialStepBoundFactor * xNorm);
}
// check orthogonality between function vector and jacobian columns
double maxCosine = 0;
if (cost != 0) {
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
double s = jacNorm[pj];
if (s != 0) {
double sum = 0;
for (int i = 0, index = pj; i <= j; ++i, index += cols) {
sum += jacobian[index] * residuals[i];
}
maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost));
}
}
}
if (maxCosine <= orthoTolerance) {
return;
}
// rescale if necessary
for (int j = 0; j < cols; ++j) {
diag[j] = Math.max(diag[j], jacNorm[j]);
}
// inner loop
for (double ratio = 0; ratio < 1.0e-4;) {
// save the state
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
oldX[pj] = parameters[pj].getEstimate();
}
double previousCost = cost;
double[] tmpVec = residuals;
residuals = oldRes;
oldRes = tmpVec;
// determine the Levenberg-Marquardt parameter
determineLMParameter(oldRes, delta, diag, work1, work2, work3);
// compute the new point and the norm of the evolution direction
double lmNorm = 0;
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
lmDir[pj] = -lmDir[pj];
parameters[pj].setEstimate(oldX[pj] + lmDir[pj]);
double s = diag[pj] * lmDir[pj];
lmNorm += s * s;
}
lmNorm = Math.sqrt(lmNorm);
// on the first iteration, adjust the initial step bound.
if (firstIteration) {
delta = Math.min(delta, lmNorm);
}
// evaluate the function at x + p and calculate its norm
updateResidualsAndCost();
// compute the scaled actual reduction
double actRed = -1.0;
if (0.1 * cost < previousCost) {
double r = cost / previousCost;
actRed = 1.0 - r * r;
}
// compute the scaled predicted reduction
// and the scaled directional derivative
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
double dirJ = lmDir[pj];
work1[j] = 0;
for (int i = 0, index = pj; i <= j; ++i, index += cols) {
work1[i] += jacobian[index] * dirJ;
}
}
double coeff1 = 0;
for (int j = 0; j < solvedCols; ++j) {
coeff1 += work1[j] * work1[j];
}
double pc2 = previousCost * previousCost;
coeff1 = coeff1 / pc2;
double coeff2 = lmPar * lmNorm * lmNorm / pc2;
double preRed = coeff1 + 2 * coeff2;
double dirDer = -(coeff1 + coeff2);
// ratio of the actual to the predicted reduction
ratio = (preRed == 0) ? 0 : (actRed / preRed);
// update the step bound
if (ratio <= 0.25) {
double tmp =
(actRed < 0) ? (0.5 * dirDer / (dirDer + 0.5 * actRed)) : 0.5;
if ((0.1 * cost >= previousCost) || (tmp < 0.1)) {
tmp = 0.1;
}
delta = tmp * Math.min(delta, 10.0 * lmNorm);
lmPar /= tmp;
} else if ((lmPar == 0) || (ratio >= 0.75)) {
delta = 2 * lmNorm;
lmPar *= 0.5;
}
// test for successful iteration.
if (ratio >= 1.0e-4) {
// successful iteration, update the norm
firstIteration = false;
xNorm = 0;
for (int k = 0; k < cols; ++k) {
double xK = diag[k] * parameters[k].getEstimate();
xNorm += xK * xK;
}
xNorm = Math.sqrt(xNorm);
} else {
// failed iteration, reset the previous values
cost = previousCost;
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
parameters[pj].setEstimate(oldX[pj]);
}
tmpVec = residuals;
residuals = oldRes;
oldRes = tmpVec;
}
// tests for convergence.
if (((Math.abs(actRed) <= costRelativeTolerance) &&
(preRed <= costRelativeTolerance) &&
(ratio <= 2.0)) ||
(delta <= parRelativeTolerance * xNorm)) {
return;
}
// tests for termination and stringent tolerances
// (2.2204e-16 is the machine epsilon for IEEE754)
if ((Math.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) {
throw new EstimationException("cost relative tolerance is too small ({0})," +
" no further reduction in the" +
" sum of squares is possible",
costRelativeTolerance);
} else if (delta <= 2.2204e-16 * xNorm) {
throw new EstimationException("parameters relative tolerance is too small" +
" ({0}), no further improvement in" +
" the approximate solution is possible",
parRelativeTolerance);
} else if (maxCosine <= 2.2204e-16) {
throw new EstimationException("orthogonality tolerance is too small ({0})," +
" solution is orthogonal to the jacobian",
orthoTolerance);
}
}
}
}
/**
* Determine the Levenberg-Marquardt parameter.
* <p>This implementation is a translation in Java of the MINPACK
* <a href="http://www.netlib.org/minpack/lmpar.f">lmpar</a>
* routine.</p>
* <p>This method sets the lmPar and lmDir attributes.</p>
* <p>The authors of the original fortran function are:</p>
* <ul>
* <li>Argonne National Laboratory. MINPACK project. March 1980</li>
* <li>Burton S. Garbow</li>
* <li>Kenneth E. Hillstrom</li>
* <li>Jorge J. More</li>
* </ul>
* <p>Luc Maisonobe did the Java translation.</p>
*
* @param qy array containing qTy
* @param delta upper bound on the euclidean norm of diagR * lmDir
* @param diag diagonal matrix
* @param work1 work array
* @param work2 work array
* @param work3 work array
*/
private void determineLMParameter(double[] qy, double delta, double[] diag,
double[] work1, double[] work2, double[] work3) {
// compute and store in x the gauss-newton direction, if the
// jacobian is rank-deficient, obtain a least squares solution
for (int j = 0; j < rank; ++j) {
lmDir[permutation[j]] = qy[j];
}
for (int j = rank; j < cols; ++j) {
lmDir[permutation[j]] = 0;
}
for (int k = rank - 1; k >= 0; --k) {
int pk = permutation[k];
double ypk = lmDir[pk] / diagR[pk];
for (int i = 0, index = pk; i < k; ++i, index += cols) {
lmDir[permutation[i]] -= ypk * jacobian[index];
}
lmDir[pk] = ypk;
}
// evaluate the function at the origin, and test
// for acceptance of the Gauss-Newton direction
double dxNorm = 0;
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
double s = diag[pj] * lmDir[pj];
work1[pj] = s;
dxNorm += s * s;
}
dxNorm = Math.sqrt(dxNorm);
double fp = dxNorm - delta;
if (fp <= 0.1 * delta) {
lmPar = 0;
return;
}
// if the jacobian is not rank deficient, the Newton step provides
// a lower bound, parl, for the zero of the function,
// otherwise set this bound to zero
double sum2, parl = 0;
if (rank == solvedCols) {
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
work1[pj] *= diag[pj] / dxNorm;
}
sum2 = 0;
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
double sum = 0;
for (int i = 0, index = pj; i < j; ++i, index += cols) {
sum += jacobian[index] * work1[permutation[i]];
}
double s = (work1[pj] - sum) / diagR[pj];
work1[pj] = s;
sum2 += s * s;
}
parl = fp / (delta * sum2);
}
// calculate an upper bound, paru, for the zero of the function
sum2 = 0;
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
double sum = 0;
for (int i = 0, index = pj; i <= j; ++i, index += cols) {
sum += jacobian[index] * qy[i];
}
sum /= diag[pj];
sum2 += sum * sum;
}
double gNorm = Math.sqrt(sum2);
double paru = gNorm / delta;
if (paru == 0) {
// 2.2251e-308 is the smallest positive real for IEE754
paru = 2.2251e-308 / Math.min(delta, 0.1);
}
// if the input par lies outside of the interval (parl,paru),
// set par to the closer endpoint
lmPar = Math.min(paru, Math.max(lmPar, parl));
if (lmPar == 0) {
lmPar = gNorm / dxNorm;
}
for (int countdown = 10; countdown >= 0; --countdown) {
// evaluate the function at the current value of lmPar
if (lmPar == 0) {
lmPar = Math.max(2.2251e-308, 0.001 * paru);
}
double sPar = Math.sqrt(lmPar);
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
work1[pj] = sPar * diag[pj];
}
determineLMDirection(qy, work1, work2, work3);
dxNorm = 0;
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
double s = diag[pj] * lmDir[pj];
work3[pj] = s;
dxNorm += s * s;
}
dxNorm = Math.sqrt(dxNorm);
double previousFP = fp;
fp = dxNorm - delta;
// if the function is small enough, accept the current value
// of lmPar, also test for the exceptional cases where parl is zero
if ((Math.abs(fp) <= 0.1 * delta) ||
((parl == 0) && (fp <= previousFP) && (previousFP < 0))) {
return;
}
// compute the Newton correction
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
work1[pj] = work3[pj] * diag[pj] / dxNorm;
}
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
work1[pj] /= work2[j];
double tmp = work1[pj];
for (int i = j + 1; i < solvedCols; ++i) {
work1[permutation[i]] -= jacobian[i * cols + pj] * tmp;
}
}
sum2 = 0;
for (int j = 0; j < solvedCols; ++j) {
double s = work1[permutation[j]];
sum2 += s * s;
}
double correction = fp / (delta * sum2);
// depending on the sign of the function, update parl or paru.
if (fp > 0) {
parl = Math.max(parl, lmPar);
} else if (fp < 0) {
paru = Math.min(paru, lmPar);
}
// compute an improved estimate for lmPar
lmPar = Math.max(parl, lmPar + correction);
}
}
/**
* Solve a*x = b and d*x = 0 in the least squares sense.
* <p>This implementation is a translation in Java of the MINPACK
* <a href="http://www.netlib.org/minpack/qrsolv.f">qrsolv</a>
* routine.</p>
* <p>This method sets the lmDir and lmDiag attributes.</p>
* <p>The authors of the original fortran function are:</p>
* <ul>
* <li>Argonne National Laboratory. MINPACK project. March 1980</li>
* <li>Burton S. Garbow</li>
* <li>Kenneth E. Hillstrom</li>
* <li>Jorge J. More</li>
* </ul>
* <p>Luc Maisonobe did the Java translation.</p>
*
* @param qy array containing qTy
* @param diag diagonal matrix
* @param lmDiag diagonal elements associated with lmDir
* @param work work array
*/
private void determineLMDirection(double[] qy, double[] diag,
double[] lmDiag, double[] work) {
// copy R and Qty to preserve input and initialize s
// in particular, save the diagonal elements of R in lmDir
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
for (int i = j + 1; i < solvedCols; ++i) {
jacobian[i * cols + pj] = jacobian[j * cols + permutation[i]];
}
lmDir[j] = diagR[pj];
work[j] = qy[j];
}
// eliminate the diagonal matrix d using a Givens rotation
for (int j = 0; j < solvedCols; ++j) {
// prepare the row of d to be eliminated, locating the
// diagonal element using p from the Q.R. factorization
int pj = permutation[j];
double dpj = diag[pj];
if (dpj != 0) {
Arrays.fill(lmDiag, j + 1, lmDiag.length, 0);
}
lmDiag[j] = dpj;
// the transformations to eliminate the row of d
// modify only a single element of Qty
// beyond the first n, which is initially zero.
double qtbpj = 0;
for (int k = j; k < solvedCols; ++k) {
int pk = permutation[k];
// determine a Givens rotation which eliminates the
// appropriate element in the current row of d
if (lmDiag[k] != 0) {
double sin, cos;
double rkk = jacobian[k * cols + pk];
if (Math.abs(rkk) < Math.abs(lmDiag[k])) {
double cotan = rkk / lmDiag[k];
sin = 1.0 / Math.sqrt(1.0 + cotan * cotan);
cos = sin * cotan;
} else {
double tan = lmDiag[k] / rkk;
cos = 1.0 / Math.sqrt(1.0 + tan * tan);
sin = cos * tan;
}
// compute the modified diagonal element of R and
// the modified element of (Qty,0)
jacobian[k * cols + pk] = cos * rkk + sin * lmDiag[k];
double temp = cos * work[k] + sin * qtbpj;
qtbpj = -sin * work[k] + cos * qtbpj;
work[k] = temp;
// accumulate the tranformation in the row of s
for (int i = k + 1; i < solvedCols; ++i) {
double rik = jacobian[i * cols + pk];
temp = cos * rik + sin * lmDiag[i];
lmDiag[i] = -sin * rik + cos * lmDiag[i];
jacobian[i * cols + pk] = temp;
}
}
}
// store the diagonal element of s and restore
// the corresponding diagonal element of R
int index = j * cols + permutation[j];
lmDiag[j] = jacobian[index];
jacobian[index] = lmDir[j];
}
// solve the triangular system for z, if the system is
// singular, then obtain a least squares solution
int nSing = solvedCols;
for (int j = 0; j < solvedCols; ++j) {
if ((lmDiag[j] == 0) && (nSing == solvedCols)) {
nSing = j;
}
if (nSing < solvedCols) {
work[j] = 0;
}
}
if (nSing > 0) {
for (int j = nSing - 1; j >= 0; --j) {
int pj = permutation[j];
double sum = 0;
for (int i = j + 1; i < nSing; ++i) {
sum += jacobian[i * cols + pj] * work[i];
}
work[j] = (work[j] - sum) / lmDiag[j];
}
}
// permute the components of z back to components of lmDir
for (int j = 0; j < lmDir.length; ++j) {
lmDir[permutation[j]] = work[j];
}
}
/**
* Decompose a matrix A as A.P = Q.R using Householder transforms.
* <p>As suggested in the P. Lascaux and R. Theodor book
* <i>Analyse num&eacute;rique matricielle appliqu&eacute;e &agrave;
* l'art de l'ing&eacute;nieur</i> (Masson, 1986), instead of representing
* the Householder transforms with u<sub>k</sub> unit vectors such that:
* <pre>
* H<sub>k</sub> = I - 2u<sub>k</sub>.u<sub>k</sub><sup>t</sup>
* </pre>
* we use <sub>k</sub> non-unit vectors such that:
* <pre>
* H<sub>k</sub> = I - beta<sub>k</sub>v<sub>k</sub>.v<sub>k</sub><sup>t</sup>
* </pre>
* where v<sub>k</sub> = a<sub>k</sub> - alpha<sub>k</sub> e<sub>k</sub>.
* The beta<sub>k</sub> coefficients are provided upon exit as recomputing
* them from the v<sub>k</sub> vectors would be costly.</p>
* <p>This decomposition handles rank deficient cases since the tranformations
* are performed in non-increasing columns norms order thanks to columns
* pivoting. The diagonal elements of the R matrix are therefore also in
* non-increasing absolute values order.</p>
* @exception EstimationException if the decomposition cannot be performed
*/
private void qrDecomposition() throws EstimationException {
// initializations
for (int k = 0; k < cols; ++k) {
permutation[k] = k;
double norm2 = 0;
for (int index = k; index < jacobian.length; index += cols) {
double akk = jacobian[index];
norm2 += akk * akk;
}
jacNorm[k] = Math.sqrt(norm2);
}
// transform the matrix column after column
for (int k = 0; k < cols; ++k) {
// select the column with the greatest norm on active components
int nextColumn = -1;
double ak2 = Double.NEGATIVE_INFINITY;
for (int i = k; i < cols; ++i) {
double norm2 = 0;
int iDiag = k * cols + permutation[i];
for (int index = iDiag; index < jacobian.length; index += cols) {
double aki = jacobian[index];
norm2 += aki * aki;
}
if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
throw new EstimationException(
"unable to perform Q.R decomposition on the {0}x{1} jacobian matrix",
rows, cols);
}
if (norm2 > ak2) {
nextColumn = i;
ak2 = norm2;
}
}
if (ak2 == 0) {
rank = k;
return;
}
int pk = permutation[nextColumn];
permutation[nextColumn] = permutation[k];
permutation[k] = pk;
// choose alpha such that Hk.u = alpha ek
int kDiag = k * cols + pk;
double akk = jacobian[kDiag];
double alpha = (akk > 0) ? -Math.sqrt(ak2) : Math.sqrt(ak2);
double betak = 1.0 / (ak2 - akk * alpha);
beta[pk] = betak;
// transform the current column
diagR[pk] = alpha;
jacobian[kDiag] -= alpha;
// transform the remaining columns
for (int dk = cols - 1 - k; dk > 0; --dk) {
int dkp = permutation[k + dk] - pk;
double gamma = 0;
for (int index = kDiag; index < jacobian.length; index += cols) {
gamma += jacobian[index] * jacobian[index + dkp];
}
gamma *= betak;
for (int index = kDiag; index < jacobian.length; index += cols) {
jacobian[index + dkp] -= gamma * jacobian[index];
}
}
}
rank = solvedCols;
}
/**
* Compute the product Qt.y for some Q.R. decomposition.
*
* @param y vector to multiply (will be overwritten with the result)
*/
private void qTy(double[] y) {
for (int k = 0; k < cols; ++k) {
int pk = permutation[k];
int kDiag = k * cols + pk;
double gamma = 0;
for (int i = k, index = kDiag; i < rows; ++i, index += cols) {
gamma += jacobian[index] * y[i];
}
gamma *= beta[pk];
for (int i = k, index = kDiag; i < rows; ++i, index += cols) {
y[i] -= gamma * jacobian[index];
}
}
}
/** Number of solved variables. */
private int solvedCols;
/** Diagonal elements of the R matrix in the Q.R. decomposition. */
private double[] diagR;
/** Norms of the columns of the jacobian matrix. */
private double[] jacNorm;
/** Coefficients of the Householder transforms vectors. */
private double[] beta;
/** Columns permutation array. */
private int[] permutation;
/** Rank of the jacobian matrix. */
private int rank;
/** Levenberg-Marquardt parameter. */
private double lmPar;
/** Parameters evolution direction associated with lmPar. */
private double[] lmDir;
/** Positive input variable used in determining the initial step bound. */
private double initialStepBoundFactor;
/** Desired relative error in the sum of squares. */
private double costRelativeTolerance;
/** Desired relative error in the approximate solution parameters. */
private double parRelativeTolerance;
/** Desired max cosine on the orthogonality between the function vector
* and the columns of the jacobian. */
private double orthoTolerance;
/** Serializable version identifier */
private static final long serialVersionUID = -5705952631533171019L;
}

111
src/java/org/apache/commons/math/estimation/SimpleEstimationProblem.java Normal file
View File

@ -0,0 +1,111 @@
/*
* 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.math.estimation;
import java.util.ArrayList;
import java.util.List;
/**
* Simple implementation of the {@link EstimationProblem
* EstimationProblem} interface for boilerplate data handling.
* <p>This class <em>only</em> handles parameters and measurements
* storage and unbound parameters filtering. It does not compute
* anything by itself. It should either be used with measurements
* implementation that are smart enough to know about the
* various parameters in order to compute the partial derivatives
* appropriately. Since the problem-specific logic is mainly related to
* the various measurements models, the simplest way to use this class
* is by extending it and using one internal class extending
* {@link WeightedMeasurement WeightedMeasurement} for each measurement
* type. The instances of the internal classes would have access to the
* various parameters and their current estimate.</p>
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*/
@Deprecated
public class SimpleEstimationProblem implements EstimationProblem {
/**
* Build an empty instance without parameters nor measurements.
*/
public SimpleEstimationProblem() {
parameters = new ArrayList<EstimatedParameter>();
measurements = new ArrayList<WeightedMeasurement>();
}
/**
* Get all the parameters of the problem.
* @return parameters
*/
public EstimatedParameter[] getAllParameters() {
return (EstimatedParameter[]) parameters.toArray(new EstimatedParameter[parameters.size()]);
}
/**
* Get the unbound parameters of the problem.
* @return unbound parameters
*/
public EstimatedParameter[] getUnboundParameters() {
// filter the unbound parameters
List<EstimatedParameter> unbound = new ArrayList<EstimatedParameter>(parameters.size());
for (EstimatedParameter p : parameters) {
if (! p.isBound()) {
unbound.add(p);
}
}
// convert to an array
return (EstimatedParameter[]) unbound.toArray(new EstimatedParameter[unbound.size()]);
}
/**
* Get the measurements of an estimation problem.
* @return measurements
*/
public WeightedMeasurement[] getMeasurements() {
return (WeightedMeasurement[]) measurements.toArray(new WeightedMeasurement[measurements.size()]);
}
/** Add a parameter to the problem.
* @param p parameter to add
*/
protected void addParameter(EstimatedParameter p) {
parameters.add(p);
}
/**
* Add a new measurement to the set.
* @param m measurement to add
*/
protected void addMeasurement(WeightedMeasurement m) {
measurements.add(m);
}
/** Estimated parameters. */
private final List<EstimatedParameter> parameters;
/** Measurements. */
private final List<WeightedMeasurement> measurements;
}

172
src/java/org/apache/commons/math/estimation/WeightedMeasurement.java Normal file
View File

@ -0,0 +1,172 @@
/*
* 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.math.estimation;
import java.io.Serializable;
/**
* This class represents measurements in estimation problems.
*
* <p>This abstract class implements all the methods needed to handle
* measurements in a general way. It defines neither the {@link
* #getTheoreticalValue getTheoreticalValue} nor the {@link
* #getPartial getPartial} methods, which should be defined by
* sub-classes according to the specific problem.</p>
*
* <p>The {@link #getTheoreticalValue getTheoreticalValue} and {@link
* #getPartial getPartial} methods must always use the current
* estimate of the parameters set by the solver in the problem. These
* parameters can be retrieved through the {@link
* EstimationProblem#getAllParameters
* EstimationProblem.getAllParameters} method if the measurements are
* independent of the problem, or directly if they are implemented as
* inner classes of the problem.</p>
*
* <p>The instances for which the <code>ignored</code> flag is set
* through the {@link #setIgnored setIgnored} method are ignored by the
* solvers. This can be used to reject wrong measurements at some
* steps of the estimation.</p>
*
* @see EstimationProblem
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*/
@Deprecated
public abstract class WeightedMeasurement implements Serializable {
/** Serializable version identifier. */
private static final long serialVersionUID = 4360046376796901941L;
/**
* Simple constructor.
* Build a measurement with the given parameters, and set its ignore
* flag to false.
* @param weight weight of the measurement in the least squares problem
* (two common choices are either to use 1.0 for all measurements, or to
* use a value proportional to the inverse of the variance of the measurement
* type)
*
* @param measuredValue measured value
*/
public WeightedMeasurement(double weight, double measuredValue) {
this.weight = weight;
this.measuredValue = measuredValue;
ignored = false;
}
/** Simple constructor.
*
* Build a measurement with the given parameters
*
* @param weight weight of the measurement in the least squares problem
* @param measuredValue measured value
* @param ignored true if the measurement should be ignored
*/
public WeightedMeasurement(double weight, double measuredValue,
boolean ignored) {
this.weight = weight;
this.measuredValue = measuredValue;
this.ignored = ignored;
}
/**
* Get the weight of the measurement in the least squares problem
*
* @return weight
*/
public double getWeight() {
return weight;
}
/**
* Get the measured value
*
* @return measured value
*/
public double getMeasuredValue() {
return measuredValue;
}
/**
* Get the residual for this measurement
* The residual is the measured value minus the theoretical value.
*
* @return residual
*/
public double getResidual() {
return measuredValue - getTheoreticalValue();
}
/**
* Get the theoretical value expected for this measurement
* <p>The theoretical value is the value expected for this measurement
* if the model and its parameter were all perfectly known.</p>
* <p>The value must be computed using the current estimate of the parameters
* set by the solver in the problem.</p>
*
* @return theoretical value
*/
public abstract double getTheoreticalValue();
/**
* Get the partial derivative of the {@link #getTheoreticalValue
* theoretical value} according to the parameter.
* <p>The value must be computed using the current estimate of the parameters
* set by the solver in the problem.</p>
*
* @param parameter parameter against which the partial derivative
* should be computed
* @return partial derivative of the {@link #getTheoreticalValue
* theoretical value}
*/
public abstract double getPartial(EstimatedParameter parameter);
/**
* Set the ignore flag to the specified value
* Setting the ignore flag to true allow to reject wrong
* measurements, which sometimes can be detected only rather late.
*
* @param ignored value for the ignore flag
*/
public void setIgnored(boolean ignored) {
this.ignored = ignored;
}
/**
* Check if this measurement should be ignored
*
* @return true if the measurement should be ignored
*/
public boolean isIgnored() {
return ignored;
}
/** Measurement weight. */
private final double weight;
/** Value of the measurements. */
private final double measuredValue;
/** Ignore measurement indicator. */
private boolean ignored;
}

25
src/java/org/apache/commons/math/estimation/package.html Normal file
View File

@ -0,0 +1,25 @@
<html>
<!--
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.
-->
<!-- $Revision$ -->
<body>
This package provided classes to solve estimation problems, it is deprecated since 2.0.
<p>This package has been deprecated as of 2.0. It is replaced by the optimization.general package.</p>
</body>
</html>

77
src/test/org/apache/commons/math/estimation/EstimatedParameterTest.java Normal file
View File

@ -0,0 +1,77 @@
/*
* 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.math.estimation;
import org.apache.commons.math.estimation.EstimatedParameter;
import junit.framework.*;
@Deprecated
public class EstimatedParameterTest
extends TestCase {
public EstimatedParameterTest(String name) {
super(name);
}
public void testConstruction() {
EstimatedParameter p1 = new EstimatedParameter("p1", 1.0);
assertTrue(p1.getName().equals("p1"));
checkValue(p1.getEstimate(), 1.0);
assertTrue(! p1.isBound());
EstimatedParameter p2 = new EstimatedParameter("p2", 2.0, true);
assertTrue(p2.getName().equals("p2"));
checkValue(p2.getEstimate(), 2.0);
assertTrue(p2.isBound());
}
public void testBound() {
EstimatedParameter p = new EstimatedParameter("p", 0.0);
assertTrue(! p.isBound());
p.setBound(true);
assertTrue(p.isBound());
p.setBound(false);
assertTrue(! p.isBound());
}
public void testEstimate() {
EstimatedParameter p = new EstimatedParameter("p", 0.0);
checkValue(p.getEstimate(), 0.0);
for (double e = 0.0; e < 10.0; e += 0.5) {
p.setEstimate(e);
checkValue(p.getEstimate(), e);
}
}
public static Test suite() {
return new TestSuite(EstimatedParameterTest.class);
}
private void checkValue(double value, double expected) {
assertTrue(Math.abs(value - expected) < 1.0e-10);
}
}

732
src/test/org/apache/commons/math/estimation/GaussNewtonEstimatorTest.java Normal file
View File

@ -0,0 +1,732 @@
/*
* 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.math.estimation;
import java.util.ArrayList;
import java.util.HashSet;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
/**
* <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>, for
* convenience, it is reproduced below.</p>
* <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
* <tr><td>
* Minpack Copyright Notice (1999) University of Chicago.
* All rights reserved
* </td></tr>
* <tr><td>
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* <ol>
* <li>Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.</li>
* <li>Redistributions in binary form must reproduce the above
* copyright notice, this list of conditions and the following
* disclaimer in the documentation and/or other materials provided
* with the distribution.</li>
* <li>The end-user documentation included with the redistribution, if any,
* must include the following acknowledgment:
* <code>This product includes software developed by the University of
* Chicago, as Operator of Argonne National Laboratory.</code>
* Alternately, this acknowledgment may appear in the software itself,
* if and wherever such third-party acknowledgments normally appear.</li>
* <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
* WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
* UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
* THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
* OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
* OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
* USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
* THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
* DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
* UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
* BE CORRECTED.</strong></li>
* <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
* HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
* ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
* INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
* ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
* PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
* SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
* (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
* EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
* POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
* <ol></td></tr>
* </table>
* @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
* @author Burton S. Garbow (original fortran minpack tests)
* @author Kenneth E. Hillstrom (original fortran minpack tests)
* @author Jorge J. More (original fortran minpack tests)
* @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
*/
@Deprecated
public class GaussNewtonEstimatorTest
extends TestCase {
public GaussNewtonEstimatorTest(String name) {
super(name);
}
public void testTrivial() throws EstimationException {
LinearProblem problem =
new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] {2},
new EstimatedParameter[] {
new EstimatedParameter("p0", 0)
}, 3.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(1.5,
problem.getUnboundParameters()[0].getEstimate(),
1.0e-10);
}
public void testQRColumnsPermutation() throws EstimationException {
EstimatedParameter[] x = {
new EstimatedParameter("p0", 0), new EstimatedParameter("p1", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { x[0], x[1] },
4.0),
new LinearMeasurement(new double[] { 2.0 },
new EstimatedParameter[] { x[1] },
6.0),
new LinearMeasurement(new double[] { 1.0, -2.0 },
new EstimatedParameter[] { x[0], x[1] },
1.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(7.0, x[0].getEstimate(), 1.0e-10);
assertEquals(3.0, x[1].getEstimate(), 1.0e-10);
}
public void testNoDependency() throws EstimationException {
EstimatedParameter[] p = new EstimatedParameter[] {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0),
new EstimatedParameter("p3", 0),
new EstimatedParameter("p4", 0),
new EstimatedParameter("p5", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[0] }, 0.0),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[1] }, 1.1),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[2] }, 2.2),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[3] }, 3.3),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[4] }, 4.4),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[5] }, 5.5)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
for (int i = 0; i < p.length; ++i) {
assertEquals(0.55 * i, p[i].getEstimate(), 1.0e-10);
}
}
public void testOneSet() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0 },
new EstimatedParameter[] { p[0] },
1.0),
new LinearMeasurement(new double[] { -1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { -1.0, 1.0 },
new EstimatedParameter[] { p[1], p[2] },
1.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
assertEquals(2.0, p[1].getEstimate(), 1.0e-10);
assertEquals(3.0, p[2].getEstimate(), 1.0e-10);
}
public void testTwoSets() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 1),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 3),
new EstimatedParameter("p4", 4),
new EstimatedParameter("p5", 5)
};
double epsilon = 1.0e-7;
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
// 4 elements sub-problem
new LinearMeasurement(new double[] { 2.0, 1.0, 4.0 },
new EstimatedParameter[] { p[0], p[1], p[3] },
2.0),
new LinearMeasurement(new double[] { -4.0, -2.0, 3.0, -7.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-9.0),
new LinearMeasurement(new double[] { 4.0, 1.0, -2.0, 8.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
2.0),
new LinearMeasurement(new double[] { -3.0, -12.0, -1.0 },
new EstimatedParameter[] { p[1], p[2], p[3] },
2.0),
// 2 elements sub-problem
new LinearMeasurement(new double[] { epsilon, 1.0 },
new EstimatedParameter[] { p[4], p[5] },
1.0 + epsilon * epsilon),
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[4], p[5] },
2.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals( 3.0, p[0].getEstimate(), 1.0e-10);
assertEquals( 4.0, p[1].getEstimate(), 1.0e-10);
assertEquals(-1.0, p[2].getEstimate(), 1.0e-10);
assertEquals(-2.0, p[3].getEstimate(), 1.0e-10);
assertEquals( 1.0 + epsilon, p[4].getEstimate(), 1.0e-10);
assertEquals( 1.0 - epsilon, p[5].getEstimate(), 1.0e-10);
}
public void testNonInversible() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0)
};
LinearMeasurement[] m = new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 2.0, -3.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
1.0),
new LinearMeasurement(new double[] { 2.0, 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
1.0),
new LinearMeasurement(new double[] { -3.0, -9.0 },
new EstimatedParameter[] { p[0], p[2] },
1.0)
};
LinearProblem problem = new LinearProblem(m);
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
try {
estimator.estimate(problem);
fail("an exception should have been caught");
} catch (EstimationException ee) {
// expected behavior
} catch (Exception e) {
fail("wrong exception type caught");
}
}
public void testIllConditioned() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 1),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 3)
};
LinearProblem problem1 = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 10.0, 7.0, 8.0, 7.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
32.0),
new LinearMeasurement(new double[] { 7.0, 5.0, 6.0, 5.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
23.0),
new LinearMeasurement(new double[] { 8.0, 6.0, 10.0, 9.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
33.0),
new LinearMeasurement(new double[] { 7.0, 5.0, 9.0, 10.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
31.0)
});
GaussNewtonEstimator estimator1 = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator1.estimate(problem1);
assertEquals(0, estimator1.getRMS(problem1), 1.0e-10);
assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
assertEquals(1.0, p[1].getEstimate(), 1.0e-10);
assertEquals(1.0, p[2].getEstimate(), 1.0e-10);
assertEquals(1.0, p[3].getEstimate(), 1.0e-10);
LinearProblem problem2 = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 10.0, 7.0, 8.1, 7.2 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
32.0),
new LinearMeasurement(new double[] { 7.08, 5.04, 6.0, 5.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
23.0),
new LinearMeasurement(new double[] { 8.0, 5.98, 9.89, 9.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
33.0),
new LinearMeasurement(new double[] { 6.99, 4.99, 9.0, 9.98 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
31.0)
});
GaussNewtonEstimator estimator2 = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator2.estimate(problem2);
assertEquals(0, estimator2.getRMS(problem2), 1.0e-10);
assertEquals(-81.0, p[0].getEstimate(), 1.0e-8);
assertEquals(137.0, p[1].getEstimate(), 1.0e-8);
assertEquals(-34.0, p[2].getEstimate(), 1.0e-8);
assertEquals( 22.0, p[3].getEstimate(), 1.0e-8);
}
public void testMoreEstimatedParametersSimple() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 7),
new EstimatedParameter("p1", 6),
new EstimatedParameter("p2", 5),
new EstimatedParameter("p3", 4)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 3.0, 2.0 },
new EstimatedParameter[] { p[0], p[1] },
7.0),
new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
new EstimatedParameter[] { p[1], p[2], p[3] },
3.0),
new LinearMeasurement(new double[] { 2.0, 1.0 },
new EstimatedParameter[] { p[0], p[2] },
5.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
try {
estimator.estimate(problem);
fail("an exception should have been caught");
} catch (EstimationException ee) {
// expected behavior
} catch (Exception e) {
fail("wrong exception type caught");
}
}
public void testMoreEstimatedParametersUnsorted() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 2),
new EstimatedParameter("p1", 2),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 2),
new EstimatedParameter("p4", 2),
new EstimatedParameter("p5", 2)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 },
new EstimatedParameter[] { p[2], p[3], p[4] },
12.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[4], p[5] },
-1.0),
new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
new EstimatedParameter[] { p[3], p[2], p[5] },
7.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[4], p[3] },
1.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
try {
estimator.estimate(problem);
fail("an exception should have been caught");
} catch (EstimationException ee) {
// expected behavior
} catch (Exception e) {
fail("wrong exception type caught");
}
}
public void testRedundantEquations() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 1),
new EstimatedParameter("p1", 1)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1] },
5.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
EstimatedParameter[] all = problem.getAllParameters();
for (int i = 0; i < all.length; ++i) {
assertEquals(all[i].getName().equals("p0") ? 2.0 : 1.0,
all[i].getEstimate(), 1.0e-10);
}
}
public void testInconsistentEquations() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 1),
new EstimatedParameter("p1", 1)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1] },
4.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertTrue(estimator.getRMS(problem) > 0.1);
}
public void testBoundParameters() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("unbound0", 2, false),
new EstimatedParameter("unbound1", 2, false),
new EstimatedParameter("bound", 2, true)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
3.0),
new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
1.0),
new LinearMeasurement(new double[] { 1.0, 3.0, 2.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
7.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertTrue(estimator.getRMS(problem) < 1.0e-10);
double[][] covariances = estimator.getCovariances(problem);
int i0 = 0, i1 = 1;
if (problem.getUnboundParameters()[0].getName().endsWith("1")) {
i0 = 1;
i1 = 0;
}
assertEquals(11.0 / 24, covariances[i0][i0], 1.0e-10);
assertEquals(-3.0 / 24, covariances[i0][i1], 1.0e-10);
assertEquals(-3.0 / 24, covariances[i1][i0], 1.0e-10);
assertEquals( 3.0 / 24, covariances[i1][i1], 1.0e-10);
double[] errors = estimator.guessParametersErrors(problem);
assertEquals(0, errors[i0], 1.0e-10);
assertEquals(0, errors[i1], 1.0e-10);
}
public void testMaxIterations() {
Circle circle = new Circle(98.680, 47.345);
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);
try {
GaussNewtonEstimator estimator = new GaussNewtonEstimator(4, 1.0e-14, 1.0e-14);
estimator.estimate(circle);
fail("an exception should have been caught");
} catch (EstimationException ee) {
// expected behavior
} catch (Exception e) {
fail("wrong exception type caught");
}
}
public void testCircleFitting() throws EstimationException {
Circle circle = new Circle(98.680, 47.345);
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);
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-10, 1.0e-10);
estimator.estimate(circle);
double rms = estimator.getRMS(circle);
assertEquals(1.768262623567235, Math.sqrt(circle.getM()) * rms, 1.0e-10);
assertEquals(69.96016176931406, circle.getRadius(), 1.0e-10);
assertEquals(96.07590211815305, circle.getX(), 1.0e-10);
assertEquals(48.13516790438953, circle.getY(), 1.0e-10);
}
public void testCircleFittingBadInit() throws EstimationException {
Circle circle = new Circle(-12, -12);
double[][] points = 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}
};
for (int i = 0; i < points.length; ++i) {
circle.addPoint(points[i][0], points[i][1]);
}
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
try {
estimator.estimate(circle);
fail("an exception should have been caught");
} catch (EstimationException ee) {
// expected behavior
} catch (Exception e) {
fail("wrong exception type caught");
}
}
private static class LinearProblem extends SimpleEstimationProblem {
public LinearProblem(LinearMeasurement[] measurements) {
HashSet<EstimatedParameter> set = new HashSet<EstimatedParameter>();
for (int i = 0; i < measurements.length; ++i) {
addMeasurement(measurements[i]);
EstimatedParameter[] parameters = measurements[i].getParameters();
for (int j = 0; j < parameters.length; ++j) {
set.add(parameters[j]);
}
}
for (EstimatedParameter p : set) {
addParameter(p);
}
}
}
private static class LinearMeasurement extends WeightedMeasurement {
public LinearMeasurement(double[] factors, EstimatedParameter[] parameters,
double setPoint) {
super(1.0, setPoint, true);
this.factors = factors;
this.parameters = parameters;
setIgnored(false);
}
public double getTheoreticalValue() {
double v = 0;
for (int i = 0; i < factors.length; ++i) {
v += factors[i] * parameters[i].getEstimate();
}
return v;
}
public double getPartial(EstimatedParameter parameter) {
for (int i = 0; i < parameters.length; ++i) {
if (parameters[i] == parameter) {
return factors[i];
}
}
return 0;
}
public EstimatedParameter[] getParameters() {
return parameters;
}
private double[] factors;
private EstimatedParameter[] parameters;
private static final long serialVersionUID = -3922448707008868580L;
}
private static class Circle implements EstimationProblem {
public Circle(double cx, double cy) {
this.cx = new EstimatedParameter("cx", cx);
this.cy = new EstimatedParameter(new EstimatedParameter("cy", cy));
points = new ArrayList<PointModel>();
}
public void addPoint(double px, double py) {
points.add(new PointModel(px, py));
}
public int getM() {
return points.size();
}
public WeightedMeasurement[] getMeasurements() {
return (WeightedMeasurement[]) points.toArray(new PointModel[points.size()]);
}
public EstimatedParameter[] getAllParameters() {
return new EstimatedParameter[] { cx, cy };
}
public EstimatedParameter[] getUnboundParameters() {
return new EstimatedParameter[] { cx, cy };
}
public double getPartialRadiusX() {
double dRdX = 0;
for (PointModel point : points) {
dRdX += point.getPartialDiX();
}
return dRdX / points.size();
}
public double getPartialRadiusY() {
double dRdY = 0;
for (PointModel point : points) {
dRdY += point.getPartialDiY();
}
return dRdY / points.size();
}
public double getRadius() {
double r = 0;
for (PointModel point : points) {
r += point.getCenterDistance();
}
return r / points.size();
}
public double getX() {
return cx.getEstimate();
}
public double getY() {
return cy.getEstimate();
}
private class PointModel extends WeightedMeasurement {
public PointModel(double px, double py) {
super(1.0, 0.0);
this.px = px;
this.py = py;
}
public double getPartial(EstimatedParameter parameter) {
if (parameter == cx) {
return getPartialDiX() - getPartialRadiusX();
} else if (parameter == cy) {
return getPartialDiY() - getPartialRadiusY();
}
return 0;
}
public double getCenterDistance() {
double dx = px - cx.getEstimate();
double dy = py - cy.getEstimate();
return Math.sqrt(dx * dx + dy * dy);
}
public double getPartialDiX() {
return (cx.getEstimate() - px) / getCenterDistance();
}
public double getPartialDiY() {
return (cy.getEstimate() - py) / getCenterDistance();
}
public double getTheoreticalValue() {
return getCenterDistance() - getRadius();
}
private double px;
private double py;
private static final long serialVersionUID = 1L;
}
private EstimatedParameter cx;
private EstimatedParameter cy;
private ArrayList<PointModel> points;
}
public static Test suite() {
return new TestSuite(GaussNewtonEstimatorTest.class);
}
}

844
src/test/org/apache/commons/math/estimation/LevenbergMarquardtEstimatorTest.java Normal file
View File

@ -0,0 +1,844 @@
/*
* 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.math.estimation;
import java.util.ArrayList;
import java.util.HashSet;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
/**
* <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>, for
* convenience, it is reproduced below.</p>
* <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
* <tr><td>
* Minpack Copyright Notice (1999) University of Chicago.
* All rights reserved
* </td></tr>
* <tr><td>
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* <ol>
* <li>Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.</li>
* <li>Redistributions in binary form must reproduce the above
* copyright notice, this list of conditions and the following
* disclaimer in the documentation and/or other materials provided
* with the distribution.</li>
* <li>The end-user documentation included with the redistribution, if any,
* must include the following acknowledgment:
* <code>This product includes software developed by the University of
* Chicago, as Operator of Argonne National Laboratory.</code>
* Alternately, this acknowledgment may appear in the software itself,
* if and wherever such third-party acknowledgments normally appear.</li>
* <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
* WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
* UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
* THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
* OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
* OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
* USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
* THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
* DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
* UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
* BE CORRECTED.</strong></li>
* <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
* HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
* ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
* INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
* ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
* PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
* SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
* (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
* EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
* POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
* <ol></td></tr>
* </table>
* @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
* @author Burton S. Garbow (original fortran minpack tests)
* @author Kenneth E. Hillstrom (original fortran minpack tests)
* @author Jorge J. More (original fortran minpack tests)
* @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
*/
@Deprecated
public class LevenbergMarquardtEstimatorTest
extends TestCase {
public LevenbergMarquardtEstimatorTest(String name) {
super(name);
}
public void testTrivial() throws EstimationException {
LinearProblem problem =
new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] {2},
new EstimatedParameter[] {
new EstimatedParameter("p0", 0)
}, 3.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
try {
estimator.guessParametersErrors(problem);
fail("an exception should have been thrown");
} catch (EstimationException ee) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
assertEquals(1.5,
problem.getUnboundParameters()[0].getEstimate(),
1.0e-10);
}
public void testQRColumnsPermutation() throws EstimationException {
EstimatedParameter[] x = {
new EstimatedParameter("p0", 0), new EstimatedParameter("p1", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { x[0], x[1] },
4.0),
new LinearMeasurement(new double[] { 2.0 },
new EstimatedParameter[] { x[1] },
6.0),
new LinearMeasurement(new double[] { 1.0, -2.0 },
new EstimatedParameter[] { x[0], x[1] },
1.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(7.0, x[0].getEstimate(), 1.0e-10);
assertEquals(3.0, x[1].getEstimate(), 1.0e-10);
}
public void testNoDependency() throws EstimationException {
EstimatedParameter[] p = new EstimatedParameter[] {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0),
new EstimatedParameter("p3", 0),
new EstimatedParameter("p4", 0),
new EstimatedParameter("p5", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[0] }, 0.0),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[1] }, 1.1),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[2] }, 2.2),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[3] }, 3.3),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[4] }, 4.4),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[5] }, 5.5)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
for (int i = 0; i < p.length; ++i) {
assertEquals(0.55 * i, p[i].getEstimate(), 1.0e-10);
}
}
public void testOneSet() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0 },
new EstimatedParameter[] { p[0] },
1.0),
new LinearMeasurement(new double[] { -1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { -1.0, 1.0 },
new EstimatedParameter[] { p[1], p[2] },
1.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
assertEquals(2.0, p[1].getEstimate(), 1.0e-10);
assertEquals(3.0, p[2].getEstimate(), 1.0e-10);
}
public void testTwoSets() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 1),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 3),
new EstimatedParameter("p4", 4),
new EstimatedParameter("p5", 5)
};
double epsilon = 1.0e-7;
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
// 4 elements sub-problem
new LinearMeasurement(new double[] { 2.0, 1.0, 4.0 },
new EstimatedParameter[] { p[0], p[1], p[3] },
2.0),
new LinearMeasurement(new double[] { -4.0, -2.0, 3.0, -7.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-9.0),
new LinearMeasurement(new double[] { 4.0, 1.0, -2.0, 8.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
2.0),
new LinearMeasurement(new double[] { -3.0, -12.0, -1.0 },
new EstimatedParameter[] { p[1], p[2], p[3] },
2.0),
// 2 elements sub-problem
new LinearMeasurement(new double[] { epsilon, 1.0 },
new EstimatedParameter[] { p[4], p[5] },
1.0 + epsilon * epsilon),
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[4], p[5] },
2.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals( 3.0, p[0].getEstimate(), 1.0e-10);
assertEquals( 4.0, p[1].getEstimate(), 1.0e-10);
assertEquals(-1.0, p[2].getEstimate(), 1.0e-10);
assertEquals(-2.0, p[3].getEstimate(), 1.0e-10);
assertEquals( 1.0 + epsilon, p[4].getEstimate(), 1.0e-10);
assertEquals( 1.0 - epsilon, p[5].getEstimate(), 1.0e-10);
}
public void testNonInversible() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0)
};
LinearMeasurement[] m = new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 2.0, -3.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
1.0),
new LinearMeasurement(new double[] { 2.0, 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
1.0),
new LinearMeasurement(new double[] { -3.0, -9.0 },
new EstimatedParameter[] { p[0], p[2] },
1.0)
};
LinearProblem problem = new LinearProblem(m);
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
double initialCost = estimator.getRMS(problem);
estimator.estimate(problem);
assertTrue(estimator.getRMS(problem) < initialCost);
assertTrue(Math.sqrt(m.length) * estimator.getRMS(problem) > 0.6);
try {
estimator.getCovariances(problem);
fail("an exception should have been thrown");
} catch (EstimationException ee) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
double dJ0 = 2 * (m[0].getResidual() * m[0].getPartial(p[0])
+ m[1].getResidual() * m[1].getPartial(p[0])
+ m[2].getResidual() * m[2].getPartial(p[0]));
double dJ1 = 2 * (m[0].getResidual() * m[0].getPartial(p[1])
+ m[1].getResidual() * m[1].getPartial(p[1]));
double dJ2 = 2 * (m[0].getResidual() * m[0].getPartial(p[2])
+ m[1].getResidual() * m[1].getPartial(p[2])
+ m[2].getResidual() * m[2].getPartial(p[2]));
assertEquals(0, dJ0, 1.0e-10);
assertEquals(0, dJ1, 1.0e-10);
assertEquals(0, dJ2, 1.0e-10);
}
public void testIllConditioned() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 1),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 3)
};
LinearProblem problem1 = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 10.0, 7.0, 8.0, 7.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
32.0),
new LinearMeasurement(new double[] { 7.0, 5.0, 6.0, 5.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
23.0),
new LinearMeasurement(new double[] { 8.0, 6.0, 10.0, 9.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
33.0),
new LinearMeasurement(new double[] { 7.0, 5.0, 9.0, 10.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
31.0)
});
LevenbergMarquardtEstimator estimator1 = new LevenbergMarquardtEstimator();
estimator1.estimate(problem1);
assertEquals(0, estimator1.getRMS(problem1), 1.0e-10);
assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
assertEquals(1.0, p[1].getEstimate(), 1.0e-10);
assertEquals(1.0, p[2].getEstimate(), 1.0e-10);
assertEquals(1.0, p[3].getEstimate(), 1.0e-10);
LinearProblem problem2 = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 10.0, 7.0, 8.1, 7.2 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
32.0),
new LinearMeasurement(new double[] { 7.08, 5.04, 6.0, 5.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
23.0),
new LinearMeasurement(new double[] { 8.0, 5.98, 9.89, 9.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
33.0),
new LinearMeasurement(new double[] { 6.99, 4.99, 9.0, 9.98 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
31.0)
});
LevenbergMarquardtEstimator estimator2 = new LevenbergMarquardtEstimator();
estimator2.estimate(problem2);
assertEquals(0, estimator2.getRMS(problem2), 1.0e-10);
assertEquals(-81.0, p[0].getEstimate(), 1.0e-8);
assertEquals(137.0, p[1].getEstimate(), 1.0e-8);
assertEquals(-34.0, p[2].getEstimate(), 1.0e-8);
assertEquals( 22.0, p[3].getEstimate(), 1.0e-8);
}
public void testMoreEstimatedParametersSimple() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 7),
new EstimatedParameter("p1", 6),
new EstimatedParameter("p2", 5),
new EstimatedParameter("p3", 4)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 3.0, 2.0 },
new EstimatedParameter[] { p[0], p[1] },
7.0),
new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
new EstimatedParameter[] { p[1], p[2], p[3] },
3.0),
new LinearMeasurement(new double[] { 2.0, 1.0 },
new EstimatedParameter[] { p[0], p[2] },
5.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
}
public void testMoreEstimatedParametersUnsorted() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 2),
new EstimatedParameter("p1", 2),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 2),
new EstimatedParameter("p4", 2),
new EstimatedParameter("p5", 2)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 },
new EstimatedParameter[] { p[2], p[3], p[4] },
12.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[4], p[5] },
-1.0),
new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
new EstimatedParameter[] { p[3], p[2], p[5] },
7.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[4], p[3] },
1.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(3.0, p[2].getEstimate(), 1.0e-10);
assertEquals(4.0, p[3].getEstimate(), 1.0e-10);
assertEquals(5.0, p[4].getEstimate(), 1.0e-10);
assertEquals(6.0, p[5].getEstimate(), 1.0e-10);
}
public void testRedundantEquations() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 1),
new EstimatedParameter("p1", 1)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1] },
5.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(2.0, p[0].getEstimate(), 1.0e-10);
assertEquals(1.0, p[1].getEstimate(), 1.0e-10);
}
public void testInconsistentEquations() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 1),
new EstimatedParameter("p1", 1)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1] },
4.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertTrue(estimator.getRMS(problem) > 0.1);
}
public void testControlParameters() throws EstimationException {
Circle circle = new Circle(98.680, 47.345);
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, 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false);
checkEstimate(circle, 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true);
checkEstimate(circle, 0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true);
circle.addPoint(300, -300);
checkEstimate(circle, 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true);
}
private void checkEstimate(EstimationProblem problem,
double initialStepBoundFactor, int maxCostEval,
double costRelativeTolerance, double parRelativeTolerance,
double orthoTolerance, boolean shouldFail) {
try {
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.setInitialStepBoundFactor(initialStepBoundFactor);
estimator.setMaxCostEval(maxCostEval);
estimator.setCostRelativeTolerance(costRelativeTolerance);
estimator.setParRelativeTolerance(parRelativeTolerance);
estimator.setOrthoTolerance(orthoTolerance);
estimator.estimate(problem);
assertTrue(! shouldFail);
} catch (EstimationException ee) {
assertTrue(shouldFail);
} catch (Exception e) {
fail("wrong exception type caught");
}
}
public void testCircleFitting() throws EstimationException {
Circle circle = new Circle(98.680, 47.345);
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);
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(circle);
assertTrue(estimator.getCostEvaluations() < 10);
assertTrue(estimator.getJacobianEvaluations() < 10);
double rms = estimator.getRMS(circle);
assertEquals(1.768262623567235, Math.sqrt(circle.getM()) * rms, 1.0e-10);
assertEquals(69.96016176931406, circle.getRadius(), 1.0e-10);
assertEquals(96.07590211815305, circle.getX(), 1.0e-10);
assertEquals(48.13516790438953, circle.getY(), 1.0e-10);
double[][] cov = estimator.getCovariances(circle);
assertEquals(1.839, cov[0][0], 0.001);
assertEquals(0.731, cov[0][1], 0.001);
assertEquals(cov[0][1], cov[1][0], 1.0e-14);
assertEquals(0.786, cov[1][1], 0.001);
double[] errors = estimator.guessParametersErrors(circle);
assertEquals(1.384, errors[0], 0.001);
assertEquals(0.905, errors[1], 0.001);
// add perfect measurements and check errors are reduced
double cx = circle.getX();
double cy = circle.getY();
double r = circle.getRadius();
for (double d= 0; d < 2 * Math.PI; d += 0.01) {
circle.addPoint(cx + r * Math.cos(d), cy + r * Math.sin(d));
}
estimator = new LevenbergMarquardtEstimator();
estimator.estimate(circle);
cov = estimator.getCovariances(circle);
assertEquals(0.004, cov[0][0], 0.001);
assertEquals(6.40e-7, cov[0][1], 1.0e-9);
assertEquals(cov[0][1], cov[1][0], 1.0e-14);
assertEquals(0.003, cov[1][1], 0.001);
errors = estimator.guessParametersErrors(circle);
assertEquals(0.004, errors[0], 0.001);
assertEquals(0.004, errors[1], 0.001);
}
public void testCircleFittingBadInit() throws EstimationException {
Circle circle = new Circle(-12, -12);
double[][] points = 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}
};
for (int i = 0; i < points.length; ++i) {
circle.addPoint(points[i][0], points[i][1]);
}
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(circle);
assertTrue(estimator.getCostEvaluations() < 15);
assertTrue(estimator.getJacobianEvaluations() < 10);
assertEquals( 0.030184491196225207, estimator.getRMS(circle), 1.0e-9);
assertEquals( 0.2922350065939634, circle.getRadius(), 1.0e-9);
assertEquals(-0.15173845023862165, circle.getX(), 1.0e-8);
assertEquals( 0.20750021499570379, circle.getY(), 1.0e-8);
}
public void testMath199() {
try {
QuadraticProblem problem = new QuadraticProblem();
problem.addPoint (0, -3.182591015485607, 0.0);
problem.addPoint (1, -2.5581184967730577, 4.4E-323);
problem.addPoint (2, -2.1488478161387325, 1.0);
problem.addPoint (3, -1.9122489313410047, 4.4E-323);
problem.addPoint (4, 1.7785661310051026, 0.0);
new LevenbergMarquardtEstimator().estimate(problem);
fail("an exception should have been thrown");
} catch (EstimationException ee) {
// expected behavior
}
}
private static class LinearProblem implements EstimationProblem {
public LinearProblem(LinearMeasurement[] measurements) {
this.measurements = measurements;
}
public WeightedMeasurement[] getMeasurements() {
return measurements;
}
public EstimatedParameter[] getUnboundParameters() {
return getAllParameters();
}
public EstimatedParameter[] getAllParameters() {
HashSet<EstimatedParameter> set = new HashSet<EstimatedParameter>();
for (int i = 0; i < measurements.length; ++i) {
EstimatedParameter[] parameters = measurements[i].getParameters();
for (int j = 0; j < parameters.length; ++j) {
set.add(parameters[j]);
}
}
return (EstimatedParameter[]) set.toArray(new EstimatedParameter[set.size()]);
}
private LinearMeasurement[] measurements;
}
private static class LinearMeasurement extends WeightedMeasurement {
public LinearMeasurement(double[] factors, EstimatedParameter[] parameters,
double setPoint) {
super(1.0, setPoint);
this.factors = factors;
this.parameters = parameters;
}
public double getTheoreticalValue() {
double v = 0;
for (int i = 0; i < factors.length; ++i) {
v += factors[i] * parameters[i].getEstimate();
}
return v;
}
public double getPartial(EstimatedParameter parameter) {
for (int i = 0; i < parameters.length; ++i) {
if (parameters[i] == parameter) {
return factors[i];
}
}
return 0;
}
public EstimatedParameter[] getParameters() {
return parameters;
}
private double[] factors;
private EstimatedParameter[] parameters;
private static final long serialVersionUID = -3922448707008868580L;
}
private static class Circle implements EstimationProblem {
public Circle(double cx, double cy) {
this.cx = new EstimatedParameter("cx", cx);
this.cy = new EstimatedParameter("cy", cy);
points = new ArrayList<PointModel>();
}
public void addPoint(double px, double py) {
points.add(new PointModel(px, py));
}
public int getM() {
return points.size();
}
public WeightedMeasurement[] getMeasurements() {
return (WeightedMeasurement[]) points.toArray(new PointModel[points.size()]);
}
public EstimatedParameter[] getAllParameters() {
return new EstimatedParameter[] { cx, cy };
}
public EstimatedParameter[] getUnboundParameters() {
return new EstimatedParameter[] { cx, cy };
}
public double getPartialRadiusX() {
double dRdX = 0;
for (PointModel point : points) {
dRdX += point.getPartialDiX();
}
return dRdX / points.size();
}
public double getPartialRadiusY() {
double dRdY = 0;
for (PointModel point : points) {
dRdY += point.getPartialDiY();
}
return dRdY / points.size();
}
public double getRadius() {
double r = 0;
for (PointModel point : points) {
r += point.getCenterDistance();
}
return r / points.size();
}
public double getX() {
return cx.getEstimate();
}
public double getY() {
return cy.getEstimate();
}
private class PointModel extends WeightedMeasurement {
public PointModel(double px, double py) {
super(1.0, 0.0);
this.px = px;
this.py = py;
}
public double getPartial(EstimatedParameter parameter) {
if (parameter == cx) {
return getPartialDiX() - getPartialRadiusX();
} else if (parameter == cy) {
return getPartialDiY() - getPartialRadiusY();
}
return 0;
}
public double getCenterDistance() {
double dx = px - cx.getEstimate();
double dy = py - cy.getEstimate();
return Math.sqrt(dx * dx + dy * dy);
}
public double getPartialDiX() {
return (cx.getEstimate() - px) / getCenterDistance();
}
public double getPartialDiY() {
return (cy.getEstimate() - py) / getCenterDistance();
}
public double getTheoreticalValue() {
return getCenterDistance() - getRadius();
}
private double px;
private double py;
private static final long serialVersionUID = 1L;
}
private EstimatedParameter cx;
private EstimatedParameter cy;
private ArrayList<PointModel> points;
}
private static class QuadraticProblem extends SimpleEstimationProblem {
private EstimatedParameter a;
private EstimatedParameter b;
private EstimatedParameter c;
public QuadraticProblem() {
a = new EstimatedParameter("a", 0.0);
b = new EstimatedParameter("b", 0.0);
c = new EstimatedParameter("c", 0.0);
addParameter(a);
addParameter(b);
addParameter(c);
}
public void addPoint(double x, double y, double w) {
addMeasurement(new LocalMeasurement(x, y, w));
}
public double getA() {
return a.getEstimate();
}
public double getB() {
return b.getEstimate();
}
public double getC() {
return c.getEstimate();
}
public double theoreticalValue(double x) {
return ( (a.getEstimate() * x + b.getEstimate() ) * x + c.getEstimate());
}
private double partial(double x, EstimatedParameter parameter) {
if (parameter == a) {
return x * x;
} else if (parameter == b) {
return x;
} else {
return 1.0;
}
}
private class LocalMeasurement extends WeightedMeasurement {
private static final long serialVersionUID = 1555043155023729130L;
private final double x;
// constructor
public LocalMeasurement(double x, double y, double w) {
super(w, y);
this.x = x;
}
public double getTheoreticalValue() {
return theoreticalValue(x);
}
public double getPartial(EstimatedParameter parameter) {
return partial(x, parameter);
}
}
}
public static Test suite() {
return new TestSuite(LevenbergMarquardtEstimatorTest.class);
}
}

1518
src/test/org/apache/commons/math/estimation/MinpackTest.java Normal file
View File

File diff suppressed because it is too large Load Diff

127
src/test/org/apache/commons/math/estimation/WeightedMeasurementTest.java Normal file
View File

@ -0,0 +1,127 @@
/*
* 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.math.estimation;
import org.apache.commons.math.estimation.EstimatedParameter;
import org.apache.commons.math.estimation.WeightedMeasurement;
import junit.framework.*;
@Deprecated
public class WeightedMeasurementTest
extends TestCase {
public WeightedMeasurementTest(String name) {
super(name);
p1 = null;
p2 = null;
}
public void testConstruction() {
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this);
checkValue(m.getWeight(), 3.0);
checkValue(m.getMeasuredValue(), theoretical() + 0.1);
}
public void testIgnored() {
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this);
assertTrue(!m.isIgnored());
m.setIgnored(true);
assertTrue(m.isIgnored());
m.setIgnored(false);
assertTrue(!m.isIgnored());
}
public void testTheory() {
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this);
checkValue(m.getTheoreticalValue(), theoretical());
checkValue(m.getResidual(), 0.1);
double oldP1 = p1.getEstimate();
p1.setEstimate(oldP1 + m.getResidual() / m.getPartial(p1));
checkValue(m.getResidual(), 0.0);
p1.setEstimate(oldP1);
checkValue(m.getResidual(), 0.1);
double oldP2 = p2.getEstimate();
p2.setEstimate(oldP2 + m.getResidual() / m.getPartial(p2));
checkValue(m.getResidual(), 0.0);
p2.setEstimate(oldP2);
checkValue(m.getResidual(), 0.1);
}
public static Test suite() {
return new TestSuite(WeightedMeasurementTest.class);
}
public void setUp() {
p1 = new EstimatedParameter("p1", 1.0);
p2 = new EstimatedParameter("p2", 2.0);
}
public void tearDown() {
p1 = null;
p2 = null;
}
private void checkValue(double value, double expected) {
assertTrue(Math.abs(value - expected) < 1.0e-10);
}
private double theoretical() {
return 3 * p1.getEstimate() - p2.getEstimate();
}
private double partial(EstimatedParameter p) {
if (p == p1) {
return 3.0;
} else if (p == p2) {
return -1.0;
} else {
return 0.0;
}
}
private static class MyMeasurement
extends WeightedMeasurement {
public MyMeasurement(double weight, double measuredValue,
WeightedMeasurementTest testInstance) {
super(weight, measuredValue);
this.testInstance = testInstance;
}
public double getTheoreticalValue() {
return testInstance.theoretical();
}
public double getPartial(EstimatedParameter p) {
return testInstance.partial(p);
}
private transient WeightedMeasurementTest testInstance;
private static final long serialVersionUID = -246712922500792332L;
}
private EstimatedParameter p1;
private EstimatedParameter p2;
}