def test_gaussian_multiple_populations_adpative_population_size(db_path, sampler):
sigma_x = 1
sigma_y = .5
y_observed = 2
def model(args):
return {"y": st.norm(args['x'], sigma_y).rvs()}
models = [model]
models = list(map(SimpleModel, models))
nr_populations = 4
population_size = AdaptivePopulationSize(600)
parameter_given_model_prior_distribution = [
Distribution(x=st.norm(0, sigma_x))]
abc = ABCSMC(models, parameter_given_model_prior_distribution,
MinMaxDistanceFunction(measures_to_use=["y"]),
population_size,
eps=MedianEpsilon(.2),
sampler=sampler)
abc.new(db_path, {"y": y_observed})
minimum_epsilon = -1
abc.do_not_stop_when_only_single_model_alive()
history = abc.run(minimum_epsilon, max_nr_populations=nr_populations)
posterior_x, posterior_weight = history.get_distribution(0, None)
posterior_x = posterior_x["x"].as_matrix()
sort_indices = sp.argsort(posterior_x)
f_empirical = sp.interpolate.interp1d(sp.hstack((-200, posterior_x[sort_indices], 200)),
sp.hstack((0, sp.cumsum(posterior_weight[sort_indices]), 1)))
sigma_x_given_y = 1 / sp.sqrt(1 / sigma_x ** 2 + 1 / sigma_y ** 2)
mu_x_given_y = sigma_x_given_y ** 2 * y_observed / sigma_y ** 2
expected_posterior_x = st.norm(mu_x_given_y, sigma_x_given_y)
x = sp.linspace(-8, 8)
max_distribution_difference = sp.absolute(f_empirical(x) - expected_posterior_x.cdf(x)).max()
assert max_distribution_difference < 0.15
assert history.max_t == nr_populations - 1
mean_emp, std_emp = mean_and_std(posterior_x, posterior_weight)
assert abs(mean_emp - mu_x_given_y) < .07
assert abs(std_emp - sigma_x_given_y) < .12
评论列表
文章目录