#!/bin/python3
"""============================================================================
sine basis 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
from numpy.polynomial import chebyshev as C
#==============================================================================
# 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
# chebyshev regression
def chebyshevRegression(x,degree):
y = TrueFunction(x)
c,_ = C.chebfit(x,y,degree, full=True)
cheb_series = C.Chebyshev(c,domain=[-1,1])
y_pred = cheb_series(x)
return c,y_pred
def chebyshevInference(x,c):
cheb_series = C.Chebyshev(c,domain=[-1,1])
y_pred = cheb_series(x)
return y_pred
# sine basis regression
def sineRegression(x,degree):
y = TrueFunction(x)
X = np.ones_like(x)
for i in range(1,degree+1):
X = np.c_[X,np.sin(np.pi*i*x).T]
beta = np.linalg.inv(X.T@X)@X.T@y
y_sample = X@beta
return beta,y_sample
def sineInference(x,degree,beta):
y = TrueFunction(x)
X = np.ones_like(x)
for i in range(1,degree+1):
X = np.c_[X,np.sin(np.pi*i*x).T]
y_pred = X@beta
return y_pred
x = np.linspace(0,1,101)
beta_5_2,y_sample_5_2 = sineRegression(ES_5,2)
y_pred_5_2 = sineInference(x,2,beta_5_2)
beta_10_2,y_sample_10_2 = sineRegression(ES_10,2)
y_pred_10_2 = sineInference(x,2,beta_10_2)
beta_15_2,y_sample_15_2 = sineRegression(ES_15,2)
y_pred_15_2 = sineInference(x,2,beta_15_2)
beta_5_3,y_sample_5_3 = sineRegression(ES_5,3)
y_pred_5_3 = sineInference(x,3,beta_5_3)
beta_10_3,y_sample_10_3 = sineRegression(ES_10,3)
y_pred_10_3 = sineInference(x,3,beta_10_3)
beta_15_3,y_sample_15_3 = sineRegression(ES_15,3)
y_pred_15_3 = sineInference(x,3,beta_15_3)
beta_5_4,y_sample_5_4 = sineRegression(ES_5,4)
y_pred_5_4 = sineInference(x,4,beta_5_4)
beta_10_4,y_sample_10_4 = sineRegression(ES_10,4)
y_pred_10_4 = sineInference(x,4,beta_10_4)
beta_15_4,y_sample_15_4 = sineRegression(ES_15,4)
y_pred_15_4 = sineInference(x,4,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")
# plt.savefig("sine_2.png",dpi=150,bbox_inches="tight")
fig.suptitle(r"(1,sin($\pi$,x),sin(2$\pi$,x))")
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")
# plt.savefig("sine_3.png",dpi=150,bbox_inches="tight")
fig.suptitle(r"(1,sin($\pi$,x),sin(2$\pi$,x),sin(3$\pi$,x))")
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") # check point
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,sin($\pi$,x),sin(2$\pi$,x),sin(3$\pi$,x),sin(4$\pi$,x))")
# plt.savefig("sine_4.png",dpi=150,bbox_inches="tight")
plt.show()
#==============================================================================