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 }