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    }