[playground] Julia Turing.jl : Bayesian Cognitive Modeling - Comparing binomial rates
Github: https://github.com/quangtiencs/bayesian-cognitive-modeling-with-turing.jl
Bayesian Cognitive Modeling is one of the classical books for Bayesian Inference. The old version used WinBUGS/JAG software as the main implementation. You can find other implementations, such as Stan and PyMC, in the below link. I reimplemented these source codes with Julia Programming Language & Turing library in this tutorial.
- WinBUGS/JAGS (official) & Stan: https://bayesmodels.com/
- PyMC: https://github.com/pymc-devs/pymc-resources/tree/main/BCM
using Pkg
using Logging
using DynamicPPL, Turing
using Zygote, ReverseDiff
using StatsPlots, Random
using LaTeXStrings
using CSV
using DataFrames
using SpecialFunctions
using LinearAlgebra
using FillArrays
using CSV, DataFrames
using LogExpFunctions
using KernelDensity
using Interpolations
using Statistics
using StatsBase
using Dierckx
format=:png
:png
Random.seed!(9)
TaskLocalRNG()
9.1 Equality of proportions
$$ s_{1} \sim \text{Binomial}(\theta_{1},n_{1}) $$
$$ s_{2} \sim \text{Binomial}(\theta_{2},n_{2}) $$
$$ \theta_{1} \sim \text{Beta}(1,1) $$
$$ \theta_{2} \sim \text{Beta}(1,1) $$
$$ \delta = \theta_{1} -\theta_{2} $$
s1 = 424
s2 = 5416
n1 = 777
n2 = 9072
@model function EqualityOfProportions(s1, n1, s2, n2)
theta1 ~ Beta(1, 1)
theta2 ~ Beta(1, 1)
s1 ~ Binomial(n1, theta1)
s2 ~ Binomial(n2, theta2)
end
iterations = 5_000
burnin = 1_000
model_equality_of_proportions = EqualityOfProportions(s1, n1, s2, n2)
chain = sample(model_equality_of_proportions, NUTS(), iterations, burnin=burnin)
┌ Info: Found initial step size
└ ϵ = 0.025
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00
Chains MCMC chain (5000×14×1 Array{Float64, 3}):
Iterations = 1001:1:6000
Number of chains = 1
Samples per chain = 5000
Wall duration = 8.34 seconds
Compute duration = 8.34 seconds
parameters = theta1, theta2
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
theta1 0.5457 0.0177 0.0003 4760.9049 3267.6844 1.0003 ⋯
theta2 0.5970 0.0051 0.0001 5112.7350 3402.6623 0.9999 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
theta1 0.5098 0.5338 0.5458 0.5573 0.5806
theta2 0.5869 0.5935 0.5969 0.6006 0.6069
plot(chain, size=(800, 500),
left_margin=10Plots.mm, bottom_margin=10Plots.mm, fmt=format)
function plot_equality_of_proportions(chain, prior_dist, zoom, size=(600, 400))
delta = (chain[:theta1] .- chain[:theta2]) |> vec
plot(zoom, [pdf(prior_dist, e) for e in zoom], size=size, fmt=format, linewidth=2, label="Prior")
k = kde(delta)
plot!(zoom, pdf(k, zoom), linewidth=2, ls=:dash, label="Posterior")
plot!([0, 0], [pdf(prior_dist, 0), pdf(k, 0)], linewidth=2, label=false)
scatter!([0, 0], [pdf(prior_dist, 0), pdf(k, 0)], label=false)
bayes_factor = pdf(k, 0) / pdf(prior_dist, 0)
println("Bayes Factor $bayes_factor")
xlabel!("Delta")
ylabel!("Density")
xlims!(min(zoom...), max(zoom...))
end
plot_equality_of_proportions (generic function with 2 methods)
plot(
plot_equality_of_proportions(chain, TriangularDist(-1, 1, 0), -1:0.001:1),
plot_equality_of_proportions(chain, TriangularDist(-1, 1, 0), -0.01:0.0001:0.01),
size=(1200, 400), left_margin=10Plots.mm, bottom_margin=10Plots.mm, fmt=format
)
Bayes Factor 0.3584135097709259
Bayes Factor 0.3584135097709259
9.2 Order-restricted equality of proportions
$$ s_{1} \sim \text{Binomial}(\theta_{1},n_{1}) $$
$$ s_{2} \sim \text{Binomial}(\theta_{2},n_{2}) $$
$$ (\theta_{1},\theta_{2}) \sim \text{Uniform}(0,1)^2, \theta_{2} \gt \theta_{1} $$
$$ \delta = \theta_{1} -\theta_{2} $$
s1 = 424
s2 = 5416
n1 = 777
n2 = 9072
phi(x) = cdf(Normal(0, 1), x)
function normcdf1(thetap1, thetap2)
angle = 45 * Float64(pi) / 180
return phi((cos(angle) * thetap1) - (sin(angle) * abs(thetap2)))
end
function normcdf2(thetap1, thetap2)
angle = 45 * Float64(pi) / 180
return phi((sin(angle) * thetap1) + (cos(angle) * abs(thetap2)))
end
@model function OrderRestrictedEqualityOfProportions()
theta2a ~ Uniform(0, 1)
theta1a ~ Uniform(0, theta2a)
deltaa = theta1a - theta2a
thetap ~ MvNormal([0, 0], [1, 1])
theta1e = normcdf1(thetap[1], thetap[2])
theta2e = normcdf2(thetap[1], thetap[2])
deltae = theta1e - theta2e
return (deltaa=deltaa, deltae=deltae, theta1e=theta1e, theta2e=theta2e)
end
iterations = 30_000
model_order_restricted_equality_of_proportions = OrderRestrictedEqualityOfProportions()
chain = sample(model_order_restricted_equality_of_proportions, SMC(), iterations)
chains_params = Turing.MCMCChains.get_sections(chain, :parameters)
quantities = generated_quantities(model_order_restricted_equality_of_proportions, chains_params)
chain
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:01
Chains MCMC chain (30000×6×1 Array{Float64, 3}):
Log evidence = 0.0
Iterations = 1:1:30000
Number of chains = 1
Samples per chain = 30000
Wall duration = 43.98 seconds
Compute duration = 43.98 seconds
parameters = theta2a, theta1a, thetap[1], thetap[2]
internals = lp, weight
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
theta2a 0.4995 0.2877 0.0017 29819.6834 29409.0199 1.0000 ⋯
theta1a 0.2493 0.2208 0.0013 29429.5913 29059.7470 1.0001 ⋯
thetap[1] 0.0054 0.9919 0.0058 29523.0793 29504.0565 1.0002 ⋯
thetap[2] 0.0031 1.0001 0.0058 29401.8339 28564.3645 1.0001 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
theta2a 0.0254 0.2495 0.4993 0.7489 0.9751
theta1a 0.0038 0.0675 0.1845 0.3808 0.7909
thetap[1] -1.9402 -0.6612 0.0024 0.6752 1.9533
thetap[2] -1.9516 -0.6748 -0.0002 0.6744 1.9698
# theta1 = chain[:theta1a] |> vec |> collect
# theta2 = chain[:theta2a] |> vec |> collect;
alpha = 0.03
theta1e = map(x -> x[:theta1e], quantities)
theta2e = map(x -> x[:theta2e], quantities);
deltaa = map(x -> x[:deltaa], quantities);
deltae = map(x -> x[:deltae], quantities);
p1 = scatter(chain[:theta1a], chain[:theta2a],
markersize=2, alpha=alpha, size=(300,300), legend=false)
xlabel!(L"\theta_1")
ylabel!(L"\theta_2")
p2 = scatter(theta1e, theta2e,
markersize=2, alpha=alpha, size=(300,300), legend=false)
xlabel!(L"\theta_1")
ylabel!(L"\theta_2")
bins = 15
h = fit(Histogram, deltaa |> vec, nbins=bins)
x_range = h.edges[1] |> collect
cen_x = [(x_range[i] + x_range[i+1])/2 for i in 1:(length(x_range)-1)]
cen_y = h.weights
spl = Spline1D(cen_x, cen_y; s=length(cen_x))
xq = minimum(cen_x):0.001:maximum(cen_x)
p3 = plot(xq, spl(xq), label="Approximate", size=(300,300),
color="green", linewidth=2, ls=:dash, yticks=false)
h = fit(Histogram, deltae |> vec, nbins=bins)
x_range = h.edges[1] |> collect
cen_x = [(x_range[i] + x_range[i+1])/2 for i in 1:(length(x_range)-1)]
cen_y = h.weights
spl = Spline1D(cen_x, cen_y; s=length(cen_x))
xq = minimum(cen_x):0.001:maximum(cen_x)
plot!(xq, spl(xq), label="Extract", color="red", linewidth=2)
xlims!(-1, 0)
xlabel!(L"\delta = \theta_1 - \theta_2")
ylabel!("Density")
plot(p1, p2, p3,
size=(900, 300),
fmt=format, layout=(1,3),bottom_margin=5Plots.mm, left_margin=5Plots.mm
)
s1 = 424
s2 = 5416
n1 = 777
n2 = 9072
phi(x) = cdf(Normal(0, 1), x)
function normcdf1(thetap1, thetap2)
angle = 45 * Float64(pi) / 180
return phi((cos(angle) * thetap1) - (sin(angle) * abs(thetap2)))
end
function normcdf2(thetap1, thetap2)
angle = 45 * Float64(pi) / 180
return phi((sin(angle) * thetap1) + (cos(angle) * abs(thetap2)))
end
@model function OrderRestrictedEqualityOfProportions(s1, n1, s2, n2)
thetap ~ MvNormal([0, 0], [1, 1])
theta1e = normcdf1(thetap[1], thetap[2])
theta2e = normcdf2(thetap[1], thetap[2])
deltae = theta1e - theta2e
s1 ~ Binomial(n1, theta1e)
s2 ~ Binomial(n2, theta2e)
return (deltae=deltae, theta1e=theta1e, theta2e=theta2e)
end
iterations = 5_000
burnin = 2_500
model_order_restricted_equality_of_proportions = OrderRestrictedEqualityOfProportions(s1, n1, s2, n2)
chain = sample(model_order_restricted_equality_of_proportions, NUTS(), iterations, burnin=burnin)
chains_params = Turing.MCMCChains.get_sections(chain, :parameters)
quantities = generated_quantities(model_order_restricted_equality_of_proportions, chains_params)
chain
┌ Info: Found initial step size
└ ϵ = 0.00625
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00
Chains MCMC chain (5000×14×1 Array{Float64, 3}):
Iterations = 1001:1:6000
Number of chains = 1
Samples per chain = 5000
Wall duration = 2.52 seconds
Compute duration = 2.52 seconds
parameters = thetap[1], thetap[2]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
thetap[1] 0.2547 0.0341 0.0009 1330.4698 1648.0405 1.0003 ⋯
thetap[2] -0.0062 0.0986 0.0216 28.2434 468.4895 1.0375 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
thetap[1] 0.1880 0.2316 0.2544 0.2776 0.3208
thetap[2] -0.1470 -0.0955 -0.0415 0.0902 0.1487
plot(chain, fmt=format)
function plot_order_restricted_equality_of_proportions(quantities, prior_dist, zoom, size=(600, 400))
delta = map(x -> x[:deltae], quantities) |> collect |> vec
plot(zoom, [pdf(prior_dist, e) for e in zoom], size=size, fmt=format, linewidth=2, label="Prior")
k = kde(delta)
plot!(zoom, pdf(k, zoom), linewidth=2, ls=:dash, label="Posterior")
plot!([0, 0], [pdf(prior_dist, 0), pdf(k, 0)], linewidth=2, label=false)
scatter!([0, 0], [pdf(prior_dist, 0), pdf(k, 0)], label=false)
bayes_factor = pdf(k, 0) / pdf(prior_dist, 0)
println("Bayes Factor $bayes_factor")
xlabel!("Delta")
ylabel!("Density")
xlims!(min(zoom...), max(zoom...))
end
plot_order_restricted_equality_of_proportions (generic function with 2 methods)
plot(
plot_order_restricted_equality_of_proportions(quantities,TriangularDist(-1, 0, 0), -1:0.001:0.02),
plot_order_restricted_equality_of_proportions(quantities,TriangularDist(-1, 0, 0), -0.02:0.0001:0.0001),
size=(1200, 400), left_margin=10Plots.mm, bottom_margin=10Plots.mm, fmt=format
)
Bayes Factor 0.19015084480044653
Bayes Factor 0.19015084480044653
delta = map(x -> x[:deltae], quantities) |> collect |> vec
k = kde(delta)
bayes_factor_2 = pdf(k, 0) / pdf(TriangularDist(-1, 0, 0), 0)
0.19015084480044653
9.3 Comparing within-subject proportions
$$ \mu \sim \text{Gaussian}(0,1)_{\mathcal I(0,∞)} $$
$$ \sigma \sim \text{Uniform}(0,10) $$
$$ \delta \sim \text{Gaussian}(0,1)_{\mathcal I(0,∞)} $$
$$ \sigma_{\alpha} \sim \text{Uniform}(0,10) $$
$$ \mu_{\alpha} = \delta\sigma_{\alpha} $$
$$ \alpha_{i} \sim \text{Gaussian}(\mu_{\alpha},1/\sigma^2_{\alpha}) $$
$$ \phi^n_{i} \sim \text{Gaussian}(\mu,1/\sigma^2) $$
$$ \phi^b_{i} = \phi^n_{i}+\alpha_{i} $$
$$ \theta^n_{i} = \Phi(\phi^n_{i}) $$
$$ \theta^b_{i} = \Phi(\phi^b_{i}) $$
$$ s^n_{i} \sim \text{Binomial}(\theta^n_{i},n^n) $$
$$ s^b_{i} \sim \text{Binomial}(\theta^b_{i},n^b) $$
# fmt: off
### Zeelenberg data:
# Study Both:
sb = [
15, 11, 15, 14, 15, 18, 16, 16, 18, 16, 15, 13, 18, 12, 11, 13, 17, 18, 16, 11, 17, 18,
12, 18, 18, 14, 21, 18, 17, 10, 11, 12, 16, 18, 17, 15, 19, 12, 21, 15, 16, 20, 15, 19,
16, 16, 14, 18, 16, 19, 17, 11, 19, 18, 16, 16, 11, 19, 18, 12, 15, 18, 20, 8, 12, 19,
16, 16, 16, 12, 18, 17, 11, 20
]
nb = 21
# Study Neither:
sn = [
15, 12, 14, 15, 13, 14, 10, 17, 13, 16, 16, 10, 15, 15, 10, 14, 17, 18, 19, 12, 19, 18,
10, 18, 16, 13, 15, 20, 13, 15, 13, 14, 19, 19, 19, 18, 13, 12, 19, 16, 14, 17, 15, 16,
15, 16, 13, 15, 14, 19, 12, 11, 17, 13, 18, 13, 13, 19, 18, 13, 13, 16, 18, 14, 14, 17,
12, 12, 16, 14, 16, 18, 13, 13
]
nn = 21
ns = length(sb)
74
phi(x) = cdf(Normal(0, 1), x)
@model function ComparingWithinSubjectProportions(sb, sn)
mu ~ truncated(Normal(0, 1), lower=0)
sigma ~ Uniform(0, 10)
delta ~ truncated(Normal(0, 1), lower=0)
sigma_alpha ~ Uniform(0, 10)
mu_alpha = delta * sigma_alpha
alpha_n ~ filldist(Normal(mu_alpha, sigma_alpha), length(sb))
phi_n ~ filldist(Normal(mu, sigma), length(sb))
phi_b = phi_n + alpha_n
theta_n = phi.(phi_n)
theta_b = phi.(phi_b)
for i in eachindex(sb)
sb[i] ~ Binomial(nb, theta_b[i])
sn[i] ~ Binomial(nn, theta_n[i])
end
return (theta_n=theta_n, theta_b=theta_b)
end
iterations = 10_000
burnin = 5_000
model_comparing_within_subject_proportions = ComparingWithinSubjectProportions(sb, sn)
chain = sample(model_comparing_within_subject_proportions, NUTS(), iterations, burnin=burnin)
chains_params = Turing.MCMCChains.get_sections(chain, :parameters)
quantities = generated_quantities(model_comparing_within_subject_proportions, chains_params)
chain
┌ Warning: failed to find valid initial parameters in 10 tries; consider providing explicit initial parameters using the `init_params` keyword
└ @ Turing.Inference ~/.julia/packages/Turing/FSBW7/src/inference/hmc.jl:172
┌ Info: Found initial step size
└ ϵ = 0.4
Sampling: 100%|█████████████████████████████████████████| Time: 0:01:49
Chains MCMC chain (10000×164×1 Array{Float64, 3}):
Iterations = 1001:1:11000
Number of chains = 1
Samples per chain = 10000
Wall duration = 117.99 seconds
Compute duration = 117.99 seconds
parameters = mu, sigma, delta, sigma_alpha, alpha_n[1], alpha_n[2], alpha_n[3], alpha_n[4], alpha_n[5], alpha_n[6], alpha_n[7], alpha_n[8], alpha_n[9], alpha_n[10], alpha_n[11], alpha_n[12], alpha_n[13], alpha_n[14], alpha_n[15], alpha_n[16], alpha_n[17], alpha_n[18], alpha_n[19], alpha_n[20], alpha_n[21], alpha_n[22], alpha_n[23], alpha_n[24], alpha_n[25], alpha_n[26], alpha_n[27], alpha_n[28], alpha_n[29], alpha_n[30], alpha_n[31], alpha_n[32], alpha_n[33], alpha_n[34], alpha_n[35], alpha_n[36], alpha_n[37], alpha_n[38], alpha_n[39], alpha_n[40], alpha_n[41], alpha_n[42], alpha_n[43], alpha_n[44], alpha_n[45], alpha_n[46], alpha_n[47], alpha_n[48], alpha_n[49], alpha_n[50], alpha_n[51], alpha_n[52], alpha_n[53], alpha_n[54], alpha_n[55], alpha_n[56], alpha_n[57], alpha_n[58], alpha_n[59], alpha_n[60], alpha_n[61], alpha_n[62], alpha_n[63], alpha_n[64], alpha_n[65], alpha_n[66], alpha_n[67], alpha_n[68], alpha_n[69], alpha_n[70], alpha_n[71], alpha_n[72], alpha_n[73], alpha_n[74], phi_n[1], phi_n[2], phi_n[3], phi_n[4], phi_n[5], phi_n[6], phi_n[7], phi_n[8], phi_n[9], phi_n[10], phi_n[11], phi_n[12], phi_n[13], phi_n[14], phi_n[15], phi_n[16], phi_n[17], phi_n[18], phi_n[19], phi_n[20], phi_n[21], phi_n[22], phi_n[23], phi_n[24], phi_n[25], phi_n[26], phi_n[27], phi_n[28], phi_n[29], phi_n[30], phi_n[31], phi_n[32], phi_n[33], phi_n[34], phi_n[35], phi_n[36], phi_n[37], phi_n[38], phi_n[39], phi_n[40], phi_n[41], phi_n[42], phi_n[43], phi_n[44], phi_n[45], phi_n[46], phi_n[47], phi_n[48], phi_n[49], phi_n[50], phi_n[51], phi_n[52], phi_n[53], phi_n[54], phi_n[55], phi_n[56], phi_n[57], phi_n[58], phi_n[59], phi_n[60], phi_n[61], phi_n[62], phi_n[63], phi_n[64], phi_n[65], phi_n[66], phi_n[67], phi_n[68], phi_n[69], phi_n[70], phi_n[71], phi_n[72], phi_n[73], phi_n[74]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
mu 0.5973 0.0496 0.0021 572.9467 673.8812 1.0002 ⋯
sigma 0.2942 0.0425 0.0026 295.9773 222.7146 1.0057 ⋯
delta 0.7468 0.4114 0.0290 239.6704 192.7663 1.0029 ⋯
sigma_alpha 0.1374 0.0672 0.0071 77.2335 81.1119 1.0153 ⋯
alpha_n[1] 0.0729 0.1402 0.0028 2451.4641 2172.0902 1.0032 ⋯
alpha_n[2] 0.0236 0.1395 0.0030 2222.5928 2500.1767 1.0025 ⋯
alpha_n[3] 0.0820 0.1413 0.0025 2676.7747 2543.1883 1.0019 ⋯
alpha_n[4] 0.0526 0.1374 0.0023 2976.3506 2765.4212 1.0027 ⋯
alpha_n[5] 0.0942 0.1410 0.0037 1527.6959 1913.7432 1.0057 ⋯
alpha_n[6] 0.1459 0.1546 0.0072 509.7721 2006.5145 1.0045 ⋯
alpha_n[7] 0.1401 0.1509 0.0070 484.8464 2013.3353 1.0049 ⋯
alpha_n[8] 0.0743 0.1410 0.0022 3811.5153 2706.9510 1.0037 ⋯
alpha_n[9] 0.1530 0.1532 0.0077 411.8584 2122.9226 1.0031 ⋯
alpha_n[10] 0.0858 0.1393 0.0028 2310.9980 3019.9447 1.0030 ⋯
alpha_n[11] 0.0646 0.1399 0.0024 2966.3679 3203.7794 1.0021 ⋯
alpha_n[12] 0.0790 0.1391 0.0027 2427.1029 2194.5175 1.0042 ⋯
alpha_n[13] 0.1423 0.1531 0.0064 623.3686 2052.0431 1.0011 ⋯
alpha_n[14] 0.0198 0.1430 0.0044 1036.2783 1623.8591 1.0039 ⋯
alpha_n[15] 0.0480 0.1411 0.0027 2582.9114 2656.8088 1.0068 ⋯
alpha_n[16] 0.0430 0.1388 0.0027 2340.8110 2686.8980 1.0017 ⋯
alpha_n[17] 0.0957 0.1460 0.0032 2175.9008 2445.5079 1.0047 ⋯
alpha_n[18] 0.1123 0.1465 0.0042 1191.3904 2545.0541 1.0010 ⋯
alpha_n[19] 0.0555 0.1390 0.0026 2397.0830 2819.8775 1.0040 ⋯
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱
1 column and 129 rows omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
mu 0.4980 0.5650 0.5976 0.6302 0.6949
sigma 0.2172 0.2646 0.2922 0.3213 0.3855
delta 0.0803 0.4403 0.6930 1.0081 1.6793
sigma_alpha 0.0331 0.0835 0.1335 0.1821 0.2778
alpha_n[1] -0.2166 -0.0060 0.0641 0.1531 0.3670
alpha_n[2] -0.2976 -0.0478 0.0325 0.1060 0.2963
alpha_n[3] -0.2084 0.0039 0.0722 0.1615 0.3828
alpha_n[4] -0.2359 -0.0204 0.0515 0.1304 0.3290
alpha_n[5] -0.1774 0.0122 0.0834 0.1721 0.4091
alpha_n[6] -0.1097 0.0407 0.1209 0.2350 0.5008
alpha_n[7] -0.1171 0.0396 0.1164 0.2259 0.4853
alpha_n[8] -0.2110 -0.0036 0.0650 0.1529 0.3780
alpha_n[9] -0.0951 0.0474 0.1292 0.2402 0.5100
alpha_n[10] -0.1838 0.0032 0.0759 0.1660 0.3806
alpha_n[11] -0.2221 -0.0123 0.0592 0.1443 0.3494
alpha_n[12] -0.2036 -0.0025 0.0711 0.1579 0.3756
alpha_n[13] -0.1145 0.0423 0.1180 0.2298 0.4902
alpha_n[14] -0.3064 -0.0517 0.0277 0.1068 0.2891
alpha_n[15] -0.2610 -0.0254 0.0496 0.1304 0.3223
alpha_n[16] -0.2654 -0.0292 0.0440 0.1248 0.3178
alpha_n[17] -0.1837 0.0101 0.0814 0.1757 0.4157
alpha_n[18] -0.1608 0.0203 0.0957 0.1948 0.4408
alpha_n[19] -0.2334 -0.0221 0.0532 0.1356 0.3367
⋮ ⋮ ⋮ ⋮ ⋮ ⋮
129 rows omitted
function plot_comparing_within_subject_proportions(chain, prior_dist, zoom, size=(600, 400))
delta = chain[:delta] |> vec
plot(zoom, [pdf(prior_dist, e) for e in zoom],
size=size, fmt=format, linewidth=2, label="Prior")
k = kde(delta)
plot!(zoom, pdf(k, zoom), linewidth=2, ls=:dash, label="Posterior")
plot!([0, 0], [pdf(prior_dist, 0), pdf(k, 0)], linewidth=2, label=false)
scatter!([0, 0], [pdf(prior_dist, 0), pdf(k, 0)], label=false)
bayes_factor = pdf(k, 0) / pdf(prior_dist, 0)
println("Bayes Factor $bayes_factor")
xlims!(min(zoom...), max(zoom...))
end
plot_comparing_within_subject_proportions (generic function with 2 methods)
plot_comparing_within_subject_proportions(chain,
truncated(Normal(0, 1), lower=0), -0.1:0.01:4)
Bayes Factor 0.20233647781390138
9.4 Comparing between-subject proportions
$$ \delta \sim \text{Gaussian}(0,1)$$ $$ \mu \sim \text{Gaussian}(0,1)$$ $$ \sigma \sim \text{Uniform}(0,10)$$ $$ \alpha = \delta\sigma$$
$$ \phi^c_{i} \sim \text{Gaussian}(\mu,+\alpha/2,1/\sigma^2)$$ $$ \phi^a_{j} \sim \text{Gaussian}(\mu,-\alpha/2,1/\sigma^2)$$
$$ \theta^c_{i} = \Phi(\phi^c_{i})$$ $$ \theta^a_{j} = \Phi(\phi^a_{j})$$
$$ s^c_{i} \sim \text{Binomial}(\theta^c_{i},n^c_{i})$$ $$ s^a_{j} \sim \text{Binomial}(\theta^a_{j},n^a_{j})$$
# fmt: off
### Geurts data:
# Normal Controls:
numerrors1 = [
15, 10, 61, 11, 60, 44, 63, 70, 57, 11, 67, 21, 89, 12, 63, 11, 96, 10, 37, 19, 44,
18, 78, 27, 60, 14
]
nc = [
89, 74, 128, 87, 128, 121, 128, 128, 128, 78, 128, 106, 128, 83, 128, 100, 128,
73, 128, 86, 128, 86, 128, 100, 128, 79
]
kc = nc - numerrors1
nsc = length(kc)
# ADHD:
numerrors2 = [
88, 50, 58, 17, 40, 18, 21, 50, 21, 69, 19, 29, 11, 76, 46, 36, 37, 72, 27, 92, 13,
39, 53, 31, 49, 57, 17, 10, 12, 21, 39, 43, 49, 17, 39, 13, 68, 24, 21, 27, 48, 54,
41, 75, 38, 76, 21, 41, 61, 24, 28, 21
]
na = [
128, 128, 128, 86, 128, 117, 89, 128, 110, 128, 93, 107, 87, 128, 128, 113, 128,
128, 98, 128, 93, 116, 128, 116, 128, 128, 93, 86, 86, 96, 128, 128, 128, 86, 128,
78, 128, 111, 100, 95, 128, 128, 128, 128, 128, 128, 98, 127, 128, 93, 110, 96
]
ka = na - numerrors2
nsa = length(ka)
52
phi(x) = cdf(Normal(0, 1), x)
@model function ComparingBetweenSubjectProportions(ka, na, kc, nc)
mu ~ Normal(0, 1)
sigma ~ Uniform(0, 10)
delta ~ Normal(0, 1)
alpha = delta * sigma
phi_c ~ filldist(Normal(mu + alpha / 2, sigma), length(kc))
phi_a ~ filldist(Normal(mu - alpha / 2, sigma), length(ka))
theta_c = phi.(phi_c)
theta_a = phi.(phi_a)
for i in eachindex(kc)
kc[i] ~ Binomial(nc[i], theta_c[i])
end
for i in eachindex(ka)
ka[i] ~ Binomial(na[i], theta_a[i])
end
return (theta_c=theta_c, theta_a=theta_a)
end
iterations = 10_000
burnin = 5_000
model_comparing_between_subject_proportions = ComparingBetweenSubjectProportions(ka, na, kc, nc)
chain = sample(model_comparing_between_subject_proportions, NUTS(), iterations, burnin=burnin)
chains_params = Turing.MCMCChains.get_sections(chain, :parameters)
quantities = generated_quantities(model_comparing_between_subject_proportions, chains_params)
chain
┌ Info: Found initial step size
└ ϵ = 0.4
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:17
Chains MCMC chain (10000×93×1 Array{Float64, 3}):
Iterations = 1001:1:11000
Number of chains = 1
Samples per chain = 10000
Wall duration = 24.1 seconds
Compute duration = 24.1 seconds
parameters = mu, sigma, delta, phi_c[1], phi_c[2], phi_c[3], phi_c[4], phi_c[5], phi_c[6], phi_c[7], phi_c[8], phi_c[9], phi_c[10], phi_c[11], phi_c[12], phi_c[13], phi_c[14], phi_c[15], phi_c[16], phi_c[17], phi_c[18], phi_c[19], phi_c[20], phi_c[21], phi_c[22], phi_c[23], phi_c[24], phi_c[25], phi_c[26], phi_a[1], phi_a[2], phi_a[3], phi_a[4], phi_a[5], phi_a[6], phi_a[7], phi_a[8], phi_a[9], phi_a[10], phi_a[11], phi_a[12], phi_a[13], phi_a[14], phi_a[15], phi_a[16], phi_a[17], phi_a[18], phi_a[19], phi_a[20], phi_a[21], phi_a[22], phi_a[23], phi_a[24], phi_a[25], phi_a[26], phi_a[27], phi_a[28], phi_a[29], phi_a[30], phi_a[31], phi_a[32], phi_a[33], phi_a[34], phi_a[35], phi_a[36], phi_a[37], phi_a[38], phi_a[39], phi_a[40], phi_a[41], phi_a[42], phi_a[43], phi_a[44], phi_a[45], phi_a[46], phi_a[47], phi_a[48], phi_a[49], phi_a[50], phi_a[51], phi_a[52]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
mu 0.4491 0.0566 0.0005 15336.5752 7682.3289 1.0000 ⋯
sigma 0.4555 0.0411 0.0003 14459.6877 7892.1440 1.0001 ⋯
delta -0.0625 0.2396 0.0020 14827.7631 8287.6669 0.9999 ⋯
phi_c[1] 0.9098 0.1452 0.0010 22685.1207 7152.0471 1.0001 ⋯
phi_c[2] 1.0183 0.1655 0.0012 20182.5693 6811.8721 1.0001 ⋯
phi_c[3] 0.0802 0.1104 0.0008 21627.7966 6682.0873 1.0010 ⋯
phi_c[4] 1.0620 0.1584 0.0011 19051.9682 6810.8421 0.9999 ⋯
phi_c[5] 0.0996 0.1083 0.0007 21239.4743 6932.1835 1.0011 ⋯
phi_c[6] 0.3558 0.1143 0.0008 19799.4003 6512.6977 1.0002 ⋯
phi_c[7] 0.0428 0.1076 0.0007 23057.4906 7492.5669 0.9999 ⋯
phi_c[8] -0.0854 0.1080 0.0008 19389.8145 6662.9344 0.9999 ⋯
phi_c[9] 0.1543 0.1081 0.0007 22027.6426 6310.8263 1.0002 ⋯
phi_c[10] 0.9997 0.1627 0.0012 18576.4177 6908.7128 0.9999 ⋯
phi_c[11] -0.0293 0.1091 0.0008 19808.7096 6918.0579 1.0003 ⋯
phi_c[12] 0.8166 0.1315 0.0009 20979.7148 6735.9554 1.0003 ⋯
phi_c[13] -0.4533 0.1096 0.0008 20534.6909 6801.2302 0.9999 ⋯
phi_c[14] 0.9895 0.1575 0.0011 21462.2113 6667.1350 0.9999 ⋯
phi_c[15] 0.0450 0.1072 0.0007 23768.8986 7307.9527 1.0005 ⋯
phi_c[16] 1.1401 0.1541 0.0012 15683.2323 6845.6433 0.9999 ⋯
phi_c[17] -0.6034 0.1149 0.0008 19118.4181 7445.9207 0.9999 ⋯
phi_c[18] 1.0100 0.1667 0.0012 18105.5386 7301.9226 1.0003 ⋯
phi_c[19] 0.5493 0.1130 0.0009 17479.9917 6935.5551 1.0005 ⋯
phi_c[20] 0.7398 0.1438 0.0010 18738.2446 6726.0646 1.0000 ⋯
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱
1 column and 58 rows omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
mu 0.3373 0.4116 0.4493 0.4872 0.5587
sigma 0.3815 0.4270 0.4529 0.4809 0.5429
delta -0.5330 -0.2227 -0.0638 0.1009 0.3985
phi_c[1] 0.6269 0.8123 0.9079 1.0056 1.1952
phi_c[2] 0.7027 0.9059 1.0138 1.1272 1.3559
phi_c[3] -0.1355 0.0065 0.0811 0.1530 0.2969
phi_c[4] 0.7519 0.9549 1.0599 1.1697 1.3735
phi_c[5] -0.1118 0.0264 0.0994 0.1722 0.3143
phi_c[6] 0.1346 0.2796 0.3550 0.4327 0.5834
phi_c[7] -0.1666 -0.0282 0.0415 0.1143 0.2550
phi_c[8] -0.2989 -0.1586 -0.0844 -0.0126 0.1252
phi_c[9] -0.0572 0.0824 0.1528 0.2270 0.3676
phi_c[10] 0.6855 0.8885 0.9975 1.1074 1.3264
phi_c[11] -0.2462 -0.1024 -0.0280 0.0450 0.1832
phi_c[12] 0.5601 0.7273 0.8145 0.9035 1.0815
phi_c[13] -0.6700 -0.5248 -0.4530 -0.3804 -0.2372
phi_c[14] 0.6841 0.8843 0.9877 1.0948 1.3049
phi_c[15] -0.1695 -0.0269 0.0457 0.1182 0.2552
phi_c[16] 0.8486 1.0347 1.1364 1.2422 1.4538
phi_c[17] -0.8275 -0.6806 -0.6039 -0.5259 -0.3783
phi_c[18] 0.6884 0.8971 1.0071 1.1204 1.3428
phi_c[19] 0.3277 0.4729 0.5494 0.6252 0.7728
phi_c[20] 0.4603 0.6420 0.7392 0.8363 1.0227
⋮ ⋮ ⋮ ⋮ ⋮ ⋮
58 rows omitted
plot(chain[[:delta]], fmt=format)
function plot_comparing_between_subject_proportions(chain, prior_dist, zoom, size=(600, 400))
delta = chain[:delta] |> vec
plot(zoom, [pdf(prior_dist, e) for e in zoom], size=size, fmt=format, linewidth=2, label="Prior")
k = kde(delta)
plot!(zoom, pdf(k, zoom), linewidth=2, ls=:dash, label="Posterior")
plot!([0, 0], [pdf(prior_dist, 0), pdf(k, 0)], linewidth=2, label=false)
println(pdf(k, 0) / pdf(prior_dist, 0))
scatter!([0, 0], [pdf(prior_dist, 0), pdf(k, 0)], label=false)
bayes_factor = pdf(k, 0) / pdf(prior_dist, 0)
println("Bayes Factor $bayes_factor")
xlabel!("Delta")
ylabel!("Density")
xlims!(min(zoom...), max(zoom...))
end
plot_comparing_between_subject_proportions (generic function with 2 methods)
plot_comparing_between_subject_proportions(chain,
Normal(0,1), -3:0.01:3)
3.920024719768321
Bayes Factor 3.920024719768321
9.5 Order-restricted between-subject proportions
$$ \delta \sim \text{Gaussian}(0,1)_{\mathcal I(0,∞)}$$ $$ \mu \sim \text{Gaussian}(0,1)$$ $$ \sigma \sim \text{Uniform}(0,10)$$ $$ \alpha = \delta\sigma$$
$$ \phi^c_{i} \sim \text{Gaussian}(\mu,+\alpha/2,1/\sigma^2)$$ $$ \phi^a_{j} \sim \text{Gaussian}(\mu,-\alpha/2,1/\sigma^2)$$
$$ \theta^c_{i} = \Phi(\phi^c_{i})$$ $$ \theta^a_{j} = \Phi(\phi^a_{j})$$
$$ s^c_{i} \sim \text{Binomial}(\theta^c_{i},n^c_{i})$$ $$ s^a_{j} \sim \text{Binomial}(\theta^a_{j},n^a_{j})$$
# fmt: off
### Geurts data:
# Normal Controls:
numerrors1 = [
15, 10, 61, 11, 60, 44, 63, 70, 57, 11, 67, 21, 89, 12, 63, 11, 96, 10, 37, 19, 44,
18, 78, 27, 60, 14
]
nc = [
89, 74, 128, 87, 128, 121, 128, 128, 128, 78, 128, 106, 128, 83, 128, 100, 128,
73, 128, 86, 128, 86, 128, 100, 128, 79
]
kc = nc - numerrors1
nsc = length(kc)
# ADHD:
numerrors2 = [
88, 50, 58, 17, 40, 18, 21, 50, 21, 69, 19, 29, 11, 76, 46, 36, 37, 72, 27, 92, 13,
39, 53, 31, 49, 57, 17, 10, 12, 21, 39, 43, 49, 17, 39, 13, 68, 24, 21, 27, 48, 54,
41, 75, 38, 76, 21, 41, 61, 24, 28, 21
]
na = [
128, 128, 128, 86, 128, 117, 89, 128, 110, 128, 93, 107, 87, 128, 128, 113, 128,
128, 98, 128, 93, 116, 128, 116, 128, 128, 93, 86, 86, 96, 128, 128, 128, 86, 128,
78, 128, 111, 100, 95, 128, 128, 128, 128, 128, 128, 98, 127, 128, 93, 110, 96
]
ka = na - numerrors2
nsa = length(ka)
52
phi(x) = cdf(Normal(0, 1), x)
@model function OrderRestrictedBetweenSubjectProportions(ka, na, kc, nc)
mu ~ Normal(0, 1)
sigma ~ Uniform(0, 10)
delta ~ truncated(Normal(0, 1), lower=0)
alpha = delta * sigma
phi_c ~ filldist(Normal(mu + alpha / 2, sigma), length(kc))
phi_a ~ filldist(Normal(mu - alpha / 2, sigma), length(ka))
theta_c = phi.(phi_c)
theta_a = phi.(phi_a)
for i in eachindex(kc)
kc[i] ~ Binomial(nc[i], theta_c[i])
end
for i in eachindex(ka)
ka[i] ~ Binomial(na[i], theta_a[i])
end
return (theta_c=theta_c, theta_a=theta_a)
end
iterations = 10_000
burnin = 5_000
model_order_restricted_between_subject_proportions = OrderRestrictedBetweenSubjectProportions(ka, na, kc, nc)
chain = sample(model_order_restricted_between_subject_proportions, NUTS(), iterations, burnin=burnin)
chains_params = Turing.MCMCChains.get_sections(chain, :parameters)
quantities = generated_quantities(model_order_restricted_between_subject_proportions, chains_params)
chain
┌ Info: Found initial step size
└ ϵ = 0.4
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:13
Chains MCMC chain (10000×93×1 Array{Float64, 3}):
Iterations = 1001:1:11000
Number of chains = 1
Samples per chain = 10000
Wall duration = 18.97 seconds
Compute duration = 18.97 seconds
parameters = mu, sigma, delta, phi_c[1], phi_c[2], phi_c[3], phi_c[4], phi_c[5], phi_c[6], phi_c[7], phi_c[8], phi_c[9], phi_c[10], phi_c[11], phi_c[12], phi_c[13], phi_c[14], phi_c[15], phi_c[16], phi_c[17], phi_c[18], phi_c[19], phi_c[20], phi_c[21], phi_c[22], phi_c[23], phi_c[24], phi_c[25], phi_c[26], phi_a[1], phi_a[2], phi_a[3], phi_a[4], phi_a[5], phi_a[6], phi_a[7], phi_a[8], phi_a[9], phi_a[10], phi_a[11], phi_a[12], phi_a[13], phi_a[14], phi_a[15], phi_a[16], phi_a[17], phi_a[18], phi_a[19], phi_a[20], phi_a[21], phi_a[22], phi_a[23], phi_a[24], phi_a[25], phi_a[26], phi_a[27], phi_a[28], phi_a[29], phi_a[30], phi_a[31], phi_a[32], phi_a[33], phi_a[34], phi_a[35], phi_a[36], phi_a[37], phi_a[38], phi_a[39], phi_a[40], phi_a[41], phi_a[42], phi_a[43], phi_a[44], phi_a[45], phi_a[46], phi_a[47], phi_a[48], phi_a[49], phi_a[50], phi_a[51], phi_a[52]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
mu 0.4682 0.0557 0.0005 11536.0478 6881.3690 0.9999 ⋯
sigma 0.4562 0.0404 0.0004 11262.5260 7597.4049 1.0001 ⋯
delta 0.1742 0.1351 0.0013 7274.7236 4971.6464 1.0000 ⋯
phi_c[1] 0.9161 0.1482 0.0013 12273.8395 6511.4078 1.0000 ⋯
phi_c[2] 1.0250 0.1642 0.0014 12873.3587 7262.7948 1.0003 ⋯
phi_c[3] 0.0854 0.1080 0.0010 12546.7066 7292.6724 1.0004 ⋯
phi_c[4] 1.0720 0.1558 0.0014 12445.3199 6612.2790 1.0002 ⋯
phi_c[5] 0.1027 0.1077 0.0009 14572.5607 6735.8726 1.0002 ⋯
phi_c[6] 0.3598 0.1143 0.0010 12565.2028 6809.1214 1.0005 ⋯
phi_c[7] 0.0470 0.1070 0.0009 13121.1828 6407.6875 1.0002 ⋯
phi_c[8] -0.0806 0.1086 0.0009 13535.3184 7100.9844 1.0003 ⋯
phi_c[9] 0.1591 0.1082 0.0009 14465.0975 6905.0753 0.9999 ⋯
phi_c[10] 1.0087 0.1596 0.0014 12996.8110 6475.8020 1.0000 ⋯
phi_c[11] -0.0268 0.1082 0.0010 11720.9801 7166.7608 1.0002 ⋯
phi_c[12] 0.8226 0.1344 0.0013 11117.5222 7349.6512 1.0002 ⋯
phi_c[13] -0.4502 0.1117 0.0010 13504.1301 6893.9841 0.9999 ⋯
phi_c[14] 0.9995 0.1573 0.0013 14473.6403 6848.9579 1.0004 ⋯
phi_c[15] 0.0484 0.1070 0.0010 11009.8813 7570.0938 1.0000 ⋯
phi_c[16] 1.1489 0.1534 0.0014 11584.7912 7257.6038 1.0002 ⋯
phi_c[17] -0.6009 0.1170 0.0012 10109.0426 7258.8666 0.9999 ⋯
phi_c[18] 1.0190 0.1711 0.0016 11424.2508 6492.8260 1.0000 ⋯
phi_c[19] 0.5552 0.1134 0.0010 13397.7681 5993.3786 1.0000 ⋯
phi_c[20] 0.7462 0.1418 0.0013 12835.9461 7301.9182 1.0001 ⋯
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱
1 column and 58 rows omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
mu 0.3606 0.4303 0.4676 0.5062 0.5790
sigma 0.3846 0.4277 0.4535 0.4820 0.5425
delta 0.0071 0.0672 0.1450 0.2506 0.4996
phi_c[1] 0.6295 0.8151 0.9141 1.0145 1.2147
phi_c[2] 0.7068 0.9167 1.0207 1.1318 1.3602
phi_c[3] -0.1268 0.0118 0.0857 0.1577 0.2951
phi_c[4] 0.7736 0.9659 1.0712 1.1749 1.3883
phi_c[5] -0.1100 0.0316 0.1022 0.1740 0.3121
phi_c[6] 0.1367 0.2836 0.3597 0.4362 0.5825
phi_c[7] -0.1609 -0.0248 0.0471 0.1184 0.2566
phi_c[8] -0.2908 -0.1556 -0.0827 -0.0079 0.1315
phi_c[9] -0.0490 0.0856 0.1585 0.2307 0.3744
phi_c[10] 0.6930 0.9023 1.0050 1.1138 1.3328
phi_c[11] -0.2385 -0.0976 -0.0264 0.0447 0.1830
phi_c[12] 0.5621 0.7308 0.8213 0.9129 1.0927
phi_c[13] -0.6726 -0.5255 -0.4494 -0.3744 -0.2336
phi_c[14] 0.6929 0.8935 0.9976 1.1047 1.3124
phi_c[15] -0.1613 -0.0243 0.0486 0.1214 0.2555
phi_c[16] 0.8554 1.0436 1.1465 1.2518 1.4599
phi_c[17] -0.8293 -0.6793 -0.5998 -0.5227 -0.3745
phi_c[18] 0.6886 0.9039 1.0161 1.1303 1.3668
phi_c[19] 0.3388 0.4786 0.5551 0.6306 0.7835
phi_c[20] 0.4682 0.6517 0.7460 0.8416 1.0298
⋮ ⋮ ⋮ ⋮ ⋮ ⋮
58 rows omitted
function plot_order_restricted_between_subject_proportions(chain, prior_dist, zoom, size=(600, 400))
delta = chain[:delta] |> vec
plot(zoom, [pdf(prior_dist, e) for e in zoom], size=size, fmt=format, linewidth=2, label="Prior")
k = kde(delta)
plot!(zoom, pdf(k, zoom), linewidth=2, ls=:dash, label="Posterior")
plot!([0, 0], [pdf(prior_dist, 0), pdf(k, 0)], linewidth=2, label=false)
println(pdf(k, 0) / pdf(prior_dist, 0))
scatter!([0, 0], [pdf(prior_dist, 0), pdf(k, 0)], label=false)
bayes_factor = pdf(k, 0) / pdf(prior_dist, 0)
println("Bayes Factor $bayes_factor")
xlabel!("Delta")
ylabel!("Density")
xlims!(min(zoom...), max(zoom...))
end
plot_order_restricted_between_subject_proportions (generic function with 2 methods)
plot_order_restricted_between_subject_proportions(chain,
truncated(Normal(0, 1), lower=0), -0.1:0.01:2)
2.3933158091453777
Bayes Factor 2.3933158091453777
Pkg.status()
Status `~/quangtiencs_projects/bayesian-cognitive-modeling-with-turing.jl/Project.toml`
[336ed68f] CSV v0.10.11
[a93c6f00] DataFrames v1.6.1
[39dd38d3] Dierckx v0.5.3
⌃ [366bfd00] DynamicPPL v0.23.0
[1a297f60] FillArrays v1.5.0
[7073ff75] IJulia v1.24.2
[a98d9a8b] Interpolations v0.14.7
[5ab0869b] KernelDensity v0.6.7
[b964fa9f] LaTeXStrings v1.3.0
[2ab3a3ac] LogExpFunctions v0.3.24
[91a5bcdd] Plots v1.38.17
[37e2e3b7] ReverseDiff v1.15.0
[276daf66] SpecialFunctions v2.3.0
[2913bbd2] StatsBase v0.34.0
[f3b207a7] StatsPlots v0.15.6
[fce5fe82] Turing v0.28.1
[e88e6eb3] Zygote v0.6.62
Info Packages marked with ⌃ have new versions available and may be upgradable.