 # Is there any way to to calculate lppd (pymc3)?

With reference to chapter 7 of the book Statistical Rethinking Second Edition by Richard Mcelreath; I am getting a different lppd score using posteriors from Pymc3 and R.

For convenience, I have Included all the codes from the book>>

`R code`
`7.1`

sppnames <-c(“afarensis”,“africanus”,“habilis”,“boisei”,
“rudolfensis”,“ergaster”,“sapiens”)
brainvolcc <-c(438,452,612,521,752,871,1350)
masskg <-c(37.0,35.5,34.5,41.5,55.5,61.0,53.5)
d <-data.frame(species=sppnames,brain=brainvolcc,mass=masskg)

Book uses R code:

`R Code`
`7.3`

m7.1 <-quap(
alist(
brain_std ~dnorm(mu,exp(log_sigma)),
mu <-a+b*mass_std,
a ~dnorm(0.5,1),
b ~dnorm(0,10),
log_sigma ~dnorm(0,1)
), data=d)

`R code`
`7.13`

set.seed(1)
lppd( m7.1,n=1e4)

We get average lppd score as below for each data point
` 0.6098668 0.6483438 0.5496093 0.6234934 0.4648143 0.4347605 -0.8444633`

And; Below this is my attempt to calculate lppd with pymc3

with pm.Model() as m_1:
b_std=pm.Data(‘b_std’,brains[‘brain_std’])
m_std=pm.Data(‘m_std’,brains[‘mass_std’])
a=pm.Normal(‘a’,0.5,1)
bM=pm.Normal(‘bM’,0,10)
mu=pm.Deterministic(‘mu’,a+bM*m_std)
log_sigma=pm.Normal(‘log_sigma’,0,1)
brain=pm.Normal(‘brain’,mu=mu,sigma=pm.math.exp(log_sigma),observed=b_std)
trace_1=pm.sample(chains=4)
posterior_1=pm.sample_posterior_predictive(trace_1,var_names=[‘mu’,‘log_sigma’])

With reference to code (7.14) which has only OLS lppd ; I wrote a similar which calculates lppd from the posteriors.

def my_lppd(pred_mu,sigma):
rows, cols = pred_mu.shape
ll=np.zeros((rows,cols))
for i in range(cols):
log_prob=stats.norm.logpdf(brains.brain_std,loc=pred_mu[:,i],scale=sigma[i])
ll[:,i]=log_prob
lppd = np.zeros(rows)
for i in range(rows):
lppd[i] = logsumexp(ll[i])-np.log(cols)
return lppd

But then I get different lppds, but the ordering is same

my_lppd(pred_mu=posterior_1[‘mu’].T,sigma=np.exp(posterior_1[‘log_sigma’]))

`array([ 0.40075526, 0.42605526, 0.3376496 , 0.41475887, 0.28507802, 0.24066785, -0.61324239])`
Am I missing something here? Also, the lppd score gets worse for higher-order regressions.
Thanks in advance!