/*
    * 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.math4.legacy.stat.regression;
  
  import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
  import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
  import org.apache.commons.math4.legacy.linear.LUDecomposition;
  import org.apache.commons.math4.legacy.linear.QRDecomposition;
  import org.apache.commons.math4.legacy.linear.RealMatrix;
  import org.apache.commons.math4.legacy.linear.RealVector;
  import org.apache.commons.math4.legacy.stat.StatUtils;
  import org.apache.commons.math4.legacy.stat.descriptive.moment.SecondMoment;
  
  /**
   * <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>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>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. */
      private QRDecomposition qr;
  
      /** Singularity threshold for QR decomposition. */
      private final double threshold;
  
      /**
       * Create an empty OLSMultipleLinearRegression instance.
       */
      public OLSMultipleLinearRegression() {
          this(0d);
      }
  
      /**
       * 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.
       * @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}
       * <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);
     }
 
     /**
      * <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());
     }
 
     /**
      * <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.
      *
      * @return residual sum of squares
      * @since 2.2
      * @throws org.apache.commons.math4.legacy.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 <div style="white-space: pre"><code>
      * R<sup>2</sup> = 1 - SSR / SSTO
      * </code></div>
      * 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.math4.legacy.linear.SingularMatrixException if the design matrix is singular
      * @since 2.2
      */
     public double calculateRSquared() {
         return 1 - calculateResidualSumOfSquares() / calculateTotalSumOfSquares();
     }
 
     /**
      * <p>Returns the adjusted R-squared statistic, defined by the formula <div style="white-space: pre"><code>
      * R<sup>2</sup><sub>adj</sub> = 1 - [SSR (n - 1)] / [SSTO (n - p)]
      * </code></div>
      * 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>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>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.math4.legacy.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}
      * <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.
      *
      * <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.math4.legacy.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());
     }
 
     /**
      * <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.math4.legacy.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());
     }
 }