SNPE_C won't draw samples from the posterior

[paarth-dudani]

Describe the bug
I have implemented SNPE_C with a uniform prior and a sine simulator with noise. However, it randomly won't sample/sample very slowly from the posterior, even after the network converges.

To Reproduce
1.
Python version: 3.9.13
SBI version: 0.23.1
2. Minimal code example:

n = 100
t = np.linspace(0,10,n)

num_sims = 100
num_rounds = 2
scale_true = torch.tensor([2.5])

class WrappedUniform(Uniform):
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

    def log_prob(*args, **kwargs):
        log_probs = Uniform.log_prob(*args, **kwargs)
        return log_probs.squeeze()

prior_mod = WrappedUniform(torch.tensor([0.1]),torch.tensor([5.0]))

## Function to introduce noise into the exponential:
def noisy_sine(scale, input, initial_cond):
    x_noise = np.zeros((len(input)),)
    for i in range(len(input)):
        value = 1/(scale)*torch.sin(scale*t[i]) + initial_cond + 0.05*np.random.normal(0,1)
        x_noise[i] = value
    return torch.tensor(x_noise).float()

simulator = lambda theta: ((1/scale_true)*torch.sin(-t/scale_true) + 0.05*torch.rand_like(theta)).float()
x_o = noisy_sine(scale_true,t,1)

inference = NPE(prior=prior_mod, density_estimator='made')
proposal = prior_mod

# Later rounds use the APT loss.
for _ in range(num_rounds):
    theta = proposal.sample((num_sims,))
    print(np.shape(theta))
    x = simulator(theta)
    x[x>10] = 10
    _ = inference.append_simulations(theta, x, proposal=proposal).train(validation_fraction=0.1,use_combined_loss = True, retrain_from_scratch = True)
    posterior = inference.build_posterior().set_default_x(x_o)
    proposal = posterior

posterior_samples_exp = posterior.sample((500,), x=x_o)

_ = pairplot(
    posterior_samples_exp, limits=[[0, 10]], figsize=(5, 5),
    labels=[r"$\theta$"],
    points=theta_true # add ground truth thetas
)

3. There is no error message but the algorithm just won't sample from the posterior in some iterations.

Example 1:

Example 2:

Expected behavior

Sampling from the posterior in a few seconds.

[michaeldeistler]

Can you try using truncated proposals?

from sbi.inference import NPE
from sbi.utils import RestrictedPrior, get_density_thresholder

inference = NPE(prior)
proposal = prior
for _ in range(num_rounds):
    theta = proposal.sample((num_sims,))
    x = simulator(theta)
    _ = inference.append_simulations(theta, x).train(force_first_round_loss=True)
    posterior = inference.build_posterior().set_default_x(x_o)

    accept_reject_fn = get_density_thresholder(posterior, quantile=1e-4)
    proposal = RestrictedPrior(prior, accept_reject_fn, sample_with="rejection")

[paarth-dudani]

I made the changes by using the truncated proposals as you have suggested.

I have a couple of questions.

1. As I added the density estimator: 'mdn'. I am getting the same issue as before:

inference = NPE(prior=prior_mod, density_estimator='mdn')

Similarly for other density estimator specifications such as 'made'.

inference = NPE(prior=prior_mod, density_estimator='made')

Why would this happen?

2. What does the user warning in the first image mean, since I have gotten it on numerous occasions?

[michaeldeistler]

1. If it works for MAF and NSF then just use those, they are generally better than MDN and MADE.

2. I am not sure and would have to have a look, my guess that this has to do with newer pytorch versions.

[manuelgloeckler]

To expand on this a bit:

1. Especially in the 1d case MAF and MADE will only be able to model Gaussian, which probably are a bad proxy for what you want to learn.

The "issue" lies within accept_reject sampling. The posterior for NPE given a prior on a constrained domain will perform rejection sampling for this domain. If the acceptance rate is too low, then the sampling progress will be slow/get stuck (there normally should be a UserWarning for that; I have to check why/if this stopped working).

The issue in your case actually is your code:

simulator = lambda theta: ((1/scale_true)*torch.sin(-t/scale_true) + 0.05*torch.rand_like(theta)).float()

Note that in your simulator, the parameters theta are just used here, torch.rand_like(theta), so do not have any effect on the simulation output, which remains constant from the global variable scale_true.

In the end this leads to the estimator learning "non-sense" because the data you give it is "non-sense". For certain input it then just predicts stuff outside the prior domains, which all gets rejected.

[manuelgloeckler]

Closing this issue, in favor of #1292