Using the do-operator to simulate data¶

This tutorial illustrates how we can use the underlying PyMC model of a given HSSM model, to forward simulate datasets that you may use for numerical studies.

In [1]:

import pymc as pm

import hssm

import pymc as pm import hssm

Simulate Data¶

We will simulate a simple dataset that contains two conditions.

In [2]:

# Simple model
data = hssm.simulate_data(model="ddm", theta=dict(v=0.5, a=1.5, z=0.5, t=0.1), size=500)

# Simple model data = hssm.simulate_data(model="ddm", theta=dict(v=0.5, a=1.5, z=0.5, t=0.1), size=500)

Specify HSSM Model¶

We will fit two separate models to the data.

In [3]:

# Model 1
hssm_model = hssm.HSSM(model="ddm", data=data)

# Model 1 hssm_model = hssm.HSSM(model="ddm", data=data)

Model initialized successfully.

Underlying PyMC model¶

In [4]:

pm.model_to_graphviz(hssm_model.pymc_model)

pm.model_to_graphviz(hssm_model.pymc_model)

Out[4]:

Using the do-operator¶

In [5]:

synth_idata, synth_model = hssm_model.sample_do(
    params={"v": 0, "a": 1.5, "z": 0.5, "t": 0.1}, draws=100, return_model=True
)

synth_idata, synth_model = hssm_model.sample_do( params={"v": 0, "a": 1.5, "z": 0.5, "t": 0.1}, draws=100, return_model=True )

/Users/yxu150/HSSM/.venv/lib/python3.13/site-packages/pytensor/link/numba/dispatch/basic.py:214: UserWarning: Numba will use object mode to run ddm_RV_rv{"(),(),(),(),()->(2)"}'s perform method. Set `pytensor.config.compiler_verbose = True` to see more details.
  warnings.warn(
Sampling: [rt,response]

Note that the process defines an new PyMC model, which sets the parameters we passed in our dictionary to Data objects.

In [6]:

pm.model_to_graphviz(synth_model)

pm.model_to_graphviz(synth_model)

Out[6]:

In [7]:

synth_idata

synth_idata

Out[7]:

<xarray.DataTree>
Group: /
├── Group: /prior
│       Dimensions:       (chain: 1, draw: 100, v_mean_dim_0: 1)
│       Coordinates:
│         * chain         (chain) int64 8B 0
│         * draw          (draw) int64 800B 0 1 2 3 4 5 6 7 ... 92 93 94 95 96 97 98 99
│         * v_mean_dim_0  (v_mean_dim_0) int64 8B 0
│       Data variables:
│           v_mean        (chain, draw, v_mean_dim_0) float64 800B 0.0 0.0 ... 0.0 0.0
│       Attributes:
│           created_at:                 2026-06-29T18:40:21.224201+00:00
│           creation_library:           ArviZ
│           creation_library_version:   1.1.0
│           creation_library_language:  Python
│           inference_library:          pymc
│           inference_library_version:  6.0.1
│           sample_dims:                ['chain', 'draw']
├── Group: /prior_predictive
│       Dimensions:          (chain: 1, draw: 100, __obs__: 500, rt,response_dim: 2)
│       Coordinates:
│         * chain            (chain) int64 8B 0
│         * draw             (draw) int64 800B 0 1 2 3 4 5 6 7 ... 93 94 95 96 97 98 99
│         * __obs__          (__obs__) int64 4kB 0 1 2 3 4 5 ... 494 495 496 497 498 499
│         * rt,response_dim  (rt,response_dim) int64 16B 0 1
│       Data variables:
│           rt,response      (chain, draw, __obs__, rt,response_dim) float64 800kB 0....
│       Attributes:
│           created_at:                 2026-06-29T18:40:21.225297+00:00
│           creation_library:           ArviZ
│           creation_library_version:   1.1.0
│           creation_library_language:  Python
│           inference_library:          pymc
│           inference_library_version:  6.0.1
│           sample_dims:                ['chain', 'draw']
├── 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 3...
│       Attributes:
│           created_at:                 2026-06-29T18:40:21.226032+00:00
│           creation_library:           ArviZ
│           creation_library_version:   1.1.0
│           creation_library_language:  Python
│           inference_library:          pymc
│           inference_library_version:  6.0.1
│           sample_dims:                []
└── Group: /constant_data
        Dimensions:  ()
        Data variables:
            v        float64 8B 0.0
            t        float64 8B 0.1
            z        float64 8B 0.5
            a        float64 8B 1.5
        Attributes:
            created_at:                 2026-06-29T18:40:21.226577+00:00
            creation_library:           ArviZ
            creation_library_version:   1.1.0
            creation_library_language:  Python
            inference_library:          pymc
            inference_library_version:  6.0.1
            sample_dims:                []

xarray.DataTree

• /prior(11)

  • Dimensions:
    • chain: 1
    • draw: 100
    • v_mean_dim_0: 1
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0

      array([0])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 6 ... 94 95 96 97 98 99

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,90, 91, 92, 93, 94, 95, 96, 97, 98, 99])

    • v_mean_dim_0
      (v_mean_dim_0)
      int64
      0

      array([0])

  • Data variables: (1)
    • v_mean
      (chain, draw, v_mean_dim_0)
      float64
      0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

      array([[[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],...[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.]]])

  • Attributes: (7)

    created_at :

    2026-06-29T18:40:21.224201+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.1.0

    creation_library_language :

    Python

    inference_library :

    pymc

    inference_library_version :

    6.0.1

    sample_dims :

    ['chain', 'draw']

  /prior_predictive(12)

  • Dimensions:
    • chain: 1
    • draw: 100
    • __obs__: 500
    • rt,response_dim: 2
  • Coordinates: (4)
    • chain
      (chain)
      int64
      0

      array([0])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 6 ... 94 95 96 97 98 99

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,90, 91, 92, 93, 94, 95, 96, 97, 98, 99])

    • __obs__
      (__obs__)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499], shape=(500,))

    • rt,response_dim
      (rt,response_dim)
      int64
      0 1

      array([0, 1])

  • Data variables: (1)
    • rt,response
      (chain, draw, __obs__, rt,response_dim)
      float64
      0.955 1.0 1.818 ... -1.0 8.824 -1.0

      array([[[[ 0.95497626,  1.        ],[ 1.81766438,  1.        ],[ 0.5413695 ,  1.        ],...,[ 0.59160769, -1.        ],[ 1.54743767, -1.        ],[ 2.24203038,  1.        ]],[[ 7.40745926,  1.        ],[ 7.04285717,  1.        ],[12.57329556,  1.        ],...,[ 3.55930996, -1.        ],[ 0.71321392,  1.        ],[ 0.70120293, -1.        ]],[[ 1.39014745,  1.        ],[ 0.53823096, -1.        ],[ 0.63175118, -1.        ],...,......,[ 3.02899003,  1.        ],[ 3.77768183, -1.        ],[ 6.13585806, -1.        ]],[[ 3.26716757, -1.        ],[ 0.81551516,  1.        ],[ 1.63665247,  1.        ],...,[ 1.80467117, -1.        ],[ 0.82886618,  1.        ],[ 0.78016767, -1.        ]],[[ 2.17258978, -1.        ],[ 2.93716359, -1.        ],[ 4.1945734 , -1.        ],...,[ 3.04802656,  1.        ],[ 4.24854136, -1.        ],[ 8.82415581, -1.        ]]]], shape=(1, 100, 500, 2))

  • Attributes: (7)

    created_at :

    2026-06-29T18:40:21.225297+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.1.0

    creation_library_language :

    Python

    inference_library :

    pymc

    inference_library_version :

    6.0.1

    sample_dims :

    ['chain', 'draw']

  /observed_data(10)

  • Dimensions:
    • __obs__: 500
    • rt,response_extra_dim_0: 2
  • Coordinates: (2)
    • __obs__
      (__obs__)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499], shape=(500,))

    • rt,response_extra_dim_0
      (rt,response_extra_dim_0)
      int64
      0 1

      array([0, 1])

  • Data variables: (1)
    • rt,response
      (__obs__, rt,response_extra_dim_0)
      float64
      3.292 1.0 0.8728 ... 1.0 1.133 1.0

      array([[ 3.29207063,  1.        ],[ 0.87280297, -1.        ],[ 3.98631644,  1.        ],[ 0.32350719,  1.        ],[ 3.95413184,  1.        ],[ 1.40061021, -1.        ],[ 3.41857004,  1.        ],[ 3.91596532,  1.        ],[ 0.48281989,  1.        ],[ 1.42482746,  1.        ],[ 0.95325053,  1.        ],[ 3.53079796, -1.        ],[ 1.30773926,  1.        ],[ 2.47188473,  1.        ],[ 2.42441487,  1.        ],[ 1.02935624,  1.        ],[ 0.40367547,  1.        ],[ 1.96804404,  1.        ],[ 2.87357235,  1.        ],[ 0.5239116 ,  1.        ],...[ 1.33189356,  1.        ],[ 1.22963524,  1.        ],[ 2.52636862, -1.        ],[ 1.09853423,  1.        ],[ 0.83103961, -1.        ],[ 0.61490166, -1.        ],[ 2.37368751,  1.        ],[ 1.3043437 ,  1.        ],[ 1.27017236,  1.        ],[ 4.64579535, -1.        ],[ 1.26372123,  1.        ],[ 2.34292579,  1.        ],[ 4.27715445,  1.        ],[ 1.57954025,  1.        ],[ 3.53849244, -1.        ],[ 0.92741507,  1.        ],[ 3.82046223, -1.        ],[ 1.8231287 ,  1.        ],[ 2.20301294,  1.        ],[ 1.13316762,  1.        ]])

  • Attributes: (7)

    created_at :

    2026-06-29T18:40:21.226032+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.1.0

    creation_library_language :

    Python

    inference_library :

    pymc

    inference_library_version :

    6.0.1

    sample_dims :

    []

  /constant_data(11)

  • Dimensions:
  • Data variables: (4)
    • v
      ()
      float64
      0.0

      array(0.)

    • t
      ()
      float64
      0.1

      array(0.1)

    • z
      ()
      float64
      0.5

      array(0.5)

    • a
      ()
      float64
      1.5

      array(1.5)

  • Attributes: (7)

    created_at :

    2026-06-29T18:40:21.226577+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.1.0

    creation_library_language :

    Python

    inference_library :

    pymc

    inference_library_version :

    6.0.1

    sample_dims :

    []

Observe two things:

1. A closer look at the prior group reveals that no parameters were actually sampled (it collects on the deterministic v_mean which, arguably, we can skip too)
2. The prior_predictive group is where we collected our synthetic data

Our synthetic data sits in the prior_predictive group, since, under the hood, were de facto called the sample_prior_predictive() method, on our new do-model.

If you prefer to get your outputs as a dataframe, you can use the predictive_dt_to_dataframe method from the hssm.utils module.

In [8]:

synth_df = hssm.utils.predictive_dt_to_dataframe(
    synth_idata, predictive_group="prior_predictive"
)

synth_df = hssm.utils.predictive_dt_to_dataframe( synth_idata, predictive_group="prior_predictive" )

In [9]:

synth_df

synth_df

Out[9]:

	chain	draw	__obs__	rt	response
0	0	0	0	0.954976	1.0
1	0	0	1	1.817664	1.0
2	0	0	2	0.541369	1.0
3	0	0	3	3.444736	1.0
4	0	0	4	3.674135	1.0
...	...	...	...	...	...
49995	0	99	495	1.368595	1.0
49996	0	99	496	5.492402	1.0
49997	0	99	497	3.048027	1.0
49998	0	99	498	4.248541	-1.0
49999	0	99	499	8.824156	-1.0

50000 rows × 5 columns

From here you may re-attach the original data used to define the hssm_model (e.g.to re-attach covariates that you want as part of the model you are testing on your synthetic data).

The sample_do() function can be extremely helpful in setting up numerical experiments to test e.g. for parameter identifiability in scenarios that you wish to control precisely.