In [1]:
import pymc3 as pm
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
%qtconsole --colors=linux
plt.style.use('ggplot')

from matplotlib import gridspec
from theano import tensor as tt
from scipy import stats

Chapter 12 - Psychophysical functions

A psychophysical function, showing the a sigmoid or S-shaped relationship between stimulus intensity and choice behavior. Important theoretical measures known as the point of subjective equality (PSE) and just noticeable difference (JND) are highlighted.

In [2]:
def logit(x):
    return 1/(1+np.exp(-x))
def invlogit(x):
    return np.log(x/(1-x))

x = np.linspace(-3.75, 3.75, 100)
fig, ax = plt.subplots(figsize=(9, 6))
x1 = invlogit(.5)
x2 = invlogit(.84)

plt.plot(x, logit(x), 'k', linewidth=2)
plt.plot([x1, x1], [0, .5], color='k', linestyle='--', linewidth=1)
plt.plot([-3.75, x1], [.5, .5], color='k', linestyle='--', linewidth=1)
plt.plot([x2, x2], [0, .84], color='k', linestyle='--', linewidth=1)
plt.plot([-3.75, x2], [.84, .84], color='k', linestyle='--', linewidth=1)

plt.scatter(x1, .5, c='k', s=75)
plt.text(x1+.5, .5-.02, "PSE", horizontalalignment='center', fontsize=20)
plt.plot([x1, x2], [.5/2+.06, .5/2+.06], color='k', linewidth=1)
plt.text((x1+x2)/2, .5/2, "JND", horizontalalignment='center', fontsize=20)

ax.set_xticks(())
ax.set_yticks((0, .5, .84, 1))
ax.set_yticklabels(('0', '0.5', '0.84', '1'), fontsize=12)
plt.xlim(-3.75, 3.75)
plt.ylim(0, 1)
plt.xlabel('Stimulus Intensity', fontsize=15)
plt.ylabel('Respone Probability', fontsize=15)
plt.show()

Data for all 8 subjects, showing the proportion of “long” responses as a function of test interval duration.

In [3]:
x = pd.read_csv("./data/data_x.txt", sep='\t', header=None)
n = pd.read_csv("./data/data_n.txt", sep='\t', header=None)
r = pd.read_csv("./data/data_r.txt", sep='\t', header=None)
rprop = pd.read_csv("./data/data_rprop.txt", sep='\t', header=None)

xmean = np.array([318.888, 311.0417, 284.4444, 301.5909, 
                  296.2000, 305.7692, 294.6429, 280.3571])
nstim = np.array([27, 24, 27, 22, 25, 26, 28, 28])
nsubjs = 8

fig = plt.figure(figsize=(16, 8))
fig.text(0.5, -0.02, 'Test Interval (ms)', ha='center', fontsize=20)
fig.text(-0.02, 0.5, 'Proportion of Long Responses', va='center', rotation='vertical', fontsize=20)
gs = gridspec.GridSpec(2, 4)

for ip in np.arange(nsubjs):
    ax = plt.subplot(gs[ip])
    xp = np.array(x.iloc[ip, :])
    yp = np.array(rprop.iloc[ip, :])
    ax.scatter(xp, yp, marker='s', alpha=.5)
    plt.axis([190, 410, -.1, 1.1])
    plt.yticks((0, .5, .84, 1))
    plt.title('Subject %s'%(nsubjs))

plt.tight_layout()
plt.show()

The psychometric function here is a logistic function with parameters $\alpha_i$ and $\beta_i$:

$$ \theta_{ij} = \frac{1}{1+\text{exp}\{-[\alpha_{i}+\beta_{i}(x_{ij}-\bar x_{i})]\}}$$


$$\text{or}$$
$$ \text{logit}(\theta_{ij}) = \alpha_{i}+\beta_{i}(x_{ij}-\bar x_{i})$$

12.1 Psychophysical functions

$$ r_{ij} \sim \text{Binomial}(\theta_{ij},n_{ij})$$$$ \text{logit}(\theta_{ij}) = \alpha_{i}+\beta_{i}(x_{ij}-\bar x_{i})$$$$ \alpha_{i} \sim \text{Gaussian}(\mu_{\alpha},\sigma_{\alpha})$$$$ \beta_{i} \sim \text{Gaussian}(\mu_{\beta},\sigma_{\beta})$$$$ \mu_{\alpha} \sim \text{Gaussian}(0,0.001)$$$$ \mu_{\beta} \sim \text{Gaussian}(0,0.001)$$$$ \sigma_{\alpha} \sim \text{Uniform}(0,1000)$$$$ \sigma_{\beta} \sim \text{Uniform}(0,1000)$$
In [4]:
xij_tmp = x.values
nij_tmp = n.values
rij_tmp = r.values
tmp,nstim2 = np.shape(xij_tmp)

xmeanvect = np.repeat(xmean, nstim2)
sbjidx = np.repeat(np.arange(nsubjs), nstim2)

# remove nans
validmask = np.isnan(xij_tmp.flatten())==False
xij2 = xij_tmp.flatten()
nij2 = nij_tmp.flatten()
rij2 = rij_tmp.flatten()

xij = xij2[validmask]
nij = nij2[validmask]
rij = rij2[validmask]
xvect = xmeanvect[validmask]
sbjid = sbjidx[validmask]

Helper function

In [5]:
def Phi(x):
    # probit transform 
    return 0.5 + 0.5 * pm.math.erf(x/pm.math.sqrt(2))
def tlogit(x):
    return 1/(1+tt.exp(-x))
In [8]:
with pm.Model() as model1:
    sigma_a = pm.Uniform('sigma_a', lower=0, upper=1000)
    sigma_b = pm.Uniform('sigma_b', lower=0, upper=1000)
    mu_a = pm.Normal('mu_a', mu=0, tau=.001)
    mu_b = pm.Normal('mu_b', mu=0, tau=.001)
    alpha = pm.Normal('alpha', mu=mu_a, sd=sigma_a, shape=nsubjs)
    beta = pm.Normal('beta', mu=mu_b, sd=sigma_b, shape=nsubjs)
    
    linerpredi = alpha[sbjid] + beta[sbjid]*(xij-xvect)
    thetaij=pm.Deterministic('thetaij', tlogit(linerpredi))
    
    rij_ = pm.Binomial('rij', p=thetaij, n=nij, observed=rij)
    trace1 = pm.sample(1e4, njobs=2, init='advi+adapt_diag')
    
pm.traceplot(trace1, varnames=['alpha', 'beta']);
Auto-assigning NUTS sampler...
Initializing NUTS using advi+adapt_diag...
Average Loss = 320.16:  37%|███▋      | 73973/200000 [00:13<00:20, 6029.83it/s]
Convergence archived at 74300
Interrupted at 74,300 [37%]: Average Loss = 370.08
100%|██████████| 10500/10500.0 [00:28<00:00, 368.88it/s]
In [9]:
fig = plt.figure(figsize=(16, 8))
fig.text(0.5, -0.02, 'Test Interval (ms)', ha='center', fontsize=20)
fig.text(-0.02, 0.5, 'Proportion of Long Responses', va='center', rotation='vertical', fontsize=20)
gs = gridspec.GridSpec(2, 4)

burnin=9000
# get MAP estimate
tmp = pm.df_summary(trace1[burnin:], varnames=['alpha', 'beta'])
alphaMAP = tmp['mean'][np.arange(nsubjs)]
betaMAP = tmp['mean'][np.arange(nsubjs)+nsubjs]

def find_nearest(array, value):
    idx = (np.abs(array-value)).argmin()
    return idx

for ip in np.arange(nsubjs):
    ax = plt.subplot(gs[ip])
    xp = np.array(x.iloc[ip, :])
    yp = np.array(rprop.iloc[ip, :])
    ax.scatter(xp, yp, marker='s', alpha=.5)
    
    xl = np.linspace(190, 410, 100)
    yl = logit(alphaMAP[ip] + betaMAP[ip]*(xl-xmean[ip]))
    x1 = xl[find_nearest(yl, .5)]
    x2 = xl[find_nearest(yl, .84)]

    plt.plot(xl, yl, 'k', linewidth=2)
    plt.plot([x1, x1],[-.1, .5], color='k', linestyle='--', linewidth=1)
    plt.plot([190, x1],[.5, .5], color='k', linestyle='--', linewidth=1)
    plt.plot([x2, x2],[-.1, .84], color='k', linestyle='--', linewidth=1)
    plt.plot([190, x2],[.84, .84], color='k', linestyle='--', linewidth=1)

    plt.axis([190, 410, -.1, 1.1])
    plt.yticks((0, .5, .84, 1))
    plt.title('Subject %s'%(ip+1))
    
plt.tight_layout()
plt.show()
In [10]:
# Posterior sample
from collections import defaultdict
alphadist = model1.alpha
betadist = model1.beta

ppcsamples = 500
ppcsize = 100
ppc = defaultdict(list)
for idx in np.random.randint(burnin, 1e4, ppcsamples):
    param = trace1[idx]
    ppc['alpha'].append(alphadist.distribution.random(point=param, size=ppcsize))
    ppc['beta'].append(betadist.distribution.random(point=param, size=ppcsize))
    
# np.asarray(ppc['alpha']).shape
alphaPPC = np.asarray(ppc['alpha']).mean(axis=1)
betaPPC = np.asarray(ppc['beta']).mean(axis=1)
In [11]:
# PLOT FOR EXERCISE 12.1.2 
fig = plt.figure(figsize=(16, 8))
fig.text(0.5, -0.02, 'Test Interval (ms)', ha='center', fontsize=20)
fig.text(-0.02, 0.5, 'Proportion of Long Responses', va='center', rotation='vertical', fontsize=20)
gs = gridspec.GridSpec(2, 4)

ppcsamples=100

for ip in np.arange(nsubjs):
    ax = plt.subplot(gs[ip])
    xp = np.array(x.iloc[ip, :])
    yp = np.array(rprop.iloc[ip, :])
    ax.scatter(xp, yp, marker='s', alpha=.5)
    
    xl = np.linspace(190, 410, 100)
    yl = logit(alphaMAP[ip] + betaMAP[ip]*(xl-xmean[ip]))

    # Posterior sample from the trace
    for ips in np.random.randint(burnin, 1e4, ppcsamples):
        param = trace1[ips]
        yl2 = logit(param['alpha'][ip] + param['beta'][ip]*(xl-xmean[ip]))
        plt.plot(xl, yl2, 'k', linewidth=2, alpha=.05)
    
    plt.plot(xl, yl, 'r', linewidth=2)
    
    plt.axis([190, 410, -.1, 1.1])
    plt.yticks((0, .5, .84, 1))
    plt.title('Subject %s'%(ip+1))
    
    
plt.tight_layout()
plt.show()