Slow Estimation of Log Likelihood

Modeling
1

Hi all,

I am currently trying to implement a grid search on my hyperparameters. I am using loo as the criterion for model fit.
To calculate loo, I am generating log likelihood estimates in the generated quantities, like so:

generated quantities {
vector[N] log_lik;
for (n in 1:N){
log_lik[n] <- bernoulli_logit_lpmf(y[n] | w0[jj[n]] + w1[jj[n]] * s[n] - beta[ii[n]]);
}

variable lengths:
N=4388
w0=20
w1=20
beta=222

chains:2
iterations=1000

Unfortunately, estimation of the log_lik is very slow. Is there a way to accelerate this, for instance by skipping iterations?

2

Hi, what interface do you use?

log_lik shape is 2000, 4388

So it might be faster to do outside Stan.

3

I currently use rstan. So you mean taking the model fit and generating log_lik within R itself could be faster?

4

Yes, because moat of the time is probably used for moving that matrix from Stan (C++) to R.

That has only 8e6 points, so I don’t think eval takes too long.

5

thanks a lot!! So it does run much faster now that I do it within R. My only trouble is that it
gives me slightly different values from what I get from the generated values block. While many values overlap, some deviate. The means are -0.5345751 [generated values] vs. -0.5897133 [R].

Any idea why that would be? Using Using Leave-one-out cross-validation for large data • loo as an example, here’s what I did:

X = list(
ii = as.factor(d$ii), # 222
jj = as.factor(d$jj), # 41
s = as.factor(d$s)) # 2
%>% as.data.frame() %>% model.matrix(~ii+jj+s:jj-1,.,contrasts.arg = lapply(., contrasts, contrasts=FALSE)) # -1: no intercept

standf = list(y=d$y,X=X) %>% as.data.frame

draws = list(
beta = extract(model)$beta,
w0 = extract(model)$w0,
w1 = extract(model)$w1
)

llfun_long_rasch ← function(data_i, draws) {
x_df = data_i[, which(grepl(colnames(data_i), pattern = “X”))]
ii_ind = which(grepl(colnames(x_df), pattern = “X.ii”))
jj_ind = which(grepl(colnames(x_df), pattern = “X.jj[0-9][0-9]?$”))
s1jj_ind = which(grepl(colnames(x_df), pattern = “X.jj[0-9][0-9]?.s1”))
x_i = as.matrix(x_df)

logit_pred ← draws$w0 * t(x_i)[jj_ind,] + draws$w1 * t(x_i)[s1jj_ind,] - draws$beta * t(x_i)[ii_ind,] # I use %*% instead of * in R
dbinom(x = data_i$y, size = 1, prob = 1/(1 + exp(-logit_pred)), log = TRUE)
}

6

There might be some different algorithms used in R, so maybe you need to use some specific functions on R. Not an expert about R so maybe e.g. @bgoodri or @paul.buerkner can help.