1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.estimation; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.Arrays; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.LocalizedFormats; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.linear.InvalidMatrixException; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.linear.LUDecompositionImpl; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.linear.MatrixUtils; 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.linear.RealMatrix; 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Base class for implementing estimators. 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>This base class handles the boilerplates methods associated to thresholds 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * settings, jacobian and error estimation.</p> 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 1.2 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * been deprecated and replaced by package org.apache.commons.math.optimization.general 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond@Deprecated 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic abstract class AbstractEstimator implements Estimator { 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Default maximal number of cost evaluations allowed. */ 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public static final int DEFAULT_MAX_COST_EVALUATIONS = 100; 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Array of measurements. */ 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected WeightedMeasurement[] measurements; 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Array of parameters. */ 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected EstimatedParameter[] parameters; 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Jacobian matrix. 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>This matrix is in canonical form just after the calls to 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #updateJacobian()}, but may be modified by the solver 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * in the derived class (the {@link LevenbergMarquardtEstimator 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Levenberg-Marquardt estimator} does this).</p> 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double[] jacobian; 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Number of columns of the jacobian matrix. */ 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected int cols; 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Number of rows of the jacobian matrix. */ 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected int rows; 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Residuals array. 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>This array is in canonical form just after the calls to 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #updateJacobian()}, but may be modified by the solver 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * in the derived class (the {@link LevenbergMarquardtEstimator 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Levenberg-Marquardt estimator} does this).</p> 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double[] residuals; 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cost value (square root of the sum of the residuals). */ 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double cost; 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Maximal allowed number of cost evaluations. */ 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private int maxCostEval; 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Number of cost evaluations. */ 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private int costEvaluations; 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Number of jacobian evaluations. */ 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private int jacobianEvaluations; 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Build an abstract estimator for least squares problems. 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The maximal number of cost evaluations allowed is set 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.</p> 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected AbstractEstimator() { 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS); 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Set the maximal number of cost evaluations allowed. 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param maxCostEval maximal number of cost evaluations allowed 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see #estimate 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public final void setMaxCostEval(int maxCostEval) { 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.maxCostEval = maxCostEval; 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Get the number of cost evaluations. 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return number of cost evaluations 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * */ 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public final int getCostEvaluations() { 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return costEvaluations; 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Get the number of jacobian evaluations. 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return number of jacobian evaluations 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * */ 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public final int getJacobianEvaluations() { 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return jacobianEvaluations; 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Update the jacobian matrix. 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected void updateJacobian() { 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond incrementJacobianEvaluationsCounter(); 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond Arrays.fill(jacobian, 0); 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int index = 0; 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < rows; i++) { 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond WeightedMeasurement wm = measurements[i]; 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double factor = -FastMath.sqrt(wm.getWeight()); 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < cols; ++j) { 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond jacobian[index++] = factor * wm.getPartial(parameters[j]); 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Increment the jacobian evaluations counter. 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected final void incrementJacobianEvaluationsCounter() { 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ++jacobianEvaluations; 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Update the residuals array and cost function value. 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception EstimationException if the number of cost evaluations 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * exceeds the maximum allowed 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected void updateResidualsAndCost() 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws EstimationException { 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (++costEvaluations > maxCostEval) { 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new EstimationException(LocalizedFormats.MAX_EVALUATIONS_EXCEEDED, 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond maxCostEval); 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cost = 0; 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int index = 0; 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < rows; i++, index += cols) { 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond WeightedMeasurement wm = measurements[i]; 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double residual = wm.getResidual(); 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond residuals[i] = FastMath.sqrt(wm.getWeight()) * residual; 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cost += wm.getWeight() * residual * residual; 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cost = FastMath.sqrt(cost); 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Get the Root Mean Square value. 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Get the Root Mean Square value, i.e. the root of the arithmetic 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * mean of the square of all weighted residuals. This is related to the 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * criterion that is minimized by the estimator as follows: if 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <em>c</em> if the criterion, and <em>n</em> is the number of 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * measurements, then the RMS is <em>sqrt (c/n)</em>. 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param problem estimation problem 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return RMS value 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getRMS(EstimationProblem problem) { 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond WeightedMeasurement[] wm = problem.getMeasurements(); 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double criterion = 0; 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < wm.length; ++i) { 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double residual = wm[i].getResidual(); 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond criterion += wm[i].getWeight() * residual * residual; 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return FastMath.sqrt(criterion / wm.length); 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Get the Chi-Square value. 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param problem estimation problem 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return chi-square value 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getChiSquare(EstimationProblem problem) { 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond WeightedMeasurement[] wm = problem.getMeasurements(); 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double chiSquare = 0; 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < wm.length; ++i) { 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double residual = wm[i].getResidual(); 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond chiSquare += residual * residual / wm[i].getWeight(); 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return chiSquare; 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Get the covariance matrix of unbound estimated parameters. 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param problem estimation problem 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return covariance matrix 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception EstimationException if the covariance matrix 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * cannot be computed (singular problem) 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[][] getCovariances(EstimationProblem problem) 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws EstimationException { 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // set up the jacobian 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond updateJacobian(); 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute transpose(J).J, avoiding building big intermediate matrices 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = problem.getMeasurements().length; 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int m = problem.getUnboundParameters().length; 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int max = m * n; 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[][] jTj = new double[m][m]; 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < m; ++i) { 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = i; j < m; ++j) { 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sum = 0; 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < max; k += m) { 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sum += jacobian[k + i] * jacobian[k + j]; 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond jTj[i][j] = sum; 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond jTj[j][i] = sum; 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond try { 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute the covariances matrix 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealMatrix inverse = 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse(); 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return inverse.getData(); 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (InvalidMatrixException ime) { 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new EstimationException(LocalizedFormats.UNABLE_TO_COMPUTE_COVARIANCE_SINGULAR_PROBLEM); 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Guess the errors in unbound estimated parameters. 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p> 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param problem estimation problem 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return errors in estimated parameters 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception EstimationException if the covariances matrix cannot be computed 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * or the number of degrees of freedom is not positive (number of measurements 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * lesser or equal to number of parameters) 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] guessParametersErrors(EstimationProblem problem) 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws EstimationException { 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int m = problem.getMeasurements().length; 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int p = problem.getUnboundParameters().length; 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (m <= p) { 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new EstimationException( 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.NO_DEGREES_OF_FREEDOM, 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond m, p); 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[] errors = new double[problem.getUnboundParameters().length]; 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double c = FastMath.sqrt(getChiSquare(problem) / (m - p)); 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[][] covar = getCovariances(problem); 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < errors.length; ++i) { 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond errors[i] = FastMath.sqrt(covar[i][i]) * c; 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return errors; 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Initialization of the common parts of the estimation. 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>This method <em>must</em> be called at the start 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * of the {@link #estimate(EstimationProblem) estimate} 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * method.</p> 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param problem estimation problem to solve 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected void initializeEstimate(EstimationProblem problem) { 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // reset counters 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond costEvaluations = 0; 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond jacobianEvaluations = 0; 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // retrieve the equations and the parameters 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond measurements = problem.getMeasurements(); 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond parameters = problem.getUnboundParameters(); 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // arrays shared with the other private methods 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond rows = measurements.length; 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cols = parameters.length; 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond jacobian = new double[rows * cols]; 295dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond residuals = new double[rows]; 296dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 297dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cost = Double.POSITIVE_INFINITY; 298dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 299dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 300dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 301dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 302dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Solve an estimation problem. 303dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 304dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The method should set the parameters of the problem to several 305dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * trial values until it reaches convergence. If this method returns 306dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * normally (i.e. without throwing an exception), then the best 307dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * estimate of the parameters can be retrieved from the problem 308dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * itself, through the {@link EstimationProblem#getAllParameters 309dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * EstimationProblem.getAllParameters} method.</p> 310dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 311dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param problem estimation problem to solve 312dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception EstimationException if the problem cannot be solved 313dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 314dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 315dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public abstract void estimate(EstimationProblem problem) 316dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws EstimationException; 317dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 318dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 319