View Javadoc

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 }