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)
In [ ]: