001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math.stat.regression;
018
019 import org.apache.commons.math.linear.LUDecompositionImpl;
020 import org.apache.commons.math.linear.RealMatrix;
021 import org.apache.commons.math.linear.Array2DRowRealMatrix;
022 import org.apache.commons.math.linear.RealVector;
023
024 /**
025 * The GLS implementation of the multiple linear regression.
026 *
027 * GLS assumes a general covariance matrix Omega of the error
028 * <pre>
029 * u ~ N(0, Omega)
030 * </pre>
031 *
032 * Estimated by GLS,
033 * <pre>
034 * b=(X' Omega^-1 X)^-1X'Omega^-1 y
035 * </pre>
036 * whose variance is
037 * <pre>
038 * Var(b)=(X' Omega^-1 X)^-1
039 * </pre>
040 * @version $Revision: 1073460 $ $Date: 2011-02-22 20:22:39 +0100 (mar. 22 f??vr. 2011) $
041 * @since 2.0
042 */
043 public class GLSMultipleLinearRegression extends AbstractMultipleLinearRegression {
044
045 /** Covariance matrix. */
046 private RealMatrix Omega;
047
048 /** Inverse of covariance matrix. */
049 private RealMatrix OmegaInverse;
050
051 /** Replace sample data, overriding any previous sample.
052 * @param y y values of the sample
053 * @param x x values of the sample
054 * @param covariance array representing the covariance matrix
055 */
056 public void newSampleData(double[] y, double[][] x, double[][] covariance) {
057 validateSampleData(x, y);
058 newYSampleData(y);
059 newXSampleData(x);
060 validateCovarianceData(x, covariance);
061 newCovarianceData(covariance);
062 }
063
064 /**
065 * Add the covariance data.
066 *
067 * @param omega the [n,n] array representing the covariance
068 */
069 protected void newCovarianceData(double[][] omega){
070 this.Omega = new Array2DRowRealMatrix(omega);
071 this.OmegaInverse = null;
072 }
073
074 /**
075 * Get the inverse of the covariance.
076 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
077 * @return inverse of the covariance
078 */
079 protected RealMatrix getOmegaInverse() {
080 if (OmegaInverse == null) {
081 OmegaInverse = new LUDecompositionImpl(Omega).getSolver().getInverse();
082 }
083 return OmegaInverse;
084 }
085
086 /**
087 * Calculates beta by GLS.
088 * <pre>
089 * b=(X' Omega^-1 X)^-1X'Omega^-1 y
090 * </pre>
091 * @return beta
092 */
093 @Override
094 protected RealVector calculateBeta() {
095 RealMatrix OI = getOmegaInverse();
096 RealMatrix XT = X.transpose();
097 RealMatrix XTOIX = XT.multiply(OI).multiply(X);
098 RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
099 return inverse.multiply(XT).multiply(OI).operate(Y);
100 }
101
102 /**
103 * Calculates the variance on the beta.
104 * <pre>
105 * Var(b)=(X' Omega^-1 X)^-1
106 * </pre>
107 * @return The beta variance matrix
108 */
109 @Override
110 protected RealMatrix calculateBetaVariance() {
111 RealMatrix OI = getOmegaInverse();
112 RealMatrix XTOIX = X.transpose().multiply(OI).multiply(X);
113 return new LUDecompositionImpl(XTOIX).getSolver().getInverse();
114 }
115
116
117 /**
118 * Calculates the estimated variance of the error term using the formula
119 * <pre>
120 * Var(u) = Tr(u' Omega^-1 u)/(n-k)
121 * </pre>
122 * where n and k are the row and column dimensions of the design
123 * matrix X.
124 *
125 * @return error variance
126 * @since 2.2
127 */
128 @Override
129 protected double calculateErrorVariance() {
130 RealVector residuals = calculateResiduals();
131 double t = residuals.dotProduct(getOmegaInverse().operate(residuals));
132 return t / (X.getRowDimension() - X.getColumnDimension());
133
134 }
135
136 }