import contextlib
import io
import logging
import os
import warnings
logging.getLogger("pytensor").setLevel(logging.ERROR)
logging.getLogger("jax._src.xla_bridge").setLevel(logging.ERROR)
logging.getLogger("pymc").setLevel(logging.ERROR)
logging.getLogger("hssm").setLevel(logging.ERROR)
import arviz as az
import marimo as mo
import numpy as np
from matplotlib import pyplot as plt
import hssm
FULL_RUN = os.environ.get("FULL_RUN", "0") == "1"
N_TRIALS = 500 if FULL_RUN else 200
N_TUNE = 750 if FULL_RUN else 100
N_DRAWS = 750 if FULL_RUN else 100
N_CHAINS = 2 if FULL_RUN else 1
N_PPC_DRAWS = 100 if FULL_RUN else 20
RANDOM_SEED = 20260830
@contextlib.contextmanager
def quiet_console():
"""Hide progress noise while replaying unique, path-free warnings."""
previous_disable_level = logging.root.manager.disable
logging.disable(logging.CRITICAL)
with warnings.catch_warnings(record=True) as caught_warnings:
warnings.simplefilter("always")
try:
with (
contextlib.redirect_stdout(io.StringIO()),
contextlib.redirect_stderr(io.StringIO()),
):
yield
finally:
logging.disable(previous_disable_level)
for message in dict.fromkeys(str(item.message) for item in caught_warnings):
print(f"Warning: {message}")
The HSSM tutorial¶
This guided tutorial starts where the Quickstart ends. Instead of fitting another intercept-only DDM, you will answer one common research question: does an experimental condition change drift rate?
We take one recommended route from data to interpretation:
simulate -> specify one regression -> sample -> diagnose -> check predictions -> interpret
You will fit one model and make one decision at each step. The optional scenic route covers alternative model families, custom priors, hierarchical variants, model comparison, and low-level extensions.
What you will learn¶
By the end you will be able to:
- express a trial-level predictor with HSSM's formula interface;
- verify that the intended parameter and predictor entered the model;
- read a compact ArviZ summary and basic chain diagnostic;
- check whether the fitted model can reproduce the observed data; and
- report the condition effect with posterior uncertainty.
For installation and the mechanics of a first fit, complete the Quickstart first. This page deliberately does not repeat those steps.
Run this tutorial¶
On Colab, install HSSM with %pip install hssm, then restart the runtime.
The published outputs use the full configuration. Routine notebook CI uses
a smaller configuration to exercise the same workflow quickly.
{
"mode": "full (published outputs)" if FULL_RUN else "quick (CI smoke check)",
"artifact_marker": (
"<!-- hssm-full-run-artifact: true -->"
if FULL_RUN
else "<!-- hssm-full-run-artifact: false -->"
),
"deterministic_init_marker": "<!-- hssm-deterministic-init: true -->",
"trials": N_TRIALS,
"chains": N_CHAINS,
"tune": N_TUNE,
"draws": N_DRAWS,
"posterior_predictive_draws": N_PPC_DRAWS,
}
{'mode': 'full (published outputs)',
'artifact_marker': '<!-- hssm-full-run-artifact: true -->',
'deterministic_init_marker': '<!-- hssm-deterministic-init: true -->',
'trials': 500,
'chains': 2,
'tune': 750,
'draws': 750,
'posterior_predictive_draws': 100}
1. Simulate the question, not just the model¶
We create two equally likely conditions. The simulated drift rate is
0.3 in the reference condition and 1.0 in the treatment condition, so
the known treatment effect is 0.7. Boundary separation, starting point,
and non-decision time stay constant.
In an applied analysis, condition would come from your experimental data;
only this simulation step would change.
rng = np.random.default_rng(RANDOM_SEED)
condition = rng.integers(0, 2, size=N_TRIALS)
true_values = {
"v_Intercept": 0.3,
"v_condition": 0.7,
"a": 1.5,
"z": 0.5,
"t": 0.25,
}
trial_v = true_values["v_Intercept"] + true_values["v_condition"] * condition
data = hssm.simulate_data(
model="ddm",
theta={
"v": trial_v,
"a": true_values["a"],
"z": true_values["z"],
"t": true_values["t"],
},
size=1,
random_state=RANDOM_SEED,
)
data["condition"] = condition
data.head()
| rt | response | condition | |
|---|---|---|---|
| 0 | 2.668714 | 1.0 | 0 |
| 1 | 5.240685 | 1.0 | 0 |
| 2 | 2.376304 | 1.0 | 1 |
| 3 | 1.463552 | 1.0 | 0 |
| 4 | 2.870995 | -1.0 | 0 |
HSSM requires rt and response; predictor columns sit beside them. Here
condition=0 is the reference level, so the intercept is its drift rate
and the condition coefficient is the change from 0 to 1.
2. Specify one parameter regression¶
We let only drift rate vary with condition. The other DDM parameters use HSSM's defaults. This is a good first specification when the scientific hypothesis concerns evidence quality rather than response caution, bias, or non-decision processes.
model = hssm.HSSM(
data=data,
model="ddm",
initval_jitter=0,
include=[
{
"name": "v",
"formula": "v ~ 1 + condition",
"link": "identity",
}
],
)
print(model)
Model initialized successfully.
Hierarchical Sequential Sampling Model
Model: ddm
Response variable: rt,response
Likelihood: analytical
Observations: 500
Parameters:
v:
Formula: v ~ 1 + condition
Priors:
v_Intercept ~ Normal(mu: 2.0, sigma: 3.0)
v_condition ~ Normal(mu: 0.0, sigma: 0.25)
Link: identity
Explicit bounds: (-inf, inf)
a:
Prior: HalfNormal(sigma: 2.0)
Explicit bounds: (0.0, inf)
z:
Prior: Uniform(lower: 0.0, upper: 1.0)
Explicit bounds: (0.0, 1.0)
t:
Prior: HalfNormal(sigma: 2.0)
Explicit bounds: (0.0, inf)
Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)
Read the model summary before sampling. It should show v_Intercept and
v_condition under drift rate, while a, z, and t remain ordinary
parameters. This check catches misspelled predictors and unintended model
structure before computation begins.
For coefficient-level priors, use Specify priors and fix parameters.
3. Sample once, then inspect a compact result¶
The full run uses two chains and enough draws for a tutorial-quality diagnostic view. The quick run is only an execution smoke test; one-chain diagnostics are not evidence of convergence.
with quiet_console():
idata = model.sample(
sampler="pymc",
chains=N_CHAINS,
cores=1,
draws=N_DRAWS,
tune=N_TUNE,
random_seed=RANDOM_SEED,
idata_kwargs={"log_likelihood": False},
progressbar=False,
)
{
"groups": tuple(idata.children),
"posterior_sizes": dict(idata.posterior.ds.sizes),
}
{'groups': ('posterior', 'sample_stats', 'observed_data'),
'posterior_sizes': {'chain': 2, 'draw': 750}}
parameter_names = ["v_Intercept", "v_condition", "a", "z", "t"]
posterior_summary = az.summary(
idata,
var_names=parameter_names,
kind="all" if FULL_RUN else "stats",
round_to=2,
)
posterior_summary
| mean | sd | eti89_lb | eti89_ub | ess_bulk | ess_tail | r_hat | mcse_mean | mcse_sd | |
|---|---|---|---|---|---|---|---|---|---|
| v_Intercept | 0.25 | 0.05 | 0.16 | 0.34 | 1539.09 | 1199.58 | 1.0 | 0.0 | 0.0 |
| v_condition | 0.67 | 0.08 | 0.55 | 0.79 | 1336.98 | 1186.11 | 1.0 | 0.0 | 0.0 |
| a | 1.48 | 0.04 | 1.42 | 1.55 | 1213.40 | 1171.91 | 1.0 | 0.0 | 0.0 |
| z | 0.53 | 0.02 | 0.50 | 0.56 | 1113.56 | 1074.06 | 1.0 | 0.0 | 0.0 |
| t | 0.28 | 0.03 | 0.24 | 0.32 | 1229.29 | 1121.21 | 1.0 | 0.0 | 0.0 |
divergences = int(idata.sample_stats.diverging.values.sum())
if FULL_RUN:
diagnostic_summary = az.summary(
idata,
var_names=parameter_names,
kind="diagnostics",
round_to="none",
)
max_rhat = float(diagnostic_summary["r_hat"].max())
min_ess = float(diagnostic_summary[["ess_bulk", "ess_tail"]].min().min())
assert divergences == 0, f"full run has {divergences} divergences"
assert max_rhat <= 1.01, f"full-run max R-hat is {max_rhat:.4f}"
assert min_ess >= 200, f"full-run minimum ESS is {min_ess:.1f}"
health = {
"scope": "full diagnostic validation",
"divergences": divergences,
"max_rhat": round(max_rhat, 4),
"min_bulk_or_tail_ess": round(min_ess, 1),
}
else:
health = {
"scope": "quick execution/specification smoke check only",
"divergences_reported_not_gated": divergences,
"rhat_and_ess": "not evaluated with one chain",
}
health
{'scope': 'full diagnostic validation',
'divergences': 0,
'max_rhat': 1.0017,
'min_bulk_or_tail_ess': 1074.1}
Start with the coefficient row: its mean and interval describe the
treatment-minus-reference change in drift. Then check r_hat (near 1),
effective sample sizes, and the trace view below. In real work, increase
draws and investigate warnings before interpretation.
See Plot posteriors and predictions for the broader diagnostics and plotting toolkit.
trace_plot = az.plot_trace_dist(idata, var_names=["v_Intercept", "v_condition"])
trace_figure = trace_plot.get_viz("figure")
plt.close(trace_figure)
trace_figure
4. Check predictions before interpreting the coefficient¶
Diagnostics tell us whether the sampler explored this model; they do not tell us whether the model reproduces the observed choices and response times. A posterior predictive check asks exactly that question.
with quiet_console():
ppc_idata = model.sample_posterior_predictive(
dt=idata,
draws=N_PPC_DRAWS,
inplace=False,
)
Warning: Numba will use object mode to run ddm_RV_rv{"(),(),(),(),()->(2)"}'s perform method. Set `pytensor.config.compiler_verbose = True` to see more details.
predictive_grid = model.plot_predictive(dt=ppc_idata, col="condition")
predictive_figure = predictive_grid.figure
# The condition facet labels carry the useful comparison; the default
# global title overlaps them in the static docs layout.
predictive_figure.suptitle("")
predictive_figure.subplots_adjust(top=0.88)
plt.close(predictive_figure)
predictive_figure
Look for broad agreement in response proportions and response-time shape.
A visible mismatch is a reason to revise the model before telling a story
about v_condition.
5. Report the condition effect¶
We summarize the posterior mean, a 94% highest-density interval, and the posterior probability that the treatment effect is positive. These are direct descriptions of uncertainty, not a binary significance test.
effect_draws = az.extract(idata, var_names=["v_condition"]).values
effect_hdi = az.hdi(effect_draws, prob=0.94)
probability_positive = float((effect_draws > 0).mean())
if FULL_RUN:
assert effect_hdi[0] <= true_values["v_condition"] <= effect_hdi[1], (
"full-run condition-effect HDI misses the known simulated value"
)
assert probability_positive >= 0.95, (
"full-run posterior does not clearly support a positive condition effect"
)
effect_report = {
"known_simulated_effect": true_values["v_condition"],
"posterior_mean": round(float(effect_draws.mean()), 3),
"94%_HDI": tuple(np.round(effect_hdi, 3)),
"P(v_condition > 0)": round(probability_positive, 3),
}
effect_report
{'known_simulated_effect': 0.7,
'posterior_mean': 0.675,
'94%_HDI': (np.float64(0.534), np.float64(0.815)),
'P(v_condition > 0)': 1.0}
Where to go next¶
You now have one reusable HSSM workflow: put predictors in the data, attach a formula to the cognitive parameter named by your hypothesis, sample once, diagnose, check predictions, and report posterior uncertainty.
Choose the next page for the question you actually have:
- add participants with Hierarchical modeling;
- refine assumptions with Specify priors and fix parameters;
- compare pre-specified candidates with Compare and interpret models;
- see an applied end-to-end analysis in A complete scientific workflow; or
- explore every alternative in The HSSM tutorial: scenic route.