1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.stat.regression; 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.linear.LUDecompositionImpl; 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.linear.RealMatrix; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.linear.Array2DRowRealMatrix; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.linear.RealVector; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The GLS implementation of the multiple linear regression. 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * GLS assumes a general covariance matrix Omega of the error 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * u ~ N(0, Omega) 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Estimated by GLS, 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * b=(X' Omega^-1 X)^-1X'Omega^-1 y 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * whose variance is 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Var(b)=(X' Omega^-1 X)^-1 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1073460 $ $Date: 2011-02-22 20:22:39 +0100 (mar. 22 févr. 2011) $ 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.0 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class GLSMultipleLinearRegression extends AbstractMultipleLinearRegression { 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Covariance matrix. */ 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix Omega; 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Inverse of covariance matrix. */ 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix OmegaInverse; 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Replace sample data, overriding any previous sample. 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y y values of the sample 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x x values of the sample 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param covariance array representing the covariance matrix 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void newSampleData(double[] y, double[][] x, double[][] covariance) { 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond validateSampleData(x, y); 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond newYSampleData(y); 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond newXSampleData(x); 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond validateCovarianceData(x, covariance); 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond newCovarianceData(covariance); 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Add the covariance data. 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param omega the [n,n] array representing the covariance 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected void newCovarianceData(double[][] omega){ 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.Omega = new Array2DRowRealMatrix(omega); 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.OmegaInverse = null; 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Get the inverse of the covariance. 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p> 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return inverse of the covariance 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected RealMatrix getOmegaInverse() { 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (OmegaInverse == null) { 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond OmegaInverse = new LUDecompositionImpl(Omega).getSolver().getInverse(); 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return OmegaInverse; 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calculates beta by GLS. 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * b=(X' Omega^-1 X)^-1X'Omega^-1 y 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return beta 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected RealVector calculateBeta() { 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealMatrix OI = getOmegaInverse(); 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealMatrix XT = X.transpose(); 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealMatrix XTOIX = XT.multiply(OI).multiply(X); 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse(); 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return inverse.multiply(XT).multiply(OI).operate(Y); 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calculates the variance on the beta. 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Var(b)=(X' Omega^-1 X)^-1 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return The beta variance matrix 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected RealMatrix calculateBetaVariance() { 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealMatrix OI = getOmegaInverse(); 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealMatrix XTOIX = X.transpose().multiply(OI).multiply(X); 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return new LUDecompositionImpl(XTOIX).getSolver().getInverse(); 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calculates the estimated variance of the error term using the formula 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Var(u) = Tr(u' Omega^-1 u)/(n-k) 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * where n and k are the row and column dimensions of the design 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * matrix X. 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return error variance 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double calculateErrorVariance() { 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond RealVector residuals = calculateResiduals(); 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double t = residuals.dotProduct(getOmegaInverse().operate(residuals)); 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return t / (X.getRowDimension() - X.getColumnDimension()); 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 137