This tutorial is a rather comprehensive introduction to HDDM with focus on the new LAN extension.
The methods behind the new HDDMnn()
, HDDMnnRegressor()
and HDDMnnStimcoding()
classes can be found in our original dedicated publication.
These are new featues. Please let us know on the HDDM forum and/or via github reports regarding bugs or other limitations and we will do our best to help as soon as we can.
Networks were trained over a fairly wide range of parameters which hopefully capture the scope of common empirical data. The networks will not accurately report likelihoods outside that range, so we explicitly limit the range of parameters that can be sampled from. If you find that your posterior samples reach and get stuck at the allowed parameter bounds (which you will see in the posterior plots), please notify us and we will do our best to provide improved networks over time.
You may encounter more print output than with standard HDDM. These are sanity checks and the verbosity will vanish progressively.
In the upper left menu click on Runtime, then Change runtime type and select GPU as hardware accelerator
# Note: Usually colab has all other packages which we may use already installed
# The basic configuration of colabs does change over time, so you may have to add
# some install commands here if imports below don't work for package xyz
!pip install scikit-learn
!pip install cython
!pip install pymc==2.3.8
Requirement already satisfied: scikit-learn in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (0.24.2) Requirement already satisfied: numpy>=1.13.3 in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (from scikit-learn) (1.19.1) Requirement already satisfied: scipy>=0.19.1 in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (from scikit-learn) (1.7.2) Requirement already satisfied: joblib>=0.11 in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (from scikit-learn) (1.0.1) Requirement already satisfied: threadpoolctl>=2.0.0 in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (from scikit-learn) (2.1.0) Requirement already satisfied: cython in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (0.29.24) Requirement already satisfied: pymc in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (2.3.8)
!pip install -U --no-deps git+https://github.com/hddm-devs/hddm
!pip install -U --no-deps git+https://github.com/hddm-devs/kabuki
Collecting git+https://github.com/hddm-devs/hddm Cloning https://github.com/hddm-devs/hddm to /private/var/folders/gx/s43vynx550qbypcxm83fv56dzq4hgg/T/pip-req-build-xzqqwrcn Running command git clone -q https://github.com/hddm-devs/hddm /private/var/folders/gx/s43vynx550qbypcxm83fv56dzq4hgg/T/pip-req-build-xzqqwrcn Running command git submodule update --init --recursive -q Building wheels for collected packages: HDDM Building wheel for HDDM (setup.py) ... |
# MODULE IMPORTS ----
# warning settings
import warnings
warnings.simplefilter(action="ignore", category=FutureWarning)
# Data management
import pandas as pd
import numpy as np
import pickle
# Plotting
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
# Stats functionality
from statsmodels.distributions.empirical_distribution import ECDF
# HDDM
import hddm
from hddm.simulators.hddm_dataset_generators import simulator_h_c
The main concern of this notebook is to present the extended capabilities of the HDDM toolbox as a result of the new HDDMnn
classes.
Primarily we are interested in the additional models we can now be fit to data. So let's take stock of the models that were added to standard HDDM.
A model with a linearly collapsing angle. Adds a parameter $\theta$, which specifies the angle of the bound.
A model that includes a collapsing bound parameterized as the scaled cdf of a Weibull distribution. This adds two parameters to the standard DDM, $\alpha$ and $\beta$.
The Levy model is essentially a standard DDM where noise is not driven by a Gaussian distribution, but the noise process is now parameterized by the new parameter $\alpha$, which interpolates between a Gausian $\alpha = 2$ and a Cauchy (heavy tailed) $\alpha = 1$.
This model implements the 2-choice LCA, which includes a an inhibition / excitation parameter $g$.
Find more details on these models in our companion paper.
The addition of 3 choice and 4 choice models, comes with slightly more limited functionality as compared to 2 choice models. Specifically, not all plot-concepts currently standard in HDDM translate immediately to models with more choice options. We are trying to align this functionality going forward.
Please find the original description in this paper.
Race models simply take out the mutual and self-inhibition of LCAs.
Implements an linearly collapsing bound as above under the respective 2 choice models
Let's first take a look at some of the useful metadata we can use to set up our models and simulators.
If we type hddm.simulators.model_config
, we get back a dictionary that stores a bunch of information
for each of the models that are currently implemented in HDDM. It lists,
doc
string that gives some information about the status of the model as it pertains to it's usability as well as some potential usage tips. Please read the doc
string before using any of the new models.params
,param_bounds
boundary
)params_default
).slice_widths
)params_trans
you can choose parameters which will be logit transformed for sampling (order as in params
)choices
determines valid choice options under the modelhddm_include
, it lists the parameters which we want to include when initializing our HDDM Model with one of the sequential sampling models available.You won't need most of these options if you are getting started, but they do provide you with useful information and a couple extra degrees of freedom when it comes to optimizing your sampler.
# List the models currently available
hddm.model_config.model_config.keys()
dict_keys(['ddm_hddm_base', 'full_ddm_hddm_base', 'ddm', 'angle', 'weibull', 'levy', 'full_ddm', 'ornstein', 'ddm_sdv', 'gamma_drift', 'gamma_drift_angle', 'ds_conflict_drift', 'ds_conflict_drift_angle', 'ddm_par2', 'ddm_par2_no_bias', 'ddm_par2_angle_no_bias', 'ddm_par2_weibull_no_bias', 'ddm_seq2', 'ddm_seq2_no_bias', 'ddm_seq2_angle_no_bias', 'ddm_seq2_weibull_no_bias', 'ddm_mic2_adj', 'ddm_mic2_adj_no_bias', 'ddm_mic2_adj_angle_no_bias', 'ddm_mic2_adj_weibull_no_bias', 'race_no_bias_3', 'race_no_bias_angle_3', 'race_no_bias_4', 'race_no_bias_angle_4', 'lca_no_bias_3', 'lca_no_bias_angle_3', 'lca_no_bias_4', 'lca_no_bias_angle_4', 'weibull_cdf', 'full_ddm2'])
You find two kinds of extra models which were not mentioned in the model listing above:
hddm_base
models are used predominantly with the basic HDDM()
classes. These models are not to be used with the HDDMnn()
classes.Now taking a closer look at the angle
model
# Metadata
model = "ddm"
n_samples = 1000
# Config for our current model
hddm.model_config.model_config[model]
{'doc': 'Basic DDM. Meant for use with the LAN extension. \nNote that the boundaries here are coded as -a, and a in line with all other models meant for the LAN extension. \nTo compare model fits between standard HDDM and HDDMnn when using the DDM model, multiply the boundary (a) parameter by 2. \nWe recommend using standard HDDM if you are interested in the basic DDM, but you might want to use this for testing.', 'params': ['v', 'a', 'z', 't'], 'params_trans': [0, 0, 1, 0], 'params_std_upper': [1.5, 1.0, None, 1.0], 'param_bounds': [[-3.0, 0.3, 0.1, 0.001], [3.0, 2.5, 0.9, 2.0]], 'boundary': <function hddm.simulators.boundary_functions.constant(t=0)>, 'params_default': [0.0, 1.0, 0.5, 0.001], 'hddm_include': ['z'], 'choices': [-1, 1], 'slice_widths': {'v': 1.5, 'v_std': 1, 'a': 1, 'a_std': 1, 'z': 0.1, 'z_trans': 0.2, 't': 0.01, 't_std': 0.15}}
# Looking at the doc string before using the model
print(hddm.model_config.model_config[model]["doc"])
Basic DDM. Meant for use with the LAN extension. Note that the boundaries here are coded as -a, and a in line with all other models meant for the LAN extension. To compare model fits between standard HDDM and HDDMnn when using the DDM model, multiply the boundary (a) parameter by 2. We recommend using standard HDDM if you are interested in the basic DDM, but you might want to use this for testing.
Let's start by generating some data from the angle
model. For this you have available the simulators
module, specifically we will start with the simulator_h_c
function.
If you are curious about all the capabilities of this function, please check the help()
function for it.
from hddm.simulators.hddm_dataset_generators import simulator_h_c
data, full_parameter_dict = simulator_h_c(
n_subjects=1,
n_trials_per_subject=n_samples,
model=model,
p_outlier=0.00,
conditions=None,
depends_on=None,
regression_models=None,
regression_covariates=None,
group_only_regressors=False,
group_only=None,
fixed_at_default=None,
)
A quick look into what the simulator spits out (you can also read about it in the docs).
We get back a tuple
of two:
First, a DataFrame which holds a rt
, a response
and a subj_idx
column as well as trial-by-trial ground truth parameters.
Second a parameter dictionary which has parameter names in accordance with HDDM()
trace names. This is useful for some of our plots.
data
rt | response | subj_idx | v | a | z | t | |
---|---|---|---|---|---|---|---|
0 | 1.956185 | 1.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
1 | 1.035191 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
2 | 1.004191 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
3 | 1.510186 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
4 | 1.164191 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
... | ... | ... | ... | ... | ... | ... | ... |
995 | 1.697184 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
996 | 1.520186 | 1.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
997 | 1.552186 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
998 | 1.038191 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
999 | 0.932191 | 1.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
1000 rows × 7 columns
# Here unspectacularly, parameter names are unchanged
# (single subject fits do not need any parameter name augmentation)
full_parameter_dict
{'v': -0.48173086489284433, 'a': 0.6556418306610691, 't': 0.8871907031605131, 'z': 0.4398408702789776}
Now that we have our simulated data, we look to visualise it. Let's look at a couple of plots that we can use for this purpose.
The HDDM.plotting
module includes the plot_from_data
function, which allows you to plot
subsets from a dataset, according to a grouping specified by the groupby
argument.
The plot creates a matplotlib.axes
object for each subset, and you can provide a function to manipulate
this axes object. Some of these axes manipulators are provided your you. Here we focus on the
_plot_func_model
axes manipulator supplied under the plot_func
argument.
Check out the arguments of plot_from_data
and _plot_func_model
using the help()
function.
You have quite some freedom in styling these plots.
We will refer to this plot as the model cartoon plot
.
hddm.plotting.plot_from_data(
df=data,
generative_model=model,
columns=1,
groupby=["subj_idx"],
figsize=(4, 3),
value_range=np.arange(0, 5, 0.1),
plot_func=hddm.plotting._plot_func_model,
**{"alpha": 1.0, "ylim": 3, "add_data_rts": True, "add_data_model": False}
)
plt.show()
subj_idx(0)
If we set add_model = True
, this will add a cartoon of the model on top of the histograms.
This model cartoon plot
will only work for 2-choice models for now.
Moreover, often useful for illustration purposes, we can include a bunch of simulations trajectories into the model plot (note the corresponding arguments). Common to all models currently included is their conceptual reliance on there particle trajectories. Reaction times and choices are simulated as boundary crossings of these particles. If you don't want to include these trajectories, just set show_trajectories = False
.
hddm.plotting.plot_from_data(
df=data,
generative_model=model,
columns=1,
groupby=["subj_idx"],
figsize=(4, 3),
value_range=np.arange(0, 5, 0.1),
plot_func=hddm.plotting._plot_func_model,
**{"alpha": 1.0, "ylim": 3, "add_data_rts": True, "add_data_model": True}
)
plt.show()
subj_idx(0)
If you are interested, you can use this plot to investigate the behavior of models across different parameters setups.
Now, we try to fit these models to data! Let's start with an simple dataset. In other words, we have one single participant who provides $n$ datatpoints (reaction times and choices) from some two alternative forced choice task paradigm.
In this demo we fit to simulated data. This serves as a template, and you can easily adapt it to your needs.
# Metadata
nmcmc = 1500
model = "angle"
n_samples = 1000
includes = hddm.model_config.model_config[model]["hddm_include"]
When defining includes
,
you can also pick only as subset of the parameters suggested under hddm.model_config.model_config
.
# Generate some simulatred data
from hddm.simulators.hddm_dataset_generators import simulator_h_c
data, full_parameter_dict = simulator_h_c(
n_subjects=1,
n_trials_per_subject=n_samples,
model=model,
p_outlier=0.00,
conditions=None,
depends_on=None,
regression_models=None,
regression_covariates=None, # need this to make initial covariate matrix from which to use dmatrix (patsy)
group_only_regressors=False,
group_only=None,
fixed_at_default=None,
)
data
rt | response | subj_idx | v | a | z | t | theta | |
---|---|---|---|---|---|---|---|---|
0 | 2.096904 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
1 | 2.154903 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
2 | 1.862907 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
3 | 1.847907 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
4 | 1.927906 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
... | ... | ... | ... | ... | ... | ... | ... | ... |
995 | 2.260902 | 1.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
996 | 1.895906 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
997 | 1.782908 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
998 | 1.864907 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
999 | 1.812907 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
1000 rows × 8 columns
# Define the HDDM model
hddmnn_model = hddm.HDDMnn(
data,
informative=False,
include=includes,
p_outlier=0.01,
w_outlier=0.1,
model=model,
)
Supplied model_config specifies params_std_upper for z as None. Changed to 10
# Sample
hddmnn_model.sample(nmcmc, burn=500)
[-----------------100%-----------------] 1500 of 1500 complete in 103.2 sec
<pymc.MCMC.MCMC at 0x141b8e410>
The plot_caterpillar()
function below displays parameterwise,
Again use the help()
function to learn more.
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(
hddm_model=hddmnn_model,
ground_truth_parameter_dict=full_parameter_dict,
figsize=(8, 5),
columns=3,
)
plt.show()
model cartoon plot
)¶Another way to examine whether or not our recovery was satisfactory is to perform posterior predictive checks. Essentially, we are looking to simulate datasets from the trace and check whether it aligns with the ground truth participant data. This answers the question of whether or not these parameters that you recovered can actually reproduce the data.
Use the plot_posterior_predictive()
function in the plotting
module for this. It is structured just like the plot_from_data()
function, but instead of providing a dataset, you supply a hddm model.
Use the help()
function to check out all the functionality.
hddm.plotting.plot_posterior_predictive(
model=hddmnn_model,
columns=1,
groupby=["subj_idx"],
figsize=(6, 4),
value_range=np.arange(0, 5, 0.1),
plot_func=hddm.plotting._plot_func_model,
parameter_recovery_mode=True,
**{"alpha": 0.01, "ylim": 3, "samples": 200}
)
plt.show()
passing
A small note on convergence:
Note that the MCMC algorithm requires the chain to converge. There are many heuristics that help you identifying problems with convergence, such as the trace plot, auto correlation plot, and marginal posterior histogram. In the trace plots, there might be a problem if you see large jumps. In the autocorrelation plot, there might be a problem if it does not drop rapidly. The HDDMnn()
classes support the computation of the Gelman-Rubin, r-hat statistic, as you would with any hddm
model. Generally, by extracting the traces, you are free to compute any convergence statistics you want of course.
# TAKING A LOOK AT THE POSTERIOR TRACES
hddmnn_model.plot_posteriors(hddm.simulators.model_config[model]["params"])
plt.show()
Plotting v Plotting a Plotting z Plotting t Plotting theta
hddm.plotting.plot_posterior_pair(
hddmnn_model, save=False, parameter_recovery_mode=True, samples=500, figsize=(6, 6)
)
The 'h' in hddm
stands for hierarchical, so let's do it! If we have data from multiple participants and we assume that the parameters of single participants are drawn from respective group or global distributions, we can model this explicitly in hddm
by specifying is_group_model = True
.
Implicitly we are fitting a model of the following kind,
$$p(\{\theta_j\}, \{\theta_g\} | \mathbf{x}) \propto \left[ \prod_j^{J} \left[ \prod_i^{N_j} p(x_i^j | \theta_j) \right] p(\theta_j | \theta_g) \right] p( \theta_g | \theta_h )$$where (let's say for the angle model),
$\theta_j = \{v_j, a_j, z_j, t_j, \theta_j \}$, are the model parameters for subject j.
$\theta_g = \{v_g^{\mu}, a_g^{\mu}, z_g^{\mu}, t_g^{\mu}, \theta_g^{\mu}, v_g^{\sigma}, a_g^{\sigma}, z_g^{\sigma}, t_g^{\sigma}, \theta_g^{\sigma} \}$ (scary, but for completeness), are the mean and variance parameters for our group level normal distributions, and $\{ \theta_h \}$ are fixed hyperparameters.
$x_i^j = \{rt_i^j, c_i^j \}$, are the choice and reaction time of subject j during trial i.
In words, the right hand side of the equation tells us that we have a global parameter distribution with certain means and variances for each parameter (we want to figure these means and variances out), from which the subject level parameters are drawn and finally subject level datapoints follow the likelihood distribution of our ddm / angle / weibull / you name it mdoels.
# Metadata
nmcmc = 1000
model = "angle"
n_trials_per_subject = 200
n_subjects = 10
# test regressors only False
# add p_outliers to the generator !
data, full_parameter_dict = simulator_h_c(
data=None,
n_subjects=n_subjects,
n_trials_per_subject=n_trials_per_subject,
model=model,
p_outlier=0.00,
conditions=None,
depends_on=None,
regression_models=None,
regression_covariates=None,
group_only_regressors=False,
group_only=None,
fixed_at_default=None,
)
hddmnn_model = hddm.HDDMnn(
data,
model=model,
informative=False,
is_group_model=True,
include=hddm.simulators.model_config[model]["hddm_include"],
p_outlier=0.0,
)
{'v': 1.5, 'v_std': 1, 'a': 1, 'a_std': 1, 'z': 0.1, 'z_trans': 0.2, 't': 0.01, 't_std': 0.15, 'theta': 0.1, 'theta_std': 0.2} Supplied model_config specifies params_std_upper for z as None. Changed to 10
hddmnn_model.sample(
nmcmc, burn=100
) # if you want to save the model specify extra arguments --> dbname='traces.db', db='pickle'. # hddmnn_model.save('test_model')
[-----------------100%-----------------] 1000 of 1000 complete in 339.0 sec
<pymc.MCMC.MCMC at 0x14bb81390>
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(
hddm_model=hddmnn_model,
ground_truth_parameter_dict=full_parameter_dict,
figsize=(8, 5),
columns=3,
)
plt.show()
hddm.plotting.plot_posterior_predictive(
model=hddmnn_model,
columns=3,
figsize=(10, 7),
groupby=["subj_idx"],
value_range=np.arange(0, 5, 0.1),
plot_func=hddm.plotting._plot_func_model,
parameter_recovery_mode=True,
**{
"alpha": 0.01,
"ylim": 3,
"add_posterior_mean_rts": True,
"add_posterior_mean_model": True,
"add_posterior_uncertainty_rts": False,
"add_posterior_uncertainty_model": False,
"samples": 200,
"legend_fontsize": 7.0,
}
)
An important aspect of these posterior analysis, is the consideration of experiment design. We may have an experiment in which subject are exposed to a variety of conditions, such as for example different degrees of difficulty of the same task
It is often reasonable to assume that all but the conceptually relevant parameters are common across conditions.
As a by-product, such experiment designs can help us with the recovery of the constant parameters, by probing those static aspects of the model across varying kinds of datasets (driven by targeted manipulation of variable aspects of the model).
Implicitly we fit the following kind of model,
$$p( \{\theta_c \}, \theta | \mathbf{x} ) \propto \left[ \prod_c^C \left[ \prod_i^{N_i} p( x_i^c | \theta_c, \theta ) \right] p(\theta_c) \right] p(\theta)$$Where $\theta_c$ is the condition dependent part of the parameter space, and $\theta$ forms the portion of parameters which remain constant across condtions.
To give a more concrete example involving the weibull model, consider a dataset for a single participant, who went through four conditions of an experiment. Think of the conditions as manipulating the payoff structure of the experiment to incentivize / disincentivize accuracy in favor of speed. We operationalize this by treating the $a$ parameter, the initial boundary separation, as affected by the manipulation, while the rest of the parameters are constant across all experiment conditions.
The resulting model would be of the form,
$$ p( {a_c}, v, z, t, \alpha, \beta | x ) \propto \left[ \prod_c^C \left[ \prod_i^{N_c} p( x_i^c | a_c, v, z, t, \alpha, \beta) \right] p(a_c) \right] p(v, z, t, \alpha, \beta)$$# Metadata
nmcmc = 1000
model = "angle"
n_trials_per_subject = 500
# We allow the boundary conditions to vary
depends_on = {"a": ["c_one"]}
# They will depend on a fictious column 'c_one' that specifies
# levels / conditions
conditions = {"c_one": ["low", "medium", "high"]}
data, full_parameter_dict = simulator_h_c(
n_subjects=1,
n_trials_per_subject=n_trials_per_subject,
model=model,
p_outlier=0.00,
conditions=conditions,
depends_on=depends_on,
regression_models=None,
regression_covariates=None,
group_only_regressors=False,
group_only=None,
fixed_at_default=None,
)
depends_on is: {'a': ['c_one']}
# Let's check the resulting parameter vector
full_parameter_dict
{'theta': 0.7406253194726012, 'v': 1.464554358239174, 'z': 0.6206249211841304, 't': 1.534252965986117, 'a(high)': 1.0519165572885651, 'a(low)': 1.2561997135872933, 'a(medium)': 0.9265856569938499}
# Make HDDM Model
hddmnn_model = hddm.HDDMnn(
data,
model=model,
informative=False,
include=hddm.simulators.model_config[model]["hddm_include"],
p_outlier=0.0,
is_group_model=False,
depends_on=depends_on,
)
# Sample
hddmnn_model.sample(nmcmc, burn=100)
[-----------------100%-----------------] 1001 of 1000 complete in 129.4 sec
<pymc.MCMC.MCMC at 0x14f44c690>
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(
hddm_model=hddmnn_model,
ground_truth_parameter_dict=full_parameter_dict,
figsize=(8, 5),
columns=3,
)
plt.show()
hddm.plotting.plot_posterior_predictive(
model=hddmnn_model,
columns=1,
groupby=["subj_idx"],
figsize=(4, 4),
value_range=np.arange(0, 5, 0.1),
plot_func=hddm.plotting._plot_func_model,
parameter_recovery_mode=True,
**{
"alpha": 0.01,
"ylim": 3,
"add_posterior_uncertainty_rts": True,
"add_posterior_uncertainty_model": True,
"samples": 200,
}
)
plt.show()
passing
passing
passing
# Metadata
nmcmc = 1500
model = "angle"
n_subjects = 5
n_trials_per_subject = 500
data, full_parameter_dict = simulator_h_c(
n_subjects=n_subjects,
n_trials_per_subject=n_trials_per_subject,
model=model,
p_outlier=0.00,
conditions={
"c_one": ["low", "medium", "high"]
}, # , 'c_three': ['low', 'medium', 'high']},
depends_on={
"v": ["c_one"]
}, # 'theta': ['c_two']}, # 'theta': ['c_two']}, #regression_models = None, #
regression_models=None, # regression_covariates = None,
regression_covariates=None, # need this to make initial covariate matrix from which to use dmatrix (patsy)
group_only_regressors=False,
group_only=None,
fixed_at_default=None,
)
depends_on is: {'v': ['c_one']}
# Make HDDM Model
hddmnn_model = hddm.HDDMnn(
data,
model=model,
informative=False,
include=hddm.simulators.model_config[model]["hddm_include"],
p_outlier=0.0,
is_group_model=True,
depends_on={"v": "c_one"},
)
hddmnn_model.sample(nmcmc, burn=100)
[-----------------100%-----------------] 1500 of 1500 complete in 919.0 sec
<pymc.MCMC.MCMC at 0x14e324a90>
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(
hddm_model=hddmnn_model,
ground_truth_parameter_dict=full_parameter_dict,
figsize=(8, 8),
columns=3,
)
plt.show()
hddm.plotting.plot_posterior_predictive(
model=hddmnn_model,
columns=2, # groupby = ['subj_idx'],
figsize=(8, 6),
value_range=np.arange(1, 2.5, 0.1),
plot_func=hddm.plotting._plot_func_model,
parameter_recovery_mode=True,
**{
"alpha": 0.01,
"ylim": 3,
"add_posterior_uncertainty_rts": True,
"add_posterior_uncertainty_model": True,
"samples": 200,
"legend_fontsize": 7,
}
)
plt.show()
passing passing passing passing passing
passing passing passing passing passing
passing passing passing passing passing
This section provides a simple working example using the Neural Networks with the Regression backend. The regression back-end allows linking parameters to trial-by-trial covariates via a (general) linear model.
# Metadata
nmcmc = 1000
model = "angle"
n_samples_by_subject = 500
from hddm.simulators.hddm_dataset_generators import simulator_h_c
data, full_parameter_dict = simulator_h_c(
n_subjects=5,
n_samples_by_subject=n_samples_by_subject,
model=model,
p_outlier=0.00,
conditions=None,
depends_on=None,
regression_models=["t ~ 1 + covariate_name", "v ~ 1 + covariate_name"],
regression_covariates={"covariate_name": {"type": "continuous", "range": (0, 1)}},
group_only_regressors=False,
group_only=None,
fixed_at_default=None,
)
# Set up the regressor a regressor:
reg_model_v = {"model": "v ~ 1 + covariate_name", "link_func": lambda x: x}
reg_model_t = {"model": "t ~ 1 + covariate_name", "link_func": lambda x: x}
reg_descr = [reg_model_t, reg_model_v]
# Make HDDM model
hddmnn_reg = hddm.HDDMnnRegressor(
data,
reg_descr,
include=hddm.simulators.model_config[model]["hddm_include"],
model=model,
informative=False,
p_outlier=0.0,
)
Supplied model_config specifies params_std_upper for z as None. Changed to 10
# Sample
hddmnn_reg.sample(nmcmc, burn=100)
[-----------------100%-----------------] 1001 of 1000 complete in 369.4 sec
<pymc.MCMC.MCMC at 0x149a07890>
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(
hddm_model=hddmnn_reg,
ground_truth_parameter_dict=full_parameter_dict,
figsize=(8, 8),
columns=3,
)
plt.show()
You can read more about stimulus coding in the documentation.
Here just an example.
# Metadata
nmcmc = 300
model = "ddm"
n_samples_by_condition = 500
split_param = "v"
sim_data_stimcoding, parameter_dict = hddm.simulators.simulator_stimcoding(
model=model, split_by=split_param, drift_criterion=0.3, n_trials_per_condition=500
)
sim_data_stimcoding
rt | response | stim | v | a | z | t | subj_idx | |
---|---|---|---|---|---|---|---|---|
0 | 3.190470 | 1.0 | 1 | 0.834704 | 2.426857 | 0.417932 | 1.507448 | none |
1 | 3.942454 | 1.0 | 1 | 0.834704 | 2.426857 | 0.417932 | 1.507448 | none |
2 | 4.186436 | 1.0 | 1 | 0.834704 | 2.426857 | 0.417932 | 1.507448 | none |
3 | 2.205442 | 1.0 | 1 | 0.834704 | 2.426857 | 0.417932 | 1.507448 | none |
4 | 4.669401 | 1.0 | 1 | 0.834704 | 2.426857 | 0.417932 | 1.507448 | none |
... | ... | ... | ... | ... | ... | ... | ... | ... |
495 | 11.207737 | 0.0 | 2 | -0.234704 | 2.426857 | 0.417932 | 1.507448 | none |
496 | 10.334385 | 1.0 | 2 | -0.234704 | 2.426857 | 0.417932 | 1.507448 | none |
497 | 7.077227 | 0.0 | 2 | -0.234704 | 2.426857 | 0.417932 | 1.507448 | none |
498 | 8.740107 | 1.0 | 2 | -0.234704 | 2.426857 | 0.417932 | 1.507448 | none |
499 | 2.621444 | 0.0 | 2 | -0.234704 | 2.426857 | 0.417932 | 1.507448 | none |
1000 rows × 8 columns
parameter_dict
{'v': -0.5347036843723503, 'a': 2.426856838254428, 'z': 0.4179319892615798, 't': 1.5074477307220377, 'dc': 0.3}
hddmnn_model = hddm.HDDMnnStimCoding(
sim_data_stimcoding,
include=hddm.simulators.model_config[model]["hddm_include"],
model=model,
stim_col="stim",
p_outlier=0.0,
split_param=split_param,
informative=False,
drift_criterion=True,
)
hddmnn_model.sample(nmcmc, burn=100)
[-----------------100%-----------------] 300 of 300 complete in 32.4 sec
<pymc.MCMC.MCMC at 0x14e56a850>
hddmnn_model.gen_stats()
mean | std | 2.5q | 25q | 50q | 75q | 97.5q | mc err | |
---|---|---|---|---|---|---|---|---|
v | -0.539954 | 0.0155719 | -0.572298 | -0.549469 | -0.540782 | -0.528279 | -0.508065 | 0.00112267 |
a | 2.49136 | 0.00895635 | 2.47002 | 2.48866 | 2.49422 | 2.49799 | 2.49988 | 0.000770613 |
z | 0.4031 | 0.0118855 | 0.37898 | 0.397112 | 0.40288 | 0.409666 | 0.431941 | 0.000958058 |
t | 1.48852 | 0.035718 | 1.41497 | 1.46917 | 1.48828 | 1.51355 | 1.56092 | 0.00286112 |
dc | 0.348321 | 0.0202259 | 0.30703 | 0.334826 | 0.349703 | 0.361927 | 0.386422 | 0.00170049 |
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(
hddm_model=hddmnn_model,
ground_truth_parameter_dict=parameter_dict,
figsize=(8, 5),
columns=3,
)
plt.show()
NOTE:
The hddm.plotting.plot_posterior_predictive()
does not yet accept stimcoding data. This will be updated as soon as possible.
A crucial exercise in statistical modeling concern model comparison.
We are going to look at model recovery, in this section: Attempt to recover which model generated a given dataset from a set of candidate models.
For the little model recovery study we conduct here, we generate data from the weibull model and fit the data once each to the weibull, angle and ddm models.
We inspect the fits visually and then use the DIC (Deviance information criterion, lower is better :)), to check if we can recover the true model.
# Metadata
model = "weibull"
n_samples = 300
# test regressors only False
# add p_outliers to the generator !
data, full_parameter_dict = simulator_h_c(
n_subjects=1,
n_samples_by_subject=n_samples,
model=model,
p_outlier=0.00,
conditions=None,
depends_on=None,
regression_models=None,
regression_covariates=None,
group_only_regressors=False,
group_only=None,
fixed_at_default=None,
)
data
rt | response | subj_idx | v | a | z | t | alpha | beta | |
---|---|---|---|---|---|---|---|---|---|
0 | 4.204582 | 0.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
1 | 4.269577 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
2 | 4.404568 | 0.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
3 | 2.960620 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
4 | 2.223596 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
95 | 2.304595 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
96 | 3.067625 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
97 | 2.379594 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
98 | 3.991597 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
99 | 2.904617 | 0.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
100 rows × 9 columns
# Now we fit for each model:
hddmnn_model_weibull = hddm.HDDMnn(
data,
informative=False,
model="weibull",
p_outlier=0.0,
include=hddm.simulators.model_config["weibull_cdf"]["hddm_include"],
is_group_model=False,
)
hddmnn_model_angle = hddm.HDDMnn(
data,
model="angle",
informative=False,
p_outlier=0.0,
include=hddm.simulators.model_config["angle"]["hddm_include"],
is_group_model=False,
)
hddmnn_model_ddm = hddm.HDDMnn(
data,
informative=False,
model="ddm",
p_outlier=0.0,
include=hddm.simulators.model_config["ddm"]["hddm_include"],
is_group_model=False,
)
nmcmc = 1000
hddmnn_model_weibull.sample(nmcmc, burn=200)
hddmnn_model_angle.sample(nmcmc, burn=200)
hddmnn_model_ddm.sample(nmcmc, burn=200)
[-----------------100%-----------------] 1000 of 1000 complete in 23.0 sec
<pymc.MCMC.MCMC at 0x14d606c90>
Posterior Predictive: Do the 'Posterior Models' also make sense?
# WEIBULL
hddm.plotting.plot_posterior_predictive(
model=hddmnn_model_weibull,
columns=1,
groupby=["subj_idx"],
figsize=(4, 4),
value_range=np.arange(0, 5, 0.1),
plot_func=hddm.plotting._plot_func_model,
parameter_recovery_mode=True,
**{
"alpha": 0.01,
"ylim": 3,
"add_posterior_uncertainty_model": True,
"add_posterior_uncertainty_rts": False,
"add_posterior_mean_rts": True,
"samples": 200,
}
)
plt.show()
passing
# ANGLE
hddm.plotting.plot_posterior_predictive(
model=hddmnn_model_angle,
columns=1,
groupby=["subj_idx"],
figsize=(4, 4),
value_range=np.arange(0, 5, 0.1),
plot_func=hddm.plotting._plot_func_model,
parameter_recovery_mode=False,
**{
"alpha": 0.01,
"ylim": 3,
"add_posterior_uncertainty_model": True,
"add_posterior_uncertainty_rts": False,
"add_posterior_mean_rts": True,
"samples": 200,
}
)
plt.show()
# DDM
hddm.plotting.plot_posterior_predictive(
model=hddmnn_model_ddm,
columns=1,
groupby=["subj_idx"],
figsize=(4, 4),
value_range=np.arange(0, 5, 0.1),
plot_func=hddm.plotting._plot_func_model,
parameter_recovery_mode=False,
**{
"alpha": 0.01,
"ylim": 3,
"add_posterior_uncertainty_model": True,
"add_posterior_uncertainty_rts": False,
"add_posterior_mean_rts": True,
"samples": 200,
}
)
plt.show()
hddmnn_model_weibull.dic
414.65114936828616
hddmnn_model_angle.dic
415.8001557922363
hddmnn_model_ddm.dic
418.04479835510256
Fingers crossed (this was a random run after all), the DIC usually gives us a result that conforms with the intuition we get from looking at the model plots.
# Metadata
nmcmc = 1000
burn = 500
model = "angle"
# Load one of the datasets shipping with HDDM
cav_data = hddm.load_csv(hddm.__path__[0] + "/examples/cavanagh_theta_nn.csv")
cav_data
subj_idx | stim | rt | response | theta | dbs | conf | |
---|---|---|---|---|---|---|---|
0 | 0 | LL | 1.210 | 1.0 | 0.656275 | 1 | HC |
1 | 0 | WL | 1.630 | 1.0 | -0.327889 | 1 | LC |
2 | 0 | WW | 1.030 | 1.0 | -0.480285 | 1 | HC |
3 | 0 | WL | 2.770 | 1.0 | 1.927427 | 1 | LC |
4 | 0 | WW | 1.140 | 0.0 | -0.213236 | 1 | HC |
... | ... | ... | ... | ... | ... | ... | ... |
3983 | 13 | LL | 1.450 | 0.0 | -1.237166 | 0 | HC |
3984 | 13 | WL | 0.711 | 1.0 | -0.377450 | 0 | LC |
3985 | 13 | WL | 0.784 | 1.0 | -0.694194 | 0 | LC |
3986 | 13 | LL | 2.350 | 0.0 | -0.546536 | 0 | HC |
3987 | 13 | WW | 1.250 | 1.0 | 0.752388 | 0 | HC |
3988 rows × 7 columns
hddmnn_model_cav = hddm.HDDMnn(
cav_data,
model=model,
informative=False,
include=hddm.simulators.model_config[model]["hddm_include"],
p_outlier=0.05,
is_group_model=False,
depends_on={"v": "stim"},
)
hddmnn_model_cav.sample(nmcmc, burn=burn)
[-----------------100%-----------------] 1000 of 1000 complete in 243.4 sec
<pymc.MCMC.MCMC at 0x154a35350>
hddm.plotting.plot_posterior_predictive(
model=hddmnn_model_cav,
columns=1,
figsize=(4, 4),
value_range=np.arange(0, 5, 0.1),
plot_func=hddm.plotting._plot_func_model,
parameter_recovery_mode=False,
**{
"alpha": 0.01,
"ylim": 3,
"add_posterior_uncertainty_model": True,
"add_posterior_uncertainty_rts": False,
"add_posterior_mean_rts": True,
"samples": 200,
}
)
plt.show()
hddmnn_model_cav = hddm.HDDMnn(
cav_data,
model=model,
informative=False,
include=hddm.simulators.model_config[model]["hddm_include"],
is_group_model=True,
p_outlier=0.05,
)
hddmnn_model_cav.sample(nmcmc, burn=burn)
[-----------------100%-----------------] 1000 of 1000 complete in 471.3 sec
<pymc.MCMC.MCMC at 0x14cca5090>
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(hddm_model=hddmnn_model_cav, figsize=(8, 8), columns=3)
plt.show()
hddm.plotting.plot_posterior_predictive(
model=hddmnn_model_cav,
columns=3,
figsize=(10, 10),
value_range=np.arange(0, 5, 0.1),
plot_func=hddm.plotting._plot_func_model,
parameter_recovery_mode=False,
**{
"alpha": 0.01,
"ylim": 3,
"add_posterior_uncertainty_model": True,
"add_posterior_uncertainty_rts": False,
"add_posterior_mean_rts": True,
"samples": 200,
"legend_fontsize": 7,
"subplots_adjust": {"top": 0.9, "hspace": 0.3, "wspace": 0.3},
}
)
plt.show()
This is just an example. The angle model might not be the best choice here, and we are moreover ignoring the supplied conditions.
The network_inspectors
module allows you to inspect the LANs directly.
You can use the hddm.network_inspectors.get_torch_mlp()
function to access network predictions.
model = "angle"
lan_angle = hddm.network_inspectors.get_torch_mlp(model=model)
Let's predict some likelihoods !
# Make some random parameter set
parameter_df = hddm.simulators.make_parameter_vectors_nn(
model=model, param_dict=None, n_parameter_vectors=1
)
parameter_matrix = np.tile(np.squeeze(parameter_df.values), (200, 1))
# Initialize network input
network_input = np.zeros(
(parameter_matrix.shape[0], parameter_matrix.shape[1] + 2)
) # Note the + 2 on the right --> we append the parameter vectors with reaction times (+1 columns) and choices (+1 columns)
# Add reaction times
network_input[:, -2] = np.linspace(0, 3, parameter_matrix.shape[0])
# Add choices
network_input[:, -1] = np.repeat(np.random.choice([-1, 1]), parameter_matrix.shape[0])
# Convert to float
network_input = network_input.astype(np.float32)
# Show example output
print(lan_angle(network_input)[:10]) # printing the first 10 outputs
print(lan_angle(network_input).shape) # original shape of output
[[-2.9323568] [ 2.078088 ] [ 0.4104141] [-0.5943402] [-1.1136726] [-1.6901499] [-2.3512228] [-3.080151 ] [-3.8215086] [-4.4257374]] (200, 1)
HDDM provides two plotting function to investigate the network outputs directly. The kde_vs_lan_likelihoods()
plot and the lan_manifold()
plot.
kde_vs_lan_likelihoods()
¶The kde_vs_lan_likelihoods()
plot allows you to check the likelihoods produced by a LAN against Kernel Density Estimates (KDEs) from model simulations.
You can supply a panda DataFrame
that holds parameter vectors as rows.
# Make some parameters
parameter_df = hddm.simulators.make_parameter_vectors_nn(
model=model, param_dict=None, n_parameter_vectors=10
)
parameter_df
v | a | z | t | theta | |
---|---|---|---|---|---|
0 | 2.149626 | 1.684902 | 0.232222 | 0.641663 | -0.070030 |
1 | 1.817911 | 0.776330 | 0.535083 | 0.006625 | 1.069452 |
2 | -0.908428 | 0.654107 | 0.301445 | 1.560911 | 0.396448 |
3 | -0.022136 | 1.140235 | 0.479664 | 0.757727 | 1.316409 |
4 | 2.281230 | 0.366558 | 0.409224 | 1.908211 | 1.059872 |
5 | 1.067632 | 1.228020 | 0.337573 | 1.447155 | 0.083665 |
6 | 2.022131 | 1.254037 | 0.262336 | 0.416462 | 0.512724 |
7 | -1.974657 | 0.793536 | 0.791707 | 0.591319 | 1.036441 |
8 | -2.002436 | 1.382722 | 0.442411 | 0.074192 | 0.356522 |
9 | -2.757462 | 0.402900 | 0.738999 | 0.755093 | 1.334423 |
hddm.network_inspectors.kde_vs_lan_likelihoods(
parameter_df=parameter_df, model=model, cols=3, n_samples=2000, n_reps=10, show=True
)
1 of 10 2 of 10 3 of 10 4 of 10 5 of 10 6 of 10 7 of 10 8 of 10 9 of 10 10 of 10
lan_manifold()
¶Lastly, you can use the lan_manifold()
plot to investigate the LAN likelihoods over a range of parameters.
The idea is to use a base parameter vector and vary one of the parameters in a prespecificed range.
This plot can be informative if you would like to understand better how a parameter affects model behavior.
# Make some parameters
parameter_df = hddm.simulators.make_parameter_vectors_nn(
model=model, param_dict=None, n_parameter_vectors=1
)
parameter_df
v | a | z | t | theta | |
---|---|---|---|---|---|
0 | -2.218164 | 0.889863 | 0.254979 | 0.707028 | 0.040745 |
# Now plotting
hddm.network_inspectors.lan_manifold(
parameter_df=parameter_df,
vary_dict={"v": np.linspace(-2, 2, 20)},
model=model,
n_rt_steps=300,
fig_scale=1.0,
max_rt=5,
save=True,
show=True,
)
Using only the first row of the supplied parameter array !
Hopefully this tutorial proves as a useful starting point for your application.