Dealing with multivariate prior distributions (Dirichlet)

[MadathilSA]

Hi,

I am trying to implement an SMCABC for with a multivariate prior. Specifically Dirichlet prior. However, the weight calculation at run gives an error. What am I doing wrong here?
Have multivariate cases implemented in pyABC?

Here is a toy example, following the univariate example in the documentation.

import os
import tempfile

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import dirichlet
import pyabc

pyabc.settings.set_figure_params('pyabc')  # for beautified plots

A toy model:

def model(parameter):
    return {"data": (parameter['mu'][0,0] - parameter['mu'][0,1])*10 }

Bivariate prior:

prior = pyabc.Distribution(mu=dirichlet(alpha=np.array([1,1])))

Distance, configuring abc and data

def distance(x, x0):
    return abs(x["data"] - x0["data"])

abc = pyabc.ABCSMC(model, prior, distance, population_size=1000)

db_path = os.path.join(tempfile.gettempdir(), "test.db")
observation = 8
abc.new("sqlite:///" + db_path, {"data": observation})

run:

history = abc.run(minimum_epsilon=0.1, max_nr_populations=10)

ABC INFO: Calibration sample t = -1.
INFO:ABC:Calibration sample t = -1.
ABC INFO: t: 0, eps: 7.90139925e+00.
INFO:ABC:t: 0, eps: 7.90139925e+00.
ABC.History INFO: Done <ABCSMC id=11, duration=0:00:04.589834, end_time=2023-06-06 01:11:29>
INFO:ABC.History:Done <ABCSMC id=11, duration=0:00:04.589834, end_time=2023-06-06 01:11:29>
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
[<ipython-input-7-427c0577e326>](https://localhost:8080/#) in <cell line: 1>()
----> 1 history = abc.run(minimum_epsilon=0.1, max_nr_populations=10)

3 frames
[/usr/local/lib/python3.10/dist-packages/pyabc/population/population.py](https://localhost:8080/#) in normalize_weights(self)
    405         total_weight_accepted = sum(p.weight for p in self.accepted_particles)
    406 
--> 407         if np.isclose(total_weight_accepted, 0.0):
    408             logger.warning(
    409                 f"The total population weight {total_weight_accepted} is "

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Thank you very much for your time.

[yannikschaelte]

Hi @MadathilSA , you are correct that pyABC's Distribution+RV concept is primarily designed for independent univariate priors. It is however possible to subclass the DistributionBase class, to implement own functionality for random variable sampling (rvs()) and density evaluation (pdf(x)). I believe that this way wrapping around a Dirichlet distribution should be possible? See here for an example: https://pyabc.readthedocs.io/en/latest/examples/custom_priors.html

[MadathilSA]

Hi @yannikschaelte, thank you very much for this lead. We will try this approach.