1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
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 * The GLS implementation of the multiple linear regression.
27 *
28 * GLS assumes a general covariance matrix Omega of the error
29 * <pre>
30 * u ~ N(0, Omega)
31 * </pre>
32 *
33 * Estimated by GLS,
34 * <pre>
35 * b=(X' Omega^-1 X)^-1X'Omega^-1 y
36 * </pre>
37 * whose variance is
38 * <pre>
39 * Var(b)=(X' Omega^-1 X)^-1
40 * </pre>
41 * @version $Revision: 783702 $ $Date: 2009-06-11 04:54:02 -0400 (Thu, 11 Jun 2009) $
42 * @since 2.0
43 */
44 public class GLSMultipleLinearRegression extends AbstractMultipleLinearRegression {
45
46 /** Covariance matrix. */
47 private RealMatrix Omega;
48
49 /** Inverse of covariance matrix. */
50 private RealMatrix OmegaInverse;
51
52 /** Replace sample data, overriding any previous sample.
53 * @param y y values of the sample
54 * @param x x values of the sample
55 * @param covariance array representing the covariance matrix
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 * Add the covariance data.
67 *
68 * @param omega the [n,n] array representing the covariance
69 */
70 protected void newCovarianceData(double[][] omega){
71 this.Omega = new Array2DRowRealMatrix(omega);
72 this.OmegaInverse = null;
73 }
74
75 /**
76 * Get the inverse of the covariance.
77 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
78 * @return inverse of the covariance
79 */
80 protected RealMatrix getOmegaInverse() {
81 if (OmegaInverse == null) {
82 OmegaInverse = new LUDecompositionImpl(Omega).getSolver().getInverse();
83 }
84 return OmegaInverse;
85 }
86
87 /**
88 * Calculates beta by GLS.
89 * <pre>
90 * b=(X' Omega^-1 X)^-1X'Omega^-1 y
91 * </pre>
92 * @return beta
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 * Calculates the variance on the beta by GLS.
105 * <pre>
106 * Var(b)=(X' Omega^-1 X)^-1
107 * </pre>
108 * @return The beta variance matrix
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 * Calculates the variance on the y by GLS.
119 * <pre>
120 * Var(y)=Tr(u' Omega^-1 u)/(n-k)
121 * </pre>
122 * @return The Y variance
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 }