Parameter Estimation
This brief tutorial explains how to performance Bayesian parameter estimation of the QPDM using Pigeons.jl. One complication in estimating the parameters of the QPDM is that the posterior distributions may have multiple modes, which leads to convergence problems with most MCMC algorithms. Pigeons.jl uses a special type of parallel tempering to overcome this challenge. An additional advantage of using Pigeons.jl is the ability to compute Bayes factors from the log marginal likelihood using the function stepping_stone
.
Load Packages
First, we will load the required packages below.
using Pigeons
using QuantumPrisonersDilemmaModel
using Random
using StatsPlots
using Turing
Generate Simulated Data
The next step is to generate some simulated data from which the parameters can be estimated. In the code block below, the utility parameter $\mu_d$ is set to one and the entanglement parameter is set to $\gamma = 2$. A total of 50 trials is generated for each of the three conditions. The resulting values represent the number of defections per condition out of 50.
Random.seed!(65)
n = 50
parms = (μd=1.0, γ=2.0)
model = QPDM(;parms...)
data = rand(model, n)
Define Turing Model
The next step is to define a Turing model with the @model
macro. For simplicity, we will fix the utility parameter $\mu_d=1$ and set the prior of the entanglement parameter to $\gamma \sim \mathrm{normal}(0,3)$.
@model function turing_model(data, parms)
γ ~ Normal(0, 3)
data ~ QPDM(;parms..., γ)
end
Estimate Parameters
To estimate the parameters, we need to pass the Turing model to pigeons
. The second command converts the output to an MCMCChain
object, which can be used for plotting
pt = pigeons(
target=TuringLogPotential(turing_model((n, data), parms)),
record=[traces])
samples = Chains(sample_array(pt), ["γ","LL"])
The trace of the pigeon
's sampler is given below:
────────────────────────────────────────────────────────────────────────────
scans Λ log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 1.58 -11.3 0.332 0.825 1 1
4 0.529 -10.1 0.522 0.941 1 1
8 1.11 -9.43 0.501 0.877 1 1
16 1.37 -9.89 0.66 0.847 1 1
32 1.48 -10.3 0.772 0.836 1 1
64 1.46 -10.1 0.735 0.837 1 1
128 1.44 -10.4 0.776 0.84 1 1
256 1.49 -10.4 0.772 0.834 0.999 1
512 1.46 -10.3 0.816 0.838 0.999 1
1.02e+03 1.48 -10.3 0.817 0.836 0.999 1
────────────────────────────────────────────────────────────────────────────
Plot Posterior Distribution
Now we can plot the posterior distribution of $\gamma$ with plot
. The posterior distribution of $\gamma$ has a primary mode around 1 and secondary modes around 2 and 3.5.
plot(samples)