In Statistical Rethinking, 2nd Edition, section 13.1, Richard McElreath says:

Why doesn’t simple quadratic approximation, using for example quap, work with multilevel models? When a prior is itself a function of parameters, there are two levels of uncertainty. This means that the probability of the data, conditional on the parameters, must average over each level. Ordinary quadratic approximation cannot handle the averaging in the likelihood, because in general it’s not possible to derive an analytical solution. That means there is no unified function for calculating the log-posterior. So your computer cannot directly find its minimum (the maximum of the posterior).

Corresponding rethinking code:

m12.2 <- ulam(
        surv ~ dbinom(density, p),
        logit(p) <- a_tank[tank],
        a_tank[tank] ~ dnorm(a, sigma),
        a ~ dnorm(0, 1),
        sigma ~ dcauchy(0, 1)
    data=d, iter=4000, chains=4

Corresponding Stan code:

    int<lower=1> N;
    int<lower=1> N_tank;
    int surv[N];
    int density[N];
    int tank[N];
    vector[N_tank] a_tank;
    real a;
    real<lower=0> sigma;
    vector[N] p;
    sigma ~ cauchy( 0 , 1 );
    a ~ normal( 0 , 1 );
    a_tank ~ normal( a , sigma );
    for ( i in 1:N ) {
        p[i] = a_tank[tank[i]];
    surv ~ binomial_logit( density , p );
generated quantities{
    vector[N] p;
    for ( i in 1:N ) {
        p[i] = a_tank[tank[i]];

I can not figure it out how "prior is itself a function of parameters" could be a problem. It seems that they're just combined into a prior, so there is nothing differ from a single-level model?

Here's my derivation of laplace approximation, maybe I'm wrong in some steps so I can't understand Richard's statement? Please correct me.

From Bayesian Data Analysis, 3rd edition, P84.

$$ \begin{align*} p(\theta|y) \approx N(\hat{\theta}, (I(\hat{\theta}))^{-1}) \\ I(\theta) = -\frac{d^2}{d\theta^2} \log p(\theta|y) \end{align*} $$


$$ I(\theta) = -\frac{d^2}{d\theta^2} \log p(\theta|y) = -\frac{d^2}{d\theta^2} (\log p(\theta,y) - \log p(y)) = -\frac{d^2}{d\theta^2} \log p(\theta,y) $$

As the code just specified joint distribution $\log p(\theta,y)$, like:

function logp_theta_y(){
    target = 0
    target += dcauchy(sigma, 0, 1, log=TRUE);
    target += dnorm(a, 0, 1, log=TRUE);
    target += dnorm(a_tank, a, sigma, log=TRUE);
    for ( i in 1:N ) {
        target += dbinom(surv[i], density[i], logistic(a_tank[tank[i]]), log=TRUE);
    return target;

I can't see any missing material leading Quadratic Approximation not working...

In fact, I run following Julia code (Sorry, I'm not that familar with R) to get Laplace approximation and found that they're close to each other:

using CSV
using DataFrames
using Distributions
using DistributionsAD
using ForwardDiff
using Optim
using LinearAlgebra

df = DataFrame(CSV.File("reedfrogs.csv")) # https://github.com/rmcelreath/rethinking/blob/master/data/reedfrogs.csv

S = df.surv
N = df.density

N_tank = length(S)

logistic(x) = 1/(1+exp(-x))

function logp_y_p(par)
    a_tank = par[1:N_tank]
    a = par[N_tank+1]
    sigma = exp(par[N_tank+2])
    target = 0
    target += logpdf(Cauchy(0, 1), sigma) + par[N_tank+2]
    target += logpdf(Normal(0, 1), a)
    target += logpdf.(Normal.(a, sigma), a_tank) |> sum
    target += logpdf.(Binomial.(N, logistic.(a_tank)), S) |> sum

f(x) = -logp_y_p(x)
x0 = randn(50)

res = optimize(f, x0, LBFGS(); autodiff = :forward)

par_est = Optim.minimizer(res)
FI = -ForwardDiff.hessian(logp_y_p, par_est)
Sigma = inv(FI)

function normal_to_lognormal(mu, sigma)
    exp(mu+sigma^2/2), sqrt((exp(sigma^2)-1) * exp(2*mu+sigma^2))

mu_sigma, sigma_sigma = normal_to_lognormal(par_est[end], sqrt(diag(Sigma)[end]))

mu_vec = copy(par_est)
mu_vec[end] = mu_sigma

sigma_vec = sqrt.(diag(Sigma))
sigma_vec[end] = sigma_sigma

[mu_vec sigma_vec]

enter image description here enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.