Psychometric Function Estimation
Warning
This example includes deprecated functions and arguments. We plan to update this page following the up-to-date usage. Sorry for the inconvenience.
Let’s start with psychometric functions as an example. The goal of the function is to figure out whether a subject can perceive a signal with varying levels of magnitude. The function has one design variable for the intensity of a stimulus, \(x\); the model has four model parameters: guess rate \(\gamma\), lapse rate \(\delta\), threshold \(\alpha\), and slope \(\beta\).
A simple diagram for the Psychometric function.
In this example, let’s use the logistic function for the model’s shape. Then, the model can compute the probability of a subject to perceive the given stimulus with the following equation:
For this example, let’s assume the true parameters as \(\gamma = 0.5\), \(\delta = 0.04\), \(\alpha = 20\), and \(\beta = 1.5\).
# Define true parameters
GR_TRUE = 0.5
LR_TRUE = 0.04
TH_TRUE = 20
SL_TRUE = 1.5
Preparing grids
To make grids for designs and parameters, you should define two dictionaries
that contain singles grids for all designs and all parameters, respectively.
In this example, we will fix the guess_rate
to 0.5 and lapse_rate
to 0.04.
import numpy as np
designs = {
'stimulus': np.linspace(20 * np.log10(.05), 20 * np.log10(400), 120)
}
params = {
'guess_rate': [0.5],
'lapse_rate': [0.04],
'threshold': np.linspace(20 * np.log10(.1), 20 * np.log10(200), 200),
'slope': np.linspace(0, 10, 200)
}
Using pre-defined classes
To use the predefined classes for specific task and models, you can use it
with adopy.tasks.<task_name>, e.g., adopy.tasks.psi
.
from adopy.tasks.psi import ModelLogistic, EnginePsi
model = ModelLogistic()
engine = EnginePsi(model=model, designs=designs, params=params)
Using compute() method of the model instance, you can compute the probability for a subject to succeed to perceive a signal.
model.compute(stimulus=10, guess_rate=0.5, lapse_rate=0.04,
threshold=10, slope=0.5)
from scipy.stats import bernoulli
p_obs = model.compute(stimulus=d_opt['stimulus'],
guess_rate=gr_true, lapse_rate=lr_true,
threshold=th_true, slope=sl_true)
y_obs = bernoulli.rvs(p_obs)
d_opt = e.get_design()
Using self-defined classes
Instead of using pre-defined classes, they can be implemented as Task
and Model
objects by the
codes below:
import numpy as np
from adopy import Task, Model
task_psi = Task(name='Psi', key='psi', design=['stimulus'])
def inv_logit(x):
return np.divide(1, 1 + np.exp(-x))
def func_logistic(stimulus, guess_rate, lapse_rate, threshold, slope):
return guess_rate + (1 - guess_rate - lapse_rate) * inv_logit(slope * (stimulus - threshold))
model_log = Model(name='Logistic', task=task_psi,
param=['guess_rate', 'lapse_rate', 'threshold', 'slope'],
func=func_logistic)
Then, you can compute the probability using compute
method in the
model object.
print(model_log.compute(stimulus=10, guess_rate=0.5, lapse_rate=0.04, threshold=10, slope=0.5))
print(model_log.compute(stimulus=15, guess_rate=0.5, lapse_rate=0.04, threshold=10, slope=0.5))
print(model_log.compute(stimulus=5, guess_rate=0.5, lapse_rate=0.04, threshold=10, slope=0.5))
Now, if you want to use an ADO engine for the task and the model,
from adopy import Engine
engine_psi = Engine(task=task_psi, model=model_log,
designs=designs, params=params, y_obs=[0, 1])
With the Engine
instance, you can get the optimal design:
d_opt = engine_psi.get_design()
Assuming \(\gamma = 0.5\), \(\delta = 0.04\), \(\alpha = 20\) and \(\beta = 1.5\), you can get the probability of perceiving the stimulus with model_log.compute.
from scipy.stats import bernoulli
p_obs = model_log.compute(stimulus=d_opt['stimulus'], guess_rate=gr_true, lapse_rate=lr_true,
threshold=th_true, slope=sl_true)
y_obs = bernoulli.rvs(p_obs)
Lastly, using the optimal design and the corresponding response, the Engine instance can update its posterior distributions on parameters.
engine_psi.update(d_opt, y_obs)