# How to sample from estimated posterior distribution?

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.

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_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_diff = attack_rate[hometeam] - defense_rate[awayteam] + \

away_diff = attack_rate[awayteam] - defense_rate[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

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> using Turing