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.math4.legacy.stat.regression;
18
19 import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
20 import org.apache.commons.math4.legacy.linear.LUDecomposition;
21 import org.apache.commons.math4.legacy.linear.RealMatrix;
22 import org.apache.commons.math4.legacy.linear.RealVector;
23
24 /**
25 * The GLS implementation of multiple linear regression.
26 *
27 * GLS assumes a general covariance matrix Omega of the error
28 * <pre>
29 * u ~ N(0, Omega)
30 * </pre>
31 *
32 * Estimated by GLS,
33 * <pre>
34 * b=(X' Omega^-1 X)^-1X'Omega^-1 y
35 * </pre>
36 * whose variance is
37 * <pre>
38 * Var(b)=(X' Omega^-1 X)^-1
39 * </pre>
40 * @since 2.0
41 */
42 public class GLSMultipleLinearRegression extends AbstractMultipleLinearRegression {
43
44 /** Covariance matrix. */
45 private RealMatrix omega;
46
47 /** Inverse of covariance matrix. */
48 private RealMatrix omegaInverse;
49
50 /** Replace sample data, overriding any previous sample.
51 * @param y y values of the sample
52 * @param x x values of the sample
53 * @param covariance array representing the covariance matrix
54 */
55 public void newSampleData(double[] y, double[][] x, double[][] covariance) {
56 validateSampleData(x, y);
57 newYSampleData(y);
58 newXSampleData(x);
59 validateCovarianceData(x, covariance);
60 newCovarianceData(covariance);
61 }
62
63 /**
64 * Add the covariance data.
65 *
66 * @param omega the [n,n] array representing the covariance
67 */
68 protected void newCovarianceData(double[][] omega){
69 this.omega = new Array2DRowRealMatrix(omega);
70 this.omegaInverse = null;
71 }
72
73 /**
74 * Get the inverse of the covariance.
75 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
76 * @return inverse of the covariance
77 */
78 protected RealMatrix getOmegaInverse() {
79 if (omegaInverse == null) {
80 omegaInverse = new LUDecomposition(omega).getSolver().getInverse();
81 }
82 return omegaInverse;
83 }
84
85 /**
86 * Calculates beta by GLS.
87 * <pre>
88 * b=(X' Omega^-1 X)^-1X'Omega^-1 y
89 * </pre>
90 * @return beta
91 */
92 @Override
93 protected RealVector calculateBeta() {
94 RealMatrix oi = getOmegaInverse();
95 RealMatrix xt = getX().transpose();
96 RealMatrix xtoix = xt.multiply(oi).multiply(getX());
97 RealMatrix inverse = new LUDecomposition(xtoix).getSolver().getInverse();
98 return inverse.multiply(xt).multiply(oi).operate(getY());
99 }
100
101 /**
102 * Calculates the variance on the beta.
103 * <pre>
104 * Var(b)=(X' Omega^-1 X)^-1
105 * </pre>
106 * @return The beta variance matrix
107 */
108 @Override
109 protected RealMatrix calculateBetaVariance() {
110 RealMatrix oi = getOmegaInverse();
111 RealMatrix xtoix = getX().transpose().multiply(oi).multiply(getX());
112 return new LUDecomposition(xtoix).getSolver().getInverse();
113 }
114
115
116 /**
117 * Calculates the estimated variance of the error term using the formula
118 * <pre>
119 * Var(u) = Tr(u' Omega^-1 u)/(n-k)
120 * </pre>
121 * where n and k are the row and column dimensions of the design
122 * matrix X.
123 *
124 * @return error variance
125 * @since 2.2
126 */
127 @Override
128 protected double calculateErrorVariance() {
129 RealVector residuals = calculateResiduals();
130 double t = residuals.dotProduct(getOmegaInverse().operate(residuals));
131 return t / (getX().getRowDimension() - getX().getColumnDimension());
132 }
133 }