Implementing Logit-Normal priors on Dirichlet-Multinomial model

Modeling
1

Hello,

I am trying to implement the model discussed in this article. I outlined it below:

Y_{n \times d} \sim Multinomial(N, p) \\ p_i \sim Dirichlet(\theta_i) \\ log(\theta_i) = \alpha_i + X_{n \times m}\beta_{i} \\ \alpha_i \sim N(0, 10^2); \\ \beta_{i, j} = \tau_j \tilde{\lambda}_{i, j} z_{i, j};\\ \tilde{\lambda}_{i, j} = \frac{exp\{\lambda_{i, j}\}}{1 + exp\{\lambda_{i, j}\}};\\ \lambda_{i, j} \sim N(0, 10^2);\\ \tau_j \sim Half-Cauchy(0, 1);\\ z_{i, j} \sim N(0, 1).\\

for i = 1, \ldots, d and j = 1, \ldots, m.

STAN implementation:

data {
 int < lower =0> n; // number of observations
 int < lower =0> d; // number of responses
 int < lower =0> m; // number of predictors
 matrix[n, m] x[d]; // outputs
 int y[n, d]; // inputs
 real < lower =1> nu_global ; // degrees of freedom for the half -t priors for tau
 real < lower =0> scale_global ; // scale for the half -t prior for tau
 real nu_local; 
 real < lower =0> scale_local ; 
 real < lower =0> scale_intercept ; 

}

parameters {
 row_vector < lower =0>[m] tau; // global shrinkage parameter
 row_vector[d] beta0; 
 row_vector[m] z[d];
 vector<lower = 0, upper = 1>[n] p[d];
 row_vector[m] lambda[d];
}

transformed parameters {
 vector[m] beta[d] ; // regression coefficients
 row_vector[m] lambda_tilde[d] ;
 vector[n]  theta[d]; // latent function values

 for (i in 1:d){
   lambda_tilde[i] = inv_logit(lambda[i]);
   beta[i] =  to_vector(tau .* lambda_tilde[i] .* z[i]);
   theta[i] = exp(x[i]*beta[i]) ;
 }
}

model {
 beta0 ~ normal(0, scale_intercept);
 tau ~ student_t(nu_global, 0, scale_global);
 
 for (i in 1:d){
   z[i] ~ normal(0, 1);
   lambda[i] ~ normal(nu_local, scale_local);
 }
 
 for (j in 1:n){
   p[j] ~ dirichlet(theta[j]);
   y[j] ~ multinomial(p[j]);
 }

In R with a reproducible example:

library(rstan)

n <- 30
d <- 10
x <- sample(c(0, 1), size = n, replace = TRUE)
y <- t(rmultinom(n = n, size = 10000, prob = rep(1/d, d)))

stan_data <- list(
  n = n,
  d = d,
  m = 1,
  y = y,
  x = rep(list(matrix(x)), d),
  scale_intercept = 10,
  scale_global = 1,
  nu_global = 1,
  scale_local = 10,
  nu_local = 0
  
)


fit <- stan(
  file = "multinomial_dirichlet.stan",  # Stan program
  data = stan_data,    # named list of data
  chains = 1,             # number of Markov chains
  warmup = 1000,       # number of warmup iterations per chain
  iter = 2000,            # total number of iterations per chain
  cores = 1,              # number of cores (could use one per chain)
  refresh = 0,             # no progress shown
  thin = 1         
)
# [1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."   
# [2] "In addition: Warning message:"                                            
# [3] "In system(paste(CXX, ARGS), ignore.stdout = TRUE, ignore.stderr = TRUE) :"
# [4] "  'C:/RBUILD~1/4.0/usr/mingw_/bin/g++' not found"                         
# [1] "error occurred during calling the sampler; sampling not done"

multinomial_dirichlet.stan (1.3 KB)

The error message is not clear for me. Any thoughts?

1 Like
2

Something in RStan is swallowing the error. Upgrading to the latest RStan might help, I dunno.
In any case CmdStan gives

  Error evaluating the log probability at the initial value.
Exception: dirichlet_lpdf: probabilities is not a valid simplex.

It’s saying that the parameter p should be declared as a simplex and not vector<lower=0,upper=1>.

parameters {
  ...
  simplex[n] p[d];
}

Next up

Exception: multinomial_lpmf: Size of number of trials variable (10) and rows of probabilities parameter (30) must match in size

Some of the dimensions are transposed, you’re mixing up n and d. Note in particular that

simplex[n] p[d];

is an array of d vectors, each of which has n elements.

3 Likes