# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: PyCharm (Project PRML)
#     language: python
#     name: pycharm-ddcd9f7a
# ---

# # PRML Chapter 1
# please follow book along with this notebook.

import numpy as np
import matplotlib.pyplot as plt


# ## Figure 1.2 - 
# generate sample data of 10 points
# due to random noise, you might not find the plot exactly same as book.

# + code_folding=[0, 3, 6, 17] pycharm={"name": "#%%\n"}
def generate_x(size):
    return np.linspace(0,1,size)

def true_function(x):
    return np.sin(np.pi * 2 * x)

def noise(mean, sigma, size):
    return np.random.normal(loc=mean, scale=sigma, size=size)

def generate_data(no_of_sample):
    """ figure 1.2 has single dimension data.
    data generation rule : y = sin(x) + gaussian noise with mean 0, variance 0.25
    """
    x = generate_x(no_of_sample)
    t = true_function(x) + noise(0, 0.25, no_of_sample)
    return x, t

def plot_fig_1_2():

    # training data for scatter plot
    x, t = generate_data(10)

    # test data for green line, the true function without noise.
    x_test = np.linspace(0,1,100)
    y_test = true_function(x_test)

    # plot
    plt.scatter(x, t, marker= "o", facecolor="none",color = 'b', label = "training data")
    plt.plot(x_test, y_test, color='g', label='true function')

    plt.xlabel("x ( feature dimension )")
    plt.ylabel("t ( labels / training label ) ")
    plt.legend()
    plt.show()

#----------------------
# plot now
#----------------------
plot_fig_1_2()


# -

# ## Figure 1.3 :
# The error function (1.2) corresponds to the sum of the squares of the displacements of each data from the funcion y(x,w)
#
# E(w) = 0.5 * sum over data points { y(x_n, w) - t_n } ** 2

# + code_folding=[0, 3] pycharm={"name": "#%%\n"}
def prediction_function(x):
    return x

def plot_fig_1_3():

    # the prediction function y(x,w)
    sample_size =5
    x = generate_x(sample_size)
    y = prediction_function(x)
    plt.plot(x, y, color="r")

    # point 1
    x = np.array([0.2,0.2])
    y = np.array([0.2,0.8])
    plt.plot(x, y, color='g')
    plt.scatter(x, y, color="b")

    # point 2
    x = np.array([0.4,0.4])
    y = np.array([0.4,0.3])
    plt.plot(x, y, color='g')
    plt.scatter(x, y, color="b")

    # point 3
    x = np.array([0.8,0.8])
    y = np.array([0.8, 1])
    plt.plot(x, y, color='g')
    plt.scatter(x, y, color="b")

    plt.text(0.8,0.7, "y(x_n, w)")
    plt.text(0.82, 1, "t_n")

    plt.xlabel("x")
    plt.ylabel("t")
    plt.show()
    plt.close()
plot_fig_1_3()


# -

# ## Figure 1.4
# a polynomial regression on multiple M's : 0, 1, 3, 9

# + code_folding=[] pycharm={"name": "#%%\n"}
def convert_to_polynomial(x, degree):
    assert x.ndim == 1, " this is written for conversion of only single dimension x"
    assert  degree >= 0, "degree must be positive"
    x = x[:,np.newaxis]

    N = x.shape[0]

    # X will be a matrix of shape N * (degree + 1)
    X = np.ones(shape=(N, 1+degree))

    for degree_dim in range(1, degree+1):
        X[:, degree_dim] = np.power ( x[:,0], degree_dim )
        
    return X 


# + code_folding=[4]
x_train, y_train = generate_data(10)
x_test = generate_x(100)
y_test = true_function(x_test)

print(f'''
        x_train : {x_train.shape}
        y_train : {y_train.shape}
        x_test  : {x_test.shape}
        y_test  : {y_test.shape}
        ''')


# + code_folding=[]
def find_optimal_weights( X, y):
    """ optimal weights via simple inversion methods. 
    not recommended for large datasets due to matrix inversion complexity. """
    if y.ndim == 1:
        y = y[:, np.newaxis ]        
    W = np.linalg.inv ( X.T.dot(X) ).dot ( X.T ).dot(y)
    return W 


# + code_folding=[]
def predict( x, W):
    """
    W : Features * 1 matrix
    X : N * features matrix
    forward : X * W => N * 1 matrix 
    """
    if x.ndim == 1:
        X = X[:, np.newaxis]
    return x.dot(W)    


# +
fig, axes = plt.subplots(2, 2, figsize=(20,20))
degrees = [0, 1, 3, 9]

for i in range(0,2):
    for j in range(0,2):
        
        degree = degrees[ 2 * i + j ]
        
        # convert to the polynomial 
        X = convert_to_polynomial(x_train, degree)

        # find the weights 
        W = find_optimal_weights(X, y_train)

        # predict 
        X = convert_to_polynomial(x_test, degree)
        y_pred = predict ( X, W ) 

        print(f"weights found are {W.squeeze()}")

        # plots
        axes[i][j].scatter(x_train, y_train, color="b", facecolor="none", label = "train data")
        axes[i][j].plot(x_test, y_test, color="g", label = "true function")
        axes[i][j].plot(x_test, y_pred.squeeze(), color="r", label = f"predicted fn of degree {degree}")
        axes[i][j].legend()
plt.show()


# + [markdown] heading_collapsed=true
# ## Figure 1.5
# generate graph of root mean square error for for training and test data for degrees 0 to 9 inclusive. 

# + hidden=true
def calculateRMS( y_predicted, y_true ):

    if y_predicted.ndim > 1 :
        y_predicted = y_predicted.squeeze() 
        y_true = y_true.squeeze()    
        
    temp = np.sqrt( np.mean( (y_predicted - y_true)**2, axis = 0) ) 
    return temp


# + hidden=true
calculateRMS(y_pred, y_test)

# + hidden=true
# fig, axes = plt.subplots(2, 2, figsize=(20,20))
degrees = list ( np.arange(0,10) )
rms_train = []
rms_test = [] 

for degree in degrees:
    # convert to the polynomial 
    X = convert_to_polynomial(x_train, degree)
    # find the weights 
    W = find_optimal_weights(X, y_train)
    y_predicted_on_train = predict(X, W)
    rms_train.append ( calculateRMS (y_predicted_on_train, y_train) )
    
    # predict 
    X = convert_to_polynomial(x_test, degree)
    y_predicted_on_test = predict ( X, W ) 
    rms_test.append( calculateRMS(y_predicted_on_test, y_test) )
    
# plots
print
plt.scatter(degrees, rms_train, color = "b", facecolor="none")
plt.plot(degrees, rms_train, color = "b", label = "train rms")

plt.scatter(degrees, rms_test, color = "r", facecolor="none")
plt.plot(degrees, rms_test, color = "r", label = "test rms")

plt.xlabel("M")
plt.ylabel("ERMS")

plt.legend()
plt.show()

# + [markdown] heading_collapsed=true
# ## Figure 1.6 
# Effect of datasize on overfitting for M = 9 

# + hidden=true
degree = 9 
x_test = generate_x(100)
y_test = true_function(x_test)

for N in [15, 100]:
    
    x_train, y_train = generate_data(N)
   
    # convert to the polynomial 
    X = convert_to_polynomial(x_train, degree)

    # find the weights 
    W = find_optimal_weights(X, y_train)

    # predict 
    X = convert_to_polynomial(x_test, degree)
    y_pred = predict ( X, W ) 

#     print(f"weights found are {W.squeeze()}")

    # plots
    plt.figure(figsize = (10,5))
    plt.scatter(x_train, y_train, color="b", facecolor="none", label = "train data")
    plt.plot(x_test, y_test, color="g", label = "true function")
    plt.plot(x_test, y_pred.squeeze(), color="r", label = f"predicted fn of degree {degree} and training N = {N}")
    
    plt.legend()
    plt.show()


# + [markdown] heading_collapsed=true
# ## Figure 1.7 and 1.8

# + hidden=true
def find_optimal_weights_with_regularisation( X, y, lambda_ ):
    """ optimal weights via simple inversion methods. 
    not recommended for large datasets due to matrix inversion complexity. """
    if y.ndim == 1:
        y = y[:, np.newaxis ]
    
    temp = X.T.dot(X)
    I = np.eye(temp.shape[0])
    
    W = np.linalg.inv ( lambda_ * I + temp ).dot ( X.T ).dot(y)
    return W 


# + [markdown] heading_collapsed=true hidden=true
# ### Figure 1.7 
# Effect of regularisation on overfitting for N=10, M=9

# + hidden=true
degree = 9 
N = 10
x_test = generate_x(100)
y_test = true_function(x_test)

for lambda_ in [np.exp(-18), 1]:
    
    x_train, y_train = generate_data(N)
   
    # convert to the polynomial 
    X = convert_to_polynomial(x_train, degree)

    # find the weights 
    W = find_optimal_weights_with_regularisation(X, y_train, lambda_)

    # predict 
    X = convert_to_polynomial(x_test, degree)
    y_pred = predict ( X, W ) 

#     print(f"weights found are {W.squeeze()}")

    # plots
    plt.figure(figsize = (10,5))
    plt.scatter(x_train, y_train, color="b", facecolor="none", label = "train data")
    plt.plot(x_test, y_test, color="g", label = "true function")
    plt.plot(x_test, y_pred.squeeze(), color="r", label = f"predicted for N={N}, degree={degree}, lambda={lambda_}")
    
    plt.legend()
    plt.show()


# + [markdown] hidden=true
# ### Figure 1.8 
# effect of lambda on ERMS for a range of lambda's 

# + hidden=true
# fig, axes = plt.subplots(2, 2, figsize=(20,20))
lambdas = [np.exp(x) for x in range(-40, 5, 5)]
degree = 9 

rms_train = []
rms_test = [] 

for lambda_ in lambdas:
    x_train, y_train = generate_data(10)
    
    # convert to the polynomial 
    X = convert_to_polynomial(x_train, degree)
    # find the weights 
    W = find_optimal_weights_with_regularisation(X, y_train, lambda_)
    
    y_predicted_on_train = predict(X, W)
    rms_train.append ( calculateRMS (y_predicted_on_train, y_train) )
    
    # predict 
    X = convert_to_polynomial(x_test, degree)
    y_predicted_on_test = predict ( X, W ) 
    rms_test.append( calculateRMS(y_predicted_on_test, y_test) )
    
# plots
print
plt.scatter(np.log(lambdas), rms_train, color = "b", facecolor="none")
plt.plot(np.log(lambdas), rms_train, color = "b", label = "train rms")

plt.scatter(np.log(lambdas), rms_test, color = "r", facecolor="none")
plt.plot(np.log(lambdas), rms_test, color = "r", label = "test rms")

plt.xlabel("ln lambda")
plt.ylabel("ERMS")

plt.legend()
plt.show()

# + [markdown] hidden=true
#

StarcoderdataJupyterScriptsDedupFiltered