Posterior_Predictive_Check_MCMC_ex.py

import pymc as pm
from pymc.step_methods import Metropolis
from pymc import Poisson
from pymc import Normal,Uniform
from pymc import sample, Model
from pymc import find_MAP
from pymc import sample_posterior_predictive
from arviz import plot_trace, plot_autocorr, summary

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import poisson


#####Synthetic Data Generation

np.random.seed(321)
X1 = np.random.randn(50).reshape(-1,1)
X2 = np.random.randn(50).reshape(-1,1)
eps = np.random.normal(0,0.0001,50)
X = np.concatenate((X1,X2),axis=1)
Beta = np.array([0.6,-0.3])
lamb = np.exp(np.matmul(X,Beta)+eps)
Y = np.zeros(50)
for i in range(len(Y)):
  Y[i] = poisson.rvs(lamb[i],size=1)[0]
data = pd.DataFrame({'X1':X1.reshape(-1),'X2':X2.reshape(-1),'Y':Y})

#####Train Test Split
from sklearn.model_selection import train_test_split

from sklearn.model_selection import train_test_split
X = data.iloc[:,:2]
y = data.iloc[:,-1]
X_tr,X_te,y_tr,y_te = train_test_split(X,y,test_size=0.3,shuffle=True, random_state=321)



#####PyMC training

with Model() as python_ex:
  beta1 = Normal('beta1',0,2)
  beta2 = Normal('beta2',0,2)
  X_1 = pm.Data('X_1',X_tr['X1'].values)
  X_2 = pm.Data('X_2',X_tr['X2'].values)
  theta = np.exp(beta1*X_1+beta2*X_2)
  Lk = Poisson('Y',mu=theta,observed=y_tr,shape=theta.shape)
with python_ex:
  results = sample(step=pm.HamiltonianMC(),draws=1500,chains=2, tune=500, random_seed=123)


plot_trace(results,compact=True,var_names=('beta1','beta2'),chain_prop={'color': ['darkcyan', 'red']})
summary(results)

#####PPC analysis(Train data)

with python_ex:
  PPC = pm.sample_posterior_predictive(trace=results,random_seed=123)


from arviz import plot_ppc


from arviz import plot_ppc
sns.set_style("darkgrid")
plot_ppc(data=PPC,kind="kde",observed_rug=True,mean=True,colors=['blue','red','orange'])


#####PPC analysis(Test data)
with python_ex:
  pm.set_data({'X_1':X_te['X1'].values})
  pm.set_data({'X_2':X_te['X2'].values})
  PPC_te = pm.sample_posterior_predictive(trace=results,random_seed=123)

sns.set_style("darkgrid")
plot_ppc(data=PPC_te,kind="kde",mean=True,observed=False)