We discussed this in Stan gathering, and @charlesm93 suggested sharing here.
1D normal prior and Poisson model. Integrating out the latent theta works with integrate_1d()
, but doesn’t work with laplace_marginal()
. The likely reason is a loss of numerical accuracy in some part.
The Stan code including both integrate_1d and laplace_marginal. Note that in the generated quantities if laplace_marginal fails, then the output of integrate_1d is also assigned as NaN. So with the difficult example need to run integrate_1d only.
// Poisson with varying intercept ("random effect")
functions {
// for laplace_marginal
real poisson_re_log_ll(vector theta, data array[] int y, vector mu) {
return poisson_log_lpmf(y | mu + theta);
}
matrix cov_fun(real sigma, data int N) {
return diag_matrix(rep_vector(sigma^2, N));
}
// for integrate_1d
real integrand(real theta, real notused, array[] real phi,
array[] real X_i, array[] int y_i) {
real sigma = phi[1];
real mu = phi[2];
real p = exp(normal_lpdf(theta | 0, sigma) + poisson_log_lpmf(y_i | theta + mu));
return p;
}
}
data {
int N;
array[N] int y;
vector[N] mu;
real sigmaz;
int run_laplace;
}
transformed data {
vector[N] theta_0 = rep_vector(0.0, N); // initial guess
real integrate_1d_reltol = 1e-6;
// real tolerance = 1e-2;
// int max_num_steps = 1000, hessian_block_size = 1, solver = 1, max_steps_line_search = 100;
}
generated quantities {
// laplace_marginal
real ll_laplace = 0;
if (run_laplace==1) {
ll_laplace = laplace_marginal(poisson_re_log_ll, (y, mu), theta_0, cov_fun, (sigmaz, N));
}
// integrate_1d
real ll_integrate_1d = 0;
for (i in 1:N) {
ll_integrate_1d += log(integrate_1d(integrand,
negative_infinity(),
positive_infinity(),
append_array({sigmaz}, {mu[i]}),
{0}, // not used, but an empty array not allowed
{y[i]},
integrate_1d_reltol));
}
}
R code
library(cmdstanr)
modt <- cmdstan_model(stan_file = "test.stan", force_recompile=TRUE)
datat <- list(N=1, y=153, mu=4.3, sigmaz=2.0, run_laplace=1)
fitt <- modt$sample(data = datat, refresh = 0, chains = 1, parallel_chains = 1, fixed_param=TRUE, iter_sampling=1)
fitt3$draws(format="df")
with output
# A draws_df: 1 iterations, 1 chains, and 2 variables
ll_laplace ll_integrate_1d
1 -6.7 -6.7
With another data
datat <- list(N=1, y=183, mu=0.5, sigmaz=2.0, run_laplace=1)
fitt3 <- modt3$sample(data = datat, refresh = 0, chains = 1, parallel_chains = 1, fixed_param=TRUE, iter_sampling=1)
There are exceptions
Running MCMC with 1 chain...
Chain 1 Exception: laplace_marginal_density: max number of iterations: 100 exceeded. (in '/tmp/RtmpCmXqxG/model-14e1311f91c94b.stan', line 36, column 4 to column 94)
Chain 1 Exception: laplace_marginal_density: max number of iterations: 100 exceeded. (in '/tmp/RtmpCmXqxG/model-14e1311f91c94b.stan', line 36, column 4 to column 94)
Chain 1 finished in 0.0 seconds.
Sometimes instead of max number of iterations, for example, with
datat <- list(N=1, y=183, mu=0.5, sigmaz=2.5, run_laplace=0)
the exception is
Chain 1 Exception: Error in laplace_marginal_density: theta contains NaN values for all values. (in '/tmp/RtmpCmXqxG/model-14e1311f91c94b.stan', line 36, column 4 to column 94)
and sometimes -nan or -inf. In all these cases, the output is NaN.
# A draws_df: 1 iterations, 1 chains, and 2 variables
ll_laplace ll_integrate_1d
1 NaN NaN
Turning off laplace_marginal we can see that integrate_1d still works
datat <- list(N=1, y=183, mu=0.5, sigmaz=2.0, run_laplace=0)
fitt3 <- modt3$sample(data = datat, refresh = 0, chains = 1, parallel_chains = 1, fixed_param=TRUE, iter_sampling=1)
fitt3$draws(format="df")
with output
# A draws_df: 1 iterations, 1 chains, and 2 variables
ll_laplace ll_integrate_1d
1 0 -9.6
As integrate_1d works, this should not be too bad for laplace_marginal, or if it is, it would be useful to at least get more informative messages