Alternatively, you can try using ADVI, which seems to give a pretty sensible local optimal:
with simple_model:
approx = pm.fit(100000, obj_n_mc=10)
z_approx = simple_model.bijection.rmap(approx.params[0].eval())
simple_model.logp(z_approx)
# ==> array(-10848.99051671)
point = simple_model.test_point
point['z'] = z_train
simple_model.logp(point)
# ==> array(-10847.51765264)