001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.estimation; 019 020 import java.util.Arrays; 021 022 import org.apache.commons.math.linear.InvalidMatrixException; 023 import org.apache.commons.math.linear.LUDecompositionImpl; 024 import org.apache.commons.math.linear.MatrixUtils; 025 import org.apache.commons.math.linear.RealMatrix; 026 027 /** 028 * Base class for implementing estimators. 029 * <p>This base class handles the boilerplates methods associated to thresholds 030 * settings, jacobian and error estimation.</p> 031 * @version $Revision: 781122 $ $Date: 2009-06-02 14:53:23 -0400 (Tue, 02 Jun 2009) $ 032 * @since 1.2 033 * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has 034 * been deprecated and replaced by package org.apache.commons.math.optimization.general 035 * 036 */ 037 @Deprecated 038 public abstract class AbstractEstimator implements Estimator { 039 040 /** Default maximal number of cost evaluations allowed. */ 041 public static final int DEFAULT_MAX_COST_EVALUATIONS = 100; 042 043 /** 044 * Build an abstract estimator for least squares problems. 045 * <p>The maximal number of cost evaluations allowed is set 046 * to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.</p> 047 */ 048 protected AbstractEstimator() { 049 setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS); 050 } 051 052 /** 053 * Set the maximal number of cost evaluations allowed. 054 * 055 * @param maxCostEval maximal number of cost evaluations allowed 056 * @see #estimate 057 */ 058 public final void setMaxCostEval(int maxCostEval) { 059 this.maxCostEval = maxCostEval; 060 } 061 062 /** 063 * Get the number of cost evaluations. 064 * 065 * @return number of cost evaluations 066 * */ 067 public final int getCostEvaluations() { 068 return costEvaluations; 069 } 070 071 /** 072 * Get the number of jacobian evaluations. 073 * 074 * @return number of jacobian evaluations 075 * */ 076 public final int getJacobianEvaluations() { 077 return jacobianEvaluations; 078 } 079 080 /** 081 * Update the jacobian matrix. 082 */ 083 protected void updateJacobian() { 084 incrementJacobianEvaluationsCounter(); 085 Arrays.fill(jacobian, 0); 086 for (int i = 0, index = 0; i < rows; i++) { 087 WeightedMeasurement wm = measurements[i]; 088 double factor = -Math.sqrt(wm.getWeight()); 089 for (int j = 0; j < cols; ++j) { 090 jacobian[index++] = factor * wm.getPartial(parameters[j]); 091 } 092 } 093 } 094 095 /** 096 * Increment the jacobian evaluations counter. 097 */ 098 protected final void incrementJacobianEvaluationsCounter() { 099 ++jacobianEvaluations; 100 } 101 102 /** 103 * Update the residuals array and cost function value. 104 * @exception EstimationException if the number of cost evaluations 105 * exceeds the maximum allowed 106 */ 107 protected void updateResidualsAndCost() 108 throws EstimationException { 109 110 if (++costEvaluations > maxCostEval) { 111 throw new EstimationException("maximal number of evaluations exceeded ({0})", 112 maxCostEval); 113 } 114 115 cost = 0; 116 for (int i = 0, index = 0; i < rows; i++, index += cols) { 117 WeightedMeasurement wm = measurements[i]; 118 double residual = wm.getResidual(); 119 residuals[i] = Math.sqrt(wm.getWeight()) * residual; 120 cost += wm.getWeight() * residual * residual; 121 } 122 cost = Math.sqrt(cost); 123 124 } 125 126 /** 127 * Get the Root Mean Square value. 128 * Get the Root Mean Square value, i.e. the root of the arithmetic 129 * mean of the square of all weighted residuals. This is related to the 130 * criterion that is minimized by the estimator as follows: if 131 * <em>c</em> if the criterion, and <em>n</em> is the number of 132 * measurements, then the RMS is <em>sqrt (c/n)</em>. 133 * 134 * @param problem estimation problem 135 * @return RMS value 136 */ 137 public double getRMS(EstimationProblem problem) { 138 WeightedMeasurement[] wm = problem.getMeasurements(); 139 double criterion = 0; 140 for (int i = 0; i < wm.length; ++i) { 141 double residual = wm[i].getResidual(); 142 criterion += wm[i].getWeight() * residual * residual; 143 } 144 return Math.sqrt(criterion / wm.length); 145 } 146 147 /** 148 * Get the Chi-Square value. 149 * @param problem estimation problem 150 * @return chi-square value 151 */ 152 public double getChiSquare(EstimationProblem problem) { 153 WeightedMeasurement[] wm = problem.getMeasurements(); 154 double chiSquare = 0; 155 for (int i = 0; i < wm.length; ++i) { 156 double residual = wm[i].getResidual(); 157 chiSquare += residual * residual / wm[i].getWeight(); 158 } 159 return chiSquare; 160 } 161 162 /** 163 * Get the covariance matrix of unbound estimated parameters. 164 * @param problem estimation problem 165 * @return covariance matrix 166 * @exception EstimationException if the covariance matrix 167 * cannot be computed (singular problem) 168 */ 169 public double[][] getCovariances(EstimationProblem problem) 170 throws EstimationException { 171 172 // set up the jacobian 173 updateJacobian(); 174 175 // compute transpose(J).J, avoiding building big intermediate matrices 176 final int rows = problem.getMeasurements().length; 177 final int cols = problem.getUnboundParameters().length; 178 final int max = cols * rows; 179 double[][] jTj = new double[cols][cols]; 180 for (int i = 0; i < cols; ++i) { 181 for (int j = i; j < cols; ++j) { 182 double sum = 0; 183 for (int k = 0; k < max; k += cols) { 184 sum += jacobian[k + i] * jacobian[k + j]; 185 } 186 jTj[i][j] = sum; 187 jTj[j][i] = sum; 188 } 189 } 190 191 try { 192 // compute the covariances matrix 193 RealMatrix inverse = 194 new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse(); 195 return inverse.getData(); 196 } catch (InvalidMatrixException ime) { 197 throw new EstimationException("unable to compute covariances: singular problem"); 198 } 199 200 } 201 202 /** 203 * Guess the errors in unbound estimated parameters. 204 * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p> 205 * @param problem estimation problem 206 * @return errors in estimated parameters 207 * @exception EstimationException if the covariances matrix cannot be computed 208 * or the number of degrees of freedom is not positive (number of measurements 209 * lesser or equal to number of parameters) 210 */ 211 public double[] guessParametersErrors(EstimationProblem problem) 212 throws EstimationException { 213 int m = problem.getMeasurements().length; 214 int p = problem.getUnboundParameters().length; 215 if (m <= p) { 216 throw new EstimationException( 217 "no degrees of freedom ({0} measurements, {1} parameters)", 218 m, p); 219 } 220 double[] errors = new double[problem.getUnboundParameters().length]; 221 final double c = Math.sqrt(getChiSquare(problem) / (m - p)); 222 double[][] covar = getCovariances(problem); 223 for (int i = 0; i < errors.length; ++i) { 224 errors[i] = Math.sqrt(covar[i][i]) * c; 225 } 226 return errors; 227 } 228 229 /** 230 * Initialization of the common parts of the estimation. 231 * <p>This method <em>must</em> be called at the start 232 * of the {@link #estimate(EstimationProblem) estimate} 233 * method.</p> 234 * @param problem estimation problem to solve 235 */ 236 protected void initializeEstimate(EstimationProblem problem) { 237 238 // reset counters 239 costEvaluations = 0; 240 jacobianEvaluations = 0; 241 242 // retrieve the equations and the parameters 243 measurements = problem.getMeasurements(); 244 parameters = problem.getUnboundParameters(); 245 246 // arrays shared with the other private methods 247 rows = measurements.length; 248 cols = parameters.length; 249 jacobian = new double[rows * cols]; 250 residuals = new double[rows]; 251 252 cost = Double.POSITIVE_INFINITY; 253 254 } 255 256 /** 257 * Solve an estimation problem. 258 * 259 * <p>The method should set the parameters of the problem to several 260 * trial values until it reaches convergence. If this method returns 261 * normally (i.e. without throwing an exception), then the best 262 * estimate of the parameters can be retrieved from the problem 263 * itself, through the {@link EstimationProblem#getAllParameters 264 * EstimationProblem.getAllParameters} method.</p> 265 * 266 * @param problem estimation problem to solve 267 * @exception EstimationException if the problem cannot be solved 268 * 269 */ 270 public abstract void estimate(EstimationProblem problem) 271 throws EstimationException; 272 273 /** Array of measurements. */ 274 protected WeightedMeasurement[] measurements; 275 276 /** Array of parameters. */ 277 protected EstimatedParameter[] parameters; 278 279 /** 280 * Jacobian matrix. 281 * <p>This matrix is in canonical form just after the calls to 282 * {@link #updateJacobian()}, but may be modified by the solver 283 * in the derived class (the {@link LevenbergMarquardtEstimator 284 * Levenberg-Marquardt estimator} does this).</p> 285 */ 286 protected double[] jacobian; 287 288 /** Number of columns of the jacobian matrix. */ 289 protected int cols; 290 291 /** Number of rows of the jacobian matrix. */ 292 protected int rows; 293 294 /** Residuals array. 295 * <p>This array is in canonical form just after the calls to 296 * {@link #updateJacobian()}, but may be modified by the solver 297 * in the derived class (the {@link LevenbergMarquardtEstimator 298 * Levenberg-Marquardt estimator} does this).</p> 299 */ 300 protected double[] residuals; 301 302 /** Cost value (square root of the sum of the residuals). */ 303 protected double cost; 304 305 /** Maximal allowed number of cost evaluations. */ 306 private int maxCostEval; 307 308 /** Number of cost evaluations. */ 309 private int costEvaluations; 310 311 /** Number of jacobian evaluations. */ 312 private int jacobianEvaluations; 313 314 }