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 */
017package org.apache.commons.math3.stat.regression;
018
019import org.apache.commons.math3.linear.LUDecomposition;
020import org.apache.commons.math3.linear.RealMatrix;
021import org.apache.commons.math3.linear.Array2DRowRealMatrix;
022import org.apache.commons.math3.linear.RealVector;
023
024/**
025 * The GLS implementation of 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 $Id: GLSMultipleLinearRegression.java 1553598 2013-12-26 22:18:02Z psteitz $
041 * @since 2.0
042 */
043public 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 LUDecomposition(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 = getX().transpose();
097        RealMatrix XTOIX = XT.multiply(OI).multiply(getX());
098        RealMatrix inverse = new LUDecomposition(XTOIX).getSolver().getInverse();
099        return inverse.multiply(XT).multiply(OI).operate(getY());
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 = getX().transpose().multiply(OI).multiply(getX());
113        return new LUDecomposition(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 / (getX().getRowDimension() - getX().getColumnDimension());
133
134    }
135
136}