<jupyter_start><jupyter_text>Simulating - single feature, single label dataset  **Without standard scaling**<jupyter_code>import numpy as np
import pandas as pd
X = np.round(np.random.uniform(0,5, size= 1000), 3)

a, b  = np.random.normal(0, 0.1), np.random.normal(0, 0.2)#y_pred = a+ bx
error = np.random.normal(0, 1, size= 1000)
f= lambda x,e : np.round(a +b*x + e, 3)
data = pd.DataFrame({'features': X, 'targets':f(X, error)})

data.head(5)
data.to_csv("reg_sim_data.csv", index=None)<jupyter_output><empty_output><jupyter_text>Training statsmodel<jupyter_code>import statsmodels.api as sm
feats, targets= data['features'].values, data['targets'].values
feats= sm.add_constant(feats, prepend=False)
#Without standard scaling of data
model_stm = sm.OLS(targets, feats, hasconst=True)
result= model_stm.fit()
print(result.summary())
#Without standard scaling of data
print('Reference values from simulated dataset weight : {} & bias: {}'.format(b,a))<jupyter_output>Reference values from simulated dataset weight : 0.10458035217582123 & bias: 0.11878560346515568<jupyter_text>Training pytorch model<jupyter_code>import torch
import torch.nn as nn
from torch.autograd import Variable
import torch.optim as optim

class toynet(nn.Module):
    def __init__(self):
        super(toynet, self).__init__()
        self.fc1 = nn.Linear(1,1)
        torch.nn.init.xavier_normal_(self.fc1.weight)
    def forward(self, x):
        x= self.fc1(x)
        return x

model = toynet()
model
optimizer = optim.Adam(model.parameters())
loss = nn.MSELoss()
def sample(df, count):
    idx= np.random.randint(df.shape[0], size=count)
    x,y = df.values[idx][:,0].reshape(-1,1), df.values[idx][:,1]
    return x,y

def eval_hessian(loss_grad, model):
    cnt = 0
    for g in loss_grad:
        g_vector = g.contiguous().view(-1) if cnt == 0 else torch.cat([g_vector, g.contiguous().view(-1)])
        cnt = 1
    l = g_vector.size(0)
    hessian = torch.zeros(l, l)
    for idx in range(l):
        grad2rd = torch.autograd.grad(g_vector[idx], model.parameters(), create_graph=True)
        cnt = 0
        for g in grad2rd:
            g2 = g.contiguous().view(-1) if cnt == 0 else torch.cat([g2, g.contiguous().view(-1)])
            cnt = 1
        hessian[idx] = g2
    return hessian.cpu().data.numpy()
wts_bias_progression= list()
all_grads= list()
losses= list()

d2loss_dg= list()
diags= list()
f= lambda x: np.sqrt(np.abs(np.diag(np.linalg.inv(x)))).tolist()

for epoch in range(4000):
    x_, y_ = sample(data, count= 64)
    #x_,y_= data.values[:,0].reshape(-1,1), data.values[:,1]
    trainx = Variable(torch.Tensor(x_).float(), requires_grad= True)
    trainy = Variable(torch.reshape(torch.Tensor(y_).float(), (-1,1)))
    optimizer.zero_grad()
    y_pred= model(trainx)
    
    l = loss(y_pred, trainy)
    losses.append(l.item())
    if epoch%500==0:
        print('epoch: ', epoch, ' loss: ', l.item());
    
    l.backward()
    optimizer.step()
    
    wts_bias = list(model.parameters())
    wts_bias_vals= torch.cat((wts_bias[0], torch.reshape(wts_bias[1], (-1,1))),1)
    wts_bias_progression.append(wts_bias_vals.tolist())
    
    wts_bias_grads = torch.cat((wts_bias[0].grad, torch.reshape(wts_bias[1].grad, (-1,1))),1)
    #i_grads= trainx.grad
    all_grads.append(wts_bias_grads.tolist())#i_grads)
    d1loss_dg = torch.autograd.grad(loss(model(trainx), trainy), wts_bias, create_graph=True)
    d2loss_dg.append(eval_hessian(d1loss_dg, model))
    diags.append(f(d2loss_dg[-1]))

all_wts_bias = np.array(wts_bias_progression).reshape(-1,wts_bias_vals.shape[1])
all_wts_bias_df= pd.DataFrame(all_wts_bias, columns=['weight_{}'.format(idx) for idx in range(1,wts_bias_vals.shape[1])]+ ['bias'])

all_grads = np.array(all_grads).reshape(-1,wts_bias_vals.shape[1])
all_grads_df= pd.DataFrame(all_grads, columns=['beta_{}'.format(idx) for idx in range(1,wts_bias_vals.shape[1]+1)])
list(model.parameters())<jupyter_output><empty_output><jupyter_text>Analytical standard error & T-value for b (slope estimator)<jupyter_code>#S_b calc
sum_sqr_err = (1/(data['features'].shape[0] - 2))*(np.sum(error**2))
sum_sqr_diffx = np.sum((data['features'].values - np.mean(data['features'].values))**2)
s_beta = np.sqrt(sum_sqr_err/sum_sqr_diffx)
s_beta
#(1/(1000-2))*(np.sum(error**2))/np.sum((data['features'].values - np.mean(data['features'].values))**2)
b_pred = list(model.parameters())[0] 
t= (b_pred - b)/s_beta
t#t-val<jupyter_output><empty_output><jupyter_text>Calculating Pytorch Standar error<jupyter_code>mini_batches= list()#contains mini_batches of gradients from 200 epochs
mean_mini_batches= list()#contains mean of each of 200 epochs each in a mini_batches of gradients from 200 epochs
mini_size= 200
for batch in range(1, (len(all_grads))//mini_size):
    mini_batches.append(all_grads[batch*mini_size:batch*mini_size+ mini_size])
    mean_mini_batches.append(np.mean(mini_batches[-1], axis=0).tolist())

grad_mini_batch_arr = np.array(mean_mini_batches).T
grad_mini_batch_arr.shape
covariant_mat= np.cov(grad_mini_batch_arr)
pytorch_stderr= np.sqrt(covariant_mat.diagonal())#**2
pytorch_stderr<jupyter_output><empty_output><jupyter_text>Stasmodel vs Pytorch comparison<jupyter_code>summary= result.summary()
statsm_df = pd.DataFrame(summary.tables[1].data[1:], columns=summary.tables[1].data[0])
statsm_df
di = dict(zip('pytorch_grads,simulated_params,pytorch_params,statsmodel_params,statsmodel_stderr'.split(','), [all_grads[-1,:], [b,a], all_wts_bias[-1,:], result.params, statsm_df["std err"].values]))#all_grads
di.update({'pytorch_stderr':pytorch_stderr})
comp_df = pd.DataFrame(di, index=['b (slope/weight)','a (intercept/bias)'])
#comp_df['pytorch_stderr']= pytorch_stderr
comp_df['pytorch_params'] = comp_df['pytorch_params'].apply(lambda x:np.round(x,3))
comp_df['statsmodel_params'] = comp_df['statsmodel_params'].apply(lambda x:np.round(x,3))
#comp_df['Simulated_params']= [b,a]
comp_df<jupyter_output><empty_output><jupyter_text>**With Standard scaling**<jupyter_code>from sklearn.preprocessing import StandardScaler
sc= StandardScaler()
data_trans = data.copy()
trans_mat = sc.fit_transform(data_trans)
data_trans['features'] = trans_mat[:,0]
data_trans['targets']= trans_mat[:,1]
data_trans.head(5)<jupyter_output><empty_output><jupyter_text>Training statsmodel<jupyter_code>#With standard scaling of data
feats, targets= data_trans['features'].values, data_trans['targets'].values
feats= sm.add_constant(feats, prepend=False)

model_stm = sm.OLS(targets, feats, hasconst=True)
result_sc= model_stm.fit()
print(result_sc.summary())
#With standard scaling of data
print('Reference values from simulated dataset weight : {} & bias: {}'.format(b,a))<jupyter_output>Reference values from simulated dataset weight : 0.10458035217582123 & bias: 0.11878560346515568<jupyter_text>Training pytorch model<jupyter_code>model_sc = toynet()

optimizer_sc = optim.Adam(model_sc.parameters())
loss = nn.MSELoss()
wts_bias_progression= list()
all_grads= list()
losses= list()

d2loss_dg= list()
diags= list()
f= lambda x: np.sqrt(np.abs(np.diag(np.linalg.inv(x)))).tolist()

for epoch in range(4000):
    x_, y_ = sample(data_trans, count= 64)
    #x_,y_= data.values[:,0].reshape(-1,1), data.values[:,1]
    trainx = Variable(torch.Tensor(x_).float(), requires_grad= True)
    trainy = Variable(torch.reshape(torch.Tensor(y_).float(), (-1,1)))
    optimizer_sc.zero_grad()
    y_pred= model_sc(trainx)
    
    l = loss(y_pred, trainy)
    losses.append(l.item())
    if epoch%500==0:
        print('epoch: ', epoch, ' loss: ', l.item());
    
    l.backward()
    optimizer_sc.step()
    
    wts_bias = list(model_sc.parameters())
    wts_bias_vals= torch.cat((wts_bias[0], torch.reshape(wts_bias[1], (-1,1))),1)
    wts_bias_progression.append(wts_bias_vals.tolist())
    
    wts_bias_grads = torch.cat((wts_bias[0].grad, torch.reshape(wts_bias[1].grad, (-1,1))),1)
    #i_grads= trainx.grad
    all_grads.append(wts_bias_grads.tolist())#i_grads)
    d1loss_dg = torch.autograd.grad(loss(model_sc(trainx), trainy), wts_bias, create_graph=True)
    d2loss_dg.append(eval_hessian(d1loss_dg, model_sc))
    diags.append(f(d2loss_dg[-1]))

all_wts_bias = np.array(wts_bias_progression).reshape(-1, wts_bias_vals.shape[1])
all_wts_bias_df= pd.DataFrame(all_wts_bias, columns=['weight_{}'.format(idx) for idx in range(1,wts_bias_vals.shape[1])]+ ['bias'])

all_grads = np.array(all_grads).reshape(-1,wts_bias_vals.shape[1])
all_grads_df= pd.DataFrame(all_grads, columns=['beta_{}'.format(idx) for idx in range(1,wts_bias_vals.shape[1]+1)])
list(model_sc.parameters())<jupyter_output><empty_output><jupyter_text>Analytical standard error & T-value for b (slope estimator)<jupyter_code>#S_b calc
sum_sqr_err = (1/(data_trans['features'].shape[0] - 2))*(np.sum(error**2))
sum_sqr_diffx = np.sum((data_trans['features'].values - np.mean(data_trans['features'].values))**2)
s_beta = np.sqrt(sum_sqr_err/sum_sqr_diffx)
s_beta
#(1/(1000-2))*(np.sum(error**2))/np.sum((data['features'].values - np.mean(data['features'].values))**2)
b_pred = list(model_sc.parameters())[0] 
t= (b_pred - b)/s_beta
t#t-val<jupyter_output><empty_output><jupyter_text>Calculating Pytorch Standar error<jupyter_code>mini_batches= list()#contains mini_batches of gradients from 200 epochs
mean_mini_batches= list()#contains mean of each of 200 epochs each in a mini_batches of gradients from 200 epochs
mini_size= 200
for batch in range(1, (len(all_grads))//mini_size):
    mini_batches.append(all_grads[batch*mini_size:batch*mini_size+ mini_size])
    mean_mini_batches.append(np.mean(mini_batches[-1], axis=0).tolist())

grad_mini_batch_arr = np.array(mean_mini_batches).T
grad_mini_batch_arr.shape
covariant_mat= np.cov(grad_mini_batch_arr)
pytorch_stderr_sc= np.sqrt(covariant_mat.diagonal())#**2
pytorch_stderr_sc<jupyter_output><empty_output><jupyter_text>Stasmodel vs Pytorch comparison<jupyter_code>summary= result_sc.summary()
statsm_df = pd.DataFrame(summary.tables[1].data[1:], columns=summary.tables[1].data[0])
statsm_df
di = dict(zip('pytorch_grads_sc,simulated_params_sc,pytorch_params_sc,statsmodel_params_sc,statsmodel_stderr_sc'.split(','), [all_grads[-1,:], [b,a], all_wts_bias[-1,:], result_sc.params, statsm_df["std err"].values]))#all_grads
di.update({'pytorch_stderr_sc':pytorch_stderr_sc})
comp_df_sc = pd.DataFrame(di, index=['b (slope/weight)','a (intercept/bias)'])
#comp_df['pytorch_stderr']= pytorch_stderr
comp_df_sc['pytorch_params_sc'] = comp_df_sc['pytorch_params_sc'].apply(lambda x:np.round(x,3))
comp_df_sc['statsmodel_params_sc'] = comp_df_sc['statsmodel_params_sc'].apply(lambda x:np.round(x,3))
#comp_df['Simulated_params']= [b,a]
comp_df_sc<jupyter_output><empty_output><jupyter_text>Scaled vs. Unscaled parameters<jupyter_code>pd.concat([comp_df, comp_df_sc], axis=1)<jupyter_output><empty_output><jupyter_text>Experiment with smaller dataset<jupyter_code>Xs = np.round(np.random.uniform(0,5, size= 20), 3)

a_s, b_s  = np.random.normal(0, 0.1), np.random.normal(0, 0.2)
error_s = np.random.normal(0, 1, size= 20)
f= lambda x,e : np.round(a_s +b_s*x + e, 3)
data_s = pd.DataFrame({'features': Xs, 'targets':f(Xs, error_s)})
sc= StandardScaler()
trans_mat_s = sc.fit_transform(data_s)

data_s['features'] = trans_mat_s[:,0]
data_s['targets']= trans_mat_s[:,1]

data_s<jupyter_output><empty_output><jupyter_text>Statsmodel<jupyter_code>import statsmodels.api as sm
feats_s, targets_s= data_s['features'].values, data_s['targets'].values
feats_s= sm.add_constant(feats_s, prepend=False)
#Without standard_scaling of data
model_stm_s = sm.OLS(targets_s, feats_s, hasconst=True)
result_s= model_stm_s.fit()
print(result_s.summary())
#Without standard_scaling of data
print('From dataset weight : {} & bias: {}'.format(b_s,a_s))<jupyter_output>From dataset weight : 0.10144210862371957 & bias: -0.2742068459149354<jupyter_text>Pytorch model<jupyter_code>model_s = toynet()
optimizer_s = optim.Adam(model_s.parameters())
loss_s = nn.MSELoss()
wts_bias_progression_s= list()
all_grads_s= list()
losses_s= list()

d2loss_dg_s= list()
diags_s= list()
f= lambda x: np.sqrt(np.abs(np.diag(np.linalg.inv(x)))).tolist()

for epoch in range(4000):
    x_, y_ = sample(data_s, count= 8)
    #x_,y_= data.values[:,0].reshape(-1,1), data.values[:,1]
    trainx = Variable(torch.Tensor(x_).float(), requires_grad= True)
    trainy = Variable(torch.reshape(torch.Tensor(y_).float(), (-1,1)))
    optimizer_s.zero_grad()
    y_pred= model_s(trainx)
    
    l = loss_s(y_pred, trainy)
    losses_s.append(l.item())
    if epoch%500==0:
        print('epoch: ', epoch, ' loss: ', l.item());
    
    l.backward()
    optimizer_s.step()
    
    wts_bias = list(model_s.parameters())
    wts_bias_vals= torch.cat((wts_bias[0], torch.reshape(wts_bias[1], (-1,1))),1)
    wts_bias_progression_s.append(wts_bias_vals.tolist())
    
    wts_bias_grads = torch.cat((wts_bias[0].grad, torch.reshape(wts_bias[1].grad, (-1,1))),1)
    #i_grads= trainx.grad
    all_grads_s.append(wts_bias_grads.tolist())#i_grads)
    d1loss_dg_s = torch.autograd.grad(loss_s(model_s(trainx), trainy), wts_bias, create_graph=True)
    d2loss_dg_s.append(eval_hessian(d1loss_dg_s, model_s))
    diags_s.append(f(d2loss_dg_s[-1]))

all_wts_bias_s = np.array(wts_bias_progression_s).reshape(-1,wts_bias_vals.shape[1])
all_wts_bias_df_s= pd.DataFrame(all_wts_bias_s, columns=['weight_{}'.format(idx) for idx in range(1,wts_bias_vals.shape[1])]+ ['bias'])

all_grads_s = np.array(all_grads_s).reshape(-1,wts_bias_vals.shape[1])
all_grads_df_s= pd.DataFrame(all_grads_s, columns=['beta_{}'.format(idx) for idx in range(1,wts_bias_vals.shape[1]+1)])<jupyter_output>epoch:  0  loss:  0.7195181250572205
epoch:  500  loss:  0.8909294009208679
epoch:  1000  loss:  0.9296377897262573
epoch:  1500  loss:  0.38164255023002625
epoch:  2000  loss:  0.6488150358200073
epoch:  2500  loss:  0.4885692000389099
epoch:  3000  loss:  0.9722121953964233
epoch:  3500  loss:  1.1143834590911865<jupyter_text>Calculating Pytorch Standar error<jupyter_code>mini_batches_s= list()#contains mini_batches of gradients from 200 epochs
mean_mini_batches_s= list()#contains mean of each of 200 epochs each in a mini_batches of gradients from 200 epochs
mini_size= 200
for batch in range(1, (len(all_grads_s))//mini_size):
    mini_batches_s.append(all_grads_s[batch*mini_size:batch*mini_size+ mini_size])
    mean_mini_batches_s.append(np.mean(mini_batches_s[-1], axis=0).tolist())

grad_mini_batch_arr_s = np.array(mean_mini_batches_s).T
grad_mini_batch_arr_s.shape
covariant_mat_s= np.cov(grad_mini_batch_arr_s)
pytorch_stderr_s= np.sqrt(covariant_mat_s.diagonal())#**2
pytorch_stderr_s
summary_s= result_s.summary()
statsm_df_s = pd.DataFrame(summary_s.tables[1].data[1:], columns=summary_s.tables[1].data[0])
statsm_df_s
#Without standard_scaling
di_s = dict(zip('pytorch_grads,simulated_params,pytorch_params,statsmodel_params,statsmodel_stderr'.split(','), [all_grads_s[-1,:], [b_s,a_s], all_wts_bias_s[-1,:], result_s.params, statsm_df_s["std err"].values]))#all_grads
di_s.update({'pytorch_stderr':pytorch_stderr_s})
comp_df_s = pd.DataFrame(di_s, index=['b (slope/weight)','a (intercept/bias)'])
comp_df_s['pytorch_params'] = comp_df_s['pytorch_params'].apply(lambda x:np.round(x,3))
comp_df_s['statsmodel_params'] = comp_df_s['statsmodel_params'].apply(lambda x:np.round(x,3))
#comp_df['Simulated_params']= [b,a]
comp_df_s<jupyter_output><empty_output><jupyter_text>Standard error & T-value for b (slope estimator)<jupyter_code>#S_b calc
sum_sqr_err_s = (1/(data_s['features'].shape[0] - 2))*(np.sum(error_s**2))
sum_sqr_diffx_s = np.sum((data_s['features'].values - np.mean(data_s['features'].values))**2)
s_beta_s = np.sqrt(sum_sqr_err_s/sum_sqr_diffx_s)
s_beta_s
#(1/(1000-2))*(np.sum(error**2))/np.sum((data['features'].values - np.mean(data['features'].values))**2)
# t-value = (b'- b)/(S_b)
b_pred_s = list(model_s.parameters())[0] 
t_s= (b_pred_s - b_s)/s_beta_s
t_s<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16186419443748434, 'fraction_numerical': 0.03614382360310699, 'mean_word_length': 6.520018841262365, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 69, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '289201', 'license': 'MIT', 'n_tokens_mistral': 6805, 'n_tokens_neox': 6500, 'n_words': 955, 'path': 'notebooks/UQ/02_prakash_pytorch_confidence_intervals.ipynb', 'repo_name': 'mlsquare/kitchen', 'size': '15964'}