Specify priors and fix parameters¶

Every HSSM model accepts custom priors — and fixed (non-sampled) parameter values — through the same small set of patterns. This page is the reference walkthrough for all of them: the include argument (dictionaries and hssm.Param), hssm.Prior, bounds, the shortcut keyword syntax, and coefficient priors in regressions.

If you have not fit a model yet, start with the Quickstart.

Run this how-to¶

On Colab, uncomment and run the installation cell below once, then restart the runtime. For local setup, GPU extras, and troubleshooting see the Installation guide.

In [1]:

# %pip install hssm

# %pip install hssm

Setup¶

We reuse the quickstart's simulated DDM dataset.

In [2]:

import arviz as az
import numpy as np

import hssm

%matplotlib inline

v_true, a_true, z_true, t_true = [0.5, 1.5, 0.5, 0.5]
dataset = hssm.simulate_data(
    model="ddm",
    theta=[v_true, a_true, z_true, t_true],
    size=1000,
)

dataset

import arviz as az import numpy as np import hssm %matplotlib inline v_true, a_true, z_true, t_true = [0.5, 1.5, 0.5, 0.5] dataset = hssm.simulate_data( model="ddm", theta=[v_true, a_true, z_true, t_true], size=1000, ) dataset

Out[2]:

	rt	response
0	2.738526	1.0
1	2.379609	1.0
2	1.780142	1.0
3	1.955865	-1.0
4	1.934673	1.0
...	...	...
995	1.741570	1.0
996	1.147650	1.0
997	1.571325	-1.0
998	1.486842	1.0
999	3.830145	1.0

1000 rows × 2 columns

The non-regression case¶

Next, let's take a look at how to specify priors in the non-regression case. In HSSM, parameter specification can be done in two ways:

1. through the include parameter, or
2. through a shortcut

Through the include parameter¶

The include parameter accepts a list of dictionaries or hssm.Param objects. Both dictionaries and hssm.Param objects are equivalent, since the content of the dictionary will be passed as parameters to the hssm.Param class during model creation, so it is more of a matter of preference. We recommend the hssm.Param object because some IDEs will be able to provide prompts for possible options of parameters. In the non-regression case, each dictionary typically looks like this:

{
    "name": "v",
    "prior": {
        "name": "Uniform",
        "lower": -5.0,
        "upper": 5.0,
    },
    "bounds": (-10.0, 10.0)
}

This is equivalent to writing:

hssm.Param(
    "v",
    prior=dict(name="Uniform", upper=-5.0, lower=5.0),
    bounds=(-10.0, 10.0)
)

The name field corresponds to the name of the parameter being specified.

The prior field specifies the distribution of the prior. There are two ways to achieve this:

1. A dictionary with the name of the distribution (typically captalized) and the parameters of the distribution that you would typically set if you were specifying a distribution in PyMC. For example, if you would like to specify pm.Normal(mu=0.0, sigma=1.0) as the prior, then in HSSM, this prior dictionary would be:

{
    "name": "Normal", ## Note it is capitalized
    "mu": 0.0,
    "sigma": 1.0,
}

or, using the dict constructor:

dict(name="Normal", mu=0.0, sigma=1.0)

2. A hssm.Prior object. This is exactly how you would specify priors using bambi (In fact, hssm.Prior is a subclass of bmb.Prior and for the most part can be used interchangeably with bmb.Prior). To specify the same normal prior as above, you would write:

hssm.Prior("Normal", mu=0.0, sigma=1.0)

The bounds field accepts a tuple of floats, indicating the lower and upper bounds for the parameter.

Fixing parameters: sometimes you might want to fix the values of a parameter. You can easily do so by specifying that value to the prior field of the dictionary. In the following example, the paramter v is fixed to 0.5.

{
    "name": "v",
    "prior": 0.5,
}

Now let's make this concrete with an example:

In [3]:

# A Normal prior for `v` without explicit bounds
param_v = {
    "name": "v",
    "prior": {
        "name": "Normal",
        "mu": 0.0,
        "sigma": 2.0,
    },
}

# A Uniform prior for `a`. Using the `dict` function
param_a = hssm.Param(
    "a",
    prior=dict(
        name="Uniform",
        lower=0.01,
        upper=5,
    ),
    bounds=(0, np.inf),
)

# A Uniform prior for `z` over (0, 1) set using hssm.Prior.
# bounds are not set, existing default bounds will be used
param_z = {"name": "z", "prior": hssm.Prior("Uniform", lower=0.0, upper=1.0)}

# A fixed value for t
param_t = {"name": "t", "prior": 0.5}

example_ddm_model = hssm.HSSM(
    data=dataset,
    model="ddm",
    include=[
        param_v,
        param_a,
        param_z,
        param_t,
    ],
)

example_ddm_model

# A Normal prior for `v` without explicit bounds param_v = { "name": "v", "prior": { "name": "Normal", "mu": 0.0, "sigma": 2.0, }, } # A Uniform prior for `a`. Using the `dict` function param_a = hssm.Param( "a", prior=dict( name="Uniform", lower=0.01, upper=5, ), bounds=(0, np.inf), ) # A Uniform prior for `z` over (0, 1) set using hssm.Prior. # bounds are not set, existing default bounds will be used param_z = {"name": "z", "prior": hssm.Prior("Uniform", lower=0.0, upper=1.0)} # A fixed value for t param_t = {"name": "t", "prior": 0.5} example_ddm_model = hssm.HSSM( data=dataset, model="ddm", include=[ param_v, param_a, param_z, param_t, ], ) example_ddm_model

Model initialized successfully.

Out[3]:

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: analytical
Observations: 1000

Parameters:

v:
    Prior: Normal(mu: 0.0, sigma: 2.0)
    Explicit bounds: (-inf, inf)

a:
    Prior: Uniform(lower: 0.01, upper: 5.0)
    Explicit bounds: (0, inf)

z:
    Prior: Uniform(lower: 0.0, upper: 1.0)
    Explicit bounds: (0.0, 1.0)

t:
    Prior: 0.5
    Explicit bounds: (0.0, inf)


Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)

In [4]:

example_ddm_model.sample()

example_ddm_model.sample()

Using default initvals. 

Initializing NUTS using adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [a, z, v]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 7 seconds.

Out[4]:

<xarray.DataTree>
Group: /
├── Group: /posterior
│       Dimensions:  (chain: 4, draw: 1000)
│       Coordinates:
│         * chain    (chain) int64 32B 0 1 2 3
│         * draw     (draw) int64 8kB 0 1 2 3 4 5 6 7 ... 993 994 995 996 997 998 999
│       Data variables:
│           v        (chain, draw) float64 32kB 0.5977 0.5554 0.5565 ... 0.5036 0.5587
│           a        (chain, draw) float64 32kB 1.529 1.478 1.498 ... 1.497 1.506 1.5
│           z        (chain, draw) float64 32kB 0.4791 0.4979 0.485 ... 0.4907 0.5159
│       Attributes:
│           created_at:                  2026-08-12T22:48:57.855420+00:00
│           creation_library:            ArviZ
│           creation_library_version:    1.3.0
│           creation_library_language:   Python
│           inference_library:           pymc
│           inference_library_version:   6.3.0
│           sample_dims:                 ['chain', 'draw']
│           sampling_time:               7.060505151748657
│           tuning_steps:                1000
│           modeling_interface:          bambi
│           modeling_interface_version:  0.20.0
├── Group: /sample_stats
│       Dimensions:                (chain: 4, draw: 1000)
│       Coordinates:
│         * chain                  (chain) int64 32B 0 1 2 3
│         * draw                   (draw) int64 8kB 0 1 2 3 4 5 ... 995 996 997 998 999
│       Data variables: (12/18)
│           largest_eigval         (chain, draw) float64 32kB nan nan nan ... nan nan
│           smallest_eigval        (chain, draw) float64 32kB nan nan nan ... nan nan
│           energy                 (chain, draw) float64 32kB 1.974e+03 ... 1.976e+03
│           energy_error           (chain, draw) float64 32kB -0.1264 ... -0.1279
│           process_time_diff      (chain, draw) float64 32kB 0.001743 ... 0.001736
│           diverging              (chain, draw) bool 4kB False False ... False False
│           ...                     ...
│           step_size_bar          (chain, draw) float64 32kB 0.7508 0.7508 ... 0.7819
│           max_energy_error       (chain, draw) float64 32kB -0.2441 ... -0.6603
│           divergences            (chain, draw) int64 32kB 0 0 0 0 0 0 ... 0 0 0 0 0 0
│           perf_counter_start     (chain, draw) float64 32kB 2.255e+05 ... 2.255e+05
│           index_in_trajectory    (chain, draw) int64 32kB 1 2 2 2 -2 -1 ... 2 2 2 5 2
│           step_size              (chain, draw) float64 32kB 0.7543 0.7543 ... 0.7689
│       Attributes:
│           created_at:                  2026-08-12T22:48:57.860696+00:00
│           creation_library:            ArviZ
│           creation_library_version:    1.3.0
│           creation_library_language:   Python
│           inference_library:           pymc
│           inference_library_version:   6.3.0
│           sample_dims:                 ['chain', 'draw']
│           sampling_time:               7.060505151748657
│           tuning_steps:                1000
│           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-12T22:48:57.862284+00:00
│           creation_library:            ArviZ
│           creation_library_version:    1.3.0
│           creation_library_language:   Python
│           inference_library:           pymc
│           inference_library_version:   6.3.0
│           sample_dims:                 []
│           modeling_interface:          bambi
│           modeling_interface_version:  0.20.0
└── Group: /log_likelihood
        Dimensions:      (chain: 4, draw: 1000, __obs__: 1000)
        Coordinates:
          * chain        (chain) int64 32B 0 1 2 3
          * draw         (draw) int64 8kB 0 1 2 3 4 5 6 ... 993 994 995 996 997 998 999
          * __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 32MB -1.767 -1.523 ... -2.605
        Attributes:
            modeling_interface:          bambi
            modeling_interface_version:  0.20.0

xarray.DataTree

• /posterior(16)

  • Dimensions:
    • chain: 4
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999], shape=(1000,))

  • Data variables: (3)
    • v
      (chain, draw)
      float64
      0.5977 0.5554 ... 0.5036 0.5587

      array([[0.59773063, 0.55544484, 0.55652921, ..., 0.57287307, 0.51571203,0.57165544],[0.58646941, 0.58684992, 0.58363237, ..., 0.61267795, 0.57356208,0.61165203],[0.56335001, 0.56573616, 0.56403897, ..., 0.54903924, 0.55142504,0.60911541],[0.54907622, 0.53640702, 0.53494452, ..., 0.5907025 , 0.50361221,0.55868833]], shape=(4, 1000))

    • a
      (chain, draw)
      float64
      1.529 1.478 1.498 ... 1.506 1.5

      array([[1.5286381 , 1.47783288, 1.49835475, ..., 1.53123384, 1.51097382,1.5027238 ],[1.47545817, 1.47851808, 1.51452509, ..., 1.50286916, 1.50732157,1.4885285 ],[1.47836183, 1.48694363, 1.49172863, ..., 1.52644198, 1.53309317,1.5348791 ],[1.53065334, 1.50420963, 1.51710089, ..., 1.49672751, 1.50593262,1.49967194]], shape=(4, 1000))

    • z
      (chain, draw)
      float64
      0.4791 0.4979 ... 0.4907 0.5159

      array([[0.47912624, 0.49793864, 0.48496443, ..., 0.51275014, 0.48688202,0.48414285],[0.48494861, 0.47915878, 0.48962155, ..., 0.4873207 , 0.47906545,0.48386143],[0.47768529, 0.49207255, 0.48505227, ..., 0.50639399, 0.50790157,0.47137813],[0.49791953, 0.49901493, 0.51247384, ..., 0.47564765, 0.4907108 ,0.51590257]], shape=(4, 1000))

  • Attributes: (11)

    created_at :

    2026-08-12T22:48:57.855420+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.3.0

    creation_library_language :

    Python

    inference_library :

    pymc

    inference_library_version :

    6.3.0

    sample_dims :

    ['chain', 'draw']

    sampling_time :

    7.060505151748657

    tuning_steps :

    1000

    modeling_interface :

    bambi

    modeling_interface_version :

    0.20.0

  /sample_stats(31)

  • Dimensions:
    • chain: 4
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999], shape=(1000,))

  • Data variables: (18)
    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],[nan, nan, nan, ..., nan, nan, nan],[nan, nan, nan, ..., nan, nan, nan],[nan, nan, nan, ..., nan, nan, nan]], shape=(4, 1000))

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],[nan, nan, nan, ..., nan, nan, nan],[nan, nan, nan, ..., nan, nan, nan],[nan, nan, nan, ..., nan, nan, nan]], shape=(4, 1000))

    • energy
      (chain, draw)
      float64
      1.974e+03 1.974e+03 ... 1.976e+03

      array([[1973.92167763, 1973.95397732, 1973.52431359, ..., 1974.94151484,1977.93508765, 1973.64142203],[1973.3508002 , 1974.70273638, 1973.82665565, ..., 1972.93997801,1972.95521581, 1972.97110423],[1973.24287554, 1973.29407194, 1972.1792358 , ..., 1973.58162142,1973.38726275, 1976.64248207],[1973.79119741, 1973.48166957, 1974.68474363, ..., 1974.06998243,1975.42994699, 1976.40157897]], shape=(4, 1000))

    • energy_error
      (chain, draw)
      float64
      -0.1264 -0.06052 ... 1.26 -0.1279

      array([[-0.12637382, -0.06051818, -0.08780733, ...,  0.64459463,0.69812635, -1.48267099],[-0.05412353, -0.11285587, -0.05071936, ...,  0.00965004,-0.03137791,  0.02395944],[ 0.10003992, -0.19921578, -0.00558297, ...,  0.04050039,0.05872448,  0.59156636],[-0.55970052, -0.1324683 ,  0.16205821, ..., -0.17905368,1.25989292, -0.12786681]], shape=(4, 1000))

    • process_time_diff
      (chain, draw)
      float64
      0.001743 0.001714 ... 0.001736

      array([[0.001743, 0.001714, 0.001706, ..., 0.000832, 0.001632, 0.001636],[0.001804, 0.001778, 0.003456, ..., 0.003292, 0.003284, 0.001687],[0.003105, 0.003535, 0.003444, ..., 0.001653, 0.001697, 0.003272],[0.000803, 0.00313 , 0.001544, ..., 0.001757, 0.003642, 0.001736]],shape=(4, 1000))

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False]], shape=(4, 1000))

    • lp
      (chain, draw)
      float64
      -1.973e+03 ... -1.975e+03

      array([[-1973.28045539, -1972.80878819, -1972.21457238, ...,-1974.66129861, -1975.06970184, -1971.97111294],[-1972.6176316 , -1972.49863828, -1971.98273481, ...,-1972.81944008, -1972.58063425, -1972.90262722],[-1972.90562515, -1972.03346132, -1972.00310675, ...,-1972.9319921 , -1973.32822889, -1974.91308849],[-1972.82416728, -1972.31463279, -1973.50364043, ...,-1972.48774521, -1975.28342646, -1975.24765476]], shape=(4, 1000))

    • perf_counter_diff
      (chain, draw)
      float64
      0.001755 0.001714 ... 0.001734

      array([[0.00175471, 0.001714  , 0.00170683, ..., 0.000832  , 0.00163225,0.00163613],[0.001807  , 0.001777  , 0.00345554, ..., 0.00329167, 0.00328146,0.00168117],[0.00310504, 0.00354254, 0.00344271, ..., 0.00166196, 0.00171296,0.00327238],[0.00080229, 0.00312971, 0.00154354, ..., 0.00176238, 0.00366263,0.00173429]], shape=(4, 1000))

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False]], shape=(4, 1000))

    • n_steps
      (chain, draw)
      float64
      3.0 3.0 3.0 3.0 ... 5.0 3.0 7.0 3.0

      array([[3., 3., 3., ..., 1., 3., 3.],[3., 3., 7., ..., 7., 7., 3.],[7., 7., 7., ..., 3., 3., 7.],[1., 7., 3., ..., 3., 7., 3.]], shape=(4, 1000))

    • tree_depth
      (chain, draw)
      int64
      2 2 2 2 2 2 2 3 ... 2 3 3 3 3 2 3 2

      array([[2, 2, 2, ..., 1, 2, 2],[2, 2, 3, ..., 3, 3, 2],[3, 3, 3, ..., 2, 2, 3],[1, 3, 2, ..., 2, 3, 2]], shape=(4, 1000))

    • acceptance_rate
      (chain, draw)
      float64
      1.0 1.0 0.7136 ... 0.968 0.5872 1.0

      array([[1.        , 1.        , 0.71361939, ..., 0.52487527, 0.6452871 ,1.        ],[1.        , 0.87249989, 0.95353672, ..., 0.99741564, 1.        ,0.99210844],[0.95716878, 0.97474164, 0.98940068, ..., 0.88692258, 0.98098885,0.46842336],[1.        , 0.98576031, 0.82348935, ..., 0.96795582, 0.58719096,1.        ]], shape=(4, 1000))

    • step_size_bar
      (chain, draw)
      float64
      0.7508 0.7508 ... 0.7819 0.7819

      array([[0.75082003, 0.75082003, 0.75082003, ..., 0.75082003, 0.75082003,0.75082003],[0.77042606, 0.77042606, 0.77042606, ..., 0.77042606, 0.77042606,0.77042606],[0.81398057, 0.81398057, 0.81398057, ..., 0.81398057, 0.81398057,0.81398057],[0.78192569, 0.78192569, 0.78192569, ..., 0.78192569, 0.78192569,0.78192569]], shape=(4, 1000))

    • max_energy_error
      (chain, draw)
      float64
      -0.2441 -0.2895 ... 1.26 -0.6603

      array([[-0.24406115, -0.28950988,  0.56173341, ...,  0.64459463,0.8247497 , -1.48267099],[-0.23070037,  0.27308823,  0.15127837, ..., -0.33820287,-0.28449098, -0.31389955],[-0.16238984, -0.19921578,  0.03830234, ...,  0.25849214,-0.09755168,  1.62669673],[-0.55970052, -0.1324683 ,  0.21497551, ..., -0.17905368,1.25989292, -0.66025069]], shape=(4, 1000))

    • divergences
      (chain, draw)
      int64
      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]], shape=(4, 1000))

    • perf_counter_start
      (chain, draw)
      float64
      2.255e+05 2.255e+05 ... 2.255e+05

      array([[225543.49866004, 225543.50044971, 225543.50219442, ...,225545.94288271, 225545.94374854, 225545.94541004],[225543.58113579, 225543.58301658, 225543.58484329, ...,225545.91809   , 225545.92141833, 225545.92855221],[225543.53827704, 225543.54141838, 225543.54502712, ...,225545.91355454, 225545.91530588, 225545.91706533],[225543.49941417, 225543.50024883, 225543.50341217, ...,225545.86372579, 225545.86553421, 225545.86925767]],shape=(4, 1000))

    • index_in_trajectory
      (chain, draw)
      int64
      1 2 2 2 -2 -1 3 ... 6 4 2 2 2 5 2

      array([[ 1,  2,  2, ...,  1,  1, -3],[-2, -3, -4, ...,  4,  2,  2],[ 1, -6, -1, ...,  2, -1, -6],[-1,  1,  3, ...,  2,  5,  2]], shape=(4, 1000))

    • step_size
      (chain, draw)
      float64
      0.7543 0.7543 ... 0.7689 0.7689

      array([[0.75430547, 0.75430547, 0.75430547, ..., 0.75430547, 0.75430547,0.75430547],[0.75930532, 0.75930532, 0.75930532, ..., 0.75930532, 0.75930532,0.75930532],[0.72762576, 0.72762576, 0.72762576, ..., 0.72762576, 0.72762576,0.72762576],[0.76892763, 0.76892763, 0.76892763, ..., 0.76892763, 0.76892763,0.76892763]], shape=(4, 1000))

  • Attributes: (11)

    created_at :

    2026-08-12T22:48:57.860696+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.3.0

    creation_library_language :

    Python

    inference_library :

    pymc

    inference_library_version :

    6.3.0

    sample_dims :

    ['chain', 'draw']

    sampling_time :

    7.060505151748657

    tuning_steps :

    1000

    modeling_interface :

    bambi

    modeling_interface_version :

    0.20.0

  /observed_data(12)

  • Dimensions:
    • __obs__: 1000
    • rt,response_extra_dim_0: 2
  • Coordinates: (2)
    • __obs__
      (__obs__)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999], shape=(1000,))

    • 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
      2.739 1.0 2.38 1.0 ... 1.0 3.83 1.0

      array([[ 2.73852587,  1.        ],[ 2.37960911,  1.        ],[ 1.78014183,  1.        ],...,[ 1.57132518, -1.        ],[ 1.4868418 ,  1.        ],[ 3.83014464,  1.        ]], shape=(1000, 2))

  • Attributes: (9)

    created_at :

    2026-08-12T22:48:57.862284+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.3.0

    creation_library_language :

    Python

    inference_library :

    pymc

    inference_library_version :

    6.3.0

    sample_dims :

    []

    modeling_interface :

    bambi

    modeling_interface_version :

    0.20.0

  /log_likelihood(6)

  • Dimensions:
    • chain: 4
    • draw: 1000
    • __obs__: 1000
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999], shape=(1000,))

    • __obs__
      (__obs__)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999], shape=(1000,))

  • Data variables: (1)
    • rt,response
      (chain, draw, __obs__)
      float64
      -1.767 -1.523 ... -0.982 -2.605

      array([[[-1.76673453, -1.52301435, -1.13749275, ..., -2.72181526,-0.98796856, -2.51492873],[-1.84438219, -1.59084294, -1.17452427, ..., -2.64573848,-0.99303288, -2.6101431 ],[-1.81001579, -1.56501292, -1.17150075, ..., -2.62559475,-1.01071106, -2.55743072],...,[-1.81622723, -1.56919342, -1.15904723, ..., -2.80679749,-0.97888086, -2.5560731 ],[-1.81595297, -1.58170385, -1.2061207 , ..., -2.56640672,-1.0546546 , -2.53075127],[-1.80110592, -1.55447052, -1.15901501, ..., -2.6577126 ,-0.99828995, -2.55402121]],[[-1.82283765, -1.56531326, -1.1491213 , ..., -2.65881705,-0.97440314, -2.6066418 ],[-1.81266659, -1.55728414, -1.14806676, ..., -2.64292858,-0.98031884, -2.59302511],[-1.79510856, -1.54778432, -1.14898066, ..., -2.71600557,-0.98505606, -2.54772865],...[-1.81666615, -1.57493351, -1.17691684, ..., -2.72528724,-1.00557678, -2.54400227],[-1.81292741, -1.57194519, -1.17479622, ..., -2.74495981,-1.0038759 , -2.53734236],[-1.74939137, -1.5070642 , -1.12906654, ..., -2.72293394,-0.98868089, -2.49811537]],[[-1.80223584, -1.56383128, -1.17649196, ..., -2.6975206 ,-1.01560834, -2.52430921],[-1.82834586, -1.58532312, -1.1877474 , ..., -2.64338393,-1.01794634, -2.56286996],[-1.83621977, -1.59331708, -1.1893731 , ..., -2.7081226 ,-1.01071672, -2.56368444],...,[-1.79114328, -1.54101985, -1.14403108, ..., -2.65851217,-0.98657231, -2.5585341 ],[-1.82875983, -1.59437664, -1.21602932, ..., -2.55054569,-1.06042646, -2.54175533],[-1.84983139, -1.5971305 , -1.17389072, ..., -2.74725445,-0.98198175, -2.60464428]]], shape=(4, 1000, 1000))

  • Attributes: (2)

    modeling_interface :

    bambi

    modeling_interface_version :

    0.20.0

In [5]:

az.summary(example_ddm_model.traces)

az.summary(example_ddm_model.traces)

Out[5]:

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
v	0.571	0.0338	0.52	0.62	2416	2678	1.00	0.00069	0.00048
a	1.5013	0.0229	1.5	1.5	2943	2533	1.00	0.00042	0.0003
z	0.4901	0.012	0.47	0.51	2230	2436	1.00	0.00025	0.00018

Specifying priors using the shortcut¶

HSSM also supports a syntax very similar to PyMC: You can directly specify priors by passing the prior to the name of the parameter in hssm.HSSM. This is convenient if the prior is simple. Below is an example almost equivalent to the above example:

In [6]:

# All ways to specify priors mentioned above are supported in the shortcut syntax
shortcut_ddm_model = hssm.HSSM(
    data=dataset,
    model="ddm",
    v={"name": "Normal", "mu": 0.0, "sigma": 2.0},
    a=dict(name="Uniform", lower=0.01, upper=5),
    z=hssm.Prior("Uniform", lower=0.01, upper=1.0),
    t=0.5,
)

shortcut_ddm_model

# All ways to specify priors mentioned above are supported in the shortcut syntax shortcut_ddm_model = hssm.HSSM( data=dataset, model="ddm", v={"name": "Normal", "mu": 0.0, "sigma": 2.0}, a=dict(name="Uniform", lower=0.01, upper=5), z=hssm.Prior("Uniform", lower=0.01, upper=1.0), t=0.5, ) shortcut_ddm_model

Model initialized successfully.

Out[6]:

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: analytical
Observations: 1000

Parameters:

v:
    Prior: Normal(mu: 0.0, sigma: 2.0)
    Explicit bounds: (-inf, inf)

a:
    Prior: Uniform(lower: 0.01, upper: 5.0)
    Explicit bounds: (0.0, inf)

z:
    Prior: Uniform(lower: 0.01, upper: 1.0)
    Explicit bounds: (0.0, 1.0)

t:
    Prior: 0.5
    Explicit bounds: (0.0, inf)


Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)

Note that the shortcut syntax also supports specifying bounds. It will be more polished in a future update. We will skip this step for now.

The regression case¶

Built on top of bambi, HSSM uses an lmer-like syntax that makes it extremely straight-forward to specify regressions.

Parameters that are targets of regressions are also specified using dictionaries in include. Below is an example for such dictionaries.

{
    "name": "v",
    "formula": "v ~ 1 + x + y",
    "prior": {
        # All ways to specify priors in the non-regression case
        # work the same way here.
        "Intercept": {"name": "Uniform", "lower": -10.0, "upper": 10.0},
        "x": dict(name="Normal", mu=0, sigma=1),
        "y": hssm.Prior("HalfNormal", sigma=0.5),
        "z": 1.0
    }
    "link": "identity",
    "bounds": (-10.0, 10.0)
}

We see that in the regression case, name and bounds are specified the exact same way as in the non-regression case. The regression formula is specified in a way that's very similar to the lmer package in R. Users that have experience with R formulas should be very familar with this syntax. In this case, the formula means that the parameter v is regressed on variables x and y, which can be found in the dataframe passed to hssm.HSSM. The 1 explicitly specifies an intercept for the regression.

In addition to the formula, users typically need to specify priors for the regression coefficients. This is done in the prior field of the dictionary. Instead of specifying priors for the parameter, the priors are now specified for the corresponding regression coefficients. If not specified, HSSM will use default priors generated in Bambi.

Users might also want to specify a link function for generalized linear models. If left unspecified, the identity link function will be used.

Now let's see an example of a regression:

In [7]:

# Generate simulated data

intercept = 1.5
x = np.random.uniform(-5.0, 5.0, size=1000)
y = np.random.uniform(-5.0, 5.0, size=1000)

v = intercept + 0.8 * x + 0.3 * y

# Generate simulated data intercept = 1.5 x = np.random.uniform(-5.0, 5.0, size=1000) y = np.random.uniform(-5.0, 5.0, size=1000) v = intercept + 0.8 * x + 0.3 * y

In [8]:

true_values = np.column_stack([v, np.repeat([[1.5, 0.5, 0.5]], axis=0, repeats=1000)])
true_values.shape

true_values = np.column_stack([v, np.repeat([[1.5, 0.5, 0.5]], axis=0, repeats=1000)]) true_values.shape

Out[8]:

(1000, 4)

In [9]:

dataset_reg_v = hssm.simulate_data(
    model="ddm",
    theta=true_values,
    size=1,  # Generate one data point for each of the 1000 set of true values
)

dataset_reg_v["x"] = x
dataset_reg_v["y"] = y

dataset_reg_v

dataset_reg_v = hssm.simulate_data( model="ddm", theta=true_values, size=1, # Generate one data point for each of the 1000 set of true values ) dataset_reg_v["x"] = x dataset_reg_v["y"] = y dataset_reg_v

Out[9]:

	rt	response	x	y
0	2.638611	-1.0	-0.340989	-4.840536
1	0.922815	1.0	2.052331	-4.828024
2	2.075011	-1.0	-3.574361	2.599690
3	3.206385	1.0	-0.338693	0.398921
4	2.437208	-1.0	-1.220605	-4.044089
...	...	...	...	...
995	2.609189	-1.0	-2.462546	0.113518
996	1.458471	1.0	-0.867655	-2.398150
997	1.162385	-1.0	-4.322686	-1.255193
998	1.447900	1.0	-2.680532	1.878721
999	1.172311	1.0	1.439159	1.969433

1000 rows × 4 columns

In [10]:

model_reg_v = hssm.HSSM(
    data=dataset_reg_v,
    include=[
        {
            "name": "v",
            "formula": "v ~ 1 + x + y",
            "prior": {
                "Intercept": {"name": "Uniform", "lower": -0.001, "upper": 0.5},
                "x": dict(name="Uniform", lower=0.0, upper=1.0),
                "y": hssm.Prior("Uniform", lower=0.0, upper=1.0),
            },
            "link": "identity",
        }
    ],
)
model_reg_v

model_reg_v = hssm.HSSM( data=dataset_reg_v, include=[ { "name": "v", "formula": "v ~ 1 + x + y", "prior": { "Intercept": {"name": "Uniform", "lower": -0.001, "upper": 0.5}, "x": dict(name="Uniform", lower=0.0, upper=1.0), "y": hssm.Prior("Uniform", lower=0.0, upper=1.0), }, "link": "identity", } ], ) model_reg_v

Model initialized successfully.

Out[10]:

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: analytical
Observations: 1000

Parameters:

v:
    Formula: v ~ 1 + x + y
    Priors:
        v_Intercept ~ Uniform(lower: -0.001, upper: 0.5)
        v_x ~ Uniform(lower: 0.0, upper: 1.0)
        v_y ~ Uniform(lower: 0.0, upper: 1.0)
    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)

In [11]:

model_reg_v.initvals

model_reg_v.initvals

Out[11]:

{'a': array(1.5),
 'z': array(0.5),
 't': array(0.025),
 'v_Intercept': array(0.),
 'v_x': array(0.5),
 'v_y': array(0.5)}

In [12]:

# Uncomment to see model graph if you have graphviz installed
# model_reg_v.graph()

# Uncomment to see model graph if you have graphviz installed # model_reg_v.graph()

In [13]:

trace_reg_v = model_reg_v.sample()

trace_reg_v = model_reg_v.sample()

Using default initvals. 

Initializing NUTS using adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [a, z, t, v_Intercept, v_x, v_y]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 14 seconds.

In [14]:

# Looks like parameter recovery was successful
az.summary(model_reg_v.traces)

# Looks like parameter recovery was successful az.summary(model_reg_v.traces)

Out[14]:

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
z	0.618	0.0085	0.6	0.63	4178	3101	1.00	0.00013	9.4e-05
v_y	0.2613	0.0144	0.24	0.28	4411	3402	1.00	0.00022	0.00016
t	0.5487	0.0075	0.54	0.56	3809	2866	1.00	0.00012	8.6e-05
a	1.269	0.0295	1.2	1.3	3330	2744	1.00	0.00051	0.00036
v_x	0.6173	0.0183	0.59	0.65	3500	2952	1.00	0.00031	0.00023
v_Intercept	0.5023	0.004	0.5	0.51	4827	3959	1.00	5.8e-05	5.7e-05

In [15]:

az.plot_trace_dist(model_reg_v.traces);

az.plot_trace_dist(model_reg_v.traces);

See also¶

• hssm.Prior API reference and hssm.Param API reference
• Setting initial values — when sampling needs a hand beyond priors
• The HSSM tutorial, which motivates these choices in context

Prior and link settings for hierarchical models¶

bambi is not designed with HSSM in mind: where priors for regression parameters are not explicitly defined, bambi's defaults are sometimes not optimal, and "identity" link functions tend not to work well for parameters that are not defined on the whole real line. HSSM therefore provides two safeguard settings, most useful for hierarchical models (the examples below use the cavanagh_theta dataset with a global_formula that gives every parameter a participant-level random intercept).

In [16]:

cav_data = hssm.load_data("cavanagh_theta")

cav_data = hssm.load_data("cavanagh_theta")

prior_settings¶

Currently we provide a "safe" strategy that uses HSSM default priors, which is turned on by default for parameters that are targets of regressions. One can compare the two models below, with safe strategy turned on and off:

In [17]:

model_safe = hssm.HSSM(
    data=cav_data,
    global_formula="y ~ 1 + (1|participant_id)",
    prior_settings="safe",
    loglik_kind="approx_differentiable",
)
model_safe

model_safe = hssm.HSSM( data=cav_data, global_formula="y ~ 1 + (1|participant_id)", prior_settings="safe", loglik_kind="approx_differentiable", ) model_safe

Model initialized successfully.

Out[17]:

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: approx_differentiable
Observations: 3988

Parameters:

v:
    Formula: v ~ 1 + (1|participant_id)
    Priors:
        v_Intercept ~ Normal(mu: 0.0, sigma: 0.25)
        v_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (-3.0, 3.0)

a:
    Formula: a ~ 1 + (1|participant_id)
    Priors:
        a_Intercept ~ Normal(mu: 1.4, sigma: 0.25)
        a_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (0.3, 2.5)

z:
    Formula: z ~ 1 + (1|participant_id)
    Priors:
        z_Intercept ~ Normal(mu: 0.5, sigma: 0.25)
        z_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (0.0, 1.0)

t:
    Formula: t ~ 1 + (1|participant_id)
    Priors:
        t_Intercept ~ Normal(mu: 1.0, sigma: 0.25)
        t_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (0.0, 2.0)


Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)

In [18]:

model_safe_off = hssm.HSSM(
    data=cav_data,
    global_formula="y ~ 1 + (1|participant_id)",
    prior_settings=None,
    loglik_kind="approx_differentiable",
)
model_safe_off

model_safe_off = hssm.HSSM( data=cav_data, global_formula="y ~ 1 + (1|participant_id)", prior_settings=None, loglik_kind="approx_differentiable", ) model_safe_off

Model initialized successfully.

Out[18]:

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: approx_differentiable
Observations: 3988

Parameters:

v:
    Formula: v ~ 1 + (1|participant_id)
    Priors:
        v_Intercept ~ Normal(mu: 0.0, sigma: 2.5)
        v_1|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 2.5))
    Link: identity
    Explicit bounds: (-3.0, 3.0)

a:
    Formula: a ~ 1 + (1|participant_id)
    Priors:
        a_Intercept ~ Normal(mu: 0.0, sigma: 1.0)
        a_1|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 1.0))
    Link: identity
    Explicit bounds: (0.3, 2.5)

z:
    Formula: z ~ 1 + (1|participant_id)
    Priors:
        z_Intercept ~ Normal(mu: 0.0, sigma: 1.0)
        z_1|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 1.0))
    Link: identity
    Explicit bounds: (0.0, 1.0)

t:
    Formula: t ~ 1 + (1|participant_id)
    Priors:
        t_Intercept ~ Normal(mu: 0.0, sigma: 1.0)
        t_1|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 1.0))
    Link: identity
    Explicit bounds: (0.0, 2.0)


Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)

link_settings¶

We also provide a link_settings switch, which changes default link functions for parameters according to their explicit bounds. See the model below with link_settings set to "log_logit":

In [19]:

model_log_logit = hssm.HSSM(
    data=cav_data,
    global_formula="y ~ 1 + (1|participant_id)",
    prior_settings=None,
    link_settings="log_logit",
)
model_log_logit

model_log_logit = hssm.HSSM( data=cav_data, global_formula="y ~ 1 + (1|participant_id)", prior_settings=None, link_settings="log_logit", ) model_log_logit

Model initialized successfully.

Out[19]:

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: analytical
Observations: 3988

Parameters:

v:
    Formula: v ~ 1 + (1|participant_id)
    Priors:
        v_Intercept ~ Normal(mu: 0.0, sigma: 2.5)
        v_1|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 2.5))
    Link: identity
    Explicit bounds: (-inf, inf)
 (ignored due to link function)
a:
    Formula: a ~ 1 + (1|participant_id)
    Priors:
        a_Intercept ~ Normal(mu: 0.0, sigma: 1.0)
        a_1|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 1.0))
    Link: log
    Explicit bounds: (0.0, inf)
 (ignored due to link function)
z:
    Formula: z ~ 1 + (1|participant_id)
    Priors:
        z_Intercept ~ Normal(mu: 0.0, sigma: 1.0)
        z_1|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 1.0))
    Link: Generalized logit link function with bounds (0.0, 1.0)
    Explicit bounds: (0.0, 1.0)
 (ignored due to link function)
t:
    Formula: t ~ 1 + (1|participant_id)
    Priors:
        t_Intercept ~ Normal(mu: 0.0, sigma: 1.0)
        t_1|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 1.0))
    Link: log
    Explicit bounds: (0.0, inf)
 (ignored due to link function)

Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)

Mixing strategies:¶

It is possible to turn on both prior_settings and link_settings:

In [20]:

model_safe_loglogit = hssm.HSSM(
    data=cav_data,
    global_formula="y ~ 1 + (1|participant_id)",
    prior_settings="safe",
    link_settings="log_logit",
)
model_safe_loglogit

model_safe_loglogit = hssm.HSSM( data=cav_data, global_formula="y ~ 1 + (1|participant_id)", prior_settings="safe", link_settings="log_logit", ) model_safe_loglogit

Model initialized successfully.

Out[20]:

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: analytical
Observations: 3988

Parameters:

v:
    Formula: v ~ 1 + (1|participant_id)
    Priors:
        v_Intercept ~ Normal(mu: 0.0, sigma: 0.25)
        v_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (-inf, inf)
 (ignored due to link function)
a:
    Formula: a ~ 1 + (1|participant_id)
    Priors:
        a_Intercept ~ Normal(mu: 0.0, sigma: 0.25)
        a_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: log
    Explicit bounds: (0.0, inf)
 (ignored due to link function)
z:
    Formula: z ~ 1 + (1|participant_id)
    Priors:
        z_Intercept ~ Normal(mu: 0.0, sigma: 0.25)
        z_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: Generalized logit link function with bounds (0.0, 1.0)
    Explicit bounds: (0.0, 1.0)
 (ignored due to link function)
t:
    Formula: t ~ 1 + (1|participant_id)
    Priors:
        t_Intercept ~ Normal(mu: 0.0, sigma: 0.25)
        t_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: log
    Explicit bounds: (0.0, inf)
 (ignored due to link function)

Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)