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 }