Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def test_gaussian_multiple_populations_crossval_kde(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 = ConstantPopulationSize(600)
parameter_given_model_prior_distribution = [Distribution(x=st.norm(0, sigma_x))]
parameter_perturbation_kernels = [GridSearchCV(MultivariateNormalTransition(),
{"scaling": np.logspace(-1, 1.5, 5)})]
abc = ABCSMC(models, parameter_given_model_prior_distribution,
MinMaxDistance(measures_to_use=["y"]),
population_size,
transitions=parameter_perturbation_kernels,
eps=MedianEpsilon(.2),
sampler=sampler)
abc.new(db_path, {"y": y_observed})
minimum_epsilon = -1
abc.do_not_stop_when_only_single_model_alive()
def test_cookie_jar(db_path, sampler):
def make_model(theta):
def model(args):
return {"result": 1 if random.random() > theta else 0}
return model
theta1 = .2
theta2 = .6
model1 = make_model(theta1)
model2 = make_model(theta2)
models = [model1, model2]
models = list(map(SimpleModel, models))
population_size = ConstantPopulationSize(1500)
parameter_given_model_prior_distribution = [Distribution(), Distribution()]
abc = ABCSMC(models, parameter_given_model_prior_distribution,
MinMaxDistance(measures_to_use=["result"]),
population_size, eps=MedianEpsilon(.1), sampler=sampler)
abc.new(db_path, {"result": 0})
minimum_epsilon = .2
history = abc.run(minimum_epsilon, max_nr_populations=1)
mp = history.get_model_probabilities(history.max_t)
expected_p1, expected_p2 = theta1 / (theta1 + theta2), theta2 / (theta1 +
theta2)
assert abs(mp.p[0] - expected_p1) + abs(mp.p[1] - expected_p2) < .05
def test_cookie_jar(db_path, sampler):
def make_model(theta):
def model(args):
return {"result": 1 if random.random() > theta else 0}
return model
theta1 = .2
theta2 = .6
model1 = make_model(theta1)
model2 = make_model(theta2)
models = [model1, model2]
models = list(map(SimpleModel, models))
model_prior = RV("randint", 0, 2)
population_size = ConstantPopulationStrategy(1500, 1)
parameter_given_model_prior_distribution = [Distribution(), Distribution()]
parameter_perturbation_kernels = [MultivariateNormalTransition() for _ in range(2)]
abc = ABCSMC(models, model_prior, ModelPerturbationKernel(2, probability_to_stay=.8),
parameter_given_model_prior_distribution,
parameter_perturbation_kernels,
MinMaxDistanceFunction(measures_to_use=["result"]),
MedianEpsilon(.1),
population_size,
sampler=sampler)
options = {'db_path': db_path}
abc.set_data({"result": 0}, 0, {}, options)
minimum_epsilon = .2
def test_two_competing_gaussians_multiple_population(db_path, sampler):
# Define a gaussian model
sigma = .5
def model(args):
return {"y": st.norm(args['x'], sigma).rvs()}
# We define two models, but they are identical so far
models = [model, model]
models = list(map(SimpleModel, models))
# However, our models' priors are not the same. Their mean differs.
mu_x_1, mu_x_2 = 0, 1
parameter_given_model_prior_distribution = [
Distribution(x=st.norm(mu_x_1, sigma)),
Distribution(x=st.norm(mu_x_2, sigma))
]
# We plug all the ABC setup together
nr_populations = 1
population_size = ConstantPopulationSize(20)
abc = ABCSMC(models, parameter_given_model_prior_distribution,
PercentileDistance(measures_to_use=["y"]),
population_size,
eps=MedianEpsilon(.2),
sampler=sampler)
def test_beta_binomial_two_identical_models(db_path, sampler):
binomial_n = 5
def model_fun(args):
return {"result": st.binom(binomial_n, args.theta).rvs()}
models = [model_fun for _ in range(2)]
models = list(map(SimpleModel, models))
model_prior = RV("randint", 0, 2)
population_size = ConstantPopulationStrategy(800, 3)
parameter_given_model_prior_distribution = [Distribution(theta=RV("beta", 1, 1))
for _ in range(2)]
parameter_perturbation_kernels = [MultivariateNormalTransition() for _ in range(2)]
abc = ABCSMC(models, model_prior, ModelPerturbationKernel(2, probability_to_stay=.8),
parameter_given_model_prior_distribution, parameter_perturbation_kernels,
MinMaxDistanceFunction(measures_to_use=["result"]), MedianEpsilon(.1),
population_size,
sampler=sampler)
options = {'db_path': db_path}
abc.set_data({"result": 2}, 0, {}, options)
minimum_epsilon = .2
history = abc.run( minimum_epsilon)
# prior
prior = pyabc.Distribution(**{key: pyabc.RV('uniform', 0, 10)
for key in theta0})
# distance
distance = pyabc.AdaptivePNormDistance()
# acceptor
acceptor = pyabc.accept_use_complete_history
# abc
db_name = "sqlite:///" + os.path.join(tempfile.gettempdir(), "tmp.db")
abc = pyabc.ABCSMC(models=pyabc.SimpleModel(data_gk),
parameter_priors=prior,
distance_function=distance,
summary_statistics=ordered_statistics_gk,
population_size=100,
acceptor=acceptor)
abc.new(db=db_name, observed_sum_stat=obs_sum_stats)
abc.run(minimum_epsilon=0, max_nr_populations=10)
# visualization
df, w = abc.history.get_distribution(m=0)
pyabc.visualization.plot_kde_matrix(df, w, limits={key: (0,10)
for key in theta0})
plt.show()
return arr
def model_2(pars):
rate = pars.rate
arr = sp.rand(4)
return arr
def distance(x, y):
return ((x - y)**2).sum()
mapper = parallel.SGE().map if parallel.sge_available() else map
abc = pyabc.ABCSMC([pyabc.SimpleModel(model_1),
pyabc.SimpleModel(model_2)],
model_prior,
pyabc.ModelPerturbationKernel(2, probability_to_stay=.8),
[rate_prior, rate_prior],
[pyabc.MultivariateNormalTransition(),
pyabc.MultivariateNormalTransition()],
distance,
pyabc.MedianEpsilon(),
population_strategy,
sampler=parallel.sampler.MappingSampler(map=mapper))
abc.stop_if_only_single_model_alive = False
options = {'db_path': "sqlite:///" + sm.output[0]}
abc.set_data(sp.rand(4), 0, {}, options)
history = abc.run(.01)
arr = sp.rand(4)
return arr
def model_2(pars):
rate = pars.rate
arr = sp.rand(4)
return arr
def distance(x, y):
return ((x - y)**2).sum()
mapper = parallel.SGE().map if parallel.sge_available() else map
abc = pyabc.ABCSMC([pyabc.SimpleModel(model_1),
pyabc.SimpleModel(model_2)],
model_prior,
pyabc.ModelPerturbationKernel(2, probability_to_stay=.8),
[rate_prior, rate_prior],
[pyabc.MultivariateNormalTransition(),
pyabc.MultivariateNormalTransition()],
distance,
pyabc.MedianEpsilon(),
population_strategy,
sampler=parallel.sampler.MappingSampler(map=mapper))
abc.stop_if_only_single_model_alive = False
options = {'db_path': "sqlite:///" + sm.output[0]}
abc.set_data(sp.rand(4), 0, {}, options)
history = abc.run(.01)