1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.legacy.stat.regression;
18
19 import org.apache.commons.math4.legacy.TestUtils;
20 import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
21 import org.apache.commons.math4.legacy.linear.RealMatrix;
22 import org.apache.commons.math4.legacy.stat.correlation.PearsonsCorrelation;
23 import org.apache.commons.math4.core.jdkmath.JdkMath;
24 import org.junit.Assert;
25 import org.junit.Test;
26
27
28
29
30 public class MillerUpdatingRegressionTest {
31
32 public MillerUpdatingRegressionTest() {
33 }
34
35
36
37 private static final double[][] airdata = {
38 new double[]{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6},
39 new double[]{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15},
40 new double[]{1140640, 1215690, 1309570, 1511530, 1676730, 1823740, 2022890, 2314760, 2639160, 3247620, 3787750, 3867750, 3996020, 4282880, 4748320, 569292, 640614, 777655, 999294, 1203970, 1358100, 1501350, 1709270, 2025400, 2548370, 3137740, 3557700, 3717740, 3962370, 4209390, 286298, 309290, 342056, 374595, 450037, 510412, 575347, 669331, 783799, 913883, 1041520, 1125800, 1096070, 1198930, 1170470, 145167, 170192, 247506, 309391, 354338, 373941, 420915, 474017, 532590, 676771, 880438, 1052020, 1193680, 1303390, 1436970, 91361, 95428, 98187, 115967, 138382, 156228, 183169, 210212, 274024, 356915, 432344, 524294, 530924, 581447, 610257, 68978, 74904, 83829, 98148, 118449, 133161, 145062, 170711, 199775, 276797, 381478, 506969, 633388, 804388, 1009500},
41 new double[]{0.952757, 0.986757, 1.09198, 1.17578, 1.16017, 1.17376, 1.29051, 1.39067, 1.61273, 1.82544, 1.54604, 1.5279, 1.6602, 1.82231, 1.93646, 0.520635, 0.534627, 0.655192, 0.791575, 0.842945, 0.852892, 0.922843, 1, 1.19845, 1.34067, 1.32624, 1.24852, 1.25432, 1.37177, 1.38974, 0.262424, 0.266433, 0.306043, 0.325586, 0.345706, 0.367517, 0.409937, 0.448023, 0.539595, 0.539382, 0.467967, 0.450544, 0.468793, 0.494397, 0.493317, 0.086393, 0.09674, 0.1415, 0.169715, 0.173805, 0.164272, 0.170906, 0.17784, 0.192248, 0.242469, 0.256505, 0.249657, 0.273923, 0.371131, 0.421411, 0.051028, 0.052646, 0.056348, 0.066953, 0.070308, 0.073961, 0.084946, 0.095474, 0.119814, 0.150046, 0.144014, 0.1693, 0.172761, 0.18667, 0.213279, 0.037682, 0.039784, 0.044331, 0.050245, 0.055046, 0.052462, 0.056977, 0.06149, 0.069027, 0.092749, 0.11264, 0.154154, 0.186461, 0.246847, 0.304013},
42 new double[]{106650, 110307, 110574, 121974, 196606, 265609, 263451, 316411, 384110, 569251, 871636, 997239, 938002, 859572, 823411, 103795, 111477, 118664, 114797, 215322, 281704, 304818, 348609, 374579, 544109, 853356, 1003200, 941977, 856533, 821361, 118788, 123798, 122882, 131274, 222037, 278721, 306564, 356073, 378311, 555267, 850322, 1015610, 954508, 886999, 844079, 114987, 120501, 121908, 127220, 209405, 263148, 316724, 363598, 389436, 547376, 850418, 1011170, 951934, 881323, 831374, 118222, 116223, 115853, 129372, 243266, 277930, 317273, 358794, 397667, 566672, 848393, 1005740, 958231, 872924, 844622, 117112, 119420, 116087, 122997, 194309, 307923, 323595, 363081, 386422, 564867, 874818, 1013170, 930477, 851676, 819476},
43 new double[]{0.534487, 0.532328, 0.547736, 0.540846, 0.591167, 0.575417, 0.594495, 0.597409, 0.638522, 0.676287, 0.605735, 0.61436, 0.633366, 0.650117, 0.625603, 0.490851, 0.473449, 0.503013, 0.512501, 0.566782, 0.558133, 0.558799, 0.57207, 0.624763, 0.628706, 0.58915, 0.532612, 0.526652, 0.540163, 0.528775, 0.524334, 0.537185, 0.582119, 0.579489, 0.606592, 0.60727, 0.582425, 0.573972, 0.654256, 0.631055, 0.56924, 0.589682, 0.587953, 0.565388, 0.577078, 0.432066, 0.439669, 0.488932, 0.484181, 0.529925, 0.532723, 0.549067, 0.55714, 0.611377, 0.645319, 0.611734, 0.580884, 0.572047, 0.59457, 0.585525, 0.442875, 0.462473, 0.519118, 0.529331, 0.557797, 0.556181, 0.569327, 0.583465, 0.631818, 0.604723, 0.587921, 0.616159, 0.605868, 0.594688, 0.635545, 0.448539, 0.475889, 0.500562, 0.500344, 0.528897, 0.495361, 0.510342, 0.518296, 0.546723, 0.554276, 0.517766, 0.580049, 0.556024, 0.537791, 0.525775}
44 };
45
46
47
48
49 @Test
50 public void testHasIntercept() {
51 MillerUpdatingRegression instance = new MillerUpdatingRegression(10, false);
52 if (instance.hasIntercept()) {
53 Assert.fail("Should not have intercept");
54 }
55 instance = new MillerUpdatingRegression(10, true);
56 if (!instance.hasIntercept()) {
57 Assert.fail("Should have intercept");
58 }
59 }
60
61
62
63
64 @Test
65 public void testAddObsGetNClear() {
66 MillerUpdatingRegression instance = new MillerUpdatingRegression(3, true);
67 double[][] xAll = new double[airdata[0].length][];
68 double[] y = new double[airdata[0].length];
69 for (int i = 0; i < airdata[0].length; i++) {
70 xAll[i] = new double[3];
71 xAll[i][0] = JdkMath.log(airdata[3][i]);
72 xAll[i][1] = JdkMath.log(airdata[4][i]);
73 xAll[i][2] = airdata[5][i];
74 y[i] = JdkMath.log(airdata[2][i]);
75 }
76 instance.addObservations(xAll, y);
77 if (instance.getN() != xAll.length) {
78 Assert.fail("Number of observations not correct in bulk addition");
79 }
80 instance.clear();
81 for (int i = 0; i < xAll.length; i++) {
82 instance.addObservation(xAll[i], y[i]);
83 }
84 if (instance.getN() != xAll.length) {
85 Assert.fail("Number of observations not correct in drip addition");
86 }
87 return;
88 }
89
90 @Test
91 public void testNegativeTestAddObs() {
92 MillerUpdatingRegression instance = new MillerUpdatingRegression(3, true);
93 try {
94 instance.addObservation(new double[]{1.0}, 0.0);
95 Assert.fail("Should throw MathIllegalArgumentException");
96 } catch (MathIllegalArgumentException iae) {
97 } catch (Exception e) {
98 Assert.fail("Should throw MathIllegalArgumentException");
99 }
100 try {
101 instance.addObservation(new double[]{1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}, 0.0);
102 Assert.fail("Should throw MathIllegalArgumentException");
103 } catch (MathIllegalArgumentException iae) {
104 } catch (Exception e) {
105 Assert.fail("Should throw MathIllegalArgumentException");
106 }
107 try {
108 instance.addObservation(new double[]{1.0, 1.0, 1.0}, 0.0);
109 } catch (Exception e) {
110 Assert.fail("Should throw MathIllegalArgumentException");
111 }
112
113
114 instance = new MillerUpdatingRegression(3, false);
115 try {
116 instance.addObservation(new double[]{1.0}, 0.0);
117 Assert.fail("Should throw MathIllegalArgumentException [NOINTERCEPT]");
118 } catch (MathIllegalArgumentException iae) {
119 } catch (Exception e) {
120 Assert.fail("Should throw MathIllegalArgumentException [NOINTERCEPT]");
121 }
122 try {
123 instance.addObservation(new double[]{1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}, 0.0);
124 Assert.fail("Should throw MathIllegalArgumentException [NOINTERCEPT]");
125 } catch (MathIllegalArgumentException iae) {
126 } catch (Exception e) {
127 Assert.fail("Should throw MathIllegalArgumentException [NOINTERCEPT]");
128 }
129 try {
130 instance.addObservation(new double[]{1.0, 1.0, 1.0}, 0.0);
131 } catch (Exception e) {
132 Assert.fail("Should throw MathIllegalArgumentException [NOINTERCEPT]");
133 }
134 }
135
136 @Test
137 public void testNegativeTestAddMultipleObs() {
138 MillerUpdatingRegression instance = new MillerUpdatingRegression(3, true);
139 try {
140 double[][] tst = {{1.0, 1.0, 1.0}, {1.20, 2.0, 2.1}};
141 double[] y = {1.0};
142 instance.addObservations(tst, y);
143
144 Assert.fail("Should throw MathIllegalArgumentException");
145 } catch (MathIllegalArgumentException iae) {
146 } catch (Exception e) {
147 Assert.fail("Should throw MathIllegalArgumentException");
148 }
149
150 try {
151 double[][] tst = {{1.0, 1.0, 1.0}, {1.20, 2.0, 2.1}};
152 double[] y = {1.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
153 instance.addObservations(tst, y);
154
155 Assert.fail("Should throw MathIllegalArgumentException");
156 } catch (MathIllegalArgumentException iae) {
157 } catch (Exception e) {
158 Assert.fail("Should throw MathIllegalArgumentException");
159 }
160 }
161
162
163
164
165 @Test
166 public void testRegressAirlineConstantExternal() {
167 MillerUpdatingRegression instance = new MillerUpdatingRegression(4, false);
168 double[][] x = new double[airdata[0].length][];
169 double[] y = new double[airdata[0].length];
170 for (int i = 0; i < airdata[0].length; i++) {
171 x[i] = new double[4];
172 x[i][0] = 1.0;
173 x[i][1] = JdkMath.log(airdata[3][i]);
174 x[i][2] = JdkMath.log(airdata[4][i]);
175 x[i][3] = airdata[5][i];
176 y[i] = JdkMath.log(airdata[2][i]);
177 }
178
179 instance.addObservations(x, y);
180 try {
181 RegressionResults result = instance.regress();
182 Assert.assertNotNull("The test case is a prototype.", result);
183 TestUtils.assertEquals(
184 new double[]{9.5169, 0.8827, 0.4540, -1.6275},
185 result.getParameterEstimates(), 1e-4);
186
187
188 TestUtils.assertEquals(
189 new double[]{.2292445, .0132545, .0203042, .345302},
190 result.getStdErrorOfEstimates(), 1.0e-4);
191
192 TestUtils.assertEquals(0.01552839, result.getMeanSquareError(), 1.0e-8);
193 } catch (Exception e) {
194 Assert.fail("Should not throw exception but does");
195 }
196 }
197
198 @Test
199 public void testRegressAirlineConstantInternal() {
200 MillerUpdatingRegression instance = new MillerUpdatingRegression(3, true);
201 double[][] x = new double[airdata[0].length][];
202 double[] y = new double[airdata[0].length];
203 for (int i = 0; i < airdata[0].length; i++) {
204 x[i] = new double[3];
205 x[i][0] = JdkMath.log(airdata[3][i]);
206 x[i][1] = JdkMath.log(airdata[4][i]);
207 x[i][2] = airdata[5][i];
208 y[i] = JdkMath.log(airdata[2][i]);
209 }
210
211 instance.addObservations(x, y);
212 try {
213 RegressionResults result = instance.regress();
214 Assert.assertNotNull("The test case is a prototype.", result);
215 TestUtils.assertEquals(
216 new double[]{9.5169, 0.8827, 0.4540, -1.6275},
217 result.getParameterEstimates(), 1e-4);
218
219
220 TestUtils.assertEquals(
221 new double[]{.2292445, .0132545, .0203042, .345302},
222 result.getStdErrorOfEstimates(), 1.0e-4);
223
224 TestUtils.assertEquals(0.9883, result.getRSquared(), 1.0e-4);
225 TestUtils.assertEquals(0.01552839, result.getMeanSquareError(), 1.0e-8);
226 } catch (Exception e) {
227 Assert.fail("Should not throw exception but does");
228 }
229 }
230
231 @Test
232 public void testFilippelli() {
233 double[] data = new double[]{
234 0.8116, -6.860120914,
235 0.9072, -4.324130045,
236 0.9052, -4.358625055,
237 0.9039, -4.358426747,
238 0.8053, -6.955852379,
239 0.8377, -6.661145254,
240 0.8667, -6.355462942,
241 0.8809, -6.118102026,
242 0.7975, -7.115148017,
243 0.8162, -6.815308569,
244 0.8515, -6.519993057,
245 0.8766, -6.204119983,
246 0.8885, -5.853871964,
247 0.8859, -6.109523091,
248 0.8959, -5.79832982,
249 0.8913, -5.482672118,
250 0.8959, -5.171791386,
251 0.8971, -4.851705903,
252 0.9021, -4.517126416,
253 0.909, -4.143573228,
254 0.9139, -3.709075441,
255 0.9199, -3.499489089,
256 0.8692, -6.300769497,
257 0.8872, -5.953504836,
258 0.89, -5.642065153,
259 0.891, -5.031376979,
260 0.8977, -4.680685696,
261 0.9035, -4.329846955,
262 0.9078, -3.928486195,
263 0.7675, -8.56735134,
264 0.7705, -8.363211311,
265 0.7713, -8.107682739,
266 0.7736, -7.823908741,
267 0.7775, -7.522878745,
268 0.7841, -7.218819279,
269 0.7971, -6.920818754,
270 0.8329, -6.628932138,
271 0.8641, -6.323946875,
272 0.8804, -5.991399828,
273 0.7668, -8.781464495,
274 0.7633, -8.663140179,
275 0.7678, -8.473531488,
276 0.7697, -8.247337057,
277 0.77, -7.971428747,
278 0.7749, -7.676129393,
279 0.7796, -7.352812702,
280 0.7897, -7.072065318,
281 0.8131, -6.774174009,
282 0.8498, -6.478861916,
283 0.8741, -6.159517513,
284 0.8061, -6.835647144,
285 0.846, -6.53165267,
286 0.8751, -6.224098421,
287 0.8856, -5.910094889,
288 0.8919, -5.598599459,
289 0.8934, -5.290645224,
290 0.894, -4.974284616,
291 0.8957, -4.64454848,
292 0.9047, -4.290560426,
293 0.9129, -3.885055584,
294 0.9209, -3.408378962,
295 0.9219, -3.13200249,
296 0.7739, -8.726767166,
297 0.7681, -8.66695597,
298 0.7665, -8.511026475,
299 0.7703, -8.165388579,
300 0.7702, -7.886056648,
301 0.7761, -7.588043762,
302 0.7809, -7.283412422,
303 0.7961, -6.995678626,
304 0.8253, -6.691862621,
305 0.8602, -6.392544977,
306 0.8809, -6.067374056,
307 0.8301, -6.684029655,
308 0.8664, -6.378719832,
309 0.8834, -6.065855188,
310 0.8898, -5.752272167,
311 0.8964, -5.132414673,
312 0.8963, -4.811352704,
313 0.9074, -4.098269308,
314 0.9119, -3.66174277,
315 0.9228, -3.2644011
316 };
317 MillerUpdatingRegression model = new MillerUpdatingRegression(10, true);
318 int off = 0;
319 double[] tmp = new double[10];
320 int nobs = 82;
321 for (int i = 0; i < nobs; i++) {
322 tmp[0] = data[off + 1];
323
324
325
326
327
328
329
330
331
332 tmp[1] = tmp[0] * tmp[0];
333 tmp[2] = tmp[0] * tmp[1];
334 tmp[3] = tmp[0] * tmp[2];
335 tmp[4] = tmp[0] * tmp[3];
336 tmp[5] = tmp[0] * tmp[4];
337 tmp[6] = tmp[0] * tmp[5];
338 tmp[7] = tmp[0] * tmp[6];
339 tmp[8] = tmp[0] * tmp[7];
340 tmp[9] = tmp[0] * tmp[8];
341 model.addObservation(tmp, data[off]);
342 off += 2;
343 }
344 RegressionResults result = model.regress();
345 double[] betaHat = result.getParameterEstimates();
346 TestUtils.assertEquals(betaHat,
347 new double[]{
348 -1467.48961422980,
349 -2772.17959193342,
350 -2316.37108160893,
351 -1127.97394098372,
352 -354.478233703349,
353 -75.1242017393757,
354 -10.8753180355343,
355 -1.06221498588947,
356 -0.670191154593408E-01,
357 -0.246781078275479E-02,
358 -0.402962525080404E-04
359 }, 1E-5);
360
361 double[] se = result.getStdErrorOfEstimates();
362 TestUtils.assertEquals(se,
363 new double[]{
364 298.084530995537,
365 559.779865474950,
366 466.477572127796,
367 227.204274477751,
368 71.6478660875927,
369 15.2897178747400,
370 2.23691159816033,
371 0.221624321934227,
372 0.142363763154724E-01,
373 0.535617408889821E-03,
374 0.896632837373868E-05
375 }, 1E-5);
376
377 TestUtils.assertEquals(0.996727416185620, result.getRSquared(), 1.0e-8);
378 TestUtils.assertEquals(0.112091743968020E-04, result.getMeanSquareError(), 1.0e-10);
379 TestUtils.assertEquals(0.795851382172941E-03, result.getErrorSumSquares(), 1.0e-10);
380 }
381
382 @Test
383 public void testWampler1() {
384 double[] data = new double[]{
385 1, 0,
386 6, 1,
387 63, 2,
388 364, 3,
389 1365, 4,
390 3906, 5,
391 9331, 6,
392 19608, 7,
393 37449, 8,
394 66430, 9,
395 111111, 10,
396 177156, 11,
397 271453, 12,
398 402234, 13,
399 579195, 14,
400 813616, 15,
401 1118481, 16,
402 1508598, 17,
403 2000719, 18,
404 2613660, 19,
405 3368421, 20};
406
407 MillerUpdatingRegression model = new MillerUpdatingRegression(5, true);
408 int off = 0;
409 double[] tmp = new double[5];
410 int nobs = 21;
411 for (int i = 0; i < nobs; i++) {
412 tmp[0] = data[off + 1];
413 tmp[1] = tmp[0] * tmp[0];
414 tmp[2] = tmp[0] * tmp[1];
415 tmp[3] = tmp[0] * tmp[2];
416 tmp[4] = tmp[0] * tmp[3];
417 model.addObservation(tmp, data[off]);
418 off += 2;
419 }
420 RegressionResults result = model.regress();
421 double[] betaHat = result.getParameterEstimates();
422 TestUtils.assertEquals(betaHat,
423 new double[]{1.0,
424 1.0, 1.0,
425 1.0, 1.0,
426 1.0}, 1E-8);
427
428 double[] se = result.getStdErrorOfEstimates();
429 TestUtils.assertEquals(se,
430 new double[]{0.0,
431 0.0, 0.0,
432 0.0, 0.0,
433 0.0}, 1E-8);
434
435 TestUtils.assertEquals(1.0, result.getRSquared(), 1.0e-10);
436 TestUtils.assertEquals(0, result.getMeanSquareError(), 1.0e-7);
437 TestUtils.assertEquals(0.00, result.getErrorSumSquares(), 1.0e-6);
438
439 return;
440 }
441
442 @Test
443 public void testWampler2() {
444 double[] data = new double[]{
445 1.00000, 0,
446 1.11111, 1,
447 1.24992, 2,
448 1.42753, 3,
449 1.65984, 4,
450 1.96875, 5,
451 2.38336, 6,
452 2.94117, 7,
453 3.68928, 8,
454 4.68559, 9,
455 6.00000, 10,
456 7.71561, 11,
457 9.92992, 12,
458 12.75603, 13,
459 16.32384, 14,
460 20.78125, 15,
461 26.29536, 16,
462 33.05367, 17,
463 41.26528, 18,
464 51.16209, 19,
465 63.00000, 20};
466
467 MillerUpdatingRegression model = new MillerUpdatingRegression(5, true);
468 int off = 0;
469 double[] tmp = new double[5];
470 int nobs = 21;
471 for (int i = 0; i < nobs; i++) {
472 tmp[0] = data[off + 1];
473 tmp[1] = tmp[0] * tmp[0];
474 tmp[2] = tmp[0] * tmp[1];
475 tmp[3] = tmp[0] * tmp[2];
476 tmp[4] = tmp[0] * tmp[3];
477 model.addObservation(tmp, data[off]);
478 off += 2;
479 }
480 RegressionResults result = model.regress();
481 double[] betaHat = result.getParameterEstimates();
482 TestUtils.assertEquals(betaHat,
483 new double[]{1.0,
484 1.0e-1, 1.0e-2,
485 1.0e-3, 1.0e-4,
486 1.0e-5}, 1E-8);
487
488 double[] se = result.getStdErrorOfEstimates();
489 TestUtils.assertEquals(se,
490 new double[]{0.0,
491 0.0, 0.0,
492 0.0, 0.0,
493 0.0}, 1E-8);
494
495 TestUtils.assertEquals(1.0, result.getRSquared(), 1.0e-10);
496 TestUtils.assertEquals(0, result.getMeanSquareError(), 1.0e-7);
497 TestUtils.assertEquals(0.00, result.getErrorSumSquares(), 1.0e-6);
498 return;
499 }
500
501 @Test
502 public void testWampler3() {
503 double[] data = new double[]{
504 760, 0,
505 -2042, 1,
506 2111, 2,
507 -1684, 3,
508 3888, 4,
509 1858, 5,
510 11379, 6,
511 17560, 7,
512 39287, 8,
513 64382, 9,
514 113159, 10,
515 175108, 11,
516 273291, 12,
517 400186, 13,
518 581243, 14,
519 811568, 15,
520 1121004, 16,
521 1506550, 17,
522 2002767, 18,
523 2611612, 19,
524 3369180, 20};
525 MillerUpdatingRegression model = new MillerUpdatingRegression(5, true);
526 int off = 0;
527 double[] tmp = new double[5];
528 int nobs = 21;
529 for (int i = 0; i < nobs; i++) {
530 tmp[0] = data[off + 1];
531 tmp[1] = tmp[0] * tmp[0];
532 tmp[2] = tmp[0] * tmp[1];
533 tmp[3] = tmp[0] * tmp[2];
534 tmp[4] = tmp[0] * tmp[3];
535 model.addObservation(tmp, data[off]);
536 off += 2;
537 }
538 RegressionResults result = model.regress();
539 double[] betaHat = result.getParameterEstimates();
540 TestUtils.assertEquals(betaHat,
541 new double[]{1.0,
542 1.0, 1.0,
543 1.0, 1.0,
544 1.0}, 1E-8);
545 double[] se = result.getStdErrorOfEstimates();
546 TestUtils.assertEquals(se,
547 new double[]{2152.32624678170,
548 2363.55173469681, 779.343524331583,
549 101.475507550350, 5.64566512170752,
550 0.112324854679312}, 1E-8);
551
552 TestUtils.assertEquals(.999995559025820, result.getRSquared(), 1.0e-10);
553 TestUtils.assertEquals(5570284.53333333, result.getMeanSquareError(), 1.0e-7);
554 TestUtils.assertEquals(83554268.0000000, result.getErrorSumSquares(), 1.0e-6);
555 return;
556 }
557
558
559 public void testWampler4() {
560 double[] data = new double[]{
561 75901, 0,
562 -204794, 1,
563 204863, 2,
564 -204436, 3,
565 253665, 4,
566 -200894, 5,
567 214131, 6,
568 -185192, 7,
569 221249, 8,
570 -138370, 9,
571 315911, 10,
572 -27644, 11,
573 455253, 12,
574 197434, 13,
575 783995, 14,
576 608816, 15,
577 1370781, 16,
578 1303798, 17,
579 2205519, 18,
580 2408860, 19,
581 3444321, 20};
582 MillerUpdatingRegression model = new MillerUpdatingRegression(5, true);
583 int off = 0;
584 double[] tmp = new double[5];
585 int nobs = 21;
586 for (int i = 0; i < nobs; i++) {
587 tmp[0] = data[off + 1];
588 tmp[1] = tmp[0] * tmp[0];
589 tmp[2] = tmp[0] * tmp[1];
590 tmp[3] = tmp[0] * tmp[2];
591 tmp[4] = tmp[0] * tmp[3];
592 model.addObservation(tmp, data[off]);
593 off += 2;
594 }
595 RegressionResults result = model.regress();
596 double[] betaHat = result.getParameterEstimates();
597 TestUtils.assertEquals(betaHat,
598 new double[]{1.0,
599 1.0, 1.0,
600 1.0, 1.0,
601 1.0}, 1E-8);
602
603 double[] se = result.getStdErrorOfEstimates();
604 TestUtils.assertEquals(se,
605 new double[]{215232.624678170,
606 236355.173469681, 77934.3524331583,
607 10147.5507550350, 564.566512170752,
608 11.2324854679312}, 1E-8);
609
610 TestUtils.assertEquals(.957478440825662, result.getRSquared(), 1.0e-10);
611 TestUtils.assertEquals(55702845333.3333, result.getMeanSquareError(), 1.0e-4);
612 TestUtils.assertEquals(835542680000.000, result.getErrorSumSquares(), 1.0e-3);
613
614 return;
615 }
616
617
618
619
620
621
622
623
624
625
626
627 @Test
628 public void testLongley() {
629
630
631 double[] design = new double[]{
632 60323, 83.0, 234289, 2356, 1590, 107608, 1947,
633 61122, 88.5, 259426, 2325, 1456, 108632, 1948,
634 60171, 88.2, 258054, 3682, 1616, 109773, 1949,
635 61187, 89.5, 284599, 3351, 1650, 110929, 1950,
636 63221, 96.2, 328975, 2099, 3099, 112075, 1951,
637 63639, 98.1, 346999, 1932, 3594, 113270, 1952,
638 64989, 99.0, 365385, 1870, 3547, 115094, 1953,
639 63761, 100.0, 363112, 3578, 3350, 116219, 1954,
640 66019, 101.2, 397469, 2904, 3048, 117388, 1955,
641 67857, 104.6, 419180, 2822, 2857, 118734, 1956,
642 68169, 108.4, 442769, 2936, 2798, 120445, 1957,
643 66513, 110.8, 444546, 4681, 2637, 121950, 1958,
644 68655, 112.6, 482704, 3813, 2552, 123366, 1959,
645 69564, 114.2, 502601, 3931, 2514, 125368, 1960,
646 69331, 115.7, 518173, 4806, 2572, 127852, 1961,
647 70551, 116.9, 554894, 4007, 2827, 130081, 1962
648 };
649
650 final int nobs = 16;
651 final int nvars = 6;
652
653
654 MillerUpdatingRegression model = new MillerUpdatingRegression(6, true);
655 int off = 0;
656 double[] tmp = new double[6];
657 for (int i = 0; i < nobs; i++) {
658 System.arraycopy(design, off + 1, tmp, 0, nvars);
659 model.addObservation(tmp, design[off]);
660 off += nvars + 1;
661 }
662
663
664 RegressionResults result = model.regress();
665 double[] betaHat = result.getParameterEstimates();
666 TestUtils.assertEquals(betaHat,
667 new double[]{-3482258.63459582, 15.0618722713733,
668 -0.358191792925910E-01, -2.02022980381683,
669 -1.03322686717359, -0.511041056535807E-01,
670 1829.15146461355}, 1E-8);
671
672
673 double[] errors = result.getStdErrorOfEstimates();
674 TestUtils.assertEquals(new double[]{890420.383607373,
675 84.9149257747669,
676 0.334910077722432E-01,
677 0.488399681651699,
678 0.214274163161675,
679 0.226073200069370,
680 455.478499142212}, errors, 1E-6);
681
682
683 TestUtils.assertEquals(0.995479004577296, result.getRSquared(), 1E-12);
684 TestUtils.assertEquals(0.992465007628826, result.getAdjustedRSquared(), 1E-12);
685
686
687
688 model = new MillerUpdatingRegression(6, false);
689 off = 0;
690 for (int i = 0; i < nobs; i++) {
691 System.arraycopy(design, off + 1, tmp, 0, nvars);
692 model.addObservation(tmp, design[off]);
693 off += nvars + 1;
694 }
695
696 result = model.regress();
697 betaHat = result.getParameterEstimates();
698 TestUtils.assertEquals(betaHat,
699 new double[]{-52.99357013868291, 0.07107319907358,
700 -0.42346585566399, -0.57256866841929,
701 -0.41420358884978, 48.41786562001326}, 1E-11);
702
703
704 errors = result.getStdErrorOfEstimates();
705 TestUtils.assertEquals(new double[]{129.54486693117232, 0.03016640003786,
706 0.41773654056612, 0.27899087467676, 0.32128496193363,
707 17.68948737819961}, errors, 1E-11);
708
709
710
711 TestUtils.assertEquals(0.9999670130706, result.getRSquared(), 1E-12);
712 TestUtils.assertEquals(0.999947220913, result.getAdjustedRSquared(), 1E-12);
713 }
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755 @Test
756 public void testOneRedundantColumn() {
757 MillerUpdatingRegression instance = new MillerUpdatingRegression(4, false);
758 MillerUpdatingRegression instance2 = new MillerUpdatingRegression(5, false);
759 double[][] x = new double[airdata[0].length][];
760 double[][] x2 = new double[airdata[0].length][];
761 double[] y = new double[airdata[0].length];
762 for (int i = 0; i < airdata[0].length; i++) {
763 x[i] = new double[4];
764 x2[i] = new double[5];
765 x[i][0] = 1.0;
766 x[i][1] = JdkMath.log(airdata[3][i]);
767 x[i][2] = JdkMath.log(airdata[4][i]);
768 x[i][3] = airdata[5][i];
769
770 x2[i][0] = x[i][0];
771 x2[i][1] = x[i][1];
772 x2[i][2] = x[i][2];
773 x2[i][3] = x[i][3];
774 x2[i][4] = x[i][3];
775
776 y[i] = JdkMath.log(airdata[2][i]);
777 }
778
779 instance.addObservations(x, y);
780 RegressionResults result = instance.regress();
781 Assert.assertNotNull("Could not estimate initial regression", result);
782
783 instance2.addObservations(x2, y);
784 RegressionResults resultRedundant = instance2.regress();
785 Assert.assertNotNull("Could not estimate redundant regression", resultRedundant);
786 double[] beta = result.getParameterEstimates();
787 double[] betar = resultRedundant.getParameterEstimates();
788 double[] se = result.getStdErrorOfEstimates();
789 double[] ser = resultRedundant.getStdErrorOfEstimates();
790
791 for (int i = 0; i < beta.length; i++) {
792 if (JdkMath.abs(beta[i] - betar[i]) > 1.0e-8) {
793 Assert.fail("Parameters not correctly estimated");
794 }
795 if (JdkMath.abs(se[i] - ser[i]) > 1.0e-8) {
796 Assert.fail("Standard errors not correctly estimated");
797 }
798 for (int j = 0; j < i; j++) {
799 if (JdkMath.abs(result.getCovarianceOfParameters(i, j)
800 - resultRedundant.getCovarianceOfParameters(i, j)) > 1.0e-8) {
801 Assert.fail("Variance Covariance not correct");
802 }
803 }
804 }
805
806
807 TestUtils.assertEquals(result.getAdjustedRSquared(), resultRedundant.getAdjustedRSquared(), 1.0e-8);
808 TestUtils.assertEquals(result.getErrorSumSquares(), resultRedundant.getErrorSumSquares(), 1.0e-8);
809 TestUtils.assertEquals(result.getMeanSquareError(), resultRedundant.getMeanSquareError(), 1.0e-8);
810 TestUtils.assertEquals(result.getRSquared(), resultRedundant.getRSquared(), 1.0e-8);
811 return;
812 }
813
814 @Test
815 public void testThreeRedundantColumn() {
816
817 MillerUpdatingRegression instance = new MillerUpdatingRegression(4, false);
818 MillerUpdatingRegression instance2 = new MillerUpdatingRegression(7, false);
819 double[][] x = new double[airdata[0].length][];
820 double[][] x2 = new double[airdata[0].length][];
821 double[] y = new double[airdata[0].length];
822 for (int i = 0; i < airdata[0].length; i++) {
823 x[i] = new double[4];
824 x2[i] = new double[7];
825 x[i][0] = 1.0;
826 x[i][1] = JdkMath.log(airdata[3][i]);
827 x[i][2] = JdkMath.log(airdata[4][i]);
828 x[i][3] = airdata[5][i];
829
830 x2[i][0] = x[i][0];
831 x2[i][1] = x[i][0];
832 x2[i][2] = x[i][1];
833 x2[i][3] = x[i][2];
834 x2[i][4] = x[i][1];
835 x2[i][5] = x[i][3];
836 x2[i][6] = x[i][2];
837
838 y[i] = JdkMath.log(airdata[2][i]);
839 }
840
841 instance.addObservations(x, y);
842 RegressionResults result = instance.regress();
843 Assert.assertNotNull("Could not estimate initial regression", result);
844
845 instance2.addObservations(x2, y);
846 RegressionResults resultRedundant = instance2.regress();
847 Assert.assertNotNull("Could not estimate redundant regression", resultRedundant);
848 double[] beta = result.getParameterEstimates();
849 double[] betar = resultRedundant.getParameterEstimates();
850 double[] se = result.getStdErrorOfEstimates();
851 double[] ser = resultRedundant.getStdErrorOfEstimates();
852
853 if (JdkMath.abs(beta[0] - betar[0]) > 1.0e-8) {
854 Assert.fail("Parameters not correct after reorder (0,3)");
855 }
856 if (JdkMath.abs(beta[1] - betar[2]) > 1.0e-8) {
857 Assert.fail("Parameters not correct after reorder (1,2)");
858 }
859 if (JdkMath.abs(beta[2] - betar[3]) > 1.0e-8) {
860 Assert.fail("Parameters not correct after reorder (2,1)");
861 }
862 if (JdkMath.abs(beta[3] - betar[5]) > 1.0e-8) {
863 Assert.fail("Parameters not correct after reorder (3,0)");
864 }
865
866 if (JdkMath.abs(se[0] - ser[0]) > 1.0e-8) {
867 Assert.fail("Se not correct after reorder (0,3)");
868 }
869 if (JdkMath.abs(se[1] - ser[2]) > 1.0e-8) {
870 Assert.fail("Se not correct after reorder (1,2)");
871 }
872 if (JdkMath.abs(se[2] - ser[3]) > 1.0e-8) {
873 Assert.fail("Se not correct after reorder (2,1)");
874 }
875 if (JdkMath.abs(se[3] - ser[5]) > 1.0e-8) {
876 Assert.fail("Se not correct after reorder (3,0)");
877 }
878
879 if (JdkMath.abs(result.getCovarianceOfParameters(0, 0)
880 - resultRedundant.getCovarianceOfParameters(0, 0)) > 1.0e-8) {
881 Assert.fail("VCV not correct after reorder (0,0)");
882 }
883 if (JdkMath.abs(result.getCovarianceOfParameters(0, 1)
884 - resultRedundant.getCovarianceOfParameters(0, 2)) > 1.0e-8) {
885 Assert.fail("VCV not correct after reorder (0,1)<->(0,2)");
886 }
887 if (JdkMath.abs(result.getCovarianceOfParameters(0, 2)
888 - resultRedundant.getCovarianceOfParameters(0, 3)) > 1.0e-8) {
889 Assert.fail("VCV not correct after reorder (0,2)<->(0,1)");
890 }
891 if (JdkMath.abs(result.getCovarianceOfParameters(0, 3)
892 - resultRedundant.getCovarianceOfParameters(0, 5)) > 1.0e-8) {
893 Assert.fail("VCV not correct after reorder (0,3)<->(0,3)");
894 }
895 if (JdkMath.abs(result.getCovarianceOfParameters(1, 0)
896 - resultRedundant.getCovarianceOfParameters(2, 0)) > 1.0e-8) {
897 Assert.fail("VCV not correct after reorder (1,0)<->(2,0)");
898 }
899 if (JdkMath.abs(result.getCovarianceOfParameters(1, 1)
900 - resultRedundant.getCovarianceOfParameters(2, 2)) > 1.0e-8) {
901 Assert.fail("VCV not correct (1,1)<->(2,1)");
902 }
903 if (JdkMath.abs(result.getCovarianceOfParameters(1, 2)
904 - resultRedundant.getCovarianceOfParameters(2, 3)) > 1.0e-8) {
905 Assert.fail("VCV not correct (1,2)<->(2,2)");
906 }
907
908 if (JdkMath.abs(result.getCovarianceOfParameters(2, 0)
909 - resultRedundant.getCovarianceOfParameters(3, 0)) > 1.0e-8) {
910 Assert.fail("VCV not correct (2,0)<->(1,0)");
911 }
912 if (JdkMath.abs(result.getCovarianceOfParameters(2, 1)
913 - resultRedundant.getCovarianceOfParameters(3, 2)) > 1.0e-8) {
914 Assert.fail("VCV not correct (2,1)<->(1,2)");
915 }
916
917 if (JdkMath.abs(result.getCovarianceOfParameters(3, 3)
918 - resultRedundant.getCovarianceOfParameters(5, 5)) > 1.0e-8) {
919 Assert.fail("VCV not correct (3,3)<->(3,2)");
920 }
921
922 TestUtils.assertEquals(result.getAdjustedRSquared(), resultRedundant.getAdjustedRSquared(), 1.0e-8);
923 TestUtils.assertEquals(result.getErrorSumSquares(), resultRedundant.getErrorSumSquares(), 1.0e-8);
924 TestUtils.assertEquals(result.getMeanSquareError(), resultRedundant.getMeanSquareError(), 1.0e-8);
925 TestUtils.assertEquals(result.getRSquared(), resultRedundant.getRSquared(), 1.0e-8);
926 return;
927 }
928
929 @Test
930 public void testPCorr() {
931 MillerUpdatingRegression instance = new MillerUpdatingRegression(4, false);
932 double[][] x = new double[airdata[0].length][];
933 double[] y = new double[airdata[0].length];
934 double[] cp = new double[10];
935 double[] yxcorr = new double[4];
936 double[] diag = new double[4];
937 double sumysq = 0.0;
938 int off = 0;
939 for (int i = 0; i < airdata[0].length; i++) {
940 x[i] = new double[4];
941 x[i][0] = 1.0;
942 x[i][1] = JdkMath.log(airdata[3][i]);
943 x[i][2] = JdkMath.log(airdata[4][i]);
944 x[i][3] = airdata[5][i];
945 y[i] = JdkMath.log(airdata[2][i]);
946 off = 0;
947 for (int j = 0; j < 4; j++) {
948 double tmp = x[i][j];
949 for (int k = 0; k <= j; k++, off++) {
950 cp[off] += tmp * x[i][k];
951 }
952 yxcorr[j] += tmp * y[i];
953 }
954 sumysq += y[i] * y[i];
955 }
956 PearsonsCorrelation pearson = new PearsonsCorrelation(x);
957 RealMatrix corr = pearson.getCorrelationMatrix();
958 off = 0;
959 for (int i = 0; i < 4; i++, off += i + 1) {
960 diag[i] = JdkMath.sqrt(cp[off]);
961 }
962
963 instance.addObservations(x, y);
964 double[] pc = instance.getPartialCorrelations(0);
965 int idx = 0;
966 off = 0;
967 int off2 = 6;
968 for (int i = 0; i < 4; i++) {
969 for (int j = 0; j < i; j++) {
970 if (JdkMath.abs(pc[idx] - cp[off] / (diag[i] * diag[j])) > 1.0e-8) {
971 Assert.fail("Failed cross products... i = " + i + " j = " + j);
972 }
973 ++idx;
974 ++off;
975 }
976 ++off;
977 if (JdkMath.abs(pc[i+off2] - yxcorr[ i] / (JdkMath.sqrt(sumysq) * diag[i])) > 1.0e-8) {
978 Assert.fail("Assert.failed cross product i = " + i + " y");
979 }
980 }
981 double[] pc2 = instance.getPartialCorrelations(1);
982
983 idx = 0;
984
985 for (int i = 1; i < 4; i++) {
986 for (int j = 1; j < i; j++) {
987 if (JdkMath.abs(pc2[idx] - corr.getEntry(j, i)) > 1.0e-8) {
988 Assert.fail("Failed cross products... i = " + i + " j = " + j);
989 }
990 ++idx;
991 }
992 }
993 double[] pc3 = instance.getPartialCorrelations(2);
994 if (pc3 == null) {
995 Assert.fail("Should not be null");
996 }
997 return;
998 }
999
1000 @Test
1001 public void testHdiag() {
1002 MillerUpdatingRegression instance = new MillerUpdatingRegression(4, false);
1003 double[][] x = new double[airdata[0].length][];
1004 double[] y = new double[airdata[0].length];
1005 for (int i = 0; i < airdata[0].length; i++) {
1006 x[i] = new double[4];
1007 x[i][0] = 1.0;
1008 x[i][1] = JdkMath.log(airdata[3][i]);
1009 x[i][2] = JdkMath.log(airdata[4][i]);
1010 x[i][3] = airdata[5][i];
1011 y[i] = JdkMath.log(airdata[2][i]);
1012 }
1013 instance.addObservations(x, y);
1014 OLSMultipleLinearRegression ols = new OLSMultipleLinearRegression();
1015 ols.setNoIntercept(true);
1016 ols.newSampleData(y, x);
1017
1018 RealMatrix rm = ols.calculateHat();
1019 for (int i = 0; i < x.length; i++) {
1020 TestUtils.assertEquals(instance.getDiagonalOfHatMatrix(x[i]), rm.getEntry(i, i), 1.0e-8);
1021 }
1022 return;
1023 }
1024 @Test
1025 public void testHdiagConstant() {
1026 MillerUpdatingRegression instance = new MillerUpdatingRegression(3, true);
1027 double[][] x = new double[airdata[0].length][];
1028 double[] y = new double[airdata[0].length];
1029 for (int i = 0; i < airdata[0].length; i++) {
1030 x[i] = new double[3];
1031 x[i][0] = JdkMath.log(airdata[3][i]);
1032 x[i][1] = JdkMath.log(airdata[4][i]);
1033 x[i][2] = airdata[5][i];
1034 y[i] = JdkMath.log(airdata[2][i]);
1035 }
1036 instance.addObservations(x, y);
1037 OLSMultipleLinearRegression ols = new OLSMultipleLinearRegression();
1038 ols.setNoIntercept(false);
1039 ols.newSampleData(y, x);
1040
1041 RealMatrix rm = ols.calculateHat();
1042 for (int i = 0; i < x.length; i++) {
1043 TestUtils.assertEquals(instance.getDiagonalOfHatMatrix(x[i]), rm.getEntry(i, i), 1.0e-8);
1044 }
1045 return;
1046 }
1047
1048
1049 private void subsetRegression(int iExclude, boolean constant){
1050 int[] indices = new int[2];
1051 int j = 0;
1052 for (int i = 0; i < 3; i++){
1053 if (i != iExclude){
1054 indices[j] = i;
1055 j++;
1056 }
1057 }
1058 int i0 = indices[0];
1059 int i1 = indices[1];
1060 MillerUpdatingRegression instance = new MillerUpdatingRegression(3, constant);
1061 MillerUpdatingRegression redRegression = new MillerUpdatingRegression(2, constant);
1062 double[][] x = new double[airdata[0].length][];
1063 double[][] xReduced = new double[airdata[0].length][];
1064 double[] y = new double[airdata[0].length];
1065 for (int i = 0; i < airdata[0].length; i++) {
1066 x[i] = new double[3];
1067 x[i][i0] = JdkMath.log(airdata[3][i]);
1068 x[i][i1] = JdkMath.log(airdata[4][i]);
1069 x[i][iExclude] = airdata[5][i];
1070
1071 xReduced[i] = new double[2];
1072 xReduced[i][0] = JdkMath.log(airdata[3][i]);
1073 xReduced[i][1] = JdkMath.log(airdata[4][i]);
1074
1075 y[i] = JdkMath.log(airdata[2][i]);
1076 }
1077
1078 instance.addObservations(x, y);
1079 redRegression.addObservations(xReduced, y);
1080
1081 int includedIndices[];
1082 if (constant){
1083 includedIndices = new int[3];
1084 includedIndices[0] = 0;
1085 includedIndices[1] = i0 + 1;
1086 includedIndices[2] = i1 + 1;
1087 } else {
1088 includedIndices = indices;
1089 }
1090
1091 RegressionResults resultsInstance = instance.regress( includedIndices );
1092 RegressionResults resultsReduced = redRegression.regress();
1093
1094 TestUtils.assertEquals(resultsInstance.getParameterEstimates(), resultsReduced.getParameterEstimates(), 1.0e-12);
1095 TestUtils.assertEquals(resultsInstance.getStdErrorOfEstimates(), resultsReduced.getStdErrorOfEstimates(), 1.0e-12);
1096 }
1097
1098
1099 @Test
1100 public void testSubsetRegression() {
1101 for (int i=0; i < 3; i++){
1102 subsetRegression(i, true);
1103 subsetRegression(i, false);
1104 }
1105 }
1106
1107 }