/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package org.statcato.statistics.inferential;

import java.util.Vector;
import Jama.Matrix;
import org.statcato.statistics.FProbabilityDistribution;

/**
 * Multiple regression.
 * 
 * @author Margaret Yau
 * @version %I%, %G%
 * @since 1.0
 */
public class MultipleRegression {
    private int n;  // sample size
    private int k;  // number of independent variables
    private Matrix A;   // n x (k + 1) matrix - 
                // [1 x_11 x_21 ... x_k1
                //  1 x_12 x_22 ... x_k2 
                //    ...
                //  1 x_1n x_2n ... x_kn]  - independent variables in columns
    private Matrix x;   // k x 1 matrix - [b_0 b_1 ... b_k]' - regression equations coefficients
    private Matrix b;   // n x 1 matrix - [y_1 y_2 ... y_n]' - dependent variable
    
    /**
     * Constructor.
     * 
     * @param IndependentVars a vector of vectors of double, where each 
     *          vector is an independent variable with the same number of values
     * @param DependentVar a vector of double, which has the same number of 
     *          values as the independent variables
     */
    public MultipleRegression(Vector<Vector<Double>> IndependentVars,
            Vector<Double> DependentVar) {
        n = DependentVar.size();    
        k = IndependentVars.size();
        double[][] matrix = new double[n][1 + k];
        for (int i = 0; i < n; ++i)
            matrix[i][0] = 1;
        
        for (int col = 0; col < k; ++col) {
            for (int row = 0; row < n; ++row) {
                double a = (Double)IndependentVars.elementAt(col).elementAt(row).doubleValue();
                matrix[row][col + 1] = a;
            }
        }
        
        A = new Matrix(matrix);        
        
        double[][] array = new double[n][1];
        for (int i = 0; i < DependentVar.size(); ++i)
            array[i][0] = DependentVar.elementAt(i);
        
        b = new Matrix(array);
        try {
            x = A.solve(b);
        }
        catch (RuntimeException e) {
            // rank deficient matrix
            throw new RuntimeException("Multiple regression cannot be " +
                    "performed on a rank deficient matrix.");
        }
    }
    
    /**
     * Returns the coefficients of the regression equation 
     * y = b_0 + b_1 * x_1 + ... + b_k * x_k as a matrix of
     * dimension k+1 by 1: [b_0 b_1 ... b_k]'.
     * 
     * @return a k+1 by 1 matrix containing the coefficients of the 
     * regression equation
     */
    public Matrix RegressionEqCoefficients() {
        return x;
    }
    
    /**
     * Returns the predicted y value given a vector of values of the
     * independent variables using the regression equation.
     * 
     * @param var a vector of double that has the same number of values as
     * the number of independent variables.
     * 
     * @return predicted y value
     */
    public double YPredicted(Matrix var) {
        double y = x.get(0, 0);
        for (int i = 0; i < k; ++i) {
            y += (Double)var.get(0, i) * x.get(i+1, 0); 
        }
        return y;
    }
    
    /**
     * Returns the average value of the dependent variable.
     * @return average y value
     */
    private double YAverage() {
        double sum = 0;
        for (int i = 0; i < n; ++i) {
            sum += b.get(i, 0);
        }
        return sum / n;
    }
    
    /**
     * Returns the total variation (the sum of squared differences 
     * between the y values and the average y value).
     * 
     * @return total variation
     */
    public double TotalVariation() {
       double yAverage = YAverage();
       double total = 0;
       for (int i = 0; i < n; ++i) {
           total += Math.pow(b.get(i, 0) - yAverage, 2);
       }
       return total;
    }
    
    /**
     * Returns the i th value of all the independent variables as a 1 by k matrix.
     * 
     * @param i index
     * @return matrix
     */
    public Matrix XVar(int i) {
        return A.getMatrix(i, i, 1, k);
    }
    
    /**
     * Returns the explained variation (the sum of squared differences 
     * between the predicted y value and the average y value).
     * 
     * @return explained variation
     */
    public double ExplainedVariation() {
       double yAverage = YAverage();
       double total = 0;
       double yPredicted;
       for (int i = 0; i < n; ++i) {
           yPredicted = YPredicted(XVar(i));    // i th predicted y value
           total += Math.pow(yPredicted - yAverage, 2);
       }
       return total;
    }
    
    /**
     * Returns the unexplained variation (the sum of squared differences 
     * between the predicted y value and the y value).
     * 
     * @return unexplained variation
     */
    public double UnexplainedVariation() {
       double total = 0;
       double yPredicted;
       for (int i = 0; i < n; ++i) {
           yPredicted = YPredicted(XVar(i));
           total += Math.pow(yPredicted - b.get(i,0), 2);
       }
       return total;
    }

    /**
     * Returns the coefficient of determination r^2, the amount of the variation 
     * in y that is explained by the regression line.
     * 
     * @return r^2
     */
    public double CoefficientOfDetermination() {
        return ExplainedVariation() / TotalVariation();
    }

    /**
     * Returns the adjusted coefficient of determination.
     * 
     * @return r^2
     */
    public double AdjustedCoefficientOfDetermination() {
        double r2 = CoefficientOfDetermination();
        return 1 - (n - 1) * (1 - r2) / (n - k - 1);
    }
    
    /**
     * Returns the standard error of estimate, 
     * sqrt(unexplained variation / (n-2)).
     * 
     * @return s
     */
    public double StandardError() {
        return Math.sqrt(UnexplainedVariation() / (n - k - 1));
    }    
    
    /**
     * Returns the test statistics F.
     * F = (explained variation) / (unexplained variation) * (n - k - 1) / k
     * 
     * @return test statistics F
     */
    public double TestStatistics() {
        double unexplained = UnexplainedVariation();
        if (unexplained == 0)
            return Double.POSITIVE_INFINITY;
        return ExplainedVariation() / UnexplainedVariation() * (n - k - 1) / k;
    }
    
    /**
     * Returns the p-Value of test statistics.
     * 
     * @return p-Value
     */
    public double PValue() {
        double ts = TestStatistics();
        if (ts == Double.POSITIVE_INFINITY)
            return 0;
        
        FProbabilityDistribution dist = new FProbabilityDistribution( 
                k, n - k - 1);
        
        double area = dist.cumulativeProbability(ts);
        return 1 - area;    // area to the right
    }
    
    /**
     * Returns the sample size.
     * 
     * @return n
     */
    public int SampleSize() {
        return n;
    }
    
    /** 
     * Returns the number of independent variables.
     * 
     * @return k
     */
    public int NumIndepVar() {
        return k;
    }
}
