Replicating bernoulli_logit_lpmf in R

Modeling
1

I am just moving some of my ex-post analyses from the generated quantities block in Stan to R.

In Stan, I use the following code:

data{
  int<lower = 1> N; 
  int<lower = 0, upper = 1> y[N];
  int<lower = 0> K;
  matrix[N, K] X;
...
parameters {
  vector[K] alpha
...
}

generated quantities {
vector[N] y_pred;
for (t in 1:N) {
y_pred[t] = bernoulli_logit_lpmf(y[t] | X[t]*alpha)
  }
}

In R, I translated this into the following formula:

ext_mod = rstan::extract(out)
N_draws = nrow(ext_mod$alpha)
y_pred = matrix(nrow = N_draws, ncol = N)

for(d in 1:N_draws){
  for (t in 1:N) {
dbinom(y[t], size = 1,
                   prob = exp(X[t, ]%*%ext_mod$alpha[d, ]) /
                              (1+exp(X[t, ]%*%ext_mod$alpha[d, ])))
  }
}

I was wondering if what I have specified in R using dbinom() corresponds to the log Bernoulli probability mass of y (bernoulli_logit_lpmf) I have specified in Stan. It is my best guess based on this page.

Background of the question is that I notice a discrepancy between the analyses performed in Stan in generated quantities and the analyses in R (…even though I exported the posterior output and then went through all the draws in R as well).

Thanks in advance.

2

The _lpmf functions return the log of the PMF for the distribution, so you would need to call dbinom with the log=TRUE option

1 Like
3

Thank you, that seems to work.