How to calculate Log Predictive Probability Desnity?

Question

I am trying to follow example in the book Statistical Rethinking (Page 210, code 7.13 and 7.14) for calculating LPPD of a simple model.

The given example code in the book is:

set.seed(1)
logprob <- sim( m7.1 , ll=TRUE , n=1e4 ) 
n <- ncol(logprob)
ns <- nrow(logprob)
f <- function( i ) log_sum_exp( logprob[,i] ) - log(ns)
( lppd <- sapply( 1:n , f ) )

I assumed ll=TRUE provides the logprobability, I tried calculating it as (in julia)

 begin
    m=10
    # Take random samples from posterior of parameters
    quap_sample = sample_from_model() 

    # evaluate mean from sampled parameters
    sim_samples = quap_sample.a .+ quap_sample.b * mass_std' 

    # log of samples, I think this part am doing wrong as it is the log of mean I am evaluating and not logprob?
    sim_samples = log.(sim_samples)

    nr,nc = size(sim_samples)
    lppd_sample = zeros(nc)
    for i = 1:nc
        lppd_sample[i] = logsumexp(sim_samples[:,i]) - log(nr)
    end
    lppd_sample
end

Of course results are not matching, Am not totally sure what is ll=TRUE doing (source code accessible here), and previous discussion were also not much helpful

Answer 1

I figured it out. ll=TRUE returns the log probabilities. I took "hints" from PyMC3 version of Statistical rethink book. Which implements it bit more clearly (it can be found here).

Corresponding Julia code would be:

 begin
    m=10000
    
    # draw m samples from posterior.of model q71.
    # model :
    #   lsigma ~ Normal(0,1)
    #   b ~ Normal(0,10)
    #   a ~ Normal(0.5,1)
    #   μ = a .+ b .* mass_std
    #   brain_std ~ MvNormal(μ, exp(lsigma))

    quap_sample = DataFrame(rand(q71.distr,m)',q71.params)

    # calculate the mean from drawn samples of a, and b
    sim_samples = quap_sample.a .+ quap_sample.b * mass_std'

    # calculate log probabilities of the actual data, brain_std 
    # using the sampled mean and lsigma
    ll = logpdf.(Normal.(sim_samples, exp.(quap_sample.lsigma)), brain_std')

    # calculate lppd = \sum_i log 1/s \sum_s p(y_i|theta_s)
    nr,nc = size(sim_samples)
    lppd_sample = zeros(nc)
    for i = 1:nc
        lppd_sample[i] = logsumexp(ll[:,i]) - log(nr)
    end
end
```