Custom Tutorial 1: Creating a custom model version of the simple DDM model

If you have not done so already, install PyBEAM. Go to the PyBEAM github (https://github.com/murrowma/pybeam) for instructions.

Once PyBEAM is installed, we can create our first model. First, download the folder “custom_model_template.zip” and unzip it. Then, create a new folder anywhere on your device and copy the files located in folder “custom_model_template” to it. There will be three files you need to copy: “setup.py”, “functions.pyx”, and “model.pyx”. The only file that needs to be modified is “model.pyx”. DO NOT TOUCH THE OTHER TWO FILES.

Un-edited, “model.pyx” will have the following appearance:

"""
# Import Python packages and C libraries.

"""

import cython
cimport cython

from libc.math cimport pow, exp, log, sin, cos

"""
# Model parameter functions.

"""

# total number of input parameters
DEF N_phi = 0

# function for the non-decision time
cdef double non_decision_time(double phi[N_phi]):
    return 0.0

# function for the relative start point
cdef double relative_start(double phi[N_phi]):
    return 0.0

# function for the drift rate
cdef double drift(double phi[N_phi], double x, double t):
    return 0.0

# function for the diffusion rate
cdef double diffusion(double phi[N_phi], double x, double t):
    return 0.0

# function for the upper decision threshold
cdef double upper_decision_threshold(double phi[N_phi], double t):
    return 0.0

# function for the lower decision threshold
cdef double lower_decision_threshold(double phi[N_phi], double t):
    return 0.0

# function for the contamination strength
cdef double contamination_strength(double phi[N_phi]):
    return 0.0

# function for the contamination probability
cdef double contamination_probability(double phi[N_phi], double t):
    return 0.0

"""
# Function to modify time step with unusual likelihood functions.

"""

# function to modify the time step
cdef double modify_dt(double phi[N_phi], double t):
    return 1.0

The first two lines,

import cython
cimport cython

import Cython. The next line,

from libc.math cimport pow, exp, log, sin, cos

imports useful C math functions. The first, pow, is the C verison of raising a number to a power. It’s syntax is pow(a, b), which is the python equivalent of a**b. Whenever using exponents, always use pow because it is substantially faster when compiled than the Python equivalent.

The next several functions (exp, log, sin, and cos) all act the same as their python equivalents. These functions are imported since they are commonly used in models in the literature, but other functions can be included by adding them to this list. For example, if you want a tangent function, you would add tan to the list, giving,

from libc.math cimport pow, exp, log, sin, cos, tan

Note that, if functions from other C libraries are desired, additional import lines can be added here.

The second block of code is where your model is defined. First, you must set “N_phi” equal to the number of parameters your model contains. If, for example, you were defining the simple DDM with contamination described in precoded Tutorial 1a, you would input N_phi = 8 here since there were 8 parameters: tnd, w, mu, sigma (the scaling parameter), b, g (contamination strength), gl, and gu (bounds for uniform contamination distribution).

# total number of input parameters
DEF N_phi = 8

Next we define outputs for our model functions. There are eight model related functions that can be modified to fit your needs (note that ONLY the outputs, NOT the inputs, can be changed). The first, non_decision_time, sets the model’s non-decision time. It accepts an input of phi, the list of all the parameters in your model (something we will use later in the tutorial). Since this function only accepts phi as an input, it can only output a constant value. In the case of the simple DDM, we get the following function:

# function for the non-decision time
cdef double non_decision_time(double phi[N_phi]):
    return phi[0] # tnd

Note how we call the return function. Since all parameters are located in list phi, we tell the function to reference those parameters, or some operation on them. Since the non-decision time is only a constant, we have it output a constant phi[0].

The second, relative start, sets the model’s relative start point. Like non_decision_time, it only accepts an input of phi, meaning it outputs a constant value only. In the case of the simple DDM, we receive the following:

# function for the non-decision time
cdef double non_decision_time(double phi[N_phi]):
    return phi[1] # w

The next two functions, drift and diffusion, set the model’s drift (v(x,t) in PyBEAM’s publication) and diffusion (D(x,t) in PyBEAM’s publication) rates (the scaling parameter in this case). They both accept three inputs: phi, x, and t. phi, as before, contains the parameters used by your model. The second, x, is the spatial coordinate and corresponds to the accumulator state. The last, t, is time. Thus, both the drift and diffusion rates can have time and spatial dependencies. Since both are constants, we write the functions return statemensts as follows:

# function for the drift rate
cdef double drift(double phi[N_phi], double x, double t):
    return phi[2] # mu

# function for the diffusion rate
cdef double diffusion(double phi[N_phi], double x, double t):
    return phi[3] # sigma

The next two functions, upper_decision_threshold and lower_decision_threshold, set the locations of the decision thresholds. These two functions accepts inputs of phi and t, meaning that the thresholds can be time dependent. In the case of the simple DDM, they are constant, giving us the following:

# function for the upper decision threshold
cdef double upper_decision_threshold(double phi[N_phi], double t):
    return phi[4] # b

# function for the lower decision threshold
cdef double lower_decision_threshold(double phi[N_phi], double t):
    return -phi[4] # b

Note here that the upper decision threshold must always be positive, while the lower decision threshold must always be negative.

The last two model functions are for adding a contamination model. The first of these, contamination_strength, sets the amount of contamination. If it returns 0, your model has no contamination; if it returns 1, your model is only contamination. This function only accepts inputs of phi, meaning it can only output a constant value. In our case, we do the following:

# function for the contamination strength
cdef double contamination_strength(double phi[N_phi]):
    return phi[5] # g

The second contamination function is contamination_probability. This function returns the probability of a contamination event occurring at time t. It accepts inputs of phi and t. In this case, a uniform contamination is added to the model, so a uniform distribution is required for this function. This gives us the following:

# function for the contamination probability
cdef double contamination_probability(double phi[N_phi], double t):
    cdef double gl = phi[6] # lower bound of uniform distribution
    cdef double gu = phi[7] # upper bound of uniform distribution
    cdef double pg = 0.0    # contamination probability

    if (t >= gl) and (t <= gu):
        pg = 1.0/(gu - gl)

    return pg

Note that we output the entire uniform distribution here. PyBEAM scales it in the background by the contamination strength and splits it between the upper and lower threshold likelihood functions automatically.

The last block of code provides a function for modifying the time step. In general this should not be touched, but it allows for your time step to be made altered in regions where it’s hard to accurately obtain the likelihood function. For example, if your drift/threshold functions have discontinuities in them, it can sometimes be helpful to add a time step modification locally around it. If it outputs 1, the time step is not modified; if it recieves a value of less than 1, the time step will get smaller; if it outputs a value more than 1; the time step will grow. For example, if it returns 0.5, the time step will be one-half of its normal value.

This gives us the following model.pyx function for the simple DDM:

"""
# Import Python packages and C libraries.

"""

import cython
cimport cython

from libc.math cimport pow, exp, log, sin, cos

"""
# Model parameter functions.

"""

# total number of input parameters
DEF N_phi = 8

# function for the non-decision time
cdef double non_decision_time(double phi[N_phi]):
    return phi[0]

# function for the relative start point
cdef double relative_start(double phi[N_phi]):
    return phi[1]

# function for the drift rate
cdef double drift(double phi[N_phi], double x, double t):
    return phi[2]

# function for the diffusion rate
cdef double diffusion(double phi[N_phi], double x, double t):
    return phi[3]

# function for the upper decision threshold
cdef double upper_decision_threshold(double phi[N_phi], double t):
    return phi[4]

# function for the lower decision threshold
cdef double lower_decision_threshold(double phi[N_phi], double t):
    return -phi[4]

# function for the contamination strength
cdef double contamination_strength(double phi[N_phi]):
    return phi[4]

# function for the contamination probability
cdef double contamination_probability(double phi[N_phi], double t):
    cdef double g = phi[6]
    cdef double gl = phi[7]
    cdef double gu = phi[8]
    cdef double pg = 0.0

    if (t >= gl) and (t <= gu):
        pg = 1.0/(gu - gl)

    return pg

"""
# Function to modify time step with unusual likelihood functions.

"""

# function to modify the time step
cdef double modify_dt(double phi[N_phi], double t):
    return 1.0

Once you have finished the model.pyx file, you now need to compile the function. To do so, open terminal, anaconda prompt, or command prompt (depending on OS), and navigate to the folder containing your model files. Then, run the following line of code:

python setup.py build_ext --inplace

Note that the program you have compiled is unique to that particular version of Python. For example, if you compile it with Python 3.9, you will need to re-compile it if you upgrade to Python 3.11.

And with that, your model is ready to be used with PyBEAM!