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 }