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

**URL:** <https://discourse.pymc.io/t/is-there-any-way-to-to-calculate-lppd-pymc3/7519>\
**Category:** Questions\
**Created:** [May 31, 2021, 6:59am UTC](https://discourse.pymc.io/t/is-there-any-way-to-to-calculate-lppd-pymc3/7519 "2021-05-31T06:59:27Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![roesta07](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/roesta07/32/6060_2.png) [@roesta07](https://discourse.pymc.io/u/roesta07)\
**Post date:** [May 31, 2021, 6:59am UTC](https://discourse.pymc.io/t/is-there-any-way-to-to-calculate-lppd-pymc3/7519/1 "2021-05-31T06:59:27Z")

</div>

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  
`[1] 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](https://github.com/pymc-devs/resources/blob/master/Rethinking_2/Chp_07.ipynb) (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!

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [June 8, 2021, 8:32pm UTC](https://discourse.pymc.io/t/is-there-any-way-to-to-calculate-lppd-pymc3/7519/2 "2021-06-08T20:32:02Z")

</div>

See [Differences of PyMC3 posterior compared to Numpyro and R rethinking packages - #2 by ricardoV94](https://discourse.pymc.io/t/differences-of-pymc3-posterior-compared-to-numpyro-and-r-rethinking-packages/7556/2)