How to use the pyabc.ModelPerturbationKernel function in pyabc

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))
    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)))

def test_all_in_one_model(db_path, sampler):
    models = [AllInOneModel() for _ in range(2)]
    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)
    mp = history.get_model_probabilities(history.max_t)
    assert abs(mp.p[0] - .5) + abs(mp.p[1] - .5) < .08

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)
    mp = history.get_model_probabilities(history.max_t)
    assert abs(mp.p[0] - .5) + abs(mp.p[1] - .5) < .08

def test_beta_binomial_two_identical_models_adaptive(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 = AdaptivePopulationStrategy(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)
    mp = history.get_model_probabilities(history.max_t)
    assert abs(mp.p[0] - .5) + abs(mp.p[1] - .5) < .08

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 = AdaptivePopulationStrategy(1500, 3)
    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(0),
                 population_size,
                 sampler=sampler)

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

    minimum_epsilon = -1
    history = abc.run(minimum_epsilon)
    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) < .1

def test_all_in_one_model(self):
        models = [AllInOneModel() for _ in range(2)]
        model_prior = RV("randint", 0, 2)
        mp_pool = multiprocessing.Pool(4)
        mp_sampler = MappingSampler(map=mp_pool.map)
        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=mp_sampler)

        options = {'db_path': self.db}
        abc.set_data({"result": 2}, 0, {}, options)

        minimum_epsilon = .2
        history = abc.run(minimum_epsilon)
        mp = history.get_model_probabilities(history.max_t)
        # self.assertLess(abs(p1 - .5) + abs(p2 - .5), .08)
        self.assertLess(abs(mp.p[0] - .5) + abs(mp.p[1] - .5), .08*5) # Dennis

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, 3)
    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(0),
                 population_size,
                 sampler=sampler)

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

    minimum_epsilon = -1
    history = abc.run(minimum_epsilon)



    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_beta_binomial_different_priors(db_path, sampler):
    binomial_n = 5

    def model(args):
        return {"result": st.binom(binomial_n, args['theta']).rvs()}

    models = [model for _ in range(2)]
    models = list(map(SimpleModel, models))
    model_prior = RV("randint", 0, 2)
    population_size = ConstantPopulationStrategy(800, 3)
    a1, b1 = 1, 1
    a2, b2 = 10, 1
    parameter_given_model_prior_distribution = [Distribution(theta=RV("beta", a1, b1)),
                                                Distribution(theta=RV("beta", a2, b2))]
    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}
    n1 = 2
    abc.set_data({"result": n1}, 0, {}, options)

    minimum_epsilon = .2
    history = abc.run(minimum_epsilon)
    mp = history.get_model_probabilities(history.max_t)

    def B(a, b):
        return gamma(a) * gamma(b) / gamma(a + b)

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)

]

# Particles are perturbed in a Gaussian fashion
parameter_perturbation_kernels = [
    lambda t, stat: Kernel(stat['cov']) for _ in range(2)
    ]

# Initialize distributed mapper:
my_distributed_mapper = MapWrapperDistribDemo(None, "Kilroy was here")

# We plug all the ABC setup together.
# We use "y" in the distance function as this
# was the result key defined for the model
nr_particles = 400
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), nr_particles,
             max_nr_allowed_sample_attempts_per_particle=2000,
             mapper=map, sampler=my_distributed_mapper)

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

# We run the ABC with 3 populations max