1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.math.estimation; 19 20 import java.util.Arrays; 21 22 import org.apache.commons.math.exception.util.LocalizedFormats; 23 import org.apache.commons.math.linear.InvalidMatrixException; 24 import org.apache.commons.math.linear.LUDecompositionImpl; 25 import org.apache.commons.math.linear.MatrixUtils; 26 import org.apache.commons.math.linear.RealMatrix; 27 import org.apache.commons.math.util.FastMath; 28 29 /** 30 * Base class for implementing estimators. 31 * <p>This base class handles the boilerplates methods associated to thresholds 32 * settings, jacobian and error estimation.</p> 33 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 aot 2010) $ 34 * @since 1.2 35 * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has 36 * been deprecated and replaced by package org.apache.commons.math.optimization.general 37 * 38 */ 39 @Deprecated 40 public abstract class AbstractEstimator implements Estimator { 41 42 /** Default maximal number of cost evaluations allowed. */ 43 public static final int DEFAULT_MAX_COST_EVALUATIONS = 100; 44 45 /** Array of measurements. */ 46 protected WeightedMeasurement[] measurements; 47 48 /** Array of parameters. */ 49 protected EstimatedParameter[] parameters; 50 51 /** 52 * Jacobian matrix. 53 * <p>This matrix is in canonical form just after the calls to 54 * {@link #updateJacobian()}, but may be modified by the solver 55 * in the derived class (the {@link LevenbergMarquardtEstimator 56 * Levenberg-Marquardt estimator} does this).</p> 57 */ 58 protected double[] jacobian; 59 60 /** Number of columns of the jacobian matrix. */ 61 protected int cols; 62 63 /** Number of rows of the jacobian matrix. */ 64 protected int rows; 65 66 /** Residuals array. 67 * <p>This array is in canonical form just after the calls to 68 * {@link #updateJacobian()}, but may be modified by the solver 69 * in the derived class (the {@link LevenbergMarquardtEstimator 70 * Levenberg-Marquardt estimator} does this).</p> 71 */ 72 protected double[] residuals; 73 74 /** Cost value (square root of the sum of the residuals). */ 75 protected double cost; 76 77 /** Maximal allowed number of cost evaluations. */ 78 private int maxCostEval; 79 80 /** Number of cost evaluations. */ 81 private int costEvaluations; 82 83 /** Number of jacobian evaluations. */ 84 private int jacobianEvaluations; 85 86 /** 87 * Build an abstract estimator for least squares problems. 88 * <p>The maximal number of cost evaluations allowed is set 89 * to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.</p> 90 */ 91 protected AbstractEstimator() { 92 setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS); 93 } 94 95 /** 96 * Set the maximal number of cost evaluations allowed. 97 * 98 * @param maxCostEval maximal number of cost evaluations allowed 99 * @see #estimate 100 */ 101 public final void setMaxCostEval(int maxCostEval) { 102 this.maxCostEval = maxCostEval; 103 } 104 105 /** 106 * Get the number of cost evaluations. 107 * 108 * @return number of cost evaluations 109 * */ 110 public final int getCostEvaluations() { 111 return costEvaluations; 112 } 113 114 /** 115 * Get the number of jacobian evaluations. 116 * 117 * @return number of jacobian evaluations 118 * */ 119 public final int getJacobianEvaluations() { 120 return jacobianEvaluations; 121 } 122 123 /** 124 * Update the jacobian matrix. 125 */ 126 protected void updateJacobian() { 127 incrementJacobianEvaluationsCounter(); 128 Arrays.fill(jacobian, 0); 129 int index = 0; 130 for (int i = 0; i < rows; i++) { 131 WeightedMeasurement wm = measurements[i]; 132 double factor = -FastMath.sqrt(wm.getWeight()); 133 for (int j = 0; j < cols; ++j) { 134 jacobian[index++] = factor * wm.getPartial(parameters[j]); 135 } 136 } 137 } 138 139 /** 140 * Increment the jacobian evaluations counter. 141 */ 142 protected final void incrementJacobianEvaluationsCounter() { 143 ++jacobianEvaluations; 144 } 145 146 /** 147 * Update the residuals array and cost function value. 148 * @exception EstimationException if the number of cost evaluations 149 * exceeds the maximum allowed 150 */ 151 protected void updateResidualsAndCost() 152 throws EstimationException { 153 154 if (++costEvaluations > maxCostEval) { 155 throw new EstimationException(LocalizedFormats.MAX_EVALUATIONS_EXCEEDED, 156 maxCostEval); 157 } 158 159 cost = 0; 160 int index = 0; 161 for (int i = 0; i < rows; i++, index += cols) { 162 WeightedMeasurement wm = measurements[i]; 163 double residual = wm.getResidual(); 164 residuals[i] = FastMath.sqrt(wm.getWeight()) * residual; 165 cost += wm.getWeight() * residual * residual; 166 } 167 cost = FastMath.sqrt(cost); 168 169 } 170 171 /** 172 * Get the Root Mean Square value. 173 * Get the Root Mean Square value, i.e. the root of the arithmetic 174 * mean of the square of all weighted residuals. This is related to the 175 * criterion that is minimized by the estimator as follows: if 176 * <em>c</em> if the criterion, and <em>n</em> is the number of 177 * measurements, then the RMS is <em>sqrt (c/n)</em>. 178 * 179 * @param problem estimation problem 180 * @return RMS value 181 */ 182 public double getRMS(EstimationProblem problem) { 183 WeightedMeasurement[] wm = problem.getMeasurements(); 184 double criterion = 0; 185 for (int i = 0; i < wm.length; ++i) { 186 double residual = wm[i].getResidual(); 187 criterion += wm[i].getWeight() * residual * residual; 188 } 189 return FastMath.sqrt(criterion / wm.length); 190 } 191 192 /** 193 * Get the Chi-Square value. 194 * @param problem estimation problem 195 * @return chi-square value 196 */ 197 public double getChiSquare(EstimationProblem problem) { 198 WeightedMeasurement[] wm = problem.getMeasurements(); 199 double chiSquare = 0; 200 for (int i = 0; i < wm.length; ++i) { 201 double residual = wm[i].getResidual(); 202 chiSquare += residual * residual / wm[i].getWeight(); 203 } 204 return chiSquare; 205 } 206 207 /** 208 * Get the covariance matrix of unbound estimated parameters. 209 * @param problem estimation problem 210 * @return covariance matrix 211 * @exception EstimationException if the covariance matrix 212 * cannot be computed (singular problem) 213 */ 214 public double[][] getCovariances(EstimationProblem problem) 215 throws EstimationException { 216 217 // set up the jacobian 218 updateJacobian(); 219 220 // compute transpose(J).J, avoiding building big intermediate matrices 221 final int n = problem.getMeasurements().length; 222 final int m = problem.getUnboundParameters().length; 223 final int max = m * n; 224 double[][] jTj = new double[m][m]; 225 for (int i = 0; i < m; ++i) { 226 for (int j = i; j < m; ++j) { 227 double sum = 0; 228 for (int k = 0; k < max; k += m) { 229 sum += jacobian[k + i] * jacobian[k + j]; 230 } 231 jTj[i][j] = sum; 232 jTj[j][i] = sum; 233 } 234 } 235 236 try { 237 // compute the covariances matrix 238 RealMatrix inverse = 239 new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse(); 240 return inverse.getData(); 241 } catch (InvalidMatrixException ime) { 242 throw new EstimationException(LocalizedFormats.UNABLE_TO_COMPUTE_COVARIANCE_SINGULAR_PROBLEM); 243 } 244 245 } 246 247 /** 248 * Guess the errors in unbound estimated parameters. 249 * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p> 250 * @param problem estimation problem 251 * @return errors in estimated parameters 252 * @exception EstimationException if the covariances matrix cannot be computed 253 * or the number of degrees of freedom is not positive (number of measurements 254 * lesser or equal to number of parameters) 255 */ 256 public double[] guessParametersErrors(EstimationProblem problem) 257 throws EstimationException { 258 int m = problem.getMeasurements().length; 259 int p = problem.getUnboundParameters().length; 260 if (m <= p) { 261 throw new EstimationException( 262 LocalizedFormats.NO_DEGREES_OF_FREEDOM, 263 m, p); 264 } 265 double[] errors = new double[problem.getUnboundParameters().length]; 266 final double c = FastMath.sqrt(getChiSquare(problem) / (m - p)); 267 double[][] covar = getCovariances(problem); 268 for (int i = 0; i < errors.length; ++i) { 269 errors[i] = FastMath.sqrt(covar[i][i]) * c; 270 } 271 return errors; 272 } 273 274 /** 275 * Initialization of the common parts of the estimation. 276 * <p>This method <em>must</em> be called at the start 277 * of the {@link #estimate(EstimationProblem) estimate} 278 * method.</p> 279 * @param problem estimation problem to solve 280 */ 281 protected void initializeEstimate(EstimationProblem problem) { 282 283 // reset counters 284 costEvaluations = 0; 285 jacobianEvaluations = 0; 286 287 // retrieve the equations and the parameters 288 measurements = problem.getMeasurements(); 289 parameters = problem.getUnboundParameters(); 290 291 // arrays shared with the other private methods 292 rows = measurements.length; 293 cols = parameters.length; 294 jacobian = new double[rows * cols]; 295 residuals = new double[rows]; 296 297 cost = Double.POSITIVE_INFINITY; 298 299 } 300 301 /** 302 * Solve an estimation problem. 303 * 304 * <p>The method should set the parameters of the problem to several 305 * trial values until it reaches convergence. If this method returns 306 * normally (i.e. without throwing an exception), then the best 307 * estimate of the parameters can be retrieved from the problem 308 * itself, through the {@link EstimationProblem#getAllParameters 309 * EstimationProblem.getAllParameters} method.</p> 310 * 311 * @param problem estimation problem to solve 312 * @exception EstimationException if the problem cannot be solved 313 * 314 */ 315 public abstract void estimate(EstimationProblem problem) 316 throws EstimationException; 317 318 } 319
