How to sample from estimated posterior distribution?

New to Julia
1

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_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!

2

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

help?> predict

should give you the information you need.