#!/bin/python3
"""============================================================================
polynomial linear regression script
Ramkumar
Sun Mar 23 05:11:57 PM IST 2025
============================================================================"""
# importing needed modules
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import qmc
#==============================================================================
# function definitions for sampling
def RandomSampling(n):
return np.random.sample(n)
def EquispacedSampling(n):
return np.linspace(0,1,n)
def LatinHypercubeSampling(n):
sampler = qmc.LatinHypercube(d=1)
return sampler.random(n).flatten()
def TrueFunction(x):
return (6*x-2)**2*np.sin(12*x-4)
# sampling data
N = 5
RS_5 = RandomSampling(N)
ES_5 = EquispacedSampling(N)
LHS_5 = LatinHypercubeSampling(N)
N = 10
RS_10 = RandomSampling(N)
ES_10 = EquispacedSampling(N)
LHS_10 = LatinHypercubeSampling(N)
N = 15
RS_15 = RandomSampling(N)
ES_15 = EquispacedSampling(N)
LHS_15 = LatinHypercubeSampling(N)
# linear regression using polynomial of given degree
def linearRegression(x,degree):
X = np.vander(x,degree+1,True)
y = TrueFunction(x)
beta = np.linalg.inv(X.T@X)@X.T@y
y_sample = X@beta
return beta,y_sample
x = np.linspace(0,1,101)
beta_5_2,y_sample_5_2 = linearRegression(ES_5,2)
y_pred_5_2 = np.vander(x,2+1,True)@beta_5_2
beta_10_2,y_sample_10_2 = linearRegression(ES_10,2)
y_pred_10_2 = np.vander(x,2+1,True)@beta_10_2
beta_15_2,y_sample_15_2 = linearRegression(ES_15,2)
y_pred_15_2 = np.vander(x,2+1,True)@beta_15_2
beta_5_3,y_sample_5_3 = linearRegression(ES_5,3)
y_pred_5_3 = np.vander(x,3+1,True)@beta_5_3
beta_10_3,y_sample_10_3 = linearRegression(ES_10,3)
y_pred_10_3 = np.vander(x,3+1,True)@beta_10_3
beta_15_3,y_sample_15_3 = linearRegression(ES_15,3)
y_pred_15_3 = np.vander(x,3+1,True)@beta_15_3
beta_5_4,y_sample_5_4 = linearRegression(ES_5,4)
y_pred_5_4 = np.vander(x,4+1,True)@beta_5_4
beta_10_4,y_sample_10_4 = linearRegression(ES_10,4)
y_pred_10_4 = np.vander(x,4+1,True)@beta_10_4
beta_15_4,y_sample_15_4 = linearRegression(ES_15,4)
y_pred_15_4 = np.vander(x,4+1,True)@beta_15_4
# plotting graphs
plt.rcParams.update({"font.size":10})
fig,ax = plt.subplots(3,1,figsize=(12,6),sharex=True,sharey=True)
x = np.linspace(0,1,101)
ax = ax.flatten()
# degree 2
ax[0].plot(x,TrueFunction(x),'-k',label = "True function")
ax[0].plot(x,y_pred_5_2,'-r',label = "predicted")
ax[0].plot(ES_5,y_sample_5_2,'og',label="sample")
ax[0].grid()
# ax[0].set_xlabel("data points")
ax[0].set_ylabel("y")
ax[0].set_title("N = 5")
ax[1].plot(x,TrueFunction(x),'-k',label = "True function")
ax[1].plot(x,y_pred_10_2,'-r',label = "predicted")
ax[1].plot(ES_10,y_sample_10_2,'og',label="sample")
ax[1].grid()
ax[1].legend(loc=[1.01,0.5])
# ax[1].set_xlabel("data points")
ax[1].set_ylabel("y")
ax[1].set_title("N = 10")
ax[2].plot(x,TrueFunction(x),'-k',label = "True function")
ax[2].plot(x,y_pred_15_2,'-r',label = "predicted")
ax[2].plot(ES_15,y_sample_15_2,'og',label="sample")
ax[2].grid()
ax[2].set_xlabel("x")
ax[2].set_ylabel("y")
ax[2].set_title("N = 15")
fig.suptitle(r"$(1,x,x^2)$")
# plt.savefig("polynomial_2.png",dpi=150,bbox_inches="tight")
plt.show()
# degree 3
fig,ax = plt.subplots(3,1,figsize=(12,6),sharex=True,sharey=True)
x = np.linspace(0,1,101)
ax = ax.flatten()
ax[0].plot(x,TrueFunction(x),'-k',label = "True function")
ax[0].plot(x,y_pred_5_3,'-r',label = "predicted")
ax[0].plot(ES_5,y_sample_5_3,'og',label="sample")
ax[0].grid()
# ax[0].set_xlabel("data points")
ax[0].set_ylabel("y")
ax[0].set_title("N = 5")
ax[1].plot(x,TrueFunction(x),'-k',label = "True function")
ax[1].plot(x,y_pred_10_3,'-r',label = "predicted")
ax[1].plot(ES_10,y_sample_10_3,'og',label="sample")
ax[1].grid()
ax[1].legend(loc=[1.01,0.5])
# ax[1].set_xlabel("data points")
ax[1].set_ylabel("y")
ax[1].set_title("N = 10")
ax[2].plot(x,TrueFunction(x),'-k',label = "True function")
ax[2].plot(x,y_pred_15_3,'-r',label = "predicted")
ax[2].plot(ES_15,y_sample_15_3,'og',label="sample")
ax[2].grid()
ax[2].set_xlabel("x")
ax[2].set_ylabel("y")
ax[2].set_title("N = 15")
fig.suptitle(r"$(1,x,x^2,x^3)$")
# plt.savefig("polynomial_3.png",dpi=150,bbox_inches="tight")
plt.show()
# degree 4
fig,ax = plt.subplots(3,1,figsize=(12,6),sharex=True,sharey=True)
x = np.linspace(0,1,101)
ax = ax.flatten()
ax[0].plot(x,TrueFunction(x),'-k',label = "True function")
ax[0].plot(x,y_pred_5_4,'-r',label = "predicted")
ax[0].plot(ES_5,y_sample_5_4,'og',label="sample")
ax[0].grid()
# ax[0].set_xlabel("data points")
ax[0].set_ylabel("y")
ax[0].set_title("N = 5")
ax[1].plot(x,TrueFunction(x),'-k',label = "True function")
ax[1].plot(x,y_pred_10_4,'-r',label = "predicted")
ax[1].plot(ES_10,y_sample_10_4,'og',label="sample")
ax[1].grid()
ax[1].legend(loc=[1.01,0.5])
# ax[1].set_xlabel("data points")
ax[1].set_ylabel("y")
ax[1].set_title("N = 10")
ax[2].plot(x,TrueFunction(x),'-k',label = "True function")
ax[2].plot(x,y_pred_15_4,'-r',label = "predicted")
ax[2].plot(ES_15,y_sample_15_4,'og',label="sample")
ax[2].grid()
ax[2].set_xlabel("x")
ax[2].set_ylabel("y")
ax[2].set_title("N = 15")
fig.suptitle(r"$(1,x,x^2,x^3,x^4)$")
# plt.savefig("polynomial_4.png",dpi=150,bbox_inches="tight")
plt.show()
#==============================================================================