this PR paves the way to allow users to add simulations without a reference posterior. This way, (I hope) it becomes easier for users to test drive and benchmark sbi for their use case even.

Closes #19

Thanks for building out and maintaining this package! There was definitely a need for something like this in the ABC/Likelihood Free community.

I'm hitting a seemingly stray dimension in mcabc.py:

from sbibm.algorithms import rej_abc 
task = sbibm.get_task("two_moons")
posterior_samples, _, _ = rej_abc(task=task, num_samples=10_000, num_observation=1, num_simulations=100_000)

which is returning a stacktrace like:

ValueError                                Traceback (most recent call last)
<ipython-input-128-10fe8b131cec> in <module>
      1 from sbibm.algorithms import rej_abc
      2 task = sbibm.get_task("two_moons")
----> 3 posterior_samples, _, _ = rej_abc(task=task, num_samples=10_000, num_observation=1, num_simul
ations=100_000)

~/.pyenv/versions/miniforge3-4.9.2/lib/python3.8/site-packages/sbibm/algorithms/sbi/mcabc.py in run(t
ask, num_samples, num_simulations, num_observation, observation, num_top_samples, quantile, eps, dist
ance, batch_size, save_distances, kde_bandwidth, sass, sass_fraction, sass_feature_expansion_degree,
lra)
    118     if num_observation is not None:
    119         true_parameters = task.get_true_parameters(num_observation=num_observation)
--> 120         log_prob_true_parameters = posterior.log_prob(true_parameters)
    121         return samples, simulator.num_simulations, log_prob_true_parameters
    122     else:

~/.pyenv/versions/miniforge3-4.9.2/lib/python3.8/site-packages/pyro/distributions/empirical.py in log
_prob(self, value)
     94         if self._validate_args:
     95             if value.shape != self.batch_shape + self.event_shape:
---> 96                 raise ValueError("``value.shape`` must be {}".format(self.batch_shape + self.
event_shape))
     97         if self.batch_shape:
     98             value = value.unsqueeze(self._aggregation_dim)

ValueError: ``value.shape`` must be torch.Size([2])

A bit of digging shows that the shape of true_parameters in this is coming out at [1,2]. Changing this line to log_prob_true_parameters = posterior.log_prob(true_parameters.squeeze()) does indeed make this run.

However, I'm not sure if the correct fix involves squeezing the tensor further upstream?

Thanks for any help!

This PR attempts a solution to #23 without (yet) introducing test categories a la @pytest.mark.slow (see https://github.com/mackelab/sbi/blob/86256e02c1080965795e65062c4ab9d3a19015d2/tests/linearGaussian_snpe_test.py#L196)

Hi

As far i can tell, this package is built using the sbi package link The sbi library currently does not seem to support multiple observations, i.e the simulator output should have batch size of 1. So generating time series data shouldn't be possible.

This is enforced in function check_for_possibly_batched_x_shape in user_input_checks

In sbibm package, you have the example code with number of observations as an argument. observation = task.get_observation(num_observation=1) # 10 per task

According to the sbi package, this shouldn't be possible. Did you use some workaround or am i misinterpreting something ?

I've noticed that the prior object from task.get_prior() is not immediately usable with the sbi package since there are no .sample() or .log_prob() methods. Specifically attempting something like this fails:

prior = task.get_prior()
inference = sbi.inference.SNPE(prior=prior, ...)

~~Looking at the code, I imagine this could be implemented by having task.get_prior() return a class instead of a function. Then the class could have a __call__() method to maintain compatibility with the current API. Happy to give this a shot if you guys agree with the change.~~

Edit: it would actually just suffice to expose the prior_dist: https://github.com/sbi-benchmark/sbibm/blob/15f068a08a938383116ffd92b92de50c580810a3/sbibm/tasks/slcp/task.py#L60

The dependence on torch 1.8 makes sense because the current sbi version, which we want to use, depends on it.

With this change we get rid of the helper function get_log_abs_det_jacobian that was distinguishing between the behavior before and after torch 1.8 and was doing the summation explicitly.

Details:

• with the dependence on torch>=1.8 the output of log_prob and log_abs_det_jacobian changes: when the input has several parameter dimensions the output will keep those dimensions and we would have to sum over them by hand to get the joint log_prob over parameter dimensions.
• this can be prevented by "reinterpreting" them as batch dimensions.
• for transforms this works via the IndependentTransform as a wrapper

See also #15

when running the sbibm demo code based on commit 074e06a, I get

import sbibm

task = sbibm.get_task("two_moons")  # See sbibm.get_available_tasks() for all tasks
prior = task.get_prior()
simulator = task.get_simulator()
observation = task.get_observation(num_observation=1)  # 10 per task

# These objects can then be used for custom inference algorithms, e.g.
# we might want to generate simulations by sampling from prior:
thetas = prior(num_samples=10_000)
xs = simulator(thetas)

# Alternatively, we can import existing algorithms, e.g:
from sbibm.algorithms import rej_abc  # See help(rej_abc) for keywords
posterior_samples, _, _ = rej_abc(task=task, num_samples=10_000, num_observation=1, num_simulations=100_000)

I get

task = <sbibm.tasks.two_moons.task.TwoMoons object at 0x7ff456f40f10>, num_samples = 50, num_simulations = 500, num_observation = 1
observation = tensor([[-0.6397,  0.1623]]), num_top_samples = 100, quantile = 0.2, eps = None, distance = 'l2', batch_size = 1000, save_distances = False
kde_bandwidth = 'cv', sass = False, sass_fraction = 0.5, sass_feature_expansion_degree = 3, lra = False

    def run(
        task: Task,
        num_samples: int,
        num_simulations: int,
        num_observation: Optional[int] = None,
        observation: Optional[torch.Tensor] = None,
        num_top_samples: Optional[int] = 100,
        quantile: Optional[float] = None,
        eps: Optional[float] = None,
        distance: str = "l2",
        batch_size: int = 1000,
        save_distances: bool = False,
        kde_bandwidth: Optional[str] = "cv",
        sass: bool = False,
        sass_fraction: float = 0.5,
        sass_feature_expansion_degree: int = 3,
        lra: bool = False,
    ) -> Tuple[torch.Tensor, int, Optional[torch.Tensor]]:
        """Runs REJ-ABC from `sbi`
    
        Choose one of `num_top_samples`, `quantile`, `eps`.
    
        Args:
            task: Task instance
            num_samples: Number of samples to generate from posterior
            num_simulations: Simulation budget
            num_observation: Observation number to load, alternative to `observation`
            observation: Observation, alternative to `num_observation`
            num_top_samples: If given, will use `top=True` with num_top_samples
            quantile: Quantile to use
            eps: Epsilon threshold to use
            distance: Distance to use
            batch_size: Batch size for simulator
            save_distances: If True, stores distances of samples to disk
            kde_bandwidth: If not None, will resample using KDE when necessary, set
                e.g. to "cv" for cross-validated bandwidth selection
            sass: If True, summary statistics are learned as in
                Fearnhead & Prangle 2012.
            sass_fraction: Fraction of simulation budget to use for sass.
            sass_feature_expansion_degree: Degree of polynomial expansion of the summary
                statistics.
            lra: If True, posterior samples are adjusted with
                linear regression as in Beaumont et al. 2002.
        Returns:
            Samples from posterior, number of simulator calls, log probability of true params if computable
        """
        assert not (num_observation is None and observation is None)
        assert not (num_observation is not None and observation is not None)
    
        assert not (num_top_samples is None and quantile is None and eps is None)
    
        log = sbibm.get_logger(__name__)
        log.info(f"Running REJ-ABC")
    
        prior = task.get_prior_dist()
        simulator = task.get_simulator(max_calls=num_simulations)
        if observation is None:
            observation = task.get_observation(num_observation)
    
        if num_top_samples is not None and quantile is None:
            if sass:
                quantile = num_top_samples / (
                    num_simulations - int(sass_fraction * num_simulations)
                )
            else:
                quantile = num_top_samples / num_simulations
    
        inference_method = MCABC(
            simulator=simulator,
            prior=prior,
            simulation_batch_size=batch_size,
            distance=distance,
            show_progress_bars=True,
        )
>       posterior, distances = inference_method(
            x_o=observation,
            num_simulations=num_simulations,
            eps=eps,
            quantile=quantile,
            return_distances=True,
            lra=lra,
            sass=sass,
            sass_expansion_degree=sass_feature_expansion_degree,
            sass_fraction=sass_fraction,
        )
E       TypeError: __call__() got an unexpected keyword argument 'return_distances'

Hello, as I mentioned in my PR #3, there seems to be some UserWarnings raised when KDE is fit with a small number of samples. I put here a small chunk of code which reproduces the warning; that is using my code from #3, so using the ABCpy inference scheme. I have not tried with the other algorithms yet.

I realize there is not much you can do about this as it is due to KDE, but maybe it can be helpful to provide a more explicit warning message saying that the number of samples for KDE are small? Not sure, I realize also this is not super important.

import sbibm

task_name = "two_moons"

task = sbibm.get_task(task_name)  # See sbibm.get_available_tasks() for all tasks
prior = task.get_prior()
simulator = task.get_simulator()
observation = task.get_observation(num_observation=1)  # 10 per task

from sbibm.algorithms.abcpy.rejection_abc import (
    run as rej_abc,
)  
num_simulations = 1000
num_samples = 10000
posterior_samples, _, _ = rej_abc(
    task=task,
    num_samples=num_samples,
    num_observation=1,
    num_simulations=num_simulations,
    num_top_samples=30,
    kde_bandwidth="cv",
)

Hi

The pip install fails currently with the following error.

ERROR: Could not find a version that satisfies the requirement sbibm
ERROR: No matching distribution found for sbibm

I tried it in a conda env and also just python3 virutal env

I know it is not much, but at least it makes the procedure more clear. One could think about adding instructions for conda. But at least these instructions can be performed with a bare python.

I was running the benchmark and found that no method was producing accurate posteriors (according to C2ST) for the gaussian_mixture task. I wondered if the simulator has somehow changed, thereby introducing a different ground truth posterior for each saved observation.

Indeed, this simple check shows that there has been some drift in the simulator

task = sbibm.get_task("gaussian_mixture")
num_observation = 5
true_theta = task.get_true_parameters(num_observation)
sbibm_obs = task.get_observation(num_observation)
new_obs = task.get_simulator()(true_theta)
obss = torch.concat([task.get_simulator()(true_theta) for _ in range(100)])
print(
    (torch.linalg.norm(sbibm_obs - obss)).mean(),
    (torch.linalg.norm(new_obs - obss)).mean(),
)

This typically returns tensor(115.6793) tensor(16.9946).

To fix the issue, either the simulator can be returned to its previous state or we could generate new ground truth parameters and observations; however, this runs the issue of not being backwards compatible with previous versions of sbibm.

refactor sbi run scripts to match the new API of sbi version >=0.20.0.

• [x] depend on newest sbi version 0.20.0 to support passing TransformedDistributions as prior
• [x] run all tests.

Hello, do you have plan to re-run the benchmark for all the 1/2-dim marginals of the tasks, at least for (S)NLE and (S)NPE?

There are some works on 1/2-dim marginal-only sbi, e.g. https://arxiv.org/abs/2107.01214. However, in Fig 1 they are comparing their method trained on marginals vs other methods trained on full distributions, which is not really an apple-to-apple comparison. It'd be useful if you could also provide the baseline for marginal-only sbi. Thanks.

sbi 0.18.0 brought in tons of changes. I was wondering if there are any plans to adopt those? If so, it might be useful to reflect performance changes in the rendered results.

For example, it might be worth considering to make the sbi version an additional field switch, e.g. like the Task currently.

Just wanted to log this here, in case sbibm will make the move to be python 3.10 compatible. Currently, pyabcranger is not compatible with python 3.10, see also https://github.com/diyabc/abcranger/issues/92

Not sure I am overseeing something, but the run methods in the algorithms only return the predicted samples - nothing else.

It might be worthwhile to consider refactoring this, so that each python module in the algorithms directory offers to return the obtained posterior. This would entail in pseudo code:

def train(...):
	return trained_objects

def infer(...)
	return predicted_objects

def run(...):
	trained_objects = train(...)
	predicted_objects = infer(trained_objects, ...)
	return predicted_objects

This refactoring should/would not change the API which is used downstream. It would however allow more analyses on the obtained posterior (mean/median map estimation versus SGD based map estimation etc).