Soooooooooooooooo
@avehtari and I have been trying to get a Laplace approximation working in Stan. For simple models, all of the components are in place. So it should work. But the algebraic solver is breaking. I’m out of ideas, so anyone who can chime in is very welcome.
The basic idea is this. Assume the following model, which is Poisson data where the log-mean is modelled with an iid gaussian random effect.
The Stan code for the basic model is as follows (Here index is the imaginatively named index to the group each object is in [from 1:M]):
data {
int N;
int M;
int y[N];
int index[N];
int<lower=1, upper=M> index[N];
}
parameters {
vector[M] group_mean;
real<lower=0> sigma;
}
model {
group_mean ~ normal(0,sigma);
sigma ~ normal(0,1);
for (i in 1:N) {
y[i] ~ poisson_log(group_mean[index[i]]);
}
}
This works fine. Everything you would expect happens happens.
Now, let’s try a Laplace approximation. This essentailly replaces the poisson log-likelihood with a normal centred at the maximum of p(group_mean | sigma,y) with a variance given by the inverse of the hessian at the maximum.
The following code implements this. The two derivatives are analytical
functions {
vector conditional_grad(vector x, vector sigma, real[] number_of_samples, int[] sums) {
vector[dims(x)[1]] result;
result = sigma[1]^2*(to_vector(sums)-to_vector(number_of_samples).*exp(x)) - x;
return result;
}
vector conditional_neg_hessian(vector x, real sigma) {
vector[dims(x)[1]] result;
result = exp(x) + 1/sigma^2;
return result;
}
}
data {
int N;
int M;
int y[N];
int<lower=1, upper=M> index[N];
}
transformed data {
vector[M] xzero = rep_vector(0.0, M);
real number_of_samples[M];
int sums[M];
for (j in 1:M) {
sums[j] = 0;
number_of_samples[j]=0.0;
}
for (i in 1:N) {
sums[index[i]] += y[i];
number_of_samples[index[i]] +=1.0;
}
}
parameters {
//vector[M] group_mean;
real<lower=0> sigma;
}
model {
vector[M] conditional_mode;
vector[M] laplace_precisions;
vector[1] sigma_tmp;
sigma_tmp[1] = sigma;
//group_mean ~ normal(0, sigma);
sigma ~ normal(0,1);
conditional_mode = algebra_solver(conditional_grad, xzero, sigma_tmp, number_of_samples, sums );
laplace_precisions = conditional_neg_hessian(conditional_mode, sigma);
for (i in 1:N) {
target += -0.5*laplace_precisions[index[i]]*(y[i]- conditional_mode[index[i]])^2;
}
}
This is a disaster. Now I’m surprised because a poisson-log-normal model is about as nice as you can get (the conditional distribution is log-concave, so the root finding should not have any problems). Instead of this working, I get this error constantly and the sampler doesn’t move.
Iteration: 2000 / 2000 [100%] (Sampling)
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: normal_lpdf: Random variable is nan, but must not be nan! (in ‘/Users/ds921/Desktop/laplace.stan’ at line 50)If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Any suggestions?
Files here:
data.R (921 Bytes)
<a class=“attachment” href=“/uploamake_data.R (131 Bytes)
ds/mc_stan/original/2X/e/e5c7815ea7b1c9e1110ee28890c3df080cdf5347.stan”>laplace.stan (1.2 KB)
mmc.stan (253 Bytes)