# Implementing normal CDF for signal detection theory model (includes full sample code)

I’m trying to adapt some WinBUGS code from the Bayesian Cognitive Modeling book:

model {
for (i in 1:k) {
# Observed counts
h[i] ~ dbin(thetah[i],s[i])
f[i] ~ dbin(thetaf[i],n[i])

# Reparameterization using equal-variance Gaussian SDT
thetah[i] <- phi(d[i]/2-c[i])
thetaf[i] <- phi(-d[i]/2-c[i])

# Priors
d[i] ~ dnorm(0, 0.5)
c[i] ~ dnorm(0, 2)
}
}


Since there does not appear to be a phi or norm_cdf function in PyMC3 or Theano, I implemented it using the error function. Here is my sample code:

import numpy as np
import pymc3 as mc
from theano.tensor import erf

# Make simulated data
obs_hit = np.ones(100)
obs_fa = np.zeros(100)

# Make it imperfect
obs_hit[0] = 0
obs_fa[0] = 1

with mc.Model() as one_param:
d = mc.Normal('d', 0, 0.5)
c = mc.Normal('c', 0, 2)
theta_hit = mc.Deterministic('hit_rate', 0.5*(1+erf((d/2-c)/np.sqrt(2))))
theta_fa = mc.Deterministic('fa_rate', 0.5*(1+erf((-d/2-c)/np.sqrt(2))))

hits = mc.Bernoulli('hit', theta_hit, observed=obs_hit)
fa = mc.Bernoulli('fa', theta_fa, observed=obs_fa)
trace = mc.sample()


While this seems to work, I have two issues:

1. I get some Theano errors. See below for output. How can I make this go away?
2. If obs_hit is all ones (i.e., the subject gets every trial correct), then the model fails to converge. How do I fix this?

Theano error output:

Auto-assigning NUTS sampler...
ERROR (theano.gof.opt): Optimization failure due to: local_grad_log_erfc_neg
ERROR (theano.gof.opt): node: Elemwise{true_div,no_inplace}(Elemwise{mul,no_inplace}.0, Elemwise{erfc,no_inplace}.0)
ERROR (theano.gof.opt): TRACEBACK:
ERROR (theano.gof.opt): Traceback (most recent call last):
File "/home/bburan/bin/miniconda3/envs/python3/lib/python3.7/site-packages/theano/gof/opt.py", line 2034, in process_node
replacements = lopt.transform(node)
File "/home/bburan/bin/miniconda3/envs/python3/lib/python3.7/site-packages/theano/tensor/opt.py", line 6789, in local_grad_log_erfc_neg
if not exp.owner.inputs[0].owner:
AttributeError: 'NoneType' object has no attribute 'owner'


# Update

I tried switching from a Bernoulli likelihood to Binomial likelihood:

with mc.Model() as one_param:
d = mc.Normal('d', 0, 0.5)
c = mc.Normal('c', 0, 2)
theta_hit = mc.Deterministic('hit_rate', 0.5*(1+erf((d/2-c)/np.sqrt(2))))
theta_fa = mc.Deterministic('fa_rate', 0.5*(1+erf((-d/2-c)/np.sqrt(2))))

#hits = mc.Bernoulli('hit', theta_hit, observed=obs_hit)
#fa = mc.Bernoulli('fa', theta_fa, observed=obs_fa)
hits = mc.Binomial('hit', p=theta_hit, n=len(obs_hit), observed=obs_hit.sum())
fa = mc.Binomial('fa', p=theta_fa, n=len(obs_fa), observed=obs_fa.sum())
trace = mc.sample()


I still get the Theano errors, but the model now converges even if obs_hit is all ones.

I have a couple of blog posts about Bayesian SDT that haven’t made it to my new site yet. This post encourages me to re-upload them.

In the meantime, here are some pointers:

1. You don’t need to alter your data to prevent “extreme” performance (0 or 100% hit or false-alarm rates under a Bayesian model.

2. Prior standard deviations for d^\prime and c should be switched to imply uniform priors on hit and false-alarm rates.

3. There is indeed a norm cdf function in PyMC3: it’s called invprobit.

4. Your “binomial” implementation will be considerably more efficient than your “Bernoulli” implementation because the number of observed variables is reduced to 2 counts rather than n trials.

1 Like

Any chance you’ve had a chance to upload them? I’d like to read them. I tracked down your blog (from your profile) but didn’t see anything specifically pertaining to Bayesian SDT (by the way, I noticed you worked with Barb – I know a lot of her students since I’m in auditory research).