Run Bayesian t-tests on the posterior¶
In this quick tutorial we illustrate how to use posterior samples to perform Bayesian hypothesis testing.
In [1]:
Copied!
import arviz as az
import numpy as np
import pandas as pd
import hssm
import arviz as az
import numpy as np
import pandas as pd
import hssm
Example 1: Separate models¶
Simulate Data¶
We will simulate a simple dataset that contains two conditions.
In [2]:
Copied!
# Condition 1
condition_1 = hssm.simulate_data(
model="ddm", theta=dict(v=0.5, a=1.5, z=0.5, t=0.1), size=500
)
# Condition 2
condition_2 = hssm.simulate_data(
model="ddm", theta=dict(v=1.0, a=1.5, z=0.5, t=0.1), size=500
)
# Condition 1
condition_1 = hssm.simulate_data(
model="ddm", theta=dict(v=0.5, a=1.5, z=0.5, t=0.1), size=500
)
# Condition 2
condition_2 = hssm.simulate_data(
model="ddm", theta=dict(v=1.0, a=1.5, z=0.5, t=0.1), size=500
)
Specify Models¶
We will fit two separate models to the data.
In [3]:
Copied!
# Model 1
m1 = hssm.HSSM(model="ddm", data=condition_1)
m1.sample(sampler="pymc", tune=500, draws=500)
# Model 2
m2 = hssm.HSSM(model="ddm", data=condition_2)
m2.sample(sampler="pymc", tune=500, draws=500)
# Model 1
m1 = hssm.HSSM(model="ddm", data=condition_1)
m1.sample(sampler="pymc", tune=500, draws=500)
# Model 2
m2 = hssm.HSSM(model="ddm", data=condition_2)
m2.sample(sampler="pymc", tune=500, draws=500)
Model initialized successfully.
Using default initvals.
Initializing NUTS using adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [z, a, t, v]
/Users/afengler/.local/share/uv/python/cpython-3.12.13-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock. self.pid = os.fork()
/Users/afengler/.local/share/uv/python/cpython-3.12.13-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock. self.pid = os.fork()
Sampling 4 chains for 500 tune and 500 draw iterations (2_000 + 2_000 draws total) took 7 seconds.
Model initialized successfully.
Using default initvals.
Initializing NUTS using adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [z, a, t, v]
/Users/afengler/.local/share/uv/python/cpython-3.12.13-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock. self.pid = os.fork()
/Users/afengler/.local/share/uv/python/cpython-3.12.13-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock. self.pid = os.fork()
Sampling 4 chains for 500 tune and 500 draw iterations (2_000 + 2_000 draws total) took 4 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
Out[3]:
<xarray.DataTree>
Group: /
├── Group: /posterior
│ Dimensions: (chain: 4, draw: 500)
│ Coordinates:
│ * chain (chain) int64 32B 0 1 2 3
│ * draw (draw) int64 4kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
│ Data variables:
│ z (chain, draw) float64 16kB 0.557 0.5482 0.5156 ... 0.6026 0.5516
│ v (chain, draw) float64 16kB 0.8854 0.9206 1.064 ... 0.7973 0.8976
│ a (chain, draw) float64 16kB 1.6 1.628 1.59 ... 1.631 1.582 1.462
│ t (chain, draw) float64 16kB 0.1156 0.0987 0.1026 ... 0.1604 0.1575
│ Attributes:
│ created_at: 2026-07-08T22:53:36.947710+00:00
│ creation_library: ArviZ
│ creation_library_version: 1.2.0
│ creation_library_language: Python
│ inference_library: pymc
│ inference_library_version: 6.1.0
│ sample_dims: ['chain', 'draw']
│ sampling_time: 3.561091899871826
│ tuning_steps: 500
│ modeling_interface: bambi
│ modeling_interface_version: 0.18.0
├── Group: /sample_stats
│ Dimensions: (chain: 4, draw: 500)
│ Coordinates:
│ * chain (chain) int64 32B 0 1 2 3
│ * draw (draw) int64 4kB 0 1 2 3 4 5 ... 495 496 497 498 499
│ Data variables: (12/18)
│ n_steps (chain, draw) float64 16kB 7.0 15.0 ... 15.0 15.0
│ smallest_eigval (chain, draw) float64 16kB nan nan nan ... nan nan
│ energy_error (chain, draw) float64 16kB -0.5117 ... 0.01835
│ perf_counter_start (chain, draw) float64 16kB 6.29e+05 ... 6.29e+05
│ energy (chain, draw) float64 16kB 706.4 707.2 ... 710.1
│ step_size (chain, draw) float64 16kB 0.2918 0.2918 ... 0.2762
│ ... ...
│ reached_max_treedepth (chain, draw) bool 2kB False False ... False False
│ index_in_trajectory (chain, draw) int64 16kB 3 -11 -10 2 ... -3 -3 5 -4
│ divergences (chain, draw) int64 16kB 0 0 0 0 0 0 ... 0 0 0 0 0 0
│ max_energy_error (chain, draw) float64 16kB -0.9361 0.1701 ... -0.1288
│ diverging (chain, draw) bool 2kB False False ... False False
│ perf_counter_diff (chain, draw) float64 16kB 0.001463 ... 0.002735
│ Attributes:
│ created_at: 2026-07-08T22:53:36.951735+00:00
│ creation_library: ArviZ
│ creation_library_version: 1.2.0
│ creation_library_language: Python
│ inference_library: pymc
│ inference_library_version: 6.1.0
│ sample_dims: ['chain', 'draw']
│ sampling_time: 3.561091899871826
│ tuning_steps: 500
│ modeling_interface: bambi
│ modeling_interface_version: 0.18.0
├── Group: /observed_data
│ Dimensions: (__obs__: 500, rt,response_extra_dim_0: 2)
│ Coordinates:
│ * __obs__ (__obs__) int64 4kB 0 1 2 3 4 ... 496 497 498 499
│ * rt,response_extra_dim_0 (rt,response_extra_dim_0) int64 16B 0 1
│ Data variables:
│ rt,response (__obs__, rt,response_extra_dim_0) float64 8kB 0...
│ Attributes:
│ created_at: 2026-07-08T22:53:36.953094+00:00
│ creation_library: ArviZ
│ creation_library_version: 1.2.0
│ creation_library_language: Python
│ inference_library: pymc
│ inference_library_version: 6.1.0
│ sample_dims: []
│ modeling_interface: bambi
│ modeling_interface_version: 0.18.0
└── Group: /log_likelihood
Dimensions: (chain: 4, draw: 500, __obs__: 500)
Coordinates:
* chain (chain) int64 32B 0 1 2 3
* draw (draw) int64 4kB 0 1 2 3 4 5 6 ... 493 494 495 496 497 498 499
* __obs__ (__obs__) int64 4kB 0 1 2 3 4 5 6 ... 494 495 496 497 498 499
Data variables:
rt,response (chain, draw, __obs__) float64 8MB -0.5061 -2.22 ... -0.4034
Attributes:
modeling_interface: bambi
modeling_interface_version: 0.18.0Bayesian t-test¶
Now, let's ask the (informed) question:
Is the v parameter higher in condition 1 than condition 2?
The beauty of the Bayesian approach is that this questions boils down to simple counting.
Reformulating the question: Are the v samples of my posterior for condition 1 higher than the samples of my posterior for condition 2?
Let's check for that, we have everything we need in our posterior samples!
In [4]:
Copied!
# Percent of posterior samples
print(
"p(v_1 > v_2) = "
f"{np.mean(m1.traces.posterior['v'].values > m2.traces.posterior['v'].values)}"
)
# Percent of posterior samples
print(
"p(v_1 > v_2) = "
f"{np.mean(m1.traces.posterior['v'].values > m2.traces.posterior['v'].values)}"
)
p(v_1 > v_2) = 0.0
Looks like our inference indicates that there is a $0\%$ chance for v parameter of condition 1 to be higher than the v parameter of condition 2.
Plotting¶
We can also plot the posterior samples to get a sense of the uncertainty.
In [5]:
Copied!
# Overlay the two conditions' posteriors for v on a single axes
pc = az.plot_dist(
m1.traces, var_names=["v"], kind="hist", visuals={"dist": dict(color="C0")}
)
az.plot_dist(
m2.traces,
var_names=["v"],
kind="hist",
plot_collection=pc,
visuals={"dist": dict(color="C1")},
)
# Retrieve the underlying matplotlib axes to add labels and a legend
ax = pc.get_target("v", {})
ax.set_xlim(0.25, 1.25)
ax.set_title("Posterior samples for v")
ax.set_xlabel("v")
from matplotlib.lines import Line2D
ax.legend(
handles=[
Line2D([0], [0], color="C0", label="Condition 1"),
Line2D([0], [0], color="C1", label="Condition 2"),
]
)
# Overlay the two conditions' posteriors for v on a single axes
pc = az.plot_dist(
m1.traces, var_names=["v"], kind="hist", visuals={"dist": dict(color="C0")}
)
az.plot_dist(
m2.traces,
var_names=["v"],
kind="hist",
plot_collection=pc,
visuals={"dist": dict(color="C1")},
)
# Retrieve the underlying matplotlib axes to add labels and a legend
ax = pc.get_target("v", {})
ax.set_xlim(0.25, 1.25)
ax.set_title("Posterior samples for v")
ax.set_xlabel("v")
from matplotlib.lines import Line2D
ax.legend(
handles=[
Line2D([0], [0], color="C0", label="Condition 1"),
Line2D([0], [0], color="C1", label="Condition 2"),
]
)
Out[5]:
<matplotlib.legend.Legend at 0x135dacd40>
A glance a these posteriors visually corroborates our simple counting analysis. The two posteriors for v are clearly separated.
To appreciate the difference, let us also plot the posteriors of the respective a parameters.
In [6]:
Copied!
# Overlay the two conditions' posteriors for a on a single axes
pc = az.plot_dist(
m1.traces, var_names=["a"], kind="hist", visuals={"dist": dict(color="C0")}
)
az.plot_dist(
m2.traces,
var_names=["a"],
kind="hist",
plot_collection=pc,
visuals={"dist": dict(color="C1")},
)
# Retrieve the underlying matplotlib axes to add labels and a legend
ax = pc.get_target("a", {})
ax.set_xlim(1.3, 1.7)
ax.set_title("Posterior samples for a")
ax.set_xlabel("a")
from matplotlib.lines import Line2D
ax.legend(
handles=[
Line2D([0], [0], color="C0", label="Condition 1"),
Line2D([0], [0], color="C1", label="Condition 2"),
]
)
# Overlay the two conditions' posteriors for a on a single axes
pc = az.plot_dist(
m1.traces, var_names=["a"], kind="hist", visuals={"dist": dict(color="C0")}
)
az.plot_dist(
m2.traces,
var_names=["a"],
kind="hist",
plot_collection=pc,
visuals={"dist": dict(color="C1")},
)
# Retrieve the underlying matplotlib axes to add labels and a legend
ax = pc.get_target("a", {})
ax.set_xlim(1.3, 1.7)
ax.set_title("Posterior samples for a")
ax.set_xlabel("a")
from matplotlib.lines import Line2D
ax.legend(
handles=[
Line2D([0], [0], color="C0", label="Condition 1"),
Line2D([0], [0], color="C1", label="Condition 2"),
]
)
Out[6]:
<matplotlib.legend.Legend at 0x135daef90>
and correspondingly,
In [7]:
Copied!
print(
"p(a_1 > a_2) = "
f"{np.mean(m1.traces.posterior['a'].values > m2.traces.posterior['a'].values)}"
)
print(
"p(a_1 > a_2) = "
f"{np.mean(m1.traces.posterior['a'].values > m2.traces.posterior['a'].values)}"
)
p(a_1 > a_2) = 0.0735
Example 2: Combined Model¶
In [8]:
Copied!
condition_1["condition"] = "C1"
condition_2["condition"] = "C2"
data = pd.concat([condition_1, condition_2]).reset_index(drop=True)
condition_1["condition"] = "C1"
condition_2["condition"] = "C2"
data = pd.concat([condition_1, condition_2]).reset_index(drop=True)
In [9]:
Copied!
m_combined = hssm.HSSM(
model="ddm",
data=data,
include=[
{
"name": "v",
"formula": "v ~ 1 + condition",
}
],
)
idata_combined = m_combined.sample(sampler="pymc", tune=500, draws=500)
m_combined = hssm.HSSM(
model="ddm",
data=data,
include=[
{
"name": "v",
"formula": "v ~ 1 + condition",
}
],
)
idata_combined = m_combined.sample(sampler="pymc", tune=500, draws=500)
Model initialized successfully.
Using default initvals.
Initializing NUTS using adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [z, a, t, v_Intercept, v_condition]
/Users/afengler/.local/share/uv/python/cpython-3.12.13-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock. self.pid = os.fork()
/Users/afengler/.local/share/uv/python/cpython-3.12.13-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock. self.pid = os.fork()
Sampling 4 chains for 500 tune and 500 draw iterations (2_000 + 2_000 draws total) took 8 seconds.
Note, now we don't have two distinct idata objects, but instead we have a single idata object with a posterior for v_condition. Let's take a closer look.
In [10]:
Copied!
m_combined.traces.posterior
m_combined.traces.posterior
Out[10]:
<xarray.DataTree 'posterior'>
Group: /posterior
Dimensions: (chain: 4, draw: 500, v_condition_dim: 1)
Coordinates:
* chain (chain) int64 32B 0 1 2 3
* draw (draw) int64 4kB 0 1 2 3 4 5 6 ... 494 495 496 497 498 499
* v_condition_dim (v_condition_dim) <U2 8B 'C2'
Data variables:
z (chain, draw) float64 16kB 0.5305 0.5075 ... 0.5127 0.4787
t (chain, draw) float64 16kB 0.09818 0.1143 ... 0.1078
v_Intercept (chain, draw) float64 16kB 0.4813 0.478 ... 0.4181 0.5055
v_condition (chain, draw, v_condition_dim) float64 16kB 0.5644 ... 0...
a (chain, draw) float64 16kB 1.616 1.489 1.532 ... 1.54 1.473
Attributes:
created_at: 2026-07-08T22:53:46.244556+00:00
creation_library: ArviZ
creation_library_version: 1.2.0
creation_library_language: Python
inference_library: pymc
inference_library_version: 6.1.0
sample_dims: ['chain', 'draw']
sampling_time: 7.684842824935913
tuning_steps: 500
modeling_interface: bambi
modeling_interface_version: 0.18.0Under the hood, Bambi created a model with a dummy variable for the C2 condition (the C1 condition represents the Intercept).
So what do we need to test here...
We still ask the same question: Is the v_condition[C2] posterior higher than the v_condition[C1] posterior?
But now we only need to check if our v_condition[C2] variable is above 0, since it represents the offset from the Intercept directly.
Let's check:
In [11]:
Copied!
print(
"p(v_condition[C2] > v_condition[C1]) = "
f"{np.mean(idata_combined.posterior['v_condition'].values > 0)}"
)
print(
"p(v_condition[C2] > v_condition[C1]) = "
f"{np.mean(idata_combined.posterior['v_condition'].values > 0)}"
)
p(v_condition[C2] > v_condition[C1]) = 1.0
or visually,
In [12]:
Copied!
az.plot_dist(idata_combined, var_names=["v_condition"], kind="hist")
az.plot_dist(idata_combined, var_names=["v_condition"], kind="hist")
Out[12]:
<arviz_plots.plot_collection.PlotCollection at 0x12e8afce0>
You can use this approach to test any number of complex statements about your parameters. There will essentially always be a way to turn your question into a simple comparison of the posterior samples, a simple counting problem.