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 /** 026 * The GLS implementation of the multiple linear regression. 027 * 028 * GLS assumes a general covariance matrix Omega of the error 029 * <pre> 030 * u ~ N(0, Omega) 031 * </pre> 032 * 033 * Estimated by GLS, 034 * <pre> 035 * b=(X' Omega^-1 X)^-1X'Omega^-1 y 036 * </pre> 037 * whose variance is 038 * <pre> 039 * Var(b)=(X' Omega^-1 X)^-1 040 * </pre> 041 * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $ 042 * @since 2.0 043 */ 044 public class GLSMultipleLinearRegression extends AbstractMultipleLinearRegression { 045 046 /** Covariance matrix. */ 047 private RealMatrix Omega; 048 049 /** Inverse of covariance matrix. */ 050 private RealMatrix OmegaInverse; 051 052 /** Replace sample data, overriding any previous sample. 053 * @param y y values of the sample 054 * @param x x values of the sample 055 * @param covariance array representing the covariance matrix 056 */ 057 public void newSampleData(double[] y, double[][] x, double[][] covariance) { 058 validateSampleData(x, y); 059 newYSampleData(y); 060 newXSampleData(x); 061 validateCovarianceData(x, covariance); 062 newCovarianceData(covariance); 063 } 064 065 /** 066 * Add the covariance data. 067 * 068 * @param omega the [n,n] array representing the covariance 069 */ 070 protected void newCovarianceData(double[][] omega){ 071 this.Omega = new Array2DRowRealMatrix(omega); 072 this.OmegaInverse = null; 073 } 074 075 /** 076 * Get the inverse of the covariance. 077 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p> 078 * @return inverse of the covariance 079 */ 080 protected RealMatrix getOmegaInverse() { 081 if (OmegaInverse == null) { 082 OmegaInverse = new LUDecompositionImpl(Omega).getSolver().getInverse(); 083 } 084 return OmegaInverse; 085 } 086 087 /** 088 * Calculates beta by GLS. 089 * <pre> 090 * b=(X' Omega^-1 X)^-1X'Omega^-1 y 091 * </pre> 092 * @return beta 093 */ 094 @Override 095 protected RealVector calculateBeta() { 096 RealMatrix OI = getOmegaInverse(); 097 RealMatrix XT = X.transpose(); 098 RealMatrix XTOIX = XT.multiply(OI).multiply(X); 099 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 }