1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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
29
30
31
32
33
34
35
36
37 @Deprecated
38 public abstract class AbstractEstimator implements Estimator {
39
40
41 public static final int DEFAULT_MAX_COST_EVALUATIONS = 100;
42
43
44
45
46
47
48 protected AbstractEstimator() {
49 setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS);
50 }
51
52
53
54
55
56
57
58 public final void setMaxCostEval(int maxCostEval) {
59 this.maxCostEval = maxCostEval;
60 }
61
62
63
64
65
66
67 public final int getCostEvaluations() {
68 return costEvaluations;
69 }
70
71
72
73
74
75
76 public final int getJacobianEvaluations() {
77 return jacobianEvaluations;
78 }
79
80
81
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
97
98 protected final void incrementJacobianEvaluationsCounter() {
99 ++jacobianEvaluations;
100 }
101
102
103
104
105
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
128
129
130
131
132
133
134
135
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
149
150
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
164
165
166
167
168
169 public double[][] getCovariances(EstimationProblem problem)
170 throws EstimationException {
171
172
173 updateJacobian();
174
175
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
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
204
205
206
207
208
209
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
231
232
233
234
235
236 protected void initializeEstimate(EstimationProblem problem) {
237
238
239 costEvaluations = 0;
240 jacobianEvaluations = 0;
241
242
243 measurements = problem.getMeasurements();
244 parameters = problem.getUnboundParameters();
245
246
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
258
259
260
261
262
263
264
265
266
267
268
269
270 public abstract void estimate(EstimationProblem problem)
271 throws EstimationException;
272
273
274 protected WeightedMeasurement[] measurements;
275
276
277 protected EstimatedParameter[] parameters;
278
279
280
281
282
283
284
285
286 protected double[] jacobian;
287
288
289 protected int cols;
290
291
292 protected int rows;
293
294
295
296
297
298
299
300 protected double[] residuals;
301
302
303 protected double cost;
304
305
306 private int maxCostEval;
307
308
309 private int costEvaluations;
310
311
312 private int jacobianEvaluations;
313
314 }