HackerRank Day 9: Multiple Linear Regression | 10 Days of Statistics solution

In this Hackerrank Day 9: Multiple Linear Regression 10 Days of Statistics problem we have an equation y = a+ b*f +…+bm*fm. so for each of the q feature sets, we need to print the value of y on the new line.

Problem solution in Python programming.

# Enter your code here. Read input from STDIN. Print output to STDOUT
#import data
import numpy as np
m,n = [int(i) for i in input().strip().split(' ')]
X = []
Y = []
for i in range(n):
    data = input().strip().split(' ')
    X.append(data[:m])
    Y.append(data[m:])
q = int(input().strip())
X_new = []
for x in range(q):
    X_new.append(input().strip().split(' '))
X = np.array(X,float)
Y = np.array(Y,float)
X_new = np.array(X_new,float)

#center
X_R = X-np.mean(X,axis=0)
Y_R = Y-np.mean(Y)

#calculate beta
beta = np.dot(np.linalg.inv(np.dot(X_R.T,X_R)),np.dot(X_R.T,Y_R))

#predict
X_new_R = X_new-np.mean(X,axis=0)
Y_new_R = np.dot(X_new_R,beta)
Y_new = Y_new_R + np.mean(Y)

#print
for i in Y_new:
    print(round(float(i),2))

Problem solution in Java Programming.

import java.io.*;
import java.util.*;
import java.util.Scanner;
import java.util.Arrays;

public class Solution {
    public static void main(String[] args) {
        /* Read input: Create and fill X,Y arrays */
        Scanner scan = new Scanner(System.in);
        int m = scan.nextInt();
        int n = scan.nextInt();
        double [][] X = new double[n][m + 1];
        double [][] Y   = new double[n][1];
        for (int row = 0; row < n; row++) {
            X[row][0] = 1;
            for (int col = 1; col <= m; col++) {
                X[row][col] = scan.nextDouble();
            }
            Y[row][0] = scan.nextDouble();
        }

        /* Calculate B */
        double [][] xtx    = multiply(transpose(X),X);
        double [][] xtxInv = invert(xtx);
        double [][] xty    = multiply(transpose(X), Y);
        double [][] B      = multiply(xtxInv, xty);
        
        int sizeB = B.length;
        
        /* Calculate and print values for the "q" feature sets */
        int q = scan.nextInt();
        for (int i = 0; i < q; i++) {
            double result = B[0][0];
            for (int row = 1; row < sizeB; row++) {
                result += scan.nextDouble() * B[row][0];
            }
            System.out.println(result);
        }
    }
    
    /* Multiplies 2 matrices in O(n^3) time */
    public static double[][] multiply(double [][] A, double [][] B) {
        int aRows = A.length;
        int aCols = A[0].length;
        int bRows = B.length;
        int bCols = B[0].length;
        
        double [][] C = new double[aRows][bCols];
        int cRows = C.length;
        int cCols = C[0].length;
        
        for (int row = 0; row < cRows; row++) {
            for (int col = 0; col < cCols; col++) {
                for (int k = 0; k < aCols; k++) {
                    C[row][col] += A[row][k] * B[k][col];
                }
            }
        }
        return C;
    }
    
    public static double[][] transpose(double [][] matrix) {
        /* Create new array with switched dimensions */
        int originalRows = matrix.length;
        int originalCols = matrix[0].length;
        int rows = originalCols;
        int cols = originalRows;
        double [][] result = new double[rows][cols];
        
        /* Fill our new 2D array */
        for (int row = 0; row < originalRows; row++) {
            for (int col = 0; col < originalCols; col++) {
                result[col][row] = matrix[row][col];
            }
        }
        return result;
    }

    public static double[][] invert(double a[][]) 
    {
        int n = a.length;
        double x[][] = new double[n][n];
        double b[][] = new double[n][n];
        int index[] = new int[n];
        for (int i=0; i<n; ++i) 
            b[i][i] = 1;
 
         // Transform the matrix into an upper triangle
        gaussian(a, index);
 
         // Update the matrix b[i][j] with the ratios stored
        for (int i=0; i<n-1; ++i)
            for (int j=i+1; j<n; ++j)
                for (int k=0; k<n; ++k)
                    b[index[j]][k]
                            -= a[index[j]][i]*b[index[i]][k];
 
         // Perform backward substitutions
        for (int i=0; i<n; ++i) 
        {
            x[n-1][i] = b[index[n-1]][i]/a[index[n-1]][n-1];
            for (int j=n-2; j>=0; --j) 
            {
                x[j][i] = b[index[j]][i];
                for (int k=j+1; k<n; ++k) 
                {
                    x[j][i] -= a[index[j]][k]*x[k][i];
                }
                x[j][i] /= a[index[j]][j];
            }
        }
        return x;
    }
 
        // Method to carry out the partial-pivoting Gaussian
        // elimination.  Here index[] stores pivoting order.
 
    public static void gaussian(double a[][], int index[]) 
    {
        int n = index.length;
        double c[] = new double[n];
 
         // Initialize the index
        for (int i=0; i<n; ++i) 
            index[i] = i;
 
         // Find the rescaling factors, one from each row
        for (int i=0; i<n; ++i) 
        {
            double c1 = 0;
            for (int j=0; j<n; ++j) 
            {
                double c0 = Math.abs(a[i][j]);
                if (c0 > c1) c1 = c0;
            }
            c[i] = c1;
        }
 
         // Search the pivoting element from each column
        int k = 0;
        for (int j=0; j<n-1; ++j) 
        {
            double pi1 = 0;
            for (int i=j; i<n; ++i) 
            {
                double pi0 = Math.abs(a[index[i]][j]);
                pi0 /= c[index[i]];
                if (pi0 > pi1) 
                {
                    pi1 = pi0;
                    k = i;
                }
            }
 
           // Interchange rows according to the pivoting order
            int itmp = index[j];
            index[j] = index[k];
            index[k] = itmp;
            for (int i=j+1; i<n; ++i)   
            {
                double pj = a[index[i]][j]/a[index[j]][j];
 
                 // Record pivoting ratios below the diagonal
                a[index[i]][j] = pj;
 
                 // Modify other elements accordingly
                for (int l=j+1; l<n; ++l)
                    a[index[i]][l] -= pj*a[index[j]][l];
            }
        }
    }
}

Problem solution in C++ programming.

#include <cmath>
#include <cstdio>
#include <vector>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <cassert>

class Matrix {
public:
    Matrix(int x, int y) : _x(x), _y(y) {
        for (int i = 0; i < x; ++i) {
            _matrix.push_back(std::vector<double>(y, 0.0));
        }
    }

    Matrix transpose() {
        Matrix ret(_y, _x);

        for (int i = 0; i < _x; ++i) {
            for (int j = 0; j < _y; ++j) {
                ret[j][i] = _matrix[i][j];
            }
        }

        return ret;
    }

    Matrix inverse() {
        Matrix temp(_x, _y * 2);

        for (int i = 0; i < _x; ++i) {
            for (int j = 0; j < _y; ++j) {
                temp[i][j] = _matrix[i][j];
            }
        }

        for (int i = 0; i < _x; ++i) {
            for (int j = _x; j < 2 * _x; ++j) {
                if (i == (j - _x)) {
                    temp[i][j] = 1.0;
                } else {
                    temp[i][j] = 0.0;
                }
            }
        }

        for (int i = 0; i < _x; ++i){
            for (int j = 0; j < _y; ++j){
                if (i != j ){
                    double ratio = temp[j][i] / temp[i][i];
                    for (int k = 0; k < 2 * _x; ++k){
                        temp[j][k] -= ratio * temp[i][k];
                    }
                }
            }
        }

        for (int i = 0; i < _x; ++i){
            double a = temp[i][i];
            for (int j = 0; j < 2 * _x; ++j){
                temp[i][j] /= a;
            }
        }

        Matrix ret(_x, _y);

        for (int i = 0; i < _x; ++i) {
            for (int j = 0; j < _y; ++j) {
                ret[i][j] = temp[i][j + _y];
            }
        }

        return ret;
    }

    void print() {
//        for (int i = 0; i < _x; ++i) {
//            for (int j = 0; j < _y; ++j) {
//                std::cout << _matrix[i][j] << ' ';
//            }
//            std::cout << std::endl;
//        }
//        std::cout << std::endl;
    }

    std::vector<double>& operator[](std::size_t idx) {
        return _matrix[idx];
    }
    const std::vector<double>& operator[](std::size_t idx) const {
        return _matrix[idx];
    }

    Matrix operator*(const Matrix& rhs) {
        assert(this->_y == rhs._x);

        Matrix ret(this->_x, rhs._y);

        for (int i = 0; i < this->_x; ++i) {
            for (int j = 0; j < rhs._y; ++j) {
                for (int k = 0; k < this->_y; ++k) {
                    ret[i][j] += (*this)[i][k] * rhs[k][j];
                }
            }
        }

        return ret;
    }

private:
    int _x;
    int _y;
    std::vector<std::vector<double>> _matrix;
};

int main() {
    int m;
    int n;

    std::cin >> m >> n;
    Matrix X(n, m + 1);
    Matrix Y(n, 1);

    for (int i = 0; i < n; ++i) {
        X[i][0] = 1;
        for (int j = 1; j < m + 1; ++j) {
            std::cin >> X[i][j];
        }
        std::cin >> Y[i][0];
    }
    X.print();
    Y.print();

    Matrix XT = X.transpose();
    XT.print();

    Matrix XTX = XT * X;
    XTX.print();

    Matrix XTXI = XTX.inverse();
    XTXI.print();

    Matrix XTXIXT = XTXI * XT;
    XTXIXT.print();

    Matrix XTXIXTY = XTXIXT * Y;
    XTXIXTY.print();

    Matrix B(XTXIXTY);
    B.print();

    int q;
    std::cin >> q;
    Matrix Q(1, m + 1);
    Q[0][0] = 1;

    std::cout << std::fixed << std::setprecision(2);
    for (int i = 0; i < q; ++i) {
        for (int j = 1; j < m + 1; ++j) {
            std::cin >> Q[0][j];
        }
        Q.print();

        Matrix R = Q * B;
        R.print();

        std::cout << R[0][0] << std::endl;
    }

    return 0;
}

Problem solution in C programming.

#include <stdio.h>

int main() {
  int m;
  int n;
  int i;
  int j;
  int k;
  double X[100][10];
  long double XT[10][100];
  long double XT_X[10][10];
  long double XT_X_I[10][10];
  long double XT_X_I_XT[10][100];
  double Y[100];
  long double B[100];
  double f[10];
  long double AI[10][20];
  int q;
  double x;

  scanf("%d %d", &m, &n);
  for (i = 0; i < n; i++) {
    X[i][0] = 1;
    XT[0][i] = X[i][0];
    for (j = 1; j < (m + 1); j++) {
      scanf("%lf", &(X[i][j]));
      XT[j][i] = X[i][j];
    }
    scanf("%lf", &Y[i]);
  }

  for (i = 0; i < (m + 1); i++) {
    for (j = 0; j < (m + 1); j++) {
      XT_X[i][j] = 0;
      for (k = 0; k < n; k++) {
        XT_X[i][j] += XT[i][k] * X[k][j];
      }
    }
  }

  for (i = 0; i < (m + 1); i++) {
    for (j = 0; j < (m + 1); j++) {
      AI[i][j] = XT_X[i][j];
    }
    for (j = (m + 1); j < (m + 1) * 2; j++) {
      AI[i][j] = 0;
      if ((i + (m + 1)) == j) {
        AI[i][j] = 1;
      }
    }
  }

  for (i = 0; i < (m + 1); i++) {
    if (AI[i][i] != 1) {
      x = AI[i][i];
      for (j = 0; j < (m + 1) * 2; j++) {
        AI[i][j] /= x;
      }
    }

    for (j = 0; j < (m + 1); j++) {
      if (i == j) {
        continue;
      }

      if (AI[j][i] != 0) {
        x = AI[j][i] / AI[i][i];
        for (k = 0; k < (m + 1) * 2; k++) {
          AI[j][k] -= AI[i][k] * x;
        }
      }
    }
  }

  for (i = 0; i < (m + 1); i++) {
    for (j = 0; j < (m + 1); j++) {
      XT_X_I[i][j] = AI[i][j + (m + 1)];
    }
  }

  for (i = 0; i < (m + 1); i++) {
    for (j = 0; j < n; j++) {
      XT_X_I_XT[i][j] = 0;
      for (k = 0; k < (m + 1); k++) {
        XT_X_I_XT[i][j] += XT_X_I[i][k] * XT[k][j];
      }
    }
  }

  for (i = 0; i < (m + 1); i++) {
    B[i] = 0;
    for (k = 0; k < n; k++) {
      B[i] += XT_X_I_XT[i][k] * Y[k];
    }
  }

  scanf("%d", &q);
  for (i = 0; i < q; i++) {
    Y[i] = B[0];
    for (j = 1; j < (m + 1); j++) {
      scanf("%lf", &f[j]);
      Y[i] += B[j] * f[j];
    }
    Y[i] += 0.0049;
    printf("%.2fn", Y[i]);
  }

  return 0;
}

Problem solution in JavaScript programming.

function multipleRegression(x, y, q) {

  /* Gabriel Giordano's implementation from C to JavaScript with some modifications
     on the research "An Efficient and Simple Algorithm for Matrix Inversion" by Ahmad Farooq and Khan Hamid */
  function matrixInverse(a) {
    let pivot = 0

    for (let p = a.length; p-- > 0;) {
      pivot = a[p][p]

      if (Math.abs(pivot) < 1e-5)
        return 0

      for (let i = a.length; i-- > 0;)
        a[i][p] /= -pivot

      for (let i = a.length; i-- > 0;)
        if (i != p)
          for (let j = a.length; j-- > 0;)
            if (j != p)
              a[i][j] += a[p][j] * a[i][p]

      for (let j = a.length; j-- > 0;)
        a[p][j] /= pivot

      a[p][p] = 1 / pivot
    }

    return a
  }

  function matrixTranspose(a) {
    return Object.keys(a[0]).map((c) => a.map((r) => r[c]))
  }

  function matrixMultiply(a, b) {
    return a.map(x => matrixTranspose(b).map(y => dotProduct(x, y)))
  }

  function dotProduct(a, b) {
    return a.map((x, i) => a[i] * b[i]).reduce((m, n) => m + n)
  }

  let t = matrixTranspose(x),
    b = matrixMultiply(matrixMultiply(matrixInverse(matrixMultiply(t, x)), t), y)

  return dotProduct([1, ...q], b)
}

function processData(input) {
  //Enter your code here
  input = input.split('n')

  let [m, n] = input[0].split(' ').map(Number),
    x = [],
    y = []

  for (let i = n + 1; i-- > 1;) {
    let e = input[i].split(' ').map(Number)
    let xe = [1]
    for (let j = m + 1; j-- > 1;) {
      xe[j] = e[j - 1]
    }
    x.push(xe)
    y.push([e[m]])
  }
  for (let i = input.length - input[n + 1]; i < input.length; ++i) {
    let f = input[i].split(' ').map(Number)
    console.log(multipleRegression(x, y, f).toFixed(2))
  }
}

process.stdin.resume();
process.stdin.setEncoding("ascii");
_input = "";
process.stdin.on("data", function (input) {
    _input += input;
});

process.stdin.on("end", function () {
   processData(_input);
});