1

I'm attempting to implement a model in stan by directly incrementing the log probability using the 'target += bernoulli_logit_lmpf(y|alpha)' functionality. The full model is below:

data {
  int<lower=0> Ntotal;
  int y[Ntotal];
  int<lower=0> Nsubj;
  int<lower=0> s[Ntotal];
  int<lower=0> Nstages;
  int<lower=0> Nstates;
  vector[Nstages*Nstates] log_ts;
  vector[Nstages*Nstates] log_ds;
  int<lower=0> a_i[Ntotal];
  int<lower=0> b_i[Ntotal];
}  

parameters {
  real log_w1;
  real log_w3;
  real delta;
  real delta_p;
  real tau_p;
}

transformed parameters {
  vector[Nstates*Nstages] gradient;
  gradient =  exp(log_ds*delta_p + log_ts*tau_p + log_w3) + exp(log_ds*delta + 
  log_w1);
}

model {
  log_w1 ~ normal(0,5);
  log_w3 ~ normal(0,5);
  delta ~ normal(0,2);
  delta_p ~ normal(0,2);
  tau_p ~ normal(0,2);

 for(n in 1:Ntotal){
    target += bernoulli_logit_lpmf(y[n] | (gradient[a_i[n]] - 
    gradient[b_i[n]]) );
  }
}

My understanding is that this should produce the same results as the sampling statement below:

model {
  log_w1 ~ normal(0,5);
  log_w3 ~ normal(0,5);
  delta ~ normal(0,2);
  delta_p ~ normal(0,2);
  tau_p ~ normal(0,2);

  y ~ bernoulli_logit(gradient[a_i] - gradient[b_i]);
}

However, the two models produce quite different results (even with the same seed in R and Stan). I thought the problem could have been with how I implemented the direct increment, so I tried the following model as well.

model {
  real lps[Ntotal];
  log_w1 ~ normal(0,5);
  log_w3 ~ normal(0,5);
  delta ~ normal(0,2);
  delta_p ~ normal(0,2);
  tau_p ~ normal(0,2);

  for(n in 1:Ntotal){
    lps[n] = bernoulli_logit_lpmf(y[n] | (gradient[a_i[n]] - gradient[b_i[n]]));
  }
  target += log_sum_exp(lps);
}

However, all three models produce different results. Interestingly, the direct increment runs about twice as fast as the sampling statement, but appears to not be giving the correct results. I need to use the direct increment method because I will eventually transform this model into a mixture model, and as far as I know, mixture models need the log probability to be updated directly. Any help with this would be much appreciated!

Regards, Tim

Tim
  • 69
  • 6
  • 1
    After playing around with this a bit more, I've found that vectorising the direct increment with: `target += bernoulli_logit_lpmf(y|gradient[a_i] - gradient[b_i]);` produced the right behaviour. – Tim Apr 16 '17 at 04:43
  • 1
    Implementing the loop as `for(n in 1:Ntotal){ lps[n] = bernoulli_logit_lpmf(y[n] | (gradient[a_i[n]] - gradient[b_i[n]]) ); } target += sum(lps);` also works. But I'm still confused about why the `log_sum_exp()` isn't what I need in this case. – Tim Apr 16 '17 at 04:56
  • Can you provide a MWE that reproduces the error when the `target +=` is inside the loop versus outside? Submitting that to stan-dev might help figure out what the underlying problem is... – BonStats Apr 16 '17 at 05:01
  • `log_sum_exp` is not what you need in this case because that would exponentiate the log-likelihood to the likelihood, _sum_ the likelihood, and take the logarithm of the summed likelihood. But, conceptually the likelihood should be _multiplied_ (not summed) over conditionally independent observations. So, you would need to use the nonexistent `log_prod_exp` function, which does not exist because it would be numerically unstable. – Ben Goodrich Apr 16 '17 at 16:56

0 Answers0