Hi Everyone,
I’m trying to use the Importance function and psis_diagnostic function to obtain the \hat{k}. Importance and psis_diagnostic.
When I run the psis_diagnostic function I’m getting the following error with dimension
ValueError: Model and guide event_dims disagree at site ‘L’: 9 vs 2
When I used Importance function to obtain the effective sample size I’m getting the warning
‘The log_weights list is empty, effective sample size is zero’
I couldn’t figure out how to solve these two issues. I would be greatly appreciated if someone could help me to figure out.
Here is the code that I used.
import os
import timeit
import torch
import pyro
from pyro.optim import Adam
from pyro.infer import SVI, Trace_ELBO, TraceGraph_ELBO, importance, Importance
from pyro import plate
import numpy as np
import pyro.distributions as dist
import torch.distributions.constraints as constraints
import scipy
from scipy.stats import beta, invgamma
from scipy.spatial import distance
import matplotlib.pyplot as plt
N = 50
M = 10
K0 = 3
K = 6
eps = 0.001
alpha = 1
a = torch.tensor([0.5,20,2])
mu = torch.tensor([100,500,1000])
seed = 150
n_steps = 100000
def simdata(N,M,K0,alpha,a,mu,seed):
pyro.set_rng_seed(seed)
alpha1 = alpha * torch.ones(K0,M)
sp1 = a * torch.ones(N,K0)
rt1 = (a/mu) * torch.ones(N,K0)
R = dist.Dirichlet(alpha1).sample()
L = dist.Gamma(sp1,rt1).sample()
LR = torch.matmul(L,R)
X = dist.Poisson(LR).sample()
return L, R, X
L1,R1,X = simdata(N,M,K0,alpha,a,mu,seed)
def model(X,K):
N = len(X)
M = len(X[0])
#defining hyperparameters
alpha0 = 1
#sampling L and R
a = pyro.sample("a",dist.InverseGamma(2*torch.ones(K),4*torch.ones(K)).to_event())
mu = pyro.sample("mu",dist.InverseGamma(0.35*torch.ones(K),0.001*torch.ones(K)).to_event())
R = pyro.sample("R",dist.Dirichlet(alpha0 * torch.ones(K,M)).to_event())
L = pyro.sample("L",dist.Gamma(a * torch.ones(N,K),(a/mu)* torch.ones(N,K)).to_event())
LR = torch.matmul(L,R)
#pyro.sample("X", dist.Poisson(LR).to_event(), obs = X)
with plate("X_1", size = len(X[0])):
with plate("X_2",size = len(X)):
pyro.sample("X", dist.Poisson(LR), obs = X)
####################################################################
def guide(X,K):
N = len(X)
M = len(X[0])
c = 10
R1 = dist.Dirichlet(1 * torch.ones(K,M)).sample()
L1 = dist.Gamma(100* torch.ones(N,K),0.01 * torch.ones(N,K)).sample()
a1 = dist.InverseGamma(2*torch.ones(K),4*torch.ones(K)).sample()
mu1 = dist.InverseGamma(0.35*torch.ones(K),0.001*torch.ones(K)).sample()
alpha0_q = pyro.param("alpha0_q", c*R1, constraint = constraints.positive)
a_q = pyro.param("a_q", c*L1, constraint = constraints.positive)
b_q = pyro.param("b_q", (1/c)*L1, constraint = constraints.positive)
a0_q = pyro.param("a0_q", 1000*mu1,constraint = constraints.positive)
b0_q = pyro.param("b0_q", 1000*mu1,constraint = constraints.positive)
c0_q = pyro.param("c0_q", c*a1,constraint = constraints.positive)
d0_q = pyro.param("d0_q", c*a1,constraint = constraints.positive)
pyro.sample("a", dist.InverseGamma(c0_q,d0_q).to_event())
pyro.sample("mu", dist.InverseGamma(a0_q,b0_q).to_event())
R = pyro.sample("R", dist.Dirichlet(alpha0_q).to_event())
L = pyro.sample("L",dist.Gamma(a_q,b_q).to_event())
torch.matmul(L,R)
#############################
#specifying the optimizing parameters
#adam_params = {"lr": 0.01, "betas": (0.90, 0.999)}
#optimizer = Adam(adam_params)
#svi = SVI(model, guide, optimizer, loss=Trace_ELBO(num_particles=10))
IS = importance(model= model,guide = guide,num_samples = 1000)
w = IS.get_normalized_weights()
ess = IS.get_ESS()
k_hat = importance.psis_diagnostic(model,guide,X,K,num_particles = 3000)
Thank you