1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math.stat.regression;
18
19 import org.apache.commons.math.linear.LUDecompositionImpl;
20 import org.apache.commons.math.linear.RealMatrix;
21 import org.apache.commons.math.linear.Array2DRowRealMatrix;
22 import org.apache.commons.math.linear.RealVector;
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44 public class GLSMultipleLinearRegression extends AbstractMultipleLinearRegression {
45
46
47 private RealMatrix Omega;
48
49
50 private RealMatrix OmegaInverse;
51
52
53
54
55
56
57 public void newSampleData(double[] y, double[][] x, double[][] covariance) {
58 validateSampleData(x, y);
59 newYSampleData(y);
60 newXSampleData(x);
61 validateCovarianceData(x, covariance);
62 newCovarianceData(covariance);
63 }
64
65
66
67
68
69
70 protected void newCovarianceData(double[][] omega){
71 this.Omega = new Array2DRowRealMatrix(omega);
72 this.OmegaInverse = null;
73 }
74
75
76
77
78
79
80 protected RealMatrix getOmegaInverse() {
81 if (OmegaInverse == null) {
82 OmegaInverse = new LUDecompositionImpl(Omega).getSolver().getInverse();
83 }
84 return OmegaInverse;
85 }
86
87
88
89
90
91
92
93
94 @Override
95 protected RealVector calculateBeta() {
96 RealMatrix OI = getOmegaInverse();
97 RealMatrix XT = X.transpose();
98 RealMatrix XTOIX = XT.multiply(OI).multiply(X);
99 RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
100 return inverse.multiply(XT).multiply(OI).operate(Y);
101 }
102
103
104
105
106
107
108
109
110 @Override
111 protected RealMatrix calculateBetaVariance() {
112 RealMatrix OI = getOmegaInverse();
113 RealMatrix XTOIX = X.transpose().multiply(OI).multiply(X);
114 return new LUDecompositionImpl(XTOIX).getSolver().getInverse();
115 }
116
117
118
119
120
121
122
123
124 @Override
125 protected double calculateYVariance() {
126 RealVector residuals = calculateResiduals();
127 double t = residuals.dotProduct(getOmegaInverse().operate(residuals));
128 return t / (X.getRowDimension() - X.getColumnDimension());
129 }
130
131 }