| /* |
| * 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; |
| |
| } |