PyMC consistently under estimating results found in paper. Possibly not sampling enough? 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

 A: 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:

