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

Questions
1

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...
Initializing NUTS using jitter+adapt_diag...
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'
2

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.

3

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
4

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