How to use the pyabc.PercentileDistanceFunction function in pyabc

y_observed = 2

    def model(args):
        return {"y": st.norm(args['x'], sigma_y).rvs()}

    models = [model]
    models = list(map(SimpleModel, models))
    model_prior = RV("randint", 0, 1)
    nr_populations = 4
    population_size = ConstantPopulationStrategy(600, nr_populations)
    parameter_given_model_prior_distribution = [Distribution(x=RV("norm", 0, sigma_x))]
    parameter_perturbation_kernels = [GridSearchCV(MultivariateNormalTransition(),
                                      {"scaling": sp.logspace(-1, 1.5, 5)})]
    abc = ABCSMC(models, model_prior, ModelPerturbationKernel(1, probability_to_stay=1),
                 parameter_given_model_prior_distribution, parameter_perturbation_kernels,
                 PercentileDistanceFunction(measures_to_use=["y"]), MedianEpsilon(.2),
                 population_size,
                 sampler=sampler)

    options = {'db_path': db_path}
    abc.set_data({"y": y_observed}, 0, {}, options)

    minimum_epsilon = -1

    abc.do_not_stop_when_only_single_model_alive()
    history = abc.run(minimum_epsilon)
    posterior_x, posterior_weight = history.get_results_distribution(0, "x")
    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)

model_prior = RV("randint", 0, 2)

    # 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=RV("norm", mu_x_1, sigma)),
                                                Distribution(x=RV("norm", mu_x_2, sigma))]

    # Particles are perturbed in a Gaussian fashion
    parameter_perturbation_kernels = [MultivariateNormalTransition() for _ in range(2)]

    # We plug all the ABC setup together
    nr_populations = 3
    population_size = AdaptivePopulationStrategy(400, 3, mean_cv=0.05)
    abc = ABCSMC(models, model_prior, ModelPerturbationKernel(2, probability_to_stay=.7),
                 parameter_given_model_prior_distribution, parameter_perturbation_kernels,
                 PercentileDistanceFunction(measures_to_use=["y"]), MedianEpsilon(.2),
                 population_size,
                 sampler=sampler)

    # Finally we add meta data such as model names and define where to store the results
    options = {'db_path': db_path}
    # y_observed is the important piece here: our actual observation.
    y_observed = 1
    abc.set_data({"y": y_observed}, 0, {}, options)

    # We run the ABC with 3 populations max
    minimum_epsilon = .05
    history = abc.run(minimum_epsilon)

    # Evaluate the model probabililties
    mp = history.get_model_probabilities(history.max_t)

y_observed = 1

    def model(args):
        return {"y": st.norm(args['x'], sigma_y).rvs()}

    models = [model, model]
    models = list(map(SimpleModel, models))
    model_prior = RV("randint", 0, 2)
    population_size = ConstantPopulationStrategy(500, 1)
    mu_x_1, mu_x_2 = 0, 1
    parameter_given_model_prior_distribution = [Distribution(x=RV("norm", mu_x_1, sigma_x)),
                                                Distribution(x=RV("norm", mu_x_2, sigma_x))]
    parameter_perturbation_kernels = [MultivariateNormalTransition() for _ in range(2)]
    abc = ABCSMC(models, model_prior, ModelPerturbationKernel(2, probability_to_stay=.7),
                 parameter_given_model_prior_distribution, parameter_perturbation_kernels,
                 PercentileDistanceFunction(measures_to_use=["y"]), MedianEpsilon(.02),
                 population_size,
                 sampler=sampler)

    options = {'db_path': db_path}
    abc.set_data({"y": y_observed}, 0, {}, options)

    minimum_epsilon = -1
    nr_populations = 1
    abc.do_not_stop_when_only_single_model_alive()
    history = abc.run(minimum_epsilon)
    mp = history.get_model_probabilities(history.max_t)

    def p_y_given_model(mu_x_model):
        return st.norm(mu_x_model, sp.sqrt(sigma_y**2+sigma_x**2)).pdf(y_observed)

    p1_expected_unnormalized = p_y_given_model(mu_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))
    model_prior = RV("randint", 0, 1)
    nr_populations = 4
    population_size = ConstantPopulationStrategy(600, nr_populations)
    parameter_given_model_prior_distribution = [Distribution(x=RV("norm", 0, sigma_x))]
    parameter_perturbation_kernels = [MultivariateNormalTransition()]
    abc = ABCSMC(models, model_prior, ModelPerturbationKernel(1, probability_to_stay=1),
                 parameter_given_model_prior_distribution, parameter_perturbation_kernels,
                 PercentileDistanceFunction(measures_to_use=["y"]), MedianEpsilon(.2),
                 population_size,
                 sampler=sampler)

    options = {'db_path': db_path}
    abc.set_data({"y": y_observed}, 0, {}, options)

    minimum_epsilon = -1

    abc.do_not_stop_when_only_single_model_alive()
    history = abc.run(minimum_epsilon)
    posterior_x, posterior_weight = history.get_results_distribution(0, "x")
    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)

model_prior = RV("randint", 0, 2)

    # 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=RV("norm", mu_x_1, sigma)),
                                                Distribution(x=RV("norm", mu_x_2, sigma))]

    # Particles are perturbed in a Gaussian fashion
    parameter_perturbation_kernels = [MultivariateNormalTransition() for _ in range(2)]

    # We plug all the ABC setup together
    nr_populations = 3
    population_size = ConstantPopulationStrategy(400, 3)
    abc = ABCSMC(models, model_prior, ModelPerturbationKernel(2, probability_to_stay=.7),
                 parameter_given_model_prior_distribution, parameter_perturbation_kernels,
                 PercentileDistanceFunction(measures_to_use=["y"]), MedianEpsilon(.2),
                 population_size,
                 sampler=sampler)

    # Finally we add meta data such as model names and define where to store the results
    options = {'db_path': db_path}
    # y_observed is the important piece here: our actual observation.
    y_observed = 1
    abc.set_data({"y": y_observed}, 0, {}, options)

    # We run the ABC with 3 populations max
    minimum_epsilon = .05
    history = abc.run(minimum_epsilon)

    # Evaluate the model probabililties
    mp = history.get_model_probabilities(history.max_t)

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))
    model_prior = RV("randint", 0, 1)
    nr_populations = 4
    population_size = AdaptivePopulationStrategy(600, nr_populations)
    parameter_given_model_prior_distribution = [Distribution(x=RV("norm", 0, sigma_x))]
    parameter_perturbation_kernels = [MultivariateNormalTransition()]
    abc = ABCSMC(models, model_prior, ModelPerturbationKernel(1, probability_to_stay=1),
                 parameter_given_model_prior_distribution, parameter_perturbation_kernels,
                 PercentileDistanceFunction(measures_to_use=["y"]), MedianEpsilon(.2),
                 population_size,
                 sampler=sampler)

    options = {'db_path': db_path}
    abc.set_data({"y": y_observed}, 0, {}, options)

    minimum_epsilon = -1

    abc.do_not_stop_when_only_single_model_alive()
    history = abc.run(minimum_epsilon)
    posterior_x, posterior_weight = history.get_results_distribution(0, "x")
    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)

def test_continuous_non_gaussian(db_path, sampler):
    def model(args):
        return {"result": sp.rand() * args['u']}

    models = [model]
    models = list(map(SimpleModel, models))
    model_prior = RV("randint", 0, 1)
    population_size = ConstantPopulationStrategy(250, 2)
    parameter_given_model_prior_distribution = [Distribution(u=RV("uniform", 0, 1))]
    parameter_perturbation_kernels = [MultivariateNormalTransition()]
    abc = ABCSMC(models, model_prior, ModelPerturbationKernel(1, probability_to_stay=1),
                 parameter_given_model_prior_distribution, parameter_perturbation_kernels,
                 PercentileDistanceFunction(measures_to_use=["result"]), MedianEpsilon(.2),
                 population_size,
                 sampler=sampler)

    options = {'db_path': db_path}
    d_observed = .5
    abc.set_data({"result": d_observed}, 0, {}, options)
    abc.do_not_stop_when_only_single_model_alive()

    minimum_epsilon = -1
    history = abc.run(minimum_epsilon)
    posterior_x, posterior_weight = history.get_results_distribution(0, "u")
    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)))

    @sp.vectorize

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))
    model_prior = RV("randint", 0, 1)
    nr_populations = 1
    population_size = ConstantPopulationStrategy(600, nr_populations)
    parameter_given_model_prior_distribution = [Distribution(x=RV("norm", 0, sigma_prior))]
    parameter_perturbation_kernels = [MultivariateNormalTransition()]
    abc = ABCSMC(models, model_prior, ModelPerturbationKernel(1, probability_to_stay=1),
                 parameter_given_model_prior_distribution, parameter_perturbation_kernels,
                 PercentileDistanceFunction(measures_to_use=["y"]), MedianEpsilon(.1),
                 population_size,
                 sampler=sampler)

    options = {'db_path': db_path}
    abc.set_data({"y": observed_data}, 0, {}, options)

    minimum_epsilon = -1


    abc.do_not_stop_when_only_single_model_alive()
    history = abc.run(minimum_epsilon)
    posterior_x, posterior_weight = history.get_results_distribution(0, "x")
    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)))

model_prior = pyabc.RV("randint", 0, 1)
population_size = pyabc.AdaptivePopulationStrategy(500, 20,
                                                   max_population_size=10000)



mapper = parallel.SGE().map if parallel.sge_available() else map
abc = pyabc.ABCSMC([pyabc.SimpleModel(abc_model)],
                    model_prior,
                    pyabc.ModelPerturbationKernel(1, probability_to_stay=.8),
                    [ABCPrior()],
                    [pyabc.MultivariateNormalTransition()],
                    pyabc.PercentileDistanceFunction(measures_to_use=["x", "y"]),
                    pyabc.MedianEpsilon(),
                    population_size,
                    sampler=parallel.sampler.MappingSampler(map=mapper))
abc.stop_if_only_single_model_alive = False




options = {'db_path': "sqlite:///" + sm.output[0]}
abc.set_data({"x": 1, "y": 1}, 0, {}, options)




history = abc.run(.01)