Embedded Laplace numerical problem

Developers
1

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

2

Thanks Aki, I was able to replicate on my machine. Also copying @stevebronder