When calling generated_quantities there is always a huge number of warnings with the following message
Warning: the following keys were not found in `vi`, and thus `kernel!` was not applied to these: ["lp"]
└ @ DynamicPPL ~/.julia/packages/DynamicPPL/h8FWT/src/varinfo.jl:1314
Is there some way I can get around this warning messages, alternative to get the same functionality without need to run generated_quantities.
Since Turing.jl now supports the great Dirac function, I don’t have that much need to use generated_quantities, but sometimes it’s nice to be able to work with non parameters or collection of parameters.
Here’s an example where I use generated_quantities to show the probability distributions of the combination of female, recovery, and drug (using my show_var_dist_pct defined below). The summary of the chains is shown first but is then “drowned” in the huge number of warnings from generated_quantities.
using Turing
# Show distribution of a variable in a MCMCChain
# Sort the dictionary in order of decreasing occurrence (percentage)
# Examples:
# - show_var_dist_pct(chains, :n) show all entries
# - show_var_dist_pct(chains, :n, 10) show first 10 entries (e.g. for large tables)
function show_var_dist_pct(chains::Chains, var, num=0)
if var in chains.name_map.parameters
println("Distributions of variable $var (num:$num)")
len = length(vcat(chains[var]...)) # handle multiple chains
c = 0
for kv in sort(collect(make_hash(chains[var])),by=x->x[2],rev=true)
c += 1
if (num == 0) || (num > 0 && c <= num)
@printf "%-3.5f => % 7d (%2.6f)\n" kv[1] kv[2] kv[2]/len
end
end
else
# println("Variable $var is not in chains")
end
end
#
# Simple model of Simpson's "paradox"
#
@model function simpson(problem)
female ~ flip(0.5)
drug ~ female ? flip(10/40.0) : flip(30/40.0)
recovery ~ if drug == true && female==true flip(0.2)
elseif drug == true && female==false flip(0.6)
elseif drug == false && female==true flip(0.3)
elseif drug == false && female==false flip(0.7)
else flip(0.0)
end
# 6. prob(recovery|(\+ drug,\+female)): 0.7
true ~ Dirac(drug == false && female == false)
return female, recovery, drug
end
model = simpson(problem)
chains = sample(model, MH(), 10_000)
display(chains)
# This the distributions of the different combinations of
# female, recovery, group
genq = generated_quantities(model, chains)
show_var_dist_pct(genq)
Here’s an example of the probability distribution for this problem:
Distributions of variable (num:0)
(false, true, false) => 29007 (0.725175)
(false, false, false) => 10970 (0.27425)
(false, false, true) => 17 (0.000425)
(true, true, false) => 4 (0.0001)
(true, true, true) => 2 (5.0e-5)
(The full model is here: http://hakank.org/julia/turing/simpson.jl )