<jupyter_start><jupyter_text>Defining pytorch model * **Following pytorch model has one linear layer with `inputs`: 3 and `outputs`:1*** **Kernel are initialized with He initialization**<jupyter_code>import pandas as pd
import numpy as np
import torch
import torch.nn as nn
from torch.autograd import Variable
import torch.nn.functional as F
import torch.optim as optim
import matplotlib.pyplot as plt

def create_model():
    class network(nn.Module):
        def __init__(self):
            super(network, self).__init__()
            self.fc1 = nn.Linear(3, 1)#100)
            
        def forward(self, X):
            X = self.fc1(X)#F.linear(self.fc1(X))#, list(self.fc1.parameters())[0])
            return X
    model= network()
    optimizer = optim.Adam(model.parameters())
    criterion= nn.MSELoss()# CrossEntropyLoss()
    return model, optimizer, criterion
model, optimizer, loss= create_model()
model<jupyter_output><empty_output><jupyter_text>Data preparation and scaling        - Each record `X` cosist of first 3 values of 4 values          from iris dataset input features.                - 4th value serve as taget label `y` .<jupyter_code>def scale(arr):
    f= lambda x: (x-np.mean(x))/np.var(x)
    for idx in range(arr.shape[1]):
        arr[:,idx]= f(arr[:, idx])
    #return arr
def plot(v1):
    plt.plot(v1)
    #plt.plot(v2)
    plt.title('Model training loss')
    plt.xlabel('epoch')
    plt.ylabel('loss')
    #plt.legend(['train', 'validation'], loc='upper right')
    return plt.show()

def samp(x, count):#For random batch generation
    idx= np.random.randint(x.shape[0], size=count)
    res = x[idx]
    return res<jupyter_output><empty_output><jupyter_text>* **Hessian calculation**<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>Model Training & plots<jupyter_code>from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
data = load_iris()
X, Y = data.data, data.target

Y[np.where(Y==2)[0]]=1
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size= 0.3)
scale(X)
scale(X_test)

#trainx= Variable(torch.Tensor(X[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
#trainy = Variable(torch.reshape(torch.Tensor(X[:,-1]).float(), (-1,1)))

#i_grads= np.empty(trainx.shape)
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_ =samp(X, count= 32)
    trainx= Variable(torch.Tensor(x_[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
    trainy = Variable(torch.reshape(torch.Tensor(x_[:,-1]).float(), (-1,1)))

    optimizer.zero_grad()
    y_pred = model(trainx)
    #print(y_pred)
    l = loss(y_pred, trainy)
    losses.append(l.item())
    if epoch%500==0:
        print('epoch: ', epoch, ' loss: ', l.item());
    # perform a backward pass (backpropagation)
    l.backward()
    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)
    # Update the parameters
    optimizer.step()

    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,4)
all_wts_bias_df= pd.DataFrame(all_wts_bias, columns=['weight_{}'.format(idx) for idx in range(1,4)]+ ['bias'])

all_grads = np.array(all_grads).reshape(-1,4)
all_grads_df= pd.DataFrame(all_grads, columns=['beta_{}'.format(idx) for idx in range(1,5)])
plot(losses)<jupyter_output>epoch:  0  loss:  2.255631685256958
epoch:  500  loss:  0.48510271310806274
epoch:  1000  loss:  0.3253370523452759
epoch:  1500  loss:  0.294547438621521
epoch:  2000  loss:  0.3795126676559448
epoch:  2500  loss:  0.13330158591270447
epoch:  3000  loss:  0.14958198368549347
epoch:  3500  loss:  0.125545933842659<jupyter_text>* **Following are the weights and bias for the ONLY layer in model, after 128 epochs**<jupyter_code>li = list(model.parameters())
li
print('Final weights gradients : {}\nFinal bias gradients: {}'.format(li[0].grad, li[1].grad))<jupyter_output>Final weights gradients : tensor([[-0.1232, -0.0649, -0.1524]])
Final bias gradients: tensor([0.1111])<jupyter_text>Gradients across epochs* **Following are gradients(for weight- Beta_1-3 and bias-beta_4) updated after each epochs**<jupyter_code>all_grads_df
all_grads_df.describe()<jupyter_output><empty_output><jupyter_text>Defining model with statsmodel* **Evaluating same data with linear regression from statsmodel**<jupyter_code>import statsmodels.api as sm
feats, targets= X[:,:-1], X[:,-1]
feats= sm.add_constant(feats, prepend=False)
model_stm = sm.OLS(targets, feats, hasconst=True)
result= model_stm.fit()
print(result.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.938
Model:                            OLS   Adj. R-squared:                  0.937
Method:                 Least Squares   F-statistic:                     734.4
Date:                Thu, 19 Mar 2020   Prob (F-statistic):           7.83e-88
Time:                        20:26:33   Log-Likelihood:                -45.701
No. Observations:                 150   AIC:                             99.40
Df Residuals:                     146   BIC:                             111.4
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>* **Weights, biases and std_error from modelstats regression model**<jupyter_code>summary= result.summary()
result.params
statsm_df = pd.DataFrame(summary.tables[1].data[1:], columns=summary.tables[1].data[0])
statsm_df<jupyter_output><empty_output><jupyter_text>* **pytorch_params vs. statsmodel_params**<jupyter_code>di = dict(zip('pytorch_grads,pytorch_params,statsmodel_params,statsmodel_stderr'.split(','), [all_grads[-1,:], all_wts_bias[-1,:], result.params, statsm_df["std err"].values]))#all_grads
di<jupyter_output><empty_output><jupyter_text>Paramters comparison<jupyter_code>comp_df = pd.DataFrame(di)
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<jupyter_output><empty_output><jupyter_text>Plots<jupyter_code>def scatter_plot(x, y, x_label, y_label):
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    min_x = min(list(x)+list(y)) - 0.5
    max_x = max(list(x)+list(y)) + 0.5
    ax.plot([min_x, max_x], [min_x, max_x])
    
    ax.set_xlim([min_x, max_x])
    ax.set_ylim([min_x, max_x])
    
    ax.set_title('{} vs. {} correlation scatterplot'.format(x_label, y_label))
    ax.scatter(x, y, c="red")
    plt.ylabel(y_label)
    plt.xlabel(x_label)
    plt.show()
    
def hexbin_plot(x, y, x_label, y_label):
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    min_x = min(list(x)+list(y)) - 0.1
    max_x = max(list(x)+list(y)) + 0.1
    ax.plot([min_x, max_x], [min_x, max_x])
    
    ax.set_xlim([min_x, max_x])
    ax.set_ylim([min_x, max_x])
    
    ax.set_title('{} vs. {} correlation scatterplot'.format(x_label, y_label))
    hbin= ax.hexbin(x, y, gridsize=25, mincnt=2)
    cb = fig.colorbar(hbin, ax=ax)
    cb.solids.set_edgecolor('face')#change colors
    cb.set_label('occurence_density')
    plt.ylabel(y_label)
    plt.xlabel(x_label)
    plt.show()
#Plot with trained cases. When predictions are calculated of train_data
predictions= model(torch.Tensor(X[:,:-1]).float())
predictions= torch.Tensor.cpu(predictions).detach().numpy().reshape(-1,)
observed= X[:,-1]


hexbin_plot(predictions,observed, 'prediction' ,'observed')
scatter_plot(comp_df['pytorch_params'].values,comp_df['statsmodel_params'].values, 'pytorch params' ,'statsmodel params')#n1,n2, 'abc','bcd')<jupyter_output><empty_output><jupyter_text>Calculating Standard error from pytorch model<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(list(np.mean(mini_batches[-1], axis=0)))

grad_mini_batch_arr = np.array(mean_mini_batches).T
grad_mini_batch_arr.shape
covariant_mat= np.cov(grad_mini_batch_arr)
covariant_mat
pytorch_stderr= np.sqrt(covariant_mat.diagonal())#**2
pytorch_stderr
di.update({'pytorch_stderr':pytorch_stderr})
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<jupyter_output><empty_output><jupyter_text>Experiments with configurations  Model with xavier normal kernel initialization<jupyter_code>import pandas as pd
import numpy as np
import torch
import torch.nn as nn
from torch.autograd import Variable
import torch.nn.functional as F
import torch.optim as optim
import matplotlib.pyplot as plt

def create_model_2():
    class network(nn.Module):
        def __init__(self):
            super(network, self).__init__()
            self.fc1 = nn.Linear(3, 1)#100)
            torch.nn.init.xavier_normal_(self.fc1.weight)

        def forward(self, X):
            X = self.fc1(X)#F.linear(self.fc1(X))#, list(self.fc1.parameters())[0])
            return X
    model= network()
    optimizer = optim.Adam(model.parameters())
    criterion= nn.MSELoss()# CrossEntropyLoss()
    return model, optimizer, criterion

model2, optimizer2, loss2= create_model_2()
model2<jupyter_output><empty_output><jupyter_text>Model training and comparison plots<jupyter_code>#trainx= Variable(torch.Tensor(X[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
#trainy = Variable(torch.reshape(torch.Tensor(X[:,-1]).float(), (-1,1)))

#i_grads= np.empty(trainx.shape)
wts_bias_progression_2= list()
all_grads_2= list()
losses_2= list()

d2loss_dg2= list()
diags_2= list()
f= lambda x: np.sqrt(np.abs(np.diag(np.linalg.inv(x)))).tolist()
    
for epoch in range(4000):
    x_transformed =samp(X, count= 32)
    trainx= Variable(torch.Tensor(x_transformed[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
    trainy = Variable(torch.reshape(torch.Tensor(x_transformed[:,-1]).float(), (-1,1)))

    optimizer2.zero_grad()
    y_pred = model2(trainx)
    #print(y_pred)
    l = loss2(y_pred, trainy)
    losses_2.append(l.item())
    if epoch%500==0:
        print('epoch: ', epoch, ' loss: ', l.item());
    # perform a backward pass (backpropagation)
    l.backward()
    optimizer2.step()

    wts_bias = list(model2.parameters())
    wts_bias_vals= torch.cat((wts_bias[0], torch.reshape(wts_bias[1], (-1,1))),1)
    wts_bias_progression_2.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_2.append(wts_bias_grads.tolist())#i_grads)
    # Update the parameters    

    d1loss_dg2 = torch.autograd.grad(loss2(model2(trainx), trainy), wts_bias, create_graph=True)
    d2loss_dg2.append(eval_hessian(d1loss_dg2, model2))
    diags_2.append(f(d2loss_dg2[-1]))

    
all_wts_bias_2 = np.array(wts_bias_progression_2).reshape(-1,4)
all_wts_bias_df_2= pd.DataFrame(all_wts_bias_2, columns=['weight_{}'.format(idx) for idx in range(1,4)]+ ['bias'])

all_grads_2 = np.array(all_grads_2).reshape(-1,4)
all_grads_df_2= pd.DataFrame(all_grads_2, columns=['beta_{}'.format(idx) for idx in range(1,5)])
plot(losses_2)
#Plot with trained cases. When predictions are calculated of train_data
predictions= model2(torch.Tensor(X[:,:-1]).float())
predictions= torch.Tensor.cpu(predictions).detach().numpy().reshape(-1,)
observed= X[:,-1]


hexbin_plot(predictions,observed, 'prediction' ,'observed')<jupyter_output><empty_output><jupyter_text>Pytorch model Standard Error & Parameters Comparison * **Calculating standard error for pytorch model**<jupyter_code>mini_batches_2= list()#contains mini_batches of gradients from 200 epochs
mean_mini_batches_2= list()#contains mean of each of 200 epochs each in a mini_batches of gradients from 200 epochs
mini_size= 200#mini_batch size
for batch in range(1, (len(all_grads_2))//mini_size):
    mini_batches_2.append(all_grads_2[batch*mini_size:batch*mini_size+ mini_size])
    mean_mini_batches_2.append(list(np.mean(mini_batches_2[-1], axis=0)))

grad_mini_batch_arr_2 = np.array(mean_mini_batches_2).T
grad_mini_batch_arr_2.shape
covariant_mat_2= np.cov(grad_mini_batch_arr_2)
pytorch_stderr_2= np.sqrt(covariant_mat_2.diagonal())
pytorch_stderr_2
di_2 = dict(zip('pytorch_grads,pytorch_params,statsmodel_params,statsmodel_stderr'.split(','), [all_grads_2[-1,:], all_wts_bias_2[-1,:], result.params, statsm_df["std err"].values]))
di_2.update({'pytorch_stderr':pytorch_stderr_2})
comp_df_2 = pd.DataFrame(di_2)
comp_df_2['pytorch_params'] = comp_df_2['pytorch_params'].apply(lambda x:np.round(x,3))
comp_df_2['statsmodel_params'] = comp_df_2['statsmodel_params'].apply(lambda x:np.round(x,3))
comp_df_2<jupyter_output><empty_output><jupyter_text>Model with xavier_uniform kernel initialization<jupyter_code>def create_model_3():
    class network(nn.Module):
        def __init__(self):
            super(network, self).__init__()
            self.fc1 = nn.Linear(3, 1)#100)
            torch.nn.init.xavier_uniform_(self.fc1.weight)

        def forward(self, X):
            X = self.fc1(X)#F.linear(self.fc1(X))#, list(self.fc1.parameters())[0])
            return X
    model= network()
    optimizer = optim.Adam(model.parameters())
    criterion= nn.MSELoss()# CrossEntropyLoss()
    return model, optimizer, criterion

model3, optimizer3, loss3= create_model_3()
model3<jupyter_output><empty_output><jupyter_text>Model training and comparison plots<jupyter_code>#trainx= Variable(torch.Tensor(X[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
#trainy = Variable(torch.reshape(torch.Tensor(X[:,-1]).float(), (-1,1)))

#i_grads= np.empty(trainx.shape)
wts_bias_progression_3= list()
all_grads_3= list()
losses_3= list()

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

for epoch in range(4000):
    x_transformed =samp(X, count= 32)
    trainx= Variable(torch.Tensor(x_transformed[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
    trainy = Variable(torch.reshape(torch.Tensor(x_transformed[:,-1]).float(), (-1,1)))

    optimizer3.zero_grad()
    y_pred = model3(trainx)
    #print(y_pred)
    l = loss3(y_pred, trainy)
    losses_3.append(l.item())
    if epoch%500==0:
        print('epoch: ', epoch, ' loss: ', l.item());
    # perform a backward pass (backpropagation)
    l.backward()
    optimizer3.step()
    
    wts_bias = list(model3.parameters())
    wts_bias_vals= torch.cat((wts_bias[0], torch.reshape(wts_bias[1], (-1,1))),1)
    wts_bias_progression_3.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_3.append(wts_bias_grads.tolist())#i_grads)
    # Update the parameters
    
    d1loss_dg3 = torch.autograd.grad(loss3(model3(trainx), trainy), wts_bias, create_graph=True)
    d2loss_dg3.append(eval_hessian(d1loss_dg3, model3))
    diags_3.append(f(d2loss_dg3[-1]))

    
all_wts_bias_3 = np.array(wts_bias_progression_3).reshape(-1,4)
all_wts_bias_df_3= pd.DataFrame(all_wts_bias_3, columns=['weight_{}'.format(idx) for idx in range(1,4)]+ ['bias'])

all_grads_3 = np.array(all_grads_3).reshape(-1,4)
all_grads_df_3= pd.DataFrame(all_grads_3, columns=['beta_{}'.format(idx) for idx in range(1,5)])
plot(losses_3)
#Plot with trained cases. When predictions are calculated of train_data
predictions= model3(torch.Tensor(X[:,:-1]).float())
predictions= torch.Tensor.cpu(predictions).detach().numpy().reshape(-1,)
observed= X[:,-1]


hexbin_plot(predictions,observed, 'prediction' ,'observed')<jupyter_output><empty_output><jupyter_text>Pytorch model Standard Error & Parameters Comparison* **Calculating standard error for pytorch model**<jupyter_code>mini_batches_3= list()#contains mini_batches of gradients from 200 epochs
mean_mini_batches_3= list()#contains mean of each of 200 epochs each in a mini_batches of gradients from 200 epochs
mini_size= 200#mini_batch size
for batch in range(1, (len(all_grads_3))//mini_size):
    mini_batches_3.append(all_grads_3[batch*mini_size:batch*mini_size+ mini_size])
    mean_mini_batches_3.append(list(np.mean(mini_batches_3[-1], axis=0)))

grad_mini_batch_arr_3 = np.array(mean_mini_batches_3).T
grad_mini_batch_arr_3.shape
covariant_mat_3= np.cov(grad_mini_batch_arr_3)
pytorch_stderr_3= np.sqrt(covariant_mat_3.diagonal())
pytorch_stderr_3
di_3 = dict(zip('pytorch_grads,pytorch_params,statsmodel_params,statsmodel_stderr'.split(','), [all_grads_3[-1,:], all_wts_bias_3[-1,:], result.params, statsm_df["std err"].values]))
di_3.update({'pytorch_stderr':pytorch_stderr_3})
comp_df_3 = pd.DataFrame(di_3)
comp_df_3['pytorch_params'] = comp_df_3['pytorch_params'].apply(lambda x:np.round(x,3))
comp_df_3['statsmodel_params'] = comp_df_3['statsmodel_params'].apply(lambda x:np.round(x,3))
comp_df_3<jupyter_output><empty_output><jupyter_text>Model training with different data scaling (MinMax)<jupyter_code>from sklearn.preprocessing import MinMaxScaler

X, Y = data.data, data.target
sc= MinMaxScaler()
X_trans = sc.fit_transform(X)

Y[np.where(Y==2)[0]]=1
X_train, X_test, y_train, y_test = train_test_split(X_trans, Y, test_size= 0.3)<jupyter_output><empty_output><jupyter_text>Model training and comparison plots<jupyter_code>model4, optimizer4, loss4= create_model()#With default He Weights Initilizing model
model4
#trainx= Variable(torch.Tensor(X_trans[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
#trainy = Variable(torch.reshape(torch.Tensor(X_trans[:,-1]).float(), (-1,1)))

wts_bias_progression_4= list()
all_grads_4= list()
losses_4= list()
d2loss_dg4= list()
diags_4= list()
f= lambda x: np.sqrt(np.abs(np.diag(np.linalg.inv(x)))).tolist()

for epoch in range(4000):
    x_transformed =samp(X_trans, count= 32)
    trainx= Variable(torch.Tensor(x_transformed[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
    trainy = Variable(torch.reshape(torch.Tensor(x_transformed[:,-1]).float(), (-1,1)))

    optimizer4.zero_grad()
    y_pred = model4(trainx)
    #print(y_pred)
    l = loss4(y_pred, trainy)
    losses_4.append(l.item())
    if epoch%500==0:
        print('epoch: ', epoch, ' loss: ', l.item());
    # perform a backward pass (backpropagation)
    l.backward()
    optimizer4.step()
    
    wts_bias = list(model4.parameters())
    wts_bias_vals= torch.cat((wts_bias[0], torch.reshape(wts_bias[1], (-1,1))),1)
    wts_bias_progression_4.append(wts_bias_vals.tolist())
    
    wts_bias_grads = torch.cat((wts_bias[0].grad, torch.reshape(wts_bias[1].grad, (-1,1))),1)
    all_grads_4.append(wts_bias_grads.tolist())#i_grads)

    d1loss_dg4 = torch.autograd.grad(loss4(model4(trainx), trainy), wts_bias, create_graph=True)
    d2loss_dg4.append(eval_hessian(d1loss_dg4, model4))
    diags_4.append(f(d2loss_dg4[-1]))

all_wts_bias_4 = np.array(wts_bias_progression_4).reshape(-1,4)
all_wts_bias_df_4= pd.DataFrame(all_wts_bias_4, columns=['weight_{}'.format(idx) for idx in range(1,4)]+ ['bias'])

all_grads_4 = np.array(all_grads_4).reshape(-1,4)
all_grads_df_4= pd.DataFrame(all_grads_4, columns=['beta_{}'.format(idx) for idx in range(1,5)])
plot(losses_4)
#Plot with trained cases. When predictions are calculated of train_data
predictions= model4(torch.Tensor(X_trans[:,:-1]).float())
predictions= torch.Tensor.cpu(predictions).detach().numpy().reshape(-1,)
observed= X_trans[:,-1]


hexbin_plot(predictions,observed, 'prediction' ,'observed')<jupyter_output><empty_output><jupyter_text>Retraining Statsmodel model with Scaled data* **Evaluating same data with linear regression from statsmodel**Statsmodel model Standard Error & Parameters Comparison<jupyter_code>feats, targets= X_trans[:,:-1], X_trans[:,-1]
feats= sm.add_constant(feats, prepend=False)

model_stm_2 = sm.OLS(targets, feats, hasconst=True)
result_2= model_stm_2.fit()
print(result_2.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.938
Model:                            OLS   Adj. R-squared:                  0.937
Method:                 Least Squares   F-statistic:                     734.4
Date:                Thu, 19 Mar 2020   Prob (F-statistic):           7.83e-88
Time:                        20:33:09   Log-Likelihood:                 168.07
No. Observations:                 150   AIC:                            -328.1
Df Residuals:                     146   BIC:                            -316.1
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>* **Weights, biases and std_error from modelstats regression model**<jupyter_code>summary_2= result_2.summary()
result_2.params
statsm_df_2 = pd.DataFrame(summary_2.tables[1].data[1:], columns=summary_2.tables[1].data[0])
statsm_df_2<jupyter_output><empty_output><jupyter_text>Pytorch model Standard Error & Parameters Comparison* **Calculating standard error for pytorch model**<jupyter_code>mini_batches_4= list()#contains mini_batches of gradients from 200 epochs
mean_mini_batches_4= list()#contains mean of each of 200 epochs each in a mini_batches of gradients from 200 epochs
mini_size= 200#mini_batch size
for batch in range(1, (len(all_grads_4))//mini_size):
    mini_batches_4.append(all_grads_4[batch*mini_size:batch*mini_size+ mini_size])
    mean_mini_batches_4.append(list(np.mean(mini_batches_4[-1], axis=0)))

grad_mini_batch_arr_4 = np.array(mean_mini_batches_4).T
grad_mini_batch_arr_4.shape
covariant_mat_4= np.cov(grad_mini_batch_arr_4)
pytorch_stderr_4= np.sqrt(covariant_mat_4.diagonal())
pytorch_stderr_4
di_4 = dict(zip('pytorch_grads,pytorch_params,statsmodel_params,statsmodel_stderr'.split(','), [all_grads_4[-1,:], all_wts_bias_4[-1,:], result_2.params, statsm_df_2["std err"].values]))
di_4.update({'pytorch_stderr':pytorch_stderr_4})
comp_df_4 = pd.DataFrame(di_4)
comp_df_4['pytorch_params'] = comp_df_4['pytorch_params'].apply(lambda x:np.round(x,3))
comp_df_4['statsmodel_params'] = comp_df_4['statsmodel_params'].apply(lambda x:np.round(x,3))
comp_df_4<jupyter_output><empty_output><jupyter_text>Model training with reference Weights from statsmodel<jupyter_code>ref_kernel = torch.reshape(torch.Tensor(result_2.params[:-1]+ np.random.normal(0,0.1, result_2.params[:-1].shape[0])).float(), (1,-1))
ref_bias = result_2.params[-1] + 0.1#Adding noise
def create_model_5():
    class network(nn.Module):
        def __init__(self):
            super(network, self).__init__()
            self.fc1 = nn.Linear(3, 1)#100)
            self.fc1.weight= torch.nn.Parameter(ref_kernel)
            #torch.nn.init.constant_(ref_kernel,0.3)
            #self.fc1.weight.data= ref_kernel
            self.fc1.bias.data.fill_(ref_bias)
            #torch.nn.init.xavier_uniform_(self.fc1.weight)

        def forward(self, X):
            X = self.fc1(X)#F.linear(self.fc1(X))#, list(self.fc1.parameters())[0])
            return X
    model= network()
    optimizer = optim.Adam(model.parameters())#LBFGS(model.parameters())
    criterion= nn.MSELoss()# CrossEntropyLoss()
    return model, optimizer, criterion

model5, optimizer5, loss5= create_model_5()
model5
result_2.params#refrenced weights, bias from statsmodel
list(model5.parameters())
#trainx= Variable(torch.Tensor(X_trans[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
#trainy = Variable(torch.reshape(torch.Tensor(X_trans[:,-1]).float(), (-1,1)))

wts_bias_progression_5= list()
all_grads_5= list()
losses_5= list()
d2loss_dg5= list()
diags_5= list()
f= lambda x: np.sqrt(np.abs(np.diag(np.linalg.inv(x)))).tolist()


for epoch in range(4000):
    x_transformed =samp(X_trans, count= 32)
    trainx= Variable(torch.Tensor(x_transformed[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
    trainy = Variable(torch.reshape(torch.Tensor(x_transformed[:,-1]).float(), (-1,1)))

    optimizer5.zero_grad()
    y_pred = model5(trainx)
    #print(y_pred)
    l = loss5(y_pred, trainy)
    losses_5.append(l.item())
    if epoch%500==0:
        print('epoch: ', epoch, ' loss: ', l.item());
    # perform a backward pass (backpropagation)
    l.backward()
    optimizer5.step()

    wts_bias = list(model5.parameters())
    wts_bias_vals= torch.cat((wts_bias[0], torch.reshape(wts_bias[1], (-1,1))),1)
    wts_bias_progression_5.append(wts_bias_vals.tolist())
    
    wts_bias_grads = torch.cat((wts_bias[0].grad, torch.reshape(wts_bias[1].grad, (-1,1))),1)
    all_grads_5.append(wts_bias_grads.tolist())#i_grads)
        
    d1loss_dg5 = torch.autograd.grad(loss5(model5(trainx), trainy), wts_bias, create_graph=True)
    d2loss_dg5.append(eval_hessian(d1loss_dg5, model5))
    diags_5.append(f(d2loss_dg5[-1]))

all_wts_bias_5 = np.array(wts_bias_progression_5).reshape(-1,4)
all_wts_bias_df_5= pd.DataFrame(all_wts_bias_5, columns=['weight_{}'.format(idx) for idx in range(1,4)]+ ['bias'])

all_grads_5 = np.array(all_grads_5).reshape(-1,4)
all_grads_df_5= pd.DataFrame(all_grads_5, columns=['beta_{}'.format(idx) for idx in range(1,5)])
plot(losses_5)
#Plot with trained cases. When predictions are calculated of train_data
predictions= model5(torch.Tensor(X_trans[:,:-1]).float())
predictions= torch.Tensor.cpu(predictions).detach().numpy().reshape(-1,)
observed= X_trans[:,-1]


hexbin_plot(predictions,observed, 'prediction' ,'observed')<jupyter_output><empty_output><jupyter_text>Pytorch model Standard Error & Parameters Comparison* **Calculating standard error for pytorch model**<jupyter_code>mini_batches_5= list()#contains mini_batches of gradients from 200 epochs
mean_mini_batches_5= list()#contains mean of each of 200 epochs each in a mini_batches of gradients from 200 epochs
mini_size= 200#mini_batch size
for batch in range(1, (len(all_grads_5))//mini_size):
    mini_batches_5.append(all_grads_5[batch*mini_size:batch*mini_size+ mini_size])
    mean_mini_batches_5.append(list(np.mean(mini_batches_5[-1], axis=0)))

grad_mini_batch_arr_5 = np.array(mean_mini_batches_5).T
grad_mini_batch_arr_5.shape
covariant_mat_5= np.cov(grad_mini_batch_arr_5)
pytorch_stderr_5= np.sqrt(covariant_mat_5.diagonal())#**2
pytorch_stderr_5
di_5 = dict(zip('pytorch_grads,pytorch_params,statsmodel_params,statsmodel_stderr'.split(','), [all_grads_5[-1,:], all_wts_bias_5[-1,:], result_2.params, statsm_df_2["std err"].values]))
di_5.update({'pytorch_stderr':pytorch_stderr_5})
comp_df_5 = pd.DataFrame(di_5)
comp_df_5['pytorch_params'] = comp_df_5['pytorch_params'].apply(lambda x:np.round(x,3))
comp_df_5['statsmodel_params'] = comp_df_5['statsmodel_params'].apply(lambda x:np.round(x,3))
comp_df_5<jupyter_output><empty_output><jupyter_text>Model training with LBFGS optimizer and all previously tried settings  LBFGS optimizer with Reference statsmodel weights & bias (Exp 4.4)<jupyter_code>ref_kernel = torch.reshape(torch.Tensor(result_2.params[:-1]+ np.random.normal(0,0.1, result_2.params[:-1].shape[0])).float(), (1,-1))
ref_bias = result_2.params[-1] + 0.1#Adding noise
def create_model_6():
    class network(nn.Module):
        def __init__(self):
            super(network, self).__init__()
            self.fc1 = nn.Linear(3, 1)#100)
            self.fc1.weight= torch.nn.Parameter(ref_kernel)
            #torch.nn.init.constant_(ref_kernel,0.3)
            #self.fc1.weight.data= ref_kernel
            self.fc1.bias.data.fill_(ref_bias)
            #torch.nn.init.xavier_uniform_(self.fc1.weight)

        def forward(self, X):
            X = self.fc1(X)#F.linear(self.fc1(X))#, list(self.fc1.parameters())[0])
            return X
    model= network()
    optimizer = optim.LBFGS(model.parameters())
    criterion= nn.MSELoss()# CrossEntropyLoss()
    return model, optimizer, criterion

model6, optimizer6, loss6= create_model_6()
model6
list(model6.parameters())
#trainx= Variable(torch.Tensor(X_trans[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
#trainy = Variable(torch.reshape(torch.Tensor(X_trans[:,-1]).float(), (-1,1)))
def report_loss(ep, loss_any, pred, obs):
    ls = loss_any(pred, obs)
    if ep%500==0:
        print('epoch: ', ep, ' loss: ', ls.item())


wts_bias_progression_6= list()
all_grads_6= list()
losses_6= list()
d2loss_dg6= list()
diags_6= list()
f= lambda x: np.sqrt(np.abs(np.diag(np.linalg.inv(x)))).tolist()

for epoch in range(4000):
    x_transformed =samp(X_trans, count= 32)
    trainx= Variable(torch.Tensor(x_transformed[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
    trainy = Variable(torch.reshape(torch.Tensor(x_transformed[:,-1]).float(), (-1,1)))

    def closure():
        optimizer6.zero_grad()
        y_pred = model6(trainx)

        l = loss6(y_pred, trainy)
        losses_6.append(l.item())
        #if epoch%500==0:
        #    print('epoch: ', epoch, ' loss: ', l.item());
        # perform a backward pass (backpropagation)
        l.backward()
        return l
    report_loss(epoch, loss6, model6(trainx), trainy)
    optimizer6.step(closure)

    wts_bias = list(model6.parameters())
    wts_bias_vals= torch.cat((wts_bias[0], torch.reshape(wts_bias[1], (-1,1))),1)
    wts_bias_progression_6.append(wts_bias_vals.tolist())

    #print('(torch.reshape(wts_bias[1].grad, (-1,1))),1):',epoch, type(wts_bias[1].grad))
    wts_bias_grads = torch.cat((wts_bias[0].grad, torch.reshape(wts_bias[1].grad, (-1,1))),1)
    all_grads_6.append(wts_bias_grads.tolist())#i_grads)
    
    d1loss_dg6 = torch.autograd.grad(loss6(model6(trainx), trainy), wts_bias, create_graph=True)
    d2loss_dg6.append(eval_hessian(d1loss_dg6, model6))
    diags_6.append(f(d2loss_dg6[-1]))

    

all_wts_bias_6 = np.array(wts_bias_progression_6).reshape(-1,4)
all_wts_bias_df_6= pd.DataFrame(all_wts_bias_6, columns=['weight_{}'.format(idx) for idx in range(1,4)]+ ['bias'])

all_grads_6 = np.array(all_grads_6).reshape(-1,4)
all_grads_df_6= pd.DataFrame(all_grads_6, columns=['beta_{}'.format(idx) for idx in range(1,5)])
plot(losses_6)
#Plot with trained cases. When predictions are calculated of train_data
predictions= model6(torch.Tensor(X_trans[:,:-1]).float())
predictions= torch.Tensor.cpu(predictions).detach().numpy().reshape(-1,)
observed= X_trans[:,-1]

hexbin_plot(predictions,observed, 'prediction' ,'observed')<jupyter_output><empty_output><jupyter_text>Pytorch model Standard Error & Parameters Comparison* **Calculating standard error for pytorch model**<jupyter_code>mini_batches_6= list()#contains mini_batches of gradients from 200 epochs
mean_mini_batches_6= list()#contains mean of each of 200 epochs each in a mini_batches of gradients from 200 epochs
mini_size= 200#mini_batch size
for batch in range(1, (len(all_grads_6))//mini_size):
    mini_batches_6.append(all_grads_6[batch*mini_size:batch*mini_size+ mini_size])
    mean_mini_batches_6.append(list(np.mean(mini_batches_6[-1], axis=0)))

grad_mini_batch_arr_6 = np.array(mean_mini_batches_6).T
grad_mini_batch_arr_6.shape
covariant_mat_6= np.cov(grad_mini_batch_arr_6)
pytorch_stderr_6= np.sqrt(covariant_mat_6.diagonal())#**2
pytorch_stderr_6
di_6 = dict(zip('pytorch_grads,pytorch_params,statsmodel_params,statsmodel_stderr'.split(','), [all_grads_6[-1,:], all_wts_bias_6[-1,:], result_2.params, statsm_df_2["std err"].values]))
di_6.update({'pytorch_stderr':pytorch_stderr_6})
comp_df_6 = pd.DataFrame(di_6)
comp_df_6['pytorch_params'] = comp_df_6['pytorch_params'].apply(lambda x:np.round(x,3))
comp_df_6['statsmodel_params'] = comp_df_6['statsmodel_params'].apply(lambda x:np.round(x,3))
comp_df_6<jupyter_output><empty_output><jupyter_text>LBFGS optimizer with Xavier kernel init (Exp 4.1)<jupyter_code>def create_model_7():
    class network(nn.Module):
        def __init__(self):
            super(network, self).__init__()
            self.fc1 = nn.Linear(3, 1)
            torch.nn.init.xavier_uniform_(self.fc1.weight)

        def forward(self, X):
            X = self.fc1(X)#F.linear(self.fc1(X))#, list(self.fc1.parameters())[0])
            return X
    model= network()
    optimizer = optim.LBFGS(model.parameters())
    criterion= nn.MSELoss()# CrossEntropyLoss()
    return model, optimizer, criterion

model7, optimizer7, loss7= create_model_7()
model7
#trainx= Variable(torch.Tensor(X_trans[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
#trainy = Variable(torch.reshape(torch.Tensor(X_trans[:,-1]).float(), (-1,1)))
def report_loss(ep, loss_any, pred, obs):
    ls = loss_any(pred, obs)
    if ep%500==0:
        print('epoch: ', ep, ' loss: ', ls.item())

wts_bias_progression_7= list()
all_grads_7= list()
losses_7= list()
d2loss_dg7= list()
diags_7= list()
f= lambda x: np.sqrt(np.abs(np.diag(np.linalg.inv(x)))).tolist()

for epoch in range(4000):
    x_transformed =samp(X_trans, count= 32)
    trainx= Variable(torch.Tensor(x_transformed[:,:-1]).float(), requires_grad=True) #df.iloc[:, :-1].values).float()
    trainy = Variable(torch.reshape(torch.Tensor(x_transformed[:,-1]).float(), (-1,1)))

    def closure():
        optimizer7.zero_grad()
        y_pred = model7(trainx)

        l = loss7(y_pred, trainy)
        losses_7.append(l.item())
        #if epoch%500==0:
        #    print('epoch: ', epoch, ' loss: ', l.item());
        # perform a backward pass (backpropagation)
        l.backward()
        return l
    report_loss(epoch, loss7, model7(trainx), trainy)
    optimizer7.step(closure)

    wts_bias = list(model7.parameters())
    wts_bias_vals= torch.cat((wts_bias[0], torch.reshape(wts_bias[1], (-1,1))),1)
    wts_bias_progression_7.append(wts_bias_vals.tolist())

    #print('(torch.reshape(wts_bias[1].grad, (-1,1))),1):',epoch, type(wts_bias[1].grad))
    wts_bias_grads = torch.cat((wts_bias[0].grad, torch.reshape(wts_bias[1].grad, (-1,1))),1)
    all_grads_7.append(wts_bias_grads.tolist())#i_grads)

    d1loss_dg7 = torch.autograd.grad(loss7(model7(trainx), trainy), wts_bias, create_graph=True)
    d2loss_dg7.append(eval_hessian(d1loss_dg7, model7))
    diags_7.append(f(d2loss_dg7[-1]))
    
all_wts_bias_7 = np.array(wts_bias_progression_7).reshape(-1,4)
all_wts_bias_df_7= pd.DataFrame(all_wts_bias_7, columns=['weight_{}'.format(idx) for idx in range(1,4)]+ ['bias'])

all_grads_7 = np.array(all_grads_7).reshape(-1,4)
all_grads_df_7= pd.DataFrame(all_grads_7, columns=['beta_{}'.format(idx) for idx in range(1,5)])
plot(losses_7)<jupyter_output>epoch:  0  loss:  0.2515413463115692
epoch:  500  loss:  0.01371206808835268
epoch:  1000  loss:  0.008589252829551697
epoch:  1500  loss:  0.006381518207490444
epoch:  2000  loss:  0.00470005115494132
epoch:  2500  loss:  0.008985985070466995
epoch:  3000  loss:  0.005079461727291346
epoch:  3500  loss:  0.007732025347650051<jupyter_text>Pytorch model Standard Error & Parameters Comparison* **Calculating standard error for pytorch model**<jupyter_code>mini_batches_7= list()#contains mini_batches of gradients from 200 epochs
mean_mini_batches_7= list()#contains mean of each of 200 epochs each in a mini_batches of gradients from 200 epochs
mini_size= 200#mini_batch size
for batch in range(1, (len(all_grads_7))//mini_size):
    mini_batches_7.append(all_grads_7[batch*mini_size:batch*mini_size+ mini_size])
    mean_mini_batches_7.append(list(np.mean(mini_batches_7[-1], axis=0)))

grad_mini_batch_arr_7 = np.array(mean_mini_batches_7).T
grad_mini_batch_arr_7.shape
covariant_mat_7= np.cov(grad_mini_batch_arr_7)
pytorch_stderr_7= np.sqrt(covariant_mat_7.diagonal())#**2
pytorch_stderr_7
di_7 = dict(zip('pytorch_grads,pytorch_params,statsmodel_params,statsmodel_stderr'.split(','), [all_grads_7[-1,:], all_wts_bias_7[-1,:], result_2.params, statsm_df_2["std err"].values]))
di_7.update({'pytorch_stderr':pytorch_stderr_7})
comp_df_7 = pd.DataFrame(di_7)
comp_df_7['pytorch_params'] = comp_df_7['pytorch_params'].apply(lambda x:np.round(x,3))
comp_df_7['statsmodel_params'] = comp_df_7['statsmodel_params'].apply(lambda x:np.round(x,3))
comp_df_7<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15809278214877773, 'fraction_numerical': 0.0401810383390785, 'mean_word_length': 4.7529518619436875, 'pattern_counts': {'":': 0, '<': 112, '<?xml version=': 0, '>': 114, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '29', 'dir': 'jupyter-structured-clean-dedup', 'id': '554904', 'license': 'MIT', 'n_tokens_mistral': 15082, 'n_tokens_neox': 14444, 'n_words': 2331, 'path': 'notebooks/UQ/01_prakash_gradients_pytorch_statsmodel.ipynb', 'repo_name': 'mlsquare/kitchen', 'size': '38003'}