# -*- coding: utf-8 -*-
"""hw1_310555029.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1JRZxfH5avZg0896ddvxZSfm2lUtMGKSu

#Linear Regression

##(0) Load the data  -> Standardize data -> Split the dataset
"""

import numpy as np
import pandas as pd

data_x_df = pd.read_csv('data_X.csv')
data_t_df = pd.read_csv('data_T.csv')

X = np.array(data_x_df)
T = np.array(data_t_df)

# standardize 
X = (X-X.mean(axis=0))/X.std(axis=0)
#T = (T-T.mean(axis=0))/T.std(axis=0)

# 80% training dataset and 20% test dataset
X_SIZE = X.shape[0]
ratio = 8/10
trainRange = int(X_SIZE*ratio) #num of training data
x_train = X[:trainRange]
x_test = X[trainRange:] 
y_train = T[:trainRange]
y_test = T[trainRange:]

print("x_train shape {} and size {}".format(x_train.shape,x_train.size))
print("y_train shape {} and size {}".format(y_train.shape,y_train.size))
print("x_test shape {} and size {}".format(x_test.shape,x_test.size))
print("y_test shape {} and size {}".format(y_test.shape,y_test.size))

"""##(1) Feature selection

### Define model
"""

class LinearRegression:
    def __init__(self,x,M):
      self.M = M
      self.size = x.shape[0]
      self.feature_num = x.shape[1]
      self.weight = None

    def closed_form(self, x, y, _lambda=0, regularize = False):
      if not regularize:
      # Using Normal Equation to find w = (xT x)^(-1) xT y
        xTx = np.dot(x.T, x)
        inv = np.linalg.inv(np.matrix(xTx))
        w = np.dot(np.dot(inv, x.T), y)
        self.weight = w
      else:
      # w = (lambdaI + xT x)^(-1) xT y
        xTx = np.dot(x.T, x)
        inv = np.linalg.inv(np.eye(xTx.shape[0]) * _lambda + xTx)
        w = np.dot(np.dot(inv, x.T), y)
        self.weight = w
        
    def predict(self, x):
      pre = np.dot(x, self.weight) 
      return pre

    def poly_features(self,x):
      #PolynomialFeatures
      # find the phi matrix
      # array_new = [[1 x1 x1^2 ...]
      #        [1 x2 x1^2 ...]
      #          .
      #          .
      #        [1 xn xn^2 ...]] 
      if(self.M == 1):
        temp_matrix1 = np.ones([x.shape[0],1])
        temp_matrix2 = x[:x.shape[0],:]
        phi_matrix = np.hstack((temp_matrix1,temp_matrix2))
        return phi_matrix

      if(self.M == 2):
        temp_matrix1 = np.ones([x.shape[0],1])
        temp_matrix2 = x[:x.shape[0],:]
        phi_matrix = np.hstack((temp_matrix1,temp_matrix2))
        position = 1 + self.feature_num
        for i in range(x.shape[1]):
          for j in range(i,x.shape[1]):
            if i <= j:
              phi_matrix = np.insert(phi_matrix, position, (x[:,i]*x[:,j]), axis=1)
              position+=1
        return phi_matrix
  
def RMSE(y_predict,y_test):
    y_predict = np.squeeze(np.asarray(y_predict)) # numpy.matrix to numpy.ndarray
    y_test = np.squeeze(np.asarray(y_test))
    return np.sqrt(np.mean((y_predict - y_test)**2))

"""####(a) In the feature selection stage, please apply polynomials of order M = 1 and M = 2 over the dimension D = 8 input data. Please evaluate the corresponding RMS error on the training set and valid set. (15%) Code Result"""

print("M = ",1)
model_1 = LinearRegression(x_train,1)
x_new = model_1.poly_features(x_train)
model_1.closed_form(x_new,y_train) #find weight
# RMSE error
y_predict = model_1.predict(x_new)
train_rmse = RMSE(y_predict,y_train)
print("training RMSE: {}".format(train_rmse))
x_new = model_1.poly_features(x_test)
y_predict = model_1.predict(x_new)
test_rmse = RMSE(y_predict,y_test)
print("testing RMSE: {}".format(test_rmse))

print("M = ",2)
model_2 = LinearRegression(x_train,2)
x_new = model_2.poly_features(x_train)
model_2.closed_form(x_new,y_train) #find weight
# RMSE error
y_predict = model_2.predict(x_new)
train_rmse = RMSE(y_predict,y_train)
print("training RMSE: {}".format(train_rmse))
x_new = model_2.poly_features(x_test)
y_predict = model_2.predict(x_new)
test_rmse = RMSE(y_predict,y_test)
print("testing RMSE: {}\n".format(test_rmse))

"""####(b)How will you analysis the weights of polynomial model M = 1 and select the most contributive feature? Code Result, Explain (10%)

**The most contributive feature is median income**
"""

for i in range(len(features)):
  x = x_train
  print("【Without ",features[i],"】")
  x = np.delete(x, i, 1)
  model = LinearRegression(x,1)
  x_new = model_1.poly_features(x)
  model.closed_form(x_new,y_train) #find weight
  # RMSE error
  y_predict = model.predict(x_new)
  train_rmse = RMSE(y_predict,y_train)
  print("training RMSE: {}".format(train_rmse))
  x = x_test
  x = np.delete(x, i, 1)
  x_new = model.poly_features(x)
  y_predict = model.predict(x_new)
  test_rmse = RMSE(y_predict,y_test)
  print("testing RMSE: {}\n".format(test_rmse))

"""##(2) Maximum likelihood approach

####(b) Introduce the basis function you just decided in (a) to linear regression model and analyze the result you get.

使用 **sigmoid basis function** 做轉換，再進行linear regression分析
"""

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

X_sigmoid = X
T_sigmoid = T

# sigmoid 
X_sigmoid = sigmoid(X_sigmoid) 

# 80% training dataset and 20% test dataset
X_SIZE = X.shape[0]
ratio = 8/10
trainRange = int(X_SIZE*ratio) #num of training data
x_train = X_sigmoid[:trainRange]
x_test = X_sigmoid[trainRange:]
y_train = T[:trainRange]
y_test = T[trainRange:]

print("【Using sigmoid function】")
M=2
for i in range(1,M+1):
  print("M = ",i)
  model = LinearRegression(x_train,i)
  x_new = model.poly_features(x_train)
  model.closed_form(x_new,y_train) #find weight
  # RMSE error
  x_new = model.poly_features(x_train)
  y_predict = model.predict(x_new)
  train_rmse = RMSE(y_predict,y_train)
  print("training RMSE: {}".format(train_rmse))
  x_new = model.poly_features(x_test)
  y_predict = model.predict(x_new)
  test_rmse = RMSE(y_predict,y_test)
  print("testing RMSE: {}".format(test_rmse))

"""#### (c)**N-fold cross-validation** in your training stage to select at least one hyperparameter for model

"""

X_sigmoid = X
T_sigmoid = T

# sigmoid 
X_sigmoid = sigmoid(X_sigmoid) 

X_SIZE = X_sigmoid.shape[0]
N = 10
ratio = 1/N
testRange = int(X_SIZE*ratio) 
M=2
print("【N-fold cross-validation】")
for i in range(1,M+1):
  train_rmse = []
  test_rmse = []
  for r in range(N):
    # data split into 1 set(test) and N-1 sets(train)
    x_test = X_sigmoid[range((testRange*r),(testRange*(r+1))),:]
    y_test = T_sigmoid[range((testRange*r),(testRange*(r+1))),:]
    x_train = np.delete(X_sigmoid, [range((testRange*r),(testRange*(r+1)))], axis=0)
    y_train = np.delete(T_sigmoid, [range((testRange*r),(testRange*(r+1)))], axis=0)

    model = LinearRegression(x_train,i)
    x_new = model.poly_features(x_train)
    model.closed_form(x_new,y_train) #find weight
    # train RMSE error
    x_new = model.poly_features(x_train)
    y_predict = model.predict(x_new)
    train_rmse.append(RMSE(y_predict,y_train))
    # test RMSE error
    x_new = model.poly_features(x_test)
    y_predict = model.predict(x_new)
    test_rmse.append(RMSE(y_predict,y_test))

  average_train_rmse = sum(train_rmse) / N
  average_test_rmse = sum(test_rmse) / N
  print(f'M = {i}, Average_train_rmse = {average_train_rmse}, Average_test_rmse = {average_test_rmse}')

"""## (3) Maximum a *posteriori* approach

####(b) Use maximum a posteriori approach method to retest the model in 2 you designed. You could choose Gaussian distribution as a prior.

使用 Maximum a posteriori 做 linear regression , 設Lambda = 0.1
"""

_lambda = 28

X_sigmoid = X
T_sigmoid = T

# sigmoid 
X_sigmoid = sigmoid(X_sigmoid)

# 80% training dataset and 20% test dataset
X_SIZE = X_sigmoid.shape[0]
ratio = 8/10
trainRange = int(X_SIZE*ratio) #num of training data
x_train = X_sigmoid[:trainRange]
x_test = X_sigmoid[trainRange:]
y_train = T_sigmoid[:trainRange]
y_test = T_sigmoid[trainRange:]

M=2
for i in range(1,M+1):
  print("M = ",i)
  model = LinearRegression(x_train,i)
  x_new = model.poly_features(x_train)
  model.closed_form(x_new,y_train,_lambda,True) #find weight
  # RMSE error
  x_new = model.poly_features(x_train)
  y_predict = model.predict(x_new)
  train_rmse = RMSE(y_predict,y_train)
  print("training RMSE: {}".format(train_rmse))
  x_new = model.poly_features(x_test)
  y_predict = model.predict(x_new)
  test_rmse = RMSE(y_predict,y_test)
  print("testing RMSE: {}".format(test_rmse))

"""使用 Maximum a posteriori 做 N-fold cross-validation """

X_sigmoid = X
T_sigmoid = T

# sigmoid 
X_sigmoid = sigmoid(X_sigmoid)

X_SIZE = X_sigmoid.shape[0]
N = 10
ratio = 1/N
testRange = int(X_SIZE*ratio) 
M=2

for i in range(1,M+1):
  train_rmse = []
  test_rmse = []
  for r in range(N):
    # data split into 1 set(test) and N-1 sets(train)
    x_test = X_sigmoid[range((testRange*r),(testRange*(r+1))),:]
    y_test = T_sigmoid[range((testRange*r),(testRange*(r+1))),:]
    x_train = np.delete(X_sigmoid, [range((testRange*r),(testRange*(r+1)))], axis=0)
    y_train = np.delete(T_sigmoid, [range((testRange*r),(testRange*(r+1)))], axis=0)

    model = LinearRegression(x_train,i)
    x_new = model.poly_features(x_train)
    model.closed_form(x_new,y_train,_lambda,True) #find weight
    # train RMSE error
    x_new = model.poly_features(x_train)
    y_predict = model.predict(x_new)
    train_rmse.append(RMSE(y_predict,y_train))
    # test RMSE error
    x_new = model.poly_features(x_test)
    y_predict = model.predict(x_new)
    test_rmse.append(RMSE(y_predict,y_test))

  average_train_rmse = sum(train_rmse) / N
  average_test_rmse = sum(test_rmse) / N
  print(f'『M = {i}』 Average_train_rmse = {average_train_rmse}, Average_test_rmse = {average_test_rmse}')