/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.math3.stat.regression;

import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.stat.StatUtils;
import org.apache.commons.math3.stat.descriptive.moment.SecondMoment;

Implements ordinary least squares (OLS) to estimate the parameters of a
multiple linear regression model.
The regression coefficients, b, satisfy the normal equations:
 XT X b = XT y 

To solve the normal equations, this implementation uses QR decomposition
of the X matrix. (See QRDecomposition for details on the decomposition algorithm.) The X matrix, also known as the design matrix, has rows corresponding to sample observations and columns corresponding to independent variables. When the model is estimated using an intercept term (i.e. when isNoIntercept is false as it is by default), the X
matrix includes an initial column identically equal to 1.  We solve the normal equations
as follows:
 XTX b = XT y
(QR)T (QR) b = (QR)Ty
RT (QTQ) R b = RT QT y
RT R b = RT QT y
(RT)-1 RT R b = (RT)-1 RT QT y
R b = QT y 

Given Q and R, the last equation is solved by back-substitution.
Since:
2.0/**
 * <p>Implements ordinary least squares (OLS) to estimate the parameters of a
 * multiple linear regression model.</p>
 *
 * <p>The regression coefficients, <code>b</code>, satisfy the normal equations:
 * <pre><code> X<sup>T</sup> X b = X<sup>T</sup> y </code></pre></p>
 *
 * <p>To solve the normal equations, this implementation uses QR decomposition
 * of the <code>X</code> matrix. (See {@link QRDecomposition} for details on the
 * decomposition algorithm.) The <code>X</code> matrix, also known as the <i>design matrix,</i>
 * has rows corresponding to sample observations and columns corresponding to independent
 * variables.  When the model is estimated using an intercept term (i.e. when
 * {@link #isNoIntercept() isNoIntercept} is false as it is by default), the <code>X</code>
 * matrix includes an initial column identically equal to 1.  We solve the normal equations
 * as follows:
 * <pre><code> X<sup>T</sup>X b = X<sup>T</sup> y
 * (QR)<sup>T</sup> (QR) b = (QR)<sup>T</sup>y
 * R<sup>T</sup> (Q<sup>T</sup>Q) R b = R<sup>T</sup> Q<sup>T</sup> y
 * R<sup>T</sup> R b = R<sup>T</sup> Q<sup>T</sup> y
 * (R<sup>T</sup>)<sup>-1</sup> R<sup>T</sup> R b = (R<sup>T</sup>)<sup>-1</sup> R<sup>T</sup> Q<sup>T</sup> y
 * R b = Q<sup>T</sup> y </code></pre></p>
 *
 * <p>Given <code>Q</code> and <code>R</code>, the last equation is solved by back-substitution.</p>
 *
 * @since 2.0
 */
public class OLSMultipleLinearRegression extends AbstractMultipleLinearRegression {

    
Cached QR decomposition of X matrix /** Cached QR decomposition of X matrix */
    private QRDecomposition qr = null;

    
Singularity threshold for QR decomposition /** Singularity threshold for QR decomposition */
    private final double threshold;

    
Create an empty OLSMultipleLinearRegression instance.
/**
     * Create an empty OLSMultipleLinearRegression instance.
     */
    public OLSMultipleLinearRegression() {
        this(0d);
    }

    
Create an empty OLSMultipleLinearRegression instance, using the given
singularity threshold for the QR decomposition.
Params:
threshold – the singularity threshold
Since:
3.3/**
     * Create an empty OLSMultipleLinearRegression instance, using the given
     * singularity threshold for the QR decomposition.
     *
     * @param threshold the singularity threshold
     * @since 3.3
     */
    public OLSMultipleLinearRegression(final double threshold) {
        this.threshold = threshold;
    }

    
Loads model x and y sample data, overriding any previous sample.
Computes and caches QR decomposition of the X matrix.
Params:
y – the [n,1] array representing the y sample
x – the [n,k] array representing the x sample
Throws:
MathIllegalArgumentException – if the x and y array data are not
            compatible for the regression/**
     * Loads model x and y sample data, overriding any previous sample.
     *
     * Computes and caches QR decomposition of the X matrix.
     * @param y the [n,1] array representing the y sample
     * @param x the [n,k] array representing the x sample
     * @throws MathIllegalArgumentException if the x and y array data are not
     *             compatible for the regression
     */
    public void newSampleData(double[] y, double[][] x) throws MathIllegalArgumentException {
        validateSampleData(x, y);
        newYSampleData(y);
        newXSampleData(x);
    }

    
{@inheritDoc}
This implementation computes and caches the QR decomposition of the X matrix.
/**
     * {@inheritDoc}
     * <p>This implementation computes and caches the QR decomposition of the X matrix.</p>
     */
    @Override
    public void newSampleData(double[] data, int nobs, int nvars) {
        super.newSampleData(data, nobs, nvars);
        qr = new QRDecomposition(getX(), threshold);
    }

    
Compute the "hat" matrix.

The hat matrix is defined in terms of the design matrix X
 by X(XTX)-1XT

The implementation here uses the QR decomposition to compute the
hat matrix as Q IpQT where Ip is the
p-dimensional identity matrix augmented by 0's.  This computational
formula is from "The Hat Matrix in Regression and ANOVA",
David C. Hoaglin and Roy E. Welsch,
The American Statistician, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.

Data for the model must have been successfully loaded using one of the newSampleData methods before invoking this method; otherwise a NullPointerException will be thrown.
Throws:
NullPointerException – unless method newSampleData has been called beforehand.
Returns:
the hat matrix/**
     * <p>Compute the "hat" matrix.
     * </p>
     * <p>The hat matrix is defined in terms of the design matrix X
     *  by X(X<sup>T</sup>X)<sup>-1</sup>X<sup>T</sup>
     * </p>
     * <p>The implementation here uses the QR decomposition to compute the
     * hat matrix as Q I<sub>p</sub>Q<sup>T</sup> where I<sub>p</sub> is the
     * p-dimensional identity matrix augmented by 0's.  This computational
     * formula is from "The Hat Matrix in Regression and ANOVA",
     * David C. Hoaglin and Roy E. Welsch,
     * <i>The American Statistician</i>, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
     * </p>
     * <p>Data for the model must have been successfully loaded using one of
     * the {@code newSampleData} methods before invoking this method; otherwise
     * a {@code NullPointerException} will be thrown.</p>
     *
     * @return the hat matrix
     * @throws NullPointerException unless method {@code newSampleData} has been
     * called beforehand.
     */
    public RealMatrix calculateHat() {
        // Create augmented identity matrix
        RealMatrix Q = qr.getQ();
        final int p = qr.getR().getColumnDimension();
        final int n = Q.getColumnDimension();
        // No try-catch or advertised NotStrictlyPositiveException - NPE above if n < 3
        Array2DRowRealMatrix augI = new Array2DRowRealMatrix(n, n);
        double[][] augIData = augI.getDataRef();
        for (int i = 0; i < n; i++) {
            for (int j =0; j < n; j++) {
                if (i == j && i < p) {
                    augIData[i][j] = 1d;
                } else {
                    augIData[i][j] = 0d;
                }
            }
        }

        // Compute and return Hat matrix
        // No DME advertised - args valid if we get here
        return Q.multiply(augI).multiply(Q.transpose());
    }

    
Returns the sum of squared deviations of Y from its mean.
If the model has no intercept term, 0 is used for the
mean of Y - i.e., what is returned is the sum of the squared Y values.
The value returned by this method is the SSTO value used in the R-squared computation.
Throws:
NullPointerException – if the sample has not been set
See Also:
AbstractMultipleLinearRegression.isNoIntercept()
Returns:
SSTO - the total sum of squares
Since:
2.2/**
     * <p>Returns the sum of squared deviations of Y from its mean.</p>
     *
     * <p>If the model has no intercept term, <code>0</code> is used for the
     * mean of Y - i.e., what is returned is the sum of the squared Y values.</p>
     *
     * <p>The value returned by this method is the SSTO value used in
     * the {@link #calculateRSquared() R-squared} computation.</p>
     *
     * @return SSTO - the total sum of squares
     * @throws NullPointerException if the sample has not been set
     * @see #isNoIntercept()
     * @since 2.2
     */
    public double calculateTotalSumOfSquares() {
        if (isNoIntercept()) {
            return StatUtils.sumSq(getY().toArray());
        } else {
            return new SecondMoment().evaluate(getY().toArray());
        }
    }

    
Returns the sum of squared residuals.
Throws:
SingularMatrixException – if the design matrix is singular
NullPointerException – if the data for the model have not been loaded
Returns:
residual sum of squares
Since:
2.2/**
     * Returns the sum of squared residuals.
     *
     * @return residual sum of squares
     * @since 2.2
     * @throws org.apache.commons.math3.linear.SingularMatrixException if the design matrix is singular
     * @throws NullPointerException if the data for the model have not been loaded
     */
    public double calculateResidualSumOfSquares() {
        final RealVector residuals = calculateResiduals();
        // No advertised DME, args are valid
        return residuals.dotProduct(residuals);
    }

    
Returns the R-Squared statistic, defined by the formula 
R2 = 1 - SSR / SSTO
 where SSR is the sum of squared residuals and SSTO is the total sum of squares 
If there is no variance in y, i.e., SSTO = 0, NaN is returned.
Throws:
NullPointerException – if the sample has not been set
SingularMatrixException – if the design matrix is singular
Returns:
R-square statistic
Since:
2.2/**
     * Returns the R-Squared statistic, defined by the formula <pre>
     * R<sup>2</sup> = 1 - SSR / SSTO
     * </pre>
     * where SSR is the {@link #calculateResidualSumOfSquares() sum of squared residuals}
     * and SSTO is the {@link #calculateTotalSumOfSquares() total sum of squares}
     *
     * <p>If there is no variance in y, i.e., SSTO = 0, NaN is returned.</p>
     *
     * @return R-square statistic
     * @throws NullPointerException if the sample has not been set
     * @throws org.apache.commons.math3.linear.SingularMatrixException if the design matrix is singular
     * @since 2.2
     */
    public double calculateRSquared() {
        return 1 - calculateResidualSumOfSquares() / calculateTotalSumOfSquares();
    }

    
Returns the adjusted R-squared statistic, defined by the formula 
R2adj = 1 - [SSR (n - 1)] / [SSTO (n - p)]
 where SSR is the sum of squared residuals, SSTO is the total sum of squares, n is the number of observations and p is the number of parameters estimated (including the intercept).

If the regression is estimated without an intercept term, what is returned is 
 1 - (1 - calculateRSquared()) * (n / (n - p)) 

If there is no variance in y, i.e., SSTO = 0, NaN is returned.
Throws:
NullPointerException – if the sample has not been set
SingularMatrixException – if the design matrix is singular
See Also:
AbstractMultipleLinearRegression.isNoIntercept()
Returns:
adjusted R-Squared statistic
Since:
2.2/**
     * <p>Returns the adjusted R-squared statistic, defined by the formula <pre>
     * R<sup>2</sup><sub>adj</sub> = 1 - [SSR (n - 1)] / [SSTO (n - p)]
     * </pre>
     * where SSR is the {@link #calculateResidualSumOfSquares() sum of squared residuals},
     * SSTO is the {@link #calculateTotalSumOfSquares() total sum of squares}, n is the number
     * of observations and p is the number of parameters estimated (including the intercept).</p>
     *
     * <p>If the regression is estimated without an intercept term, what is returned is <pre>
     * <code> 1 - (1 - {@link #calculateRSquared()}) * (n / (n - p)) </code>
     * </pre></p>
     *
     * <p>If there is no variance in y, i.e., SSTO = 0, NaN is returned.</p>
     *
     * @return adjusted R-Squared statistic
     * @throws NullPointerException if the sample has not been set
     * @throws org.apache.commons.math3.linear.SingularMatrixException if the design matrix is singular
     * @see #isNoIntercept()
     * @since 2.2
     */
    public double calculateAdjustedRSquared() {
        final double n = getX().getRowDimension();
        if (isNoIntercept()) {
            return 1 - (1 - calculateRSquared()) * (n / (n - getX().getColumnDimension()));
        } else {
            return 1 - (calculateResidualSumOfSquares() * (n - 1)) /
                (calculateTotalSumOfSquares() * (n - getX().getColumnDimension()));
        }
    }

    
{@inheritDoc}
This implementation computes and caches the QR decomposition of the X matrix
once it is successfully loaded.
/**
     * {@inheritDoc}
     * <p>This implementation computes and caches the QR decomposition of the X matrix
     * once it is successfully loaded.</p>
     */
    @Override
    protected void newXSampleData(double[][] x) {
        super.newXSampleData(x);
        qr = new QRDecomposition(getX(), threshold);
    }

    
Calculates the regression coefficients using OLS.
Data for the model must have been successfully loaded using one of the newSampleData methods before invoking this method; otherwise a NullPointerException will be thrown.
Throws:
SingularMatrixException – if the design matrix is singular
NullPointerException – if the data for the model have not been loaded
Returns:
beta/**
     * Calculates the regression coefficients using OLS.
     *
     * <p>Data for the model must have been successfully loaded using one of
     * the {@code newSampleData} methods before invoking this method; otherwise
     * a {@code NullPointerException} will be thrown.</p>
     *
     * @return beta
     * @throws org.apache.commons.math3.linear.SingularMatrixException if the design matrix is singular
     * @throws NullPointerException if the data for the model have not been loaded
     */
    @Override
    protected RealVector calculateBeta() {
        return qr.getSolver().solve(getY());
    }

    
Calculates the variance-covariance matrix of the regression parameters.

Var(b) = (XTX)-1

Uses QR decomposition to reduce (XTX)-1
to (RTR)-1, with only the top p rows of
R included, where p = the length of the beta vector.
Data for the model must have been successfully loaded using one of the newSampleData methods before invoking this method; otherwise a NullPointerException will be thrown.
Throws:
SingularMatrixException – if the design matrix is singular
NullPointerException – if the data for the model have not been loaded
Returns:
The beta variance-covariance matrix/**
     * <p>Calculates the variance-covariance matrix of the regression parameters.
     * </p>
     * <p>Var(b) = (X<sup>T</sup>X)<sup>-1</sup>
     * </p>
     * <p>Uses QR decomposition to reduce (X<sup>T</sup>X)<sup>-1</sup>
     * to (R<sup>T</sup>R)<sup>-1</sup>, with only the top p rows of
     * R included, where p = the length of the beta vector.</p>
     *
     * <p>Data for the model must have been successfully loaded using one of
     * the {@code newSampleData} methods before invoking this method; otherwise
     * a {@code NullPointerException} will be thrown.</p>
     *
     * @return The beta variance-covariance matrix
     * @throws org.apache.commons.math3.linear.SingularMatrixException if the design matrix is singular
     * @throws NullPointerException if the data for the model have not been loaded
     */
    @Override
    protected RealMatrix calculateBetaVariance() {
        int p = getX().getColumnDimension();
        RealMatrix Raug = qr.getR().getSubMatrix(0, p - 1 , 0, p - 1);
        RealMatrix Rinv = new LUDecomposition(Raug).getSolver().getInverse();
        return Rinv.multiply(Rinv.transpose());
    }

}