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)
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
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
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 (
ftag)? The chain doesn’t have it as I (naively) expected
Thanks in advance!