#!/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
from scipy.optimize import minimize
#==============================================================================
# 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)
def f_obj_1(beta,x,y):
val = np.mean(np.square(y-(beta[0]+beta[1]*x**2*np.sin(beta[2]*x+beta[3]))))
return val
def f_obj_2(beta,x,y):
val = np.mean(np.square(y-(beta[0]*x+beta[1])**2*np.sin(beta[2]*x+beta[3])))
return val
def inferModel1(x,beta):
return beta[0]+beta[1]*x**2*np.sin(beta[2]*x+beta[3])
def inferModel2(x,beta):
return (beta[0]*x+beta[1])**2*np.sin(beta[2]*x+beta[3])
def optimizeModel1(x):
y = TrueFunction(x)
# res = minimize(f_obj_1,[1,0.24,0.5,0.1],args=(x,y))
res = minimize(f_obj_1,np.random.randint(-5,15,4),args=(x,y))
return res.x, y
def optimizeModel2(x):
y = TrueFunction(x)
# res = minimize(f_obj_2,[1,0.24,0.5,0.1],args=(x,y))
res = minimize(f_obj_2,np.random.randint(-5,15,4),args=(x,y))
return res.x, y
x = np.linspace(0,1,101)
beta_5_2,y_sample_5_2 = optimizeModel1(ES_5)
y_pred_5_2 = inferModel1(x,beta_5_2)
beta_10_2,y_sample_10_2 = optimizeModel1(ES_10)
y_pred_10_2 = inferModel1(x,beta_10_2)
beta_15_2,y_sample_15_2 = optimizeModel1(ES_15)
y_pred_15_2 = inferModel1(x,beta_15_2)
beta_5_3,y_sample_5_3 = optimizeModel2(ES_5)
y_pred_5_3 = inferModel2(x,beta_5_3)
beta_10_3,y_sample_10_3 = optimizeModel2(ES_10)
y_pred_10_3 = inferModel2(x,beta_10_3)
beta_15_3,y_sample_15_3 = optimizeModel2(ES_15)
y_pred_15_3 = inferModel2(x,beta_15_3)
# 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")
b0 = str(np.round(beta_5_2[0],2))
b1 = str(np.round(beta_5_2[1],2))
b2 = str(np.round(beta_5_2[2],2))
b3 = str(np.round(beta_5_2[3],2))
ax[0].set_title(r"N = 5, $\beta$=["+b0+", "+b1+", "+b2+", "+b3+"]")
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")
b0 = str(np.round(beta_10_2[0],2))
b1 = str(np.round(beta_10_2[1],2))
b2 = str(np.round(beta_10_2[2],2))
b3 = str(np.round(beta_10_2[3],2))
ax[1].set_title(r"N = 10, $\beta$=["+b0+", "+b1+", "+b2+", "+b3+"]")
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")
b0 = str(np.round(beta_15_2[0],2))
b1 = str(np.round(beta_15_2[1],2))
b2 = str(np.round(beta_15_2[2],2))
b3 = str(np.round(beta_15_2[3],2))
ax[2].set_title(r"N = 15, $\beta$=["+b0+", "+b1+", "+b2+", "+b3+"]")
# plt.savefig("model1.png",dpi=150,bbox_inches="tight")
fig.suptitle(r"$\hat{y} = \beta_1+\beta_2 x^2\sin(\beta_3 x+\beta_4)$")
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")
b0 = str(np.round(beta_5_3[0],2))
b1 = str(np.round(beta_5_3[1],2))
b2 = str(np.round(beta_5_3[2],2))
b3 = str(np.round(beta_5_3[3],2))
ax[0].set_title(r"N = 5, $\beta$=["+b0+", "+b1+", "+b2+", "+b3+"]")
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")
b0 = str(np.round(beta_10_3[0],2))
b1 = str(np.round(beta_10_3[1],2))
b2 = str(np.round(beta_10_3[2],2))
b3 = str(np.round(beta_10_3[3],2))
ax[1].set_title(r"N = 10, $\beta$=["+b0+", "+b1+", "+b2+", "+b3+"]")
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")
b0 = str(np.round(beta_15_3[0],2))
b1 = str(np.round(beta_15_3[1],2))
b2 = str(np.round(beta_15_3[2],2))
b3 = str(np.round(beta_15_3[3],2))
ax[2].set_title(r"N = 15, $\beta$=["+b0+", "+b1+", "+b2+", "+b3+"]")
# plt.savefig("model2.png",dpi=150,bbox_inches="tight")
fig.suptitle(r"$\hat{y} = (\beta_1x+\beta_2)^2\sin(\beta_3 x+\beta_4)$")
plt.show()
#==============================================================================