def test_gaussian_single_population(db_path, sampler):
sigma_prior = 1
sigma_ground_truth = 1
observed_data = 1
def model(args):
return {"y": st.norm(args['x'], sigma_ground_truth).rvs()}
models = [model]
models = list(map(SimpleModel, models))
nr_populations = 1
population_size = ConstantPopulationSize(600)
parameter_given_model_prior_distribution = [Distribution(x=RV("norm", 0,
sigma_prior))
]
abc = ABCSMC(models, parameter_given_model_prior_distribution,
MinMaxDistanceFunction(measures_to_use=["y"]),
population_size,
eps=MedianEpsilon(.1),
sampler=sampler)
abc.new(db_path, {"y": observed_data})
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_prior**2 + 1 / sigma_ground_truth**2)
mu_x_given_y = sigma_x_given_y**2 * observed_data / sigma_ground_truth**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.12
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) < .1
评论列表
文章目录