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