/*
* 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.exception.util.LocalizedFormats;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.util.FastMath;
/**
* Base class for implementing estimators.
* <p>This base class handles the boilerplates methods associated to thresholds
* settings, jacobian and error estimation.</p>
* @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
* @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;
/** 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;
/**
* 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);
int index = 0;
for (int i = 0; i < rows; i++) {
WeightedMeasurement wm = measurements[i];
double factor = -FastMath.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(LocalizedFormats.MAX_EVALUATIONS_EXCEEDED,
maxCostEval);
}
cost = 0;
int index = 0;
for (int i = 0; i < rows; i++, index += cols) {
WeightedMeasurement wm = measurements[i];
double residual = wm.getResidual();
residuals[i] = FastMath.sqrt(wm.getWeight()) * residual;
cost += wm.getWeight() * residual * residual;
}
cost = FastMath.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 FastMath.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 n = problem.getMeasurements().length;
final int m = problem.getUnboundParameters().length;
final int max = m * n;
double[][] jTj = new double[m][m];
for (int i = 0; i < m; ++i) {
for (int j = i; j < m; ++j) {
double sum = 0;
for (int k = 0; k < max; k += m) {
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(LocalizedFormats.UNABLE_TO_COMPUTE_COVARIANCE_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(
LocalizedFormats.NO_DEGREES_OF_FREEDOM,
m, p);
}
double[] errors = new double[problem.getUnboundParameters().length];
final double c = FastMath.sqrt(getChiSquare(problem) / (m - p));
double[][] covar = getCovariances(problem);
for (int i = 0; i < errors.length; ++i) {
errors[i] = FastMath.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;
}