Thanks for your reply!
I’ve tried to play around with an additional plate, nested plates and the to_event function without any luck. I will continue to try different combinations, but the ones that make sense in my head all fail. (When I use the to_event(1) you suggested, I get a new size error that I can’t get around.) I’ve added a new version of the code below since, the example above might miss some details that I’m not aware of.
Here is a simplified version of what I want to do:
- Each data point contains 64 pixels
- The intensity of each pixel is an integer, modeled with a Poisson variable (the pixels are independent, so they potentially have 64 different intensities)
- The prior of the intensities have a gamma distribution.
- I have a set of (178) data points, each data point is a vector of size 64
I can’t seem to get the multi-dimensional Poisson distribution to do what I want to do.
import pyro
import pyro.distributions as dist
import pyro.distributions.constraints as constraints
import torch
import logging
from sklearn import datasets
digits = datasets.load_digits()
X_train = digits.images[digits.target == 0]
X_train = torch.tensor(X_train)
print(X_train.shape)
nr_pixels = 64
X_train = X_train.reshape((-1, nr_pixels)) #
nr_samples, nr_pixels = X_train.shape
print(nr_samples, nr_pixels)
def model(data):
a = 16.0 * torch.ones(nr_pixels)
b = 1.0
theta = pyro.sample("theta", dist.Gamma(a, b))
assert len(theta) == nr_pixels
with pyro.plate("data_loop", nr_samples):
pyro.sample("obs", dist.Poisson(theta).to_event(1), obs=data)
def guide(data):
aq = pyro.param("aq", torch.ones(nr_pixels), constraint=constraints.positive)
bq = pyro.param("bq", torch.tensor([1]), constraint=constraints.positive)
pyro.sample("theta", dist.Gamma(aq, bq))
adam = pyro.optim.Adam({"lr": 0.01})
elbo = pyro.infer.Trace_ELBO()
svi = pyro.infer.SVI(model, guide, adam, elbo)
pyro.render_model(model, model_args=(X_train,), render_distributions=True)
losses = []
for step in range(100):
loss = svi.step(X_train)
losses.append(loss)
logging.info(f"ELBO loss: {loss}")
edit:
some of the configurations I’ve tried:
# Shape mismatch inside plate('data_loop') at site obs dim -1, 178 vs 64
# with pyro.plate("data_loop", nr_samples):
# pyro.sample("obs", dist.Poisson(theta), obs=data)
# Expected [], actual [64]
# with pyro.plate("data_loop", nr_samples):
# pyro.sample("obs", dist.Poisson(theta).to_event(1), obs=data)
# Expected [], actual [64]
# for s in range(nr_samples):
# for p in range(nr_pixels):
# pyro.sample(f"obs_{s},{p}", dist.Poisson(theta[p]), obs=data[s])
# Expected [], actual [64]
# with pyro.plate("data_loop", nr_samples):
# with pyro.plate("inner_loop", nr_pixels):
# pyro.sample("obs", dist.Poisson(theta).to_event(1), obs=data)
# Shape mismatch inside plate('data_loop') at site obs dim -1, 178 vs 64
with pyro.plate("data_loop", nr_samples):
with pyro.plate("inner_loop", nr_pixels):
pyro.sample("obs", dist.Poisson(theta), obs=data)
Edit 2:
I’ve tried with a MultivariateNormal as well, and I can’t get that to work either, so something with my approach is definitely wrong, I just can’t figure out what. Sometimes I can render the model, but I’ve never been able to run a step.
Using the language here this is what I want to do:
- Draw rates \theta from 64 different Gamma distributions (batch_shape = 64)
- Use each rate to configure a Poisson distribution (64 different distributions)
- Observe the 178 samples from the data (sample_shape = 178)