Reparameterization of a multinomial distribution

Modeling
1

I am trying to do a reliability simulation project, and I kept getting this ‘Chain 1: Log probability evaluates to log(0), i.e. negative infinity. Chain 1: Stan can’t start sampling from this initial value.’ error. I guess I need some reparameterization but have no clue where to begin.

Here is how I simulated the data.

set.seed(666)

alpha = matrix(c(0.001, 0.05, 0.0001, 0.08), ncol = 2, byrow = TRUE)
w = c(45, 55)
t = matrix(seq(10, 40, 10), ncol = 2, byrow = TRUE)
K = c(10, 50, 100)


lambda = function(r = 1:2, j = 1:2) {
  alpha_r0 = alpha[r, 1]
  alpha_r1 = alpha[r, 2]
  return(alpha_r0*exp(alpha_r1*w[j]))
}

p = function(i = 1:2, j = 1:2){
  lam1j = lambda(1, j)
  lam2j = lambda(2, j)
  
  p0 = exp(-(lam1j + lam2j)*t[i, j])
  p1 = (lam1j/(lam1j + lam2j))*(1 - p0)
  p2 = (lam2j/(lam1j + lam2j))*(1 - p0)
  
  return(c(p0, p1, p2))
}

Tmatrxi = t(rmultinom(10, K[1], prob = p(1, 1)))

This simulated data is then estimated using a multinomial distribution. I guess the only thing that makes this a little bit trickier is the \theta (p in Stan code) in the multinomial distribution are dependent on other parameters (a00 and a11 in this case).

data {
  int<lower=0> n;
  int DAT[n, 3];
  vector[2] W;
  vector[4] Tim;
}
parameters{
  real<lower=0> a00; 
  real<lower=0> a01;
}
transformed parameters{
  real<lower=0, upper=1> p0; 
  real<lower=0, upper=1> p1; 
  real<lower=0, upper=1> p2;
  vector[3] p;
  
  p0 = exp(-(a00*exp(a01*W[1]) + a00*exp(a01*W[2]))*Tim[1]);
  p1 = (a00*exp(a01*W[1])/(a00*exp(a01*W[1]) + a00*exp(a01*W[2])))*(1 - p0);
  p2 = 1-p0-p1;
  
  p[1] = p0; p[2] = p1; p[3] = p2;
}
model{
  for (i in 1:n)
    target += multinomial_lpmf(DAT[i,] | p);
    //DAT[i,] ~ multinomial(p);

  a00 ~ gamma(1, 10);
  a01 ~ gamma(1, 10);
}

Then I run this model with rstan using the following code.

library(rstan)
rstan_options(auto_write = TRUE)

stan_dat = list(
  n = nrow(Tmatrxi),
  DAT = Tmatrxi,
  W = w,
  Tim = seq(10, 40, 10)
)

fit <- stan(
  model_code = reliab, data = stan_dat, init_r = 0.0001, 
  warmup = 500, iter = 1000, chains = 1, cores = 1)

However, Stan keeps giving me the error Chain 1: Log probability evaluates to log(0), i.e. negative infinity. Chain 1: Stan can't start sampling from this initial value.

Specifying initial values does not work here as I have tried different init_r. Even if I narrow that down to 0.0001, Stan can not sample from this distribution. I guess the other approach could be reparameterizing the model. My prior experience of reparameterization is limited non-centered parameterization, which may not be applied here.

Is there another way to get around this error? Thanks!

2

Before looking for reparametrizations, I suggest you try to debug what is actually happening. The error message tells you that your model does not have support for the initial values. Maybe the data are problematic?

Some reasons the model might not have support:

  • The model and data are admissible, you may try putting additional constraints e.g.:
    • int<lower=0> DAT[n, 3]
    • simplex[3] p; (enforces that p sums to 1)
  • Check that elements of p are not 0 (for example by printing them)
3

If it’s helpful, a colleague and I present code for a multinomial (generalized bernoulli, aka “categorical”) model here – check the supplemental files:
https://link.springer.com/article/10.1007/s00265-017-2363-8