# 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("%.2f\n", 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) {
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);
});