# How to sample from estimated posterior distribution?
**URL:** https://discourse.julialang.org/t/how-to-sample-from-estimated-posterior-distribution/59975
**Category:** New to Julia
**Tags:** turing
**Created:** [April 25, 2021, 8:24am UTC](https://discourse.julialang.org/t/how-to-sample-from-estimated-posterior-distribution/59975 "2021-04-25T08:24:14Z")
**Posts on this page:** 2
**Page:** 1
### Author:  [@Abhishek\_Sharma](https://discourse.julialang.org/u/Abhishek_Sharma)
#### Post date: [April 25, 2021, 8:24am UTC](https://discourse.julialang.org/t/how-to-sample-from-estimated-posterior-distribution/59975/1 "2021-04-25T08:24:14Z")
Hi,
I am new to both Bayesian stuff and Julia. I’m not sure if this question is too basic or has already been answered but I have some code which I wanted to translate into turing from pymc3. It’s from this [blog post](http://gijskoot.nl/bayesian/sports/soccer/predictions/pymc3/2018/02/07/knvb-model.html).
It’s about predicting goals for soccer matches using a multi-level bayesian approach.The relevant pymc3 code for the model is this:
> with pm.Model() as model:
>
> ```
> attack_hyper = pm.Normal('attack_hyper', 0, 1)
> attack_hyper_sd = pm.Gamma('attack_hyper_sd', 2, 2)
> attack_rate_nc = pm.Normal('attack_rate_nc', 0, 1, shape = n_teams)
> attack_rate = pm.Deterministic('attack_rate', attack_hyper + attack_rate_nc * attack_hyper_sd)
>
> defense_rate_nc = pm.Normal('defense_rate_nc', 0, 1, shape = n_teams)
> defense_hyper_sd = pm.Gamma('defense_hyper_sd', 2, 2)
> defense_rate = pm.Deterministic('defense_rate', defense_rate_nc * defense_hyper_sd)
>
> home_attack_hyper_sd = pm.Gamma('hahsd', 2, 2)
> home_attack_hyper = pm.Normal('hah', 0, 1)
> home_attack_nc = pm.Normal('hah_nc', 0, 1, shape = n_teams)
> home_attack_advantage = pm.Deterministic('home_aadv', home_attack_hyper + home_attack_nc * home_attack_hyper_sd)
>
> home_defense_hyper = pm.Normal('hdh', 0, 1)
> home_defense_hyper_sd = pm.Gamma("hdhsd", 2, 2)
> home_defense_nc = pm.Normal('hdh_nc', 0, 1, shape = n_teams)
> home_defense_advantage = pm.Deterministic('home_dadv', home_defense_hyper + home_defense_nc * home_defense_hyper_sd)
>
> home_diff = attack_rate[hometeam] - defense_rate[awayteam] + \
> home_attack_advantage[hometeam]
>
> away_diff = attack_rate[awayteam] - defense_rate[hometeam] - \
> home_defense_advantage[hometeam]
>
> home_goals = pm.Poisson('home_goals', pm.math.exp(home_diff), observed = played.fthg.values)
> away_goals = pm.Poisson('away_goals', pm.math.exp(away_diff), observed = played.ftag.values)
>
> trace = pm.sample(1000, tune=500)
> ppc = pm.sample_posterior_predictive(trace, samples = 1000)
>
> ```
>
> ppc
I’m mostly unsure about one thing:
1.) And how do I sample the posterior predictions?
I’m still not sure how to get back the home\_goals and away\_goals from the model like I would in pymc3 when I call `sample_posterior_predictive`.
Here’s the Turing code that I translated:
> @model function bayesian\_model(fthg, ftag)
> n\_teams = 20
> attack\_hyper ~ Normal(0, 1)
> attack\_hyper\_sd ~ Gamma(2, 2)
> attack\_rate\_nc ~ MvNormal(n\_teams, 1)
> attack\_rate = attack\_hyper .+ attack\_rate\_nc .\* attack\_hyper\_sd
>
> ```
> defense_rate_nc ~ MvNormal(n_teams, 1)
> defense_hyper_sd ~ Gamma(2, 2)
> defense_rate = defense_rate_nc .* defense_hyper_sd
>
> home_attack_hyper_sd ~ Gamma(2, 2)
> home_attack_hyper ~ Normal(0, 1)
> home_attack_nc ~ MvNormal(n_teams, 1)
> home_attack_advantage = home_attack_hyper .+ home_attack_nc .* home_attack_hyper_sd
>
> home_defense_hyper ~ Normal(0, 1)
> home_defense_hyper_sd ~ Gamma(2, 2)
> home_defense_nc ~ MvNormal(n_teams, 1)
> home_defense_advantage = home_defense_hyper .+ home_defense_nc .* home_defense_hyper_sd
>
> home_diff = attack_rate[hometeam] .- defense_rate[awayteam] .+ home_attack_advantage[hometeam]
> away_diff = attack_rate[awayteam] .- defense_rate[hometeam] .- home_defense_advantage[hometeam]
>
> for i in 1:length(home_diff)
> fthg[i] ~ Poisson(exp(home_diff[i]))
> ftag[i] ~ Poisson(exp(away_diff[i]))
>
> end
>
> ```
>
> end
>
> ####
>
> model = bayesian\_model(fthg, ftag)
> chain = sample(model, NUTS(0.65), iterations);
TLDR: How do I get back the Poisson distributed goals estimated by the model after seeing the actual goals values (`fthg` and `ftag`)? The chain doesn’t have it as I (naively) expected
Thanks in advance!
---
### Author:  [@trahflow](https://discourse.julialang.org/u/trahflow)
#### Post date: [July 3, 2021, 12:00pm UTC](https://discourse.julialang.org/t/how-to-sample-from-estimated-posterior-distribution/59975/2 "2021-07-03T12:00:34Z")
Hi @Abhishek_Sharma ,
very late answer, so maybe you’ve figured this out in the meantime, but I’ll still post for reference:
Turing has a `predict` method, that is probably doing exactly what you’re looking for. Unfortunately I think the online documentation is missing that part. However
```julia
julia> using Turing
help?> predict
```
should give you the information you need.