View Javadoc

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.linear.InvalidMatrixException;
23  import org.apache.commons.math.linear.LUDecompositionImpl;
24  import org.apache.commons.math.linear.MatrixUtils;
25  import org.apache.commons.math.linear.RealMatrix;
26  
27  /**
28   * Base class for implementing estimators.
29   * <p>This base class handles the boilerplates methods associated to thresholds
30   * settings, jacobian and error estimation.</p>
31   * @version $Revision: 781122 $ $Date: 2009-06-02 14:53:23 -0400 (Tue, 02 Jun 2009) $
32   * @since 1.2
33   * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
34   * been deprecated and replaced by package org.apache.commons.math.optimization.general
35   *
36   */
37  @Deprecated
38  public abstract class AbstractEstimator implements Estimator {
39  
40      /** Default maximal number of cost evaluations allowed. */
41      public static final int DEFAULT_MAX_COST_EVALUATIONS = 100;
42  
43      /**
44       * Build an abstract estimator for least squares problems.
45       * <p>The maximal number of cost evaluations allowed is set
46       * to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.</p>
47       */
48      protected AbstractEstimator() {
49          setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS);
50      }
51  
52      /**
53       * Set the maximal number of cost evaluations allowed.
54       * 
55       * @param maxCostEval maximal number of cost evaluations allowed
56       * @see #estimate
57       */
58      public final void setMaxCostEval(int maxCostEval) {
59          this.maxCostEval = maxCostEval;
60      }
61  
62      /**
63       * Get the number of cost evaluations.
64       * 
65       * @return number of cost evaluations
66       * */
67      public final int getCostEvaluations() {
68          return costEvaluations;
69      }
70  
71      /** 
72       * Get the number of jacobian evaluations.
73       * 
74       * @return number of jacobian evaluations
75       * */
76      public final int getJacobianEvaluations() {
77          return jacobianEvaluations;
78      }
79  
80      /** 
81       * Update the jacobian matrix.
82       */
83      protected void updateJacobian() {
84          incrementJacobianEvaluationsCounter();
85          Arrays.fill(jacobian, 0);
86          for (int i = 0, index = 0; i < rows; i++) {
87              WeightedMeasurement wm = measurements[i];
88              double factor = -Math.sqrt(wm.getWeight());
89              for (int j = 0; j < cols; ++j) {
90                  jacobian[index++] = factor * wm.getPartial(parameters[j]);
91              }
92          }
93      }
94  
95      /**
96       * Increment the jacobian evaluations counter.
97       */
98      protected final void incrementJacobianEvaluationsCounter() {
99        ++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 }