1
$\begingroup$

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.

Possible explanations:

  • 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')[1].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[0]
        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 Artificial Data Set ACE

$\endgroup$
1
$\begingroup$

I don't know how to get the .55 value from the paper, either, but I do know a way to speed up your model so that you can make sure MCMC convergence is not the culprit.

The key is to have a computation graph that can be updated easily during the MCMC steps. Here is a change to your model factory that speeds things up:

def model_factory(list_data):
    df_data     = pd.DataFrame(list_data, columns=['z','d','y'])
    N   =   df_data.shape[0]
    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)
    d   =   np.empty(N, dtype=object)
    y   =   np.empty(N, dtype=object)
    D   =   np.empty(N, dtype=object)
    Y   =   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(name='d%i' % m)
        def dm(Z=Z[m], C=C[m]):
            return np.where( ( (C==3) | ((Z == False) & (C==2)) | ((Z== True) & (C==1)) ) , .9999, .0001 )
        d[m] = dm

        """ Equation (3) """
        @pm.deterministic(name='y%i' % m)
        def ym(d=d[m], R=R[m]):
            return np.where( ( (R==3) | ((d == .0001) & (R==2)) | ((d== .9999) & (R==1)) ) , .9999, .0001 )
        y[m] = ym

        D[m]   =   pm.Bernoulli('D%i'%m, p=d[m], value=df_data['d'].values[m].astype(bool), observed=True)
        Y[m]   =   pm.Bernoulli('Y%i'%m, p=y[m], value=df_data['y'].values[m].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()

I also like to make sure that PyMC is using an appropriate step method:

arti_model = pm.MCMC(model_factory(arti_ds))

for i in range(1000):
    arti_model.use_step_method(pm.DiscreteMetropolis, arti_model.C[i])
    arti_model.use_step_method(pm.DiscreteMetropolis, arti_model.R[i])

arti_model.sample(iter=30000, burn=10000, thin=10)

With this faster model and longer chain, I find a mean ACE of .27, after 1-2 hours of computation. Too bad it is not 0.55. Here is the summary plot, which looks nice and converged:

ACE summary plot

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.