Pyro - Likelihood function and sampling dimension

Question

I am learning Pyro and I find the dimensions confusing, despite the rich and detailed documentation

This is the sketch of my model:

DATA_SIZE = 1000

simulated_daily_demand = torch.distributions.Beta(torch.tensor(2.0), torch.tensor(2.0)).sample([DATA_SIZE,])


def model(SIMULATION_DAYS = 30):
    alpha = pyro.param("alpha", pyro.distributions.Uniform(0, 10))
    beta = pyro.param("beta", pyro.distributions.Uniform(0, 10))    
    total_demand = 0
    for t in range(0, SIMULATION_DAYS):
        daily_demand = pyro.sample("daily_demand", pyro.distributions.Beta(alpha, beta), obs=simulated_daily_demand)
        total_demand = total_demand + daily_demand
    return total_demand

model()

I set here priors on the concentrations (alpha, beta). My understanding is that calling pyro.sample with observations fits on the data - I suppose that it maximizes the likelihood of the concentrations given the data.

The output I get:

len(model())
C:\ProgramData\Anaconda3\lib\site-packages\pyro\primitives.py:85: RuntimeWarning: trying to observe a value outside of inference at daily_demand
  RuntimeWarning)
1000

The size()

The values I get seem to look good. The mean of simulated_daily_demand is roughly 0.5 and the mean of model() is ~15 which is ~30*0.5. I don't get the size of the tensor. I would have expected it to be of .size() torch.Size([1]).

I also noticed the warning. I suppose Pyro is complaining because it would like me to write a guide and run some inference on the parameters (eg SVI) before being able to sample from "daily demand". I am also wondering how I could the run the model after the inference of the latent concentrations. A little sketch of the code would really help, thanks!

---

On a hindsight, I think I may have misunderstood the use of plates. Now, if I assume that the observations as independent, I would need to set up a plate. Something like:

import torch
import pyro

NUM_RUNS = 5
DATA_SIZE = 1000

simulated_daily_demand = torch.distributions.Beta(torch.tensor(2.0), torch.tensor(2.0)).sample([DATA_SIZE,])


def model(hist_demand, START_INVENTORY = torch.tensor(100.0), SIMULATION_DAYS = 30):
    alpha = pyro.param("alpha", pyro.distributions.Uniform(0, 10))
    beta = pyro.param("beta", pyro.distributions.Uniform(0, 10))    
    total_demand = 0
    for t in range(0, SIMULATION_DAYS):
        with pyro.plate("obs_loop"):
            daily_demand = pyro.sample("daily_demand", pyro.distributions.Beta(alpha, beta), obs=simulated_daily_demand)
        total_demand = total_demand + daily_demand
    return total_demand

total_demand_runs = []
for r in range(0, NUM_RUNS):
    total_demand_runs.append(model(simulated_daily_demand))

Which returns a nested list of size (NUM_RUNS, SIMULATION_DAYS) that contains a tensor of size DATA_SIZE. The elements (daily_demand) are the same across simulation days. Probably getting closer, but no cigar.

---

import torch
import pyro
import torch.distributions.constraints as constraints

NUM_RUNS = 5
SIMULATION_DAYS = 30
DATA_SIZE = 1000


simulated_daily_demand = torch.distributions.Beta(torch.tensor(2.0), torch.tensor(2.0)).sample([DATA_SIZE, SIMULATION_DAYS])

def model(hist_demand = None, START_INVENTORY = torch.tensor(100.0), SIMULATION_DAYS = SIMULATION_DAYS):
    alpha = pyro.param("alpha", pyro.distributions.Uniform(0, 10))
    beta = pyro.param("beta", pyro.distributions.Uniform(0, 10))
    with pyro.plate("obs_loop"):
        daily_demand_vector = pyro.sample("daily_demand", pyro.distributions.Beta(
            alpha * torch.ones([SIMULATION_DAYS]), 
            beta * torch.ones([SIMULATION_DAYS])), 
            obs=hist_demand
        )
    total_demand = 0
    for t in range(0, SIMULATION_DAYS):
        total_demand = total_demand + daily_demand_vector[t]
    return total_demand

def guide(hist_demand):
    alpha = pyro.param(
        "alpha", 
        pyro.distributions.Normal(torch.tensor(2.0), torch.tensor(0.1)),
        constraint = constraints.positive
        )
    beta = pyro.param(
        "beta",
        pyro.distributions.Normal(torch.tensor(2.0), torch.tensor(0.1)),
        constraint = constraints.positive
        )
    return alpha, beta

from pyro.optim import Adam
adam_params = {"lr": 0.005, "betas": (0.95, 0.999)}
optimizer = Adam(adam_params)
svi = pyro.infer.SVI(model, guide, optimizer, loss=pyro.infer.Trace_ELBO())

n_steps = 5000
# do gradient steps
for step in range(n_steps):
    svi.step(simulated_daily_demand)

alpha_q = pyro.param("alpha").item()
beta_q = pyro.param("beta").item()

Something like this seems to make sense and seems to converge: the SVI spits out approximately right parameter values. Now, the question still holds - how can I run a simulation using the inferred alpha and beta?