ABC

In [1]:
import pymc as pm
In [132]:
obs = np.array( [0,0,1] )

uni = pm.Uniform("prop", 0, 1)
fake_obs = pm.Bernoulli("fake_obs", uni, size=3 )

@pm.deterministic
def accept(uni = uni, fake_obs = fake_obs, obs = obs):
    if np.array_equal( fake_obs, obs ):
        return uni
    else:
        return None
In [133]:
mcmc = pm.MCMC([uni, fake_obs, accept] )
In [134]:
mcmc.sample(10000)
[****************100%******************]  10000 of 10000 complete
In [135]:
samples = mcmc.trace('accept')[:]
In [136]:
hist( samples[samples > 0] )
Out[136]:
(array([  47.,  118.,  158.,  125.,  141.,   95.,   74.,   48.,   33.,   10.]),
 array([ 0.01977077,  0.11150485,  0.20323893,  0.29497301,  0.38670709,
        0.47844116,  0.57017524,  0.66190932,  0.7536434 ,  0.84537748,
        0.93711155]),
 <a list of 10 Patch objects>)
In [104]:
samples[:1000]
Out[104]:
array([None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       0.641398409499, None, None, None, None, None, None, None, None,
       None, None, None, None, 0.346642876027, None, None, 0.307153130839,
       None, None, None, None, None, None, None, None, None, 0.14974390582,
       None, None, 0.760792945516, None, None, None, None, 0.273355388346,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, 0.762049281736, None, None, None, None, None, None,
       None, None, None, None, None, None, 0.386238065957, None,
       0.474620548494, None, None, None, None, None, None, None,
       0.233549246555, None, None, None, None, None, None, None,
       0.25702449777, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, 0.219752257765, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None,
       0.30753209596, None, None, None, None, None, None, None,
       0.604963812539, None, None, None, None, None, None, None, None,
       None, None, None, 0.600129961666, None, None, None, 0.423292896498,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, 0.416366401192, None, None, None, None,
       None, None, None, 0.206467531249, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None,
       0.827333806075, 0.367881810314, None, None, None, 0.463397296568,
       None, None, None, None, None, None, None, None, None, None, None,
       0.210841925823, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, 0.540566216778, None, None, 0.448174456742,
       None, 0.392659849708, 0.210200730086, None, None, None, None, None,
       None, None, None, None, None, None, 0.583595071636, None, None,
       None, None, None, None, None, None, None, 0.263724512974, None,
       None, None, None, None, None, None, None, 0.310354018644, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, 0.707879317692, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, 0.274577981416, None, None, None, None, None,
       0.632703271583, None, None, None, None, None, None, None, None,
       0.244699027015, None, None, 0.703190955046, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, 0.339908119414, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None,
       0.451298757871, None, None, 0.653381413287, None, None, None, None,
       None, None, None, None, None, None, None, None, 0.602144732268,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, 0.445724522823, None, None,
       None, None, None, None, None, 0.315490930854, None, None, None,
       None, None, None, None, None, None, None, None, None,
       0.318855322497, None, None, None, 0.717182055351, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, 0.670501111935, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, 0.363765665679, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, 0.305447088556, None, None, None,
       None, None, None, None, None, None, None, None, None,
       0.282125872166, None, None, None, None, None, None, None, None,
       None, None, None, None, None, 0.251070061765, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       0.50442411879, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, 0.641889799058, None, None, None, 0.178913518161, None, None,
       None, None, None, None, 0.194966290843, None, None, None, None,
       None, None, None, None, None, None, None, None, 0.59455012334, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, 0.153238190081, None, None, None, None,
       None, None, None, None, None, None, None, 0.500331587088, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None,
       0.443738792346, None, None, None, None, None, 0.530030253902, None,
       None, None, None, None, None, None, None, None, None,
       0.372842415902, None, None, None, None, None, None, 0.690298120711,
       None, None, 0.750438143415, None, None, 0.280714706648, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, 0.350938621096, None, None,
       None, None, None, None, None, 0.714571072204, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, 0.897063779285, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None,
       0.319977971294, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, 0.141266113341, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, None, None, None, None, None, None, None, None, None,
       None, None, 0.0703177567736, None, None, None, None, None, None], dtype=object)