Custom models from ONNX files (blackbox)¶
Build an HSSM model directly from an onnx file, via the blackbox route. For our purposes, the onnx file-format provides nice translation layer from deep learning frameworks into a common layer from which we can then reconstruct computation graph to use through PyMC.
New to HSSM's likelihood kinds? Understanding likelihoods in HSSM explains how analytical,
approx_differentiable, and blackbox likelihoods differ.
import os
import matplotlib.pyplot as plt
import numpy as np
import hssm
Loading the network¶
# Networks
network_path = os.path.join("data", "race_3_no_bias_lan_no_batch.onnx")
The network we load here does not have dynamic input dimensions, which prevents us from batching computations.
Instead of fixing things behind the scenes and loading a fixed network, we provide a useful snippet below that shows how to rectify this situation.
⚠️ Blackbox route only. The dynamic-axes rewrite below applies only to
loglik_kind="blackbox", where the ONNX file runs inside a plain Python function viaonnxruntimeand HSSM never converts the graph. HSSM'sapprox_differentiableloader requires the exact opposite — every input dimension concrete — and will reject a file rewritten this way. See The ONNX likelihood contract.
import onnx
import onnxruntime as ort
# Load model from path
onnx_model = onnx.load(network_path)
# Change input and output dimensions to be dynamic to allow for batching
# (in case this is not already done)
for input_tensor in onnx_model.graph.input:
dim_proto = input_tensor.type.tensor_type.shape.dim[0]
if not dim_proto.dim_param == "None":
dim_proto.dim_param = "None"
for output_tensor in onnx_model.graph.output:
dim_proto = output_tensor.type.tensor_type.shape.dim[0]
if not dim_proto.dim_param == "None":
dim_proto.dim_param = "None"
input_name = onnx_model.graph.input[0].name
# Please uncomment the below line to save the adjusted model
# onnx.save(onnx_model, "data/race_3_no_bias_lan_batch.onnx")
2026-08-28 01:18:36.464080706 [W:onnxruntime:Default, device_discovery.cc:146 GetPciBusId] Skipping pci_bus_id for PCI path at "/sys/devices/LNXSYSTM:00/LNXSYBUS:00/ACPI0004:00/VMBUS:00/5620e0c7-8062-4dce-aeb7-520c7ef76171" because filename "5620e0c7-8062-4dce-aeb7-520c7ef76171" did not match expected pattern of [0-9a-f]+:[0-9a-f]+:[0-9a-f]+[.][0-9a-f]+
Armed with the corrected network, let's test inference speed on a data-batch of $1000$ trials.
# Load model batch ready model
ort_session = ort.InferenceSession("data/race_3_no_bias_lan_batch.onnx")
# Test inference speed
import time
start = time.time()
for i in range(100):
ort_session.run(
None, {input_name: np.random.uniform(size=(1000, 8)).astype(np.float32)}
)
end = time.time()
print(f"Time taken: {(end - start) / 100} seconds")
Time taken: 0.0006541347503662109 seconds
Defining the Likelihood¶
The network we loaded corresponds to a LAN, for a Race model with three choice alternatives.
This model has three drift parameters v0, v1, v2, a boundary parameter a, a starting point bias z and a non-decision-time t.
Data from this model has the usual rt, choice format.
We use this to construct a simple blackbox likelihood function below. This likelihood function takes the respective data and model parameters as arguments.
The function body shapes these input arguments into a matrix and performs a batched forward pass through the loaded network via the onnx.runtime.
def my_blackbox_race_model(data, v0, v1, v2, a, z, t):
"""Calculate log-likelihood for a 3-choice race model.
Parameters
----------
data : np.ndarray
Array of shape (n_trials, 2) containing response times in first column
and choices (0, 1, or 2) in second column
v0 : float
Drift rate for accumulator 0
v1 : float
Drift rate for accumulator 1
v2 : float
Drift rate for accumulator 2
a : float
Decision threshold/boundary
z : float
Starting point bias
t : float
Non-decision time
Notes
-----
HSSM calls blackbox likelihoods positionally, in ``list_params`` order —
keep the signature and the stacked column order aligned with
``list_params`` (here: v0, v1, v2, a, z, t), which is also the order the
network was trained on.
Returns
-------
np.ndarray
Array of log-likelihood values for each trial
"""
data_nrows = data.shape[0]
data = np.vstack(
[np.full(data_nrows, param_) for param_ in [v0, v1, v2, a, z, t]]
+ [data[:, 0], data[:, 1]]
).T.astype(np.float32)
return ort_session.run(None, {input_name: data})[0].squeeze()
Simulate example data¶
# Set parameters
v0 = 1.0
v1 = 0.5
v2 = 0.25
a = 1.5
t = 0.3
z = 0.5
# simulate some data from the model
obs_race3 = hssm.simulate_data(
theta=dict(v0=v0, v1=v1, v2=v2, a=a, t=t, z=z), model="race_no_bias_3", size=1000
)
Test Likelihood Outputs¶
# Test that outputs are reasonable
for choice in [0, 1, 2]:
rts = np.linspace(0, 20, 1000)
choices = np.repeat(choice, 1000)
data = np.vstack([rts, choices]).T
out = my_blackbox_race_model(data, v0, v1, v2, a, t, z)
plt.plot(rts, np.exp(out), label=f"choice: {choice}")
plt.legend()
plt.show()
Build HSSM Model¶
We can now build a simple HSSM model that takes in our new blackbox likelihood.
model = hssm.HSSM(
data=obs_race3,
model="race_no_bias_3", # some name for the model
model_config={
"list_params": ["v0", "v1", "v2", "a", "z", "t"],
"bounds": {
"v0": (0.0, 2.5),
"v1": (0.0, 2.5),
"v2": (0.0, 2.5),
"a": (1.0, 3.0),
"z": (0.0, 0.9),
"t": (0.001, 2),
},
}, # minimal specification of model parameters and parameter bounds
loglik_kind="blackbox", # use the blackbox loglik
loglik=my_blackbox_race_model,
choices=[0, 1, 2], # list the legal choice options
z=0.5,
p_outlier=0,
)
Model initialized successfully.
model.graph()
model.sample(draws=500, tune=200, discard_tuned_samples=False)
Using default initvals.
Multiprocess sampling (2 chains in 2 jobs)
CompoundStep
>Slice: [v2]
>Slice: [v1]
>Slice: [a]
>Slice: [t]
>Slice: [v0]
<python-install>/lib/python3.13/multiprocessing/popen_fork.py:73: 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()
<python-install>/lib/python3.13/multiprocessing/popen_fork.py:73: 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 2 chains for 200 tune and 500 draw iterations (400 + 1_000 draws total) took 32 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
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
<xarray.DataTree>
Group: /
├── Group: /warmup_posterior
│ Dimensions: (chain: 2, draw: 200, v0_mean_dim_0: 1)
│ Coordinates:
│ * chain (chain) int64 16B 0 1
│ * draw (draw) int64 2kB 0 1 2 3 4 5 6 ... 194 195 196 197 198 199
│ * v0_mean_dim_0 (v0_mean_dim_0) int64 8B 0
│ Data variables:
│ v2 (chain, draw) float64 3kB 3.345e-41 6.858e-21 ... 0.2001
│ v1 (chain, draw) float64 3kB 5.591e-116 2.514e-24 ... 0.4933
│ a (chain, draw) float64 3kB 3.0 3.0 3.0 ... 1.419 1.414 1.414
│ t (chain, draw) float64 3kB 0.001002 0.00137 ... 0.2933 0.2947
│ v0 (chain, draw) float64 3kB 1.451 1.65 1.578 ... 0.9951 1.003
│ v0_mean (chain, draw, v0_mean_dim_0) float64 3kB 1.451 1.65 ... 1.003
│ Attributes:
│ created_at: 2026-08-28T01:19:18.195355+00:00
│ creation_library: ArviZ
│ creation_library_version: 1.3.0
│ creation_library_language: Python
│ inference_library: pymc
│ inference_library_version: 6.3.1
│ sample_dims: ['chain', 'draw']
│ sampling_time: 32.14871883392334
│ tuning_steps: 200
│ modeling_interface: bambi
│ modeling_interface_version: 0.20.0
├── Group: /posterior
│ Dimensions: (chain: 2, draw: 500)
│ Coordinates:
│ * chain (chain) int64 16B 0 1
│ * draw (draw) int64 4kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
│ Data variables:
│ a (chain, draw) float64 8kB 1.327 1.342 1.331 ... 1.432 1.396 1.391
│ v0 (chain, draw) float64 8kB 0.8268 0.8183 0.7995 ... 0.9167 0.9316
│ t (chain, draw) float64 8kB 0.3049 0.3049 0.3012 ... 0.2988 0.3012
│ v2 (chain, draw) float64 8kB 0.1468 0.2687 0.1719 ... 0.3808 0.3039
│ v1 (chain, draw) float64 8kB 0.3688 0.3202 0.3361 ... 0.586 0.4373
│ Attributes:
│ created_at: 2026-08-28T01:19:18.199772+00:00
│ creation_library: ArviZ
│ creation_library_version: 1.3.0
│ creation_library_language: Python
│ inference_library: pymc
│ inference_library_version: 6.3.1
│ sample_dims: ['chain', 'draw']
│ sampling_time: 32.14871883392334
│ tuning_steps: 200
│ modeling_interface: bambi
│ modeling_interface_version: 0.20.0
├── Group: /warmup_sample_stats
│ Dimensions: (chain: 2, draw: 200, nstep_in_dim_0: 5, nstep_out_dim_0: 5)
│ Coordinates:
│ * chain (chain) int64 16B 0 1
│ * draw (draw) int64 2kB 0 1 2 3 4 5 6 ... 194 195 196 197 198 199
│ * nstep_in_dim_0 (nstep_in_dim_0) int64 40B 0 1 2 3 4
│ * nstep_out_dim_0 (nstep_out_dim_0) int64 40B 0 1 2 3 4
│ Data variables:
│ nstep_in (chain, draw, nstep_in_dim_0) int64 16kB 0 0 0 0 ... 2 2 1
│ nstep_out (chain, draw, nstep_out_dim_0) int64 16kB 379 332 ... 0 0
│ Attributes:
│ created_at: 2026-08-28T01:19:18.202599+00:00
│ creation_library: ArviZ
│ creation_library_version: 1.3.0
│ creation_library_language: Python
│ inference_library: pymc
│ inference_library_version: 6.3.1
│ sample_dims: ['chain', 'draw']
│ sampling_time: 32.14871883392334
│ tuning_steps: 200
│ modeling_interface: bambi
│ modeling_interface_version: 0.20.0
├── Group: /sample_stats
│ Dimensions: (chain: 2, draw: 500, nstep_in_dim_0: 5, nstep_out_dim_0: 5)
│ Coordinates:
│ * chain (chain) int64 16B 0 1
│ * draw (draw) int64 4kB 0 1 2 3 4 5 6 ... 494 495 496 497 498 499
│ * nstep_in_dim_0 (nstep_in_dim_0) int64 40B 0 1 2 3 4
│ * nstep_out_dim_0 (nstep_out_dim_0) int64 40B 0 1 2 3 4
│ Data variables:
│ nstep_in (chain, draw, nstep_in_dim_0) int64 40kB 0 3 4 3 ... 6 0 2
│ nstep_out (chain, draw, nstep_out_dim_0) int64 40kB 0 0 0 0 ... 0 0 0
│ Attributes:
│ created_at: 2026-08-28T01:19:18.204794+00:00
│ creation_library: ArviZ
│ creation_library_version: 1.3.0
│ creation_library_language: Python
│ inference_library: pymc
│ inference_library_version: 6.3.1
│ sample_dims: ['chain', 'draw']
│ sampling_time: 32.14871883392334
│ tuning_steps: 200
│ modeling_interface: bambi
│ modeling_interface_version: 0.20.0
├── Group: /observed_data
│ Dimensions: (__obs__: 1000, rt,response_extra_dim_0: 2)
│ Coordinates:
│ * __obs__ (__obs__) int64 8kB 0 1 2 3 4 ... 996 997 998 999
│ * 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 16kB ...
│ Attributes:
│ created_at: 2026-08-28T01:19:18.206324+00:00
│ creation_library: ArviZ
│ creation_library_version: 1.3.0
│ creation_library_language: Python
│ inference_library: pymc
│ inference_library_version: 6.3.1
│ sample_dims: []
│ modeling_interface: bambi
│ modeling_interface_version: 0.20.0
└── Group: /log_likelihood
Dimensions: (chain: 2, draw: 500, __obs__: 1000)
Coordinates:
* chain (chain) int64 16B 0 1
* draw (draw) int64 4kB 0 1 2 3 4 5 6 ... 493 494 495 496 497 498 499
* __obs__ (__obs__) int64 8kB 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
Data variables:
rt,response (chain, draw, __obs__) float64 8MB 0.2207 -0.3237 ... -2.76
Attributes:
modeling_interface: bambi
modeling_interface_version: 0.20.0import arviz as az
az.plot_trace(model.traces, var_names=["~v0_mean"])
plt.tight_layout()
<environment>/site-packages/arviz_base/utils.py:149: UserWarning: Items starting with ~: ['v0_mean'] have not been found and will be ignored warnings.warn(
az.plot_pair(model.traces, var_names=["~v0_mean"])
plt.tight_layout()
<environment>/site-packages/arviz_base/utils.py:149: UserWarning: Items starting with ~: ['v0_mean'] have not been found and will be ignored warnings.warn(