I have been trying to build confidence in (my ability to correctly use) PyMC by working examples. Namely, I have been working on Chickering and Pearl 1997, and more specifically on their 'artificial' and 'lipid' data sets.
Issue / Statement of my problem: For a symmetric data set, my model listed below estimates ACE( D -> Y ) at zero as expected. For the 'artificial' and 'lipid' data sets found in Table 1 and Table 2, my model constantly underestimates ACE( D -> Y ). For example, Chickering and Pearl insist the posterior distribution collapses around .55, in my model it seems to be collapsing closer to .35. Similar results were found for the Lipid data set, but I didn't save the results, close to .42 as I recall.
- Improperly specified model. I can't have the deterministic function d return '1 or 0' as probabilities to be passed to D, so I set them very highly(lowly) at .9999(.0001).
Incorrectly tuned sampling. I'm doing this exercise to gain familiarity with PyMC, and while I have tried adjusting the number of iteration/burn in period, this takes forever to run and so the most iterations I have run is 1000. Am I not close to the number of iterations required? I also haven't used thinning bc Chickering and Pearl declare the chain is Ergodic.
import pymc as pm import numpy as np import matplotlib.pyplot as plt from pprint import pprint import pandas as pd def make_dataset(counts): """ counts: a list of 8 integers returns: a dataset to be fed to the model_factory """ zdy = [tuple(int(i) for i in list(bin(j).split('b').zfill(3))) for j in range(8)] return [obs for group in [ (zdy[i],) * counts[i] for i in range(len(zdy))] for obs in group] def model_factory(list_data): df_data = pd.DataFrame(list_data, columns=['z','d','y']) N = df_data.shape Z = pm.Bernoulli('Z', p=.5, value=df_data['z'].values, observed=True) C = np.empty(N, dtype=object) R = np.empty(N, dtype=object) for m in range(N): C[m] = pm.DiscreteUniform('c%i' % m, lower=0, upper=3) R[m] = pm.DiscreteUniform('r%i' % m, lower=0, upper=3) """ Equation (2) """ @pm.deterministic def d(Z=Z, C=C): return np.where( ( (C==3) | ((Z == False) & (C==2)) | ((Z== True) & (C==1)) ) , .9999, .0001 ) """ Equation (3) """ @pm.deterministic def y(d=d, R=R): return np.where( ( (R==3) | ((d == .0001) & (R==2)) | ((d== .9999) & (R==1)) ) , .9999, .0001 ) D = pm.Bernoulli('D', p=d, value=df_data['d'].values.astype(bool), observed=True) Y = pm.Bernoulli('Y', p=y, value=df_data['y'].values.astype(bool), observed=True) @pm.deterministic def v_r_1(R=R,d=d): return float( sum(np.where( (R==1) , 1, 0 ) )) / N @pm.deterministic def v_r_2(R=R,d=d): return float( sum(np.where( (R==2) , 1, 0 ) )) / N """ Equation (4) """ @pm.deterministic def ACE(v_r_1=v_r_1, v_r_2=v_r_2): return (v_r_1 - v_r_2) return locals() """ Chickering and Pearl's 'artifical' dataset """ arti_counts = [275, 0, 225, 0, 225, 0, 0, 275] arti_ds = make_dataset(arti_counts) arti_model = pm.MCMC(model_factory(arti_ds)) arti_model.sample(300,100) binwidth=.005 plot_data= arti_model.trace('ACE')[:] plt.hist(plot_data, bins=np.arange(min(plot_data), max(plot_data) + binwidth, binwidth)); plt.show() """ Real Data Example: Effect of Cholestyramine on Reduced Cholesterol Data Set """ lipid_counts = [158, 14, 0, 0, 52, 12, 23, 78] lipid_ds = make_dataset(lipid_counts) lipid_model = pm.MCMC(model_factory(lipid_ds)) lipid_model.sample(100,25) binwidth=.01 plot_data= lipid_model.trace('ACE')[:] plt.hist(plot_data, bins=np.arange(min(plot_data), max(plot_data) + binwidth, binwidth)); plt.show()
Adding Diagnostic Image