Samples of a Bernoulli random variable have more than one element during MCMC

I’m new to Pyro and looking for a package that can help me simulate (and possibly do other things like inference with) stochastic processes that have deterministic or discrete elements. Pyro seems to be a good solution and I’m just trying to get a simple ‘MVP’ demonstration to work.

I initially posted a question on stackoverflow but got no responses.

Should I be posting these kind of beginner questions here or on stackoverflow?

I notice that the tag ‘Pyro’ is already being used on stack overflow for some other package.

Please advise. If this forum is the right place I will re-post the question here, otherwise, please go to the link above and help answer my question!

thanks.

Hi Bill, to my knowledge, none of Pyro devs is active on stackoverflow. We use this forum to discuss about ideas and questions related to Pyro. If you have questions, please ask here. Welcome! :slight_smile:

1 Like

Thanks @fehiepsi. I’ll post my question here then.

I’m trying to get my first stochastic process model working. I adapted the code from here to suit my example problem which is simply two Gaussians with a discrete probability of the sample coming from one or the other. Right now, I’m just trying to do something simple like estimate the probability distribution of the output of this process.

Here’s my code:

import torch
import pyro
import pyro.distributions as dist
from pyro.infer.mcmc import HMC, MCMC

# Actual data sample
observations = torch.tensor(
    [0.00528813, -0.00589001, -1.20608593, 0.00190794,
     0.89052784,  0.66690464,  0.57295968, 0.02605967]
)

# Define the process
def model(observations):
    
    a_prior = dist.Beta(2, 2)
    a = pyro.sample("a", a_prior)
    c = pyro.sample('c', dist.Bernoulli(a))
    if c.item() == 1.0:
        my_dist = dist.Normal(0.785, 1.0)
    else:
        my_dist = dist.Normal(0.0, 0.01)
    
    for i, observation in enumerate(observations):
        measurement = pyro.sample(f'obs_{i}', my_dist, obs=observation)

# Clear parameters
pyro.clear_param_store()

# Define the MCMC kernel function
my_kernel = HMC(model)

# Define the MCMC algorithm
my_mcmc = MCMC(my_kernel,
               num_samples=5000,
               warmup_steps=50)

# Run the algorithm, passing the observations 
my_mcmc.run(observations)

The exception raised is:

<ipython-input-2-a668622a0fb9> in model(observations)
     11     a = pyro.sample("a", a_prior)
     12     c = pyro.sample('c', dist.Bernoulli(a))
---> 13     if c.item() == 1.0:
     14         my_dist = dist.Normal(0.785, 1.0)
     15     else:

ValueError: only one element tensors can be converted to Python scalars
Trace Shapes:    
 Param Sites:    
Sample Sites:    
       a dist   |
        value   |
       c dist   |
        value 2 |

I had a look at c using the debugger and for some reason it has two elements the second time model() is called:

tensor([0., 1.])

What is causing this? I wanted it to be a simple scalar having the values 0 or 1.

As a further test, the condition statement works fine when taking samples in the normal way:

# Conditional switch test
a_prior = dist.Beta(2, 2)
a = pyro.sample("a", a_prior)
for i in range(5):
    c = pyro.sample('c', dist.Bernoulli(a))
    if c.item() == 1.0:
        print(1, end=' ')
    else:
        print(0, end=' ')

# 0 0 1 0 0

Hi @billtubbs, HMC only works for continuous variables, so we have to integrate out the discrete latent variable. You can write your model as follows

    c = c.long()
    my_dist = dist.Normal(torch.tensor([0.785, 0.0])[c], torch.tensor([1.0, 0.01])[c])
    
    with pyro.plate("N", observations.shape[0]):
        measurement = pyro.sample('obs', my_dist, obs=observations)
1 Like

Thanks. I will have to read that documentation link in full so I can understand how it works but it does work now!