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: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 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 /** Array of measurements. */ 044 protected WeightedMeasurement[] measurements; 045 046 /** Array of parameters. */ 047 protected EstimatedParameter[] parameters; 048 049 /** 050 * Jacobian matrix. 051 * <p>This matrix is in canonical form just after the calls to 052 * {@link #updateJacobian()}, but may be modified by the solver 053 * in the derived class (the {@link LevenbergMarquardtEstimator 054 * Levenberg-Marquardt estimator} does this).</p> 055 */ 056 protected double[] jacobian; 057 058 /** Number of columns of the jacobian matrix. */ 059 protected int cols; 060 061 /** Number of rows of the jacobian matrix. */ 062 protected int rows; 063 064 /** Residuals array. 065 * <p>This array is in canonical form just after the calls to 066 * {@link #updateJacobian()}, but may be modified by the solver 067 * in the derived class (the {@link LevenbergMarquardtEstimator 068 * Levenberg-Marquardt estimator} does this).</p> 069 */ 070 protected double[] residuals; 071 072 /** Cost value (square root of the sum of the residuals). */ 073 protected double cost; 074 075 /** Maximal allowed number of cost evaluations. */ 076 private int maxCostEval; 077 078 /** Number of cost evaluations. */ 079 private int costEvaluations; 080 081 /** Number of jacobian evaluations. */ 082 private int jacobianEvaluations; 083 084 /** 085 * Build an abstract estimator for least squares problems. 086 * <p>The maximal number of cost evaluations allowed is set 087 * to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.</p> 088 */ 089 protected AbstractEstimator() { 090 setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS); 091 } 092 093 /** 094 * Set the maximal number of cost evaluations allowed. 095 * 096 * @param maxCostEval maximal number of cost evaluations allowed 097 * @see #estimate 098 */ 099 public final void setMaxCostEval(int maxCostEval) { 100 this.maxCostEval = maxCostEval; 101 } 102 103 /** 104 * Get the number of cost evaluations. 105 * 106 * @return number of cost evaluations 107 * */ 108 public final int getCostEvaluations() { 109 return costEvaluations; 110 } 111 112 /** 113 * Get the number of jacobian evaluations. 114 * 115 * @return number of jacobian evaluations 116 * */ 117 public final int getJacobianEvaluations() { 118 return jacobianEvaluations; 119 } 120 121 /** 122 * Update the jacobian matrix. 123 */ 124 protected void updateJacobian() { 125 incrementJacobianEvaluationsCounter(); 126 Arrays.fill(jacobian, 0); 127 int index = 0; 128 for (int i = 0; i < rows; i++) { 129 WeightedMeasurement wm = measurements[i]; 130 double factor = -Math.sqrt(wm.getWeight()); 131 for (int j = 0; j < cols; ++j) { 132 jacobian[index++] = factor * wm.getPartial(parameters[j]); 133 } 134 } 135 } 136 137 /** 138 * Increment the jacobian evaluations counter. 139 */ 140 protected final void incrementJacobianEvaluationsCounter() { 141 ++jacobianEvaluations; 142 } 143 144 /** 145 * Update the residuals array and cost function value. 146 * @exception EstimationException if the number of cost evaluations 147 * exceeds the maximum allowed 148 */ 149 protected void updateResidualsAndCost() 150 throws EstimationException { 151 152 if (++costEvaluations > maxCostEval) { 153 throw new EstimationException("maximal number of evaluations exceeded ({0})", 154 maxCostEval); 155 } 156 157 cost = 0; 158 int index = 0; 159 for (int i = 0; i < rows; i++, index += cols) { 160 WeightedMeasurement wm = measurements[i]; 161 double residual = wm.getResidual(); 162 residuals[i] = Math.sqrt(wm.getWeight()) * residual; 163 cost += wm.getWeight() * residual * residual; 164 } 165 cost = Math.sqrt(cost); 166 167 } 168 169 /** 170 * Get the Root Mean Square value. 171 * Get the Root Mean Square value, i.e. the root of the arithmetic 172 * mean of the square of all weighted residuals. This is related to the 173 * criterion that is minimized by the estimator as follows: if 174 * <em>c</em> if the criterion, and <em>n</em> is the number of 175 * measurements, then the RMS is <em>sqrt (c/n)</em>. 176 * 177 * @param problem estimation problem 178 * @return RMS value 179 */ 180 public double getRMS(EstimationProblem problem) { 181 WeightedMeasurement[] wm = problem.getMeasurements(); 182 double criterion = 0; 183 for (int i = 0; i < wm.length; ++i) { 184 double residual = wm[i].getResidual(); 185 criterion += wm[i].getWeight() * residual * residual; 186 } 187 return Math.sqrt(criterion / wm.length); 188 } 189 190 /** 191 * Get the Chi-Square value. 192 * @param problem estimation problem 193 * @return chi-square value 194 */ 195 public double getChiSquare(EstimationProblem problem) { 196 WeightedMeasurement[] wm = problem.getMeasurements(); 197 double chiSquare = 0; 198 for (int i = 0; i < wm.length; ++i) { 199 double residual = wm[i].getResidual(); 200 chiSquare += residual * residual / wm[i].getWeight(); 201 } 202 return chiSquare; 203 } 204 205 /** 206 * Get the covariance matrix of unbound estimated parameters. 207 * @param problem estimation problem 208 * @return covariance matrix 209 * @exception EstimationException if the covariance matrix 210 * cannot be computed (singular problem) 211 */ 212 public double[][] getCovariances(EstimationProblem problem) 213 throws EstimationException { 214 215 // set up the jacobian 216 updateJacobian(); 217 218 // compute transpose(J).J, avoiding building big intermediate matrices 219 final int n = problem.getMeasurements().length; 220 final int m = problem.getUnboundParameters().length; 221 final int max = m * n; 222 double[][] jTj = new double[m][m]; 223 for (int i = 0; i < m; ++i) { 224 for (int j = i; j < m; ++j) { 225 double sum = 0; 226 for (int k = 0; k < max; k += m) { 227 sum += jacobian[k + i] * jacobian[k + j]; 228 } 229 jTj[i][j] = sum; 230 jTj[j][i] = sum; 231 } 232 } 233 234 try { 235 // compute the covariances matrix 236 RealMatrix inverse = 237 new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse(); 238 return inverse.getData(); 239 } catch (InvalidMatrixException ime) { 240 throw new EstimationException("unable to compute covariances: singular problem"); 241 } 242 243 } 244 245 /** 246 * Guess the errors in unbound estimated parameters. 247 * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p> 248 * @param problem estimation problem 249 * @return errors in estimated parameters 250 * @exception EstimationException if the covariances matrix cannot be computed 251 * or the number of degrees of freedom is not positive (number of measurements 252 * lesser or equal to number of parameters) 253 */ 254 public double[] guessParametersErrors(EstimationProblem problem) 255 throws EstimationException { 256 int m = problem.getMeasurements().length; 257 int p = problem.getUnboundParameters().length; 258 if (m <= p) { 259 throw new EstimationException( 260 "no degrees of freedom ({0} measurements, {1} parameters)", 261 m, p); 262 } 263 double[] errors = new double[problem.getUnboundParameters().length]; 264 final double c = Math.sqrt(getChiSquare(problem) / (m - p)); 265 double[][] covar = getCovariances(problem); 266 for (int i = 0; i < errors.length; ++i) { 267 errors[i] = Math.sqrt(covar[i][i]) * c; 268 } 269 return errors; 270 } 271 272 /** 273 * Initialization of the common parts of the estimation. 274 * <p>This method <em>must</em> be called at the start 275 * of the {@link #estimate(EstimationProblem) estimate} 276 * method.</p> 277 * @param problem estimation problem to solve 278 */ 279 protected void initializeEstimate(EstimationProblem problem) { 280 281 // reset counters 282 costEvaluations = 0; 283 jacobianEvaluations = 0; 284 285 // retrieve the equations and the parameters 286 measurements = problem.getMeasurements(); 287 parameters = problem.getUnboundParameters(); 288 289 // arrays shared with the other private methods 290 rows = measurements.length; 291 cols = parameters.length; 292 jacobian = new double[rows * cols]; 293 residuals = new double[rows]; 294 295 cost = Double.POSITIVE_INFINITY; 296 297 } 298 299 /** 300 * Solve an estimation problem. 301 * 302 * <p>The method should set the parameters of the problem to several 303 * trial values until it reaches convergence. If this method returns 304 * normally (i.e. without throwing an exception), then the best 305 * estimate of the parameters can be retrieved from the problem 306 * itself, through the {@link EstimationProblem#getAllParameters 307 * EstimationProblem.getAllParameters} method.</p> 308 * 309 * @param problem estimation problem to solve 310 * @exception EstimationException if the problem cannot be solved 311 * 312 */ 313 public abstract void estimate(EstimationProblem problem) 314 throws EstimationException; 315 316 }