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.RealMatrixImpl;
24
25
26
27
28
29
30
31
32
33 public abstract class AbstractEstimator implements Estimator {
34
35
36
37
38 protected AbstractEstimator() {
39 }
40
41
42
43
44
45
46
47 public final void setMaxCostEval(int maxCostEval) {
48 this.maxCostEval = maxCostEval;
49 }
50
51
52
53
54
55
56 public final int getCostEvaluations() {
57 return costEvaluations;
58 }
59
60
61
62
63
64
65 public final int getJacobianEvaluations() {
66 return jacobianEvaluations;
67 }
68
69
70
71
72 protected void updateJacobian() {
73 incrementJacobianEvaluationsCounter();
74 Arrays.fill(jacobian, 0);
75 for (int i = 0, index = 0; i < rows; i++) {
76 WeightedMeasurement wm = measurements[i];
77 double factor = -Math.sqrt(wm.getWeight());
78 for (int j = 0; j < cols; ++j) {
79 jacobian[index++] = factor * wm.getPartial(parameters[j]);
80 }
81 }
82 }
83
84
85
86
87 protected final void incrementJacobianEvaluationsCounter() {
88 ++jacobianEvaluations;
89 }
90
91
92
93
94
95
96 protected void updateResidualsAndCost()
97 throws EstimationException {
98
99 if (++costEvaluations > maxCostEval) {
100 throw new EstimationException("maximal number of evaluations exceeded ({0})",
101 new Object[] { new Integer(maxCostEval) });
102 }
103
104 cost = 0;
105 for (int i = 0, index = 0; i < rows; i++, index += cols) {
106 WeightedMeasurement wm = measurements[i];
107 double residual = wm.getResidual();
108 residuals[i] = Math.sqrt(wm.getWeight()) * residual;
109 cost += wm.getWeight() * residual * residual;
110 }
111 cost = Math.sqrt(cost);
112
113 }
114
115
116
117
118
119
120
121
122
123
124
125
126 public double getRMS(EstimationProblem problem) {
127 WeightedMeasurement[] wm = problem.getMeasurements();
128 double criterion = 0;
129 for (int i = 0; i < wm.length; ++i) {
130 double residual = wm[i].getResidual();
131 criterion += wm[i].getWeight() * residual * residual;
132 }
133 return Math.sqrt(criterion / wm.length);
134 }
135
136
137
138
139
140
141 public double getChiSquare(EstimationProblem problem) {
142 WeightedMeasurement[] wm = problem.getMeasurements();
143 double chiSquare = 0;
144 for (int i = 0; i < wm.length; ++i) {
145 double residual = wm[i].getResidual();
146 chiSquare += residual * residual / wm[i].getWeight();
147 }
148 return chiSquare;
149 }
150
151
152
153
154
155
156
157
158 public double[][] getCovariances(EstimationProblem problem)
159 throws EstimationException {
160
161
162 updateJacobian();
163
164
165 final int rows = problem.getMeasurements().length;
166 final int cols = problem.getAllParameters().length;
167 final int max = cols * rows;
168 double[][] jTj = new double[cols][cols];
169 for (int i = 0; i < cols; ++i) {
170 for (int j = i; j < cols; ++j) {
171 double sum = 0;
172 for (int k = 0; k < max; k += cols) {
173 sum += jacobian[k + i] * jacobian[k + j];
174 }
175 jTj[i][j] = sum;
176 jTj[j][i] = sum;
177 }
178 }
179
180 try {
181
182 return new RealMatrixImpl(jTj).inverse().getData();
183 } catch (InvalidMatrixException ime) {
184 throw new EstimationException("unable to compute covariances: singular problem",
185 new Object[0]);
186 }
187
188 }
189
190
191
192
193
194
195
196
197
198
199 public double[] guessParametersErrors(EstimationProblem problem)
200 throws EstimationException {
201 int m = problem.getMeasurements().length;
202 int p = problem.getAllParameters().length;
203 if (m <= p) {
204 throw new EstimationException("no degrees of freedom ({0} measurements, {1} parameters)",
205 new Object[] { new Integer(m), new Integer(p)});
206 }
207 double[] errors = new double[problem.getAllParameters().length];
208 final double c = Math.sqrt(getChiSquare(problem) / (m - p));
209 double[][] covar = getCovariances(problem);
210 for (int i = 0; i < errors.length; ++i) {
211 errors[i] = Math.sqrt(covar[i][i]) * c;
212 }
213 return errors;
214 }
215
216
217
218
219
220
221
222
223 protected void initializeEstimate(EstimationProblem problem) {
224
225
226 costEvaluations = 0;
227 jacobianEvaluations = 0;
228
229
230 measurements = problem.getMeasurements();
231 parameters = problem.getUnboundParameters();
232
233
234 rows = measurements.length;
235 cols = parameters.length;
236 jacobian = new double[rows * cols];
237 residuals = new double[rows];
238
239 cost = Double.POSITIVE_INFINITY;
240
241 }
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257 public abstract void estimate(EstimationProblem problem)
258 throws EstimationException;
259
260
261 protected WeightedMeasurement[] measurements;
262
263
264 protected EstimatedParameter[] parameters;
265
266
267
268
269
270
271
272
273 protected double[] jacobian;
274
275
276 protected int cols;
277
278
279 protected int rows;
280
281
282
283
284
285
286
287 protected double[] residuals;
288
289
290 protected double cost;
291
292
293 private int maxCostEval;
294
295
296 private int costEvaluations;
297
298
299 private int jacobianEvaluations;
300
301 }