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\).

../_images/Psychometricfn.svg — 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:

\[\Psi(x \mid \alpha, \beta, \gamma, \delta) = \gamma + (1 - \gamma - \delta) \; \sigma\big( \beta (x - \alpha) \big) \quad \text{where } \sigma(x) = \frac{1}{1 + e^{-x}}\]

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)