Bayesian IRT Model in Python using pymc

Question

I would like to estimate an Item Response Theory (IRT) model in Python. More specifically, take the canonical IRT example of students taking an exam. For each student we observe whether or not they gave the correct answer to the questions they answered on an exam. This gives us an observed results matrix X from which for each question we would like to estimate (1) a difficulty parameter α and (2) a discrimination parameter β such that we can also estimate each students latent ability Y as a function of whether they got a correct answer or not on each test question, i.e. α + βX. The best example I could find of how to estimate this type of IRT Bayesian model using MCMC in Python was this example. What I do not understand from this example is where the X matrix of whether or not a student got a correct answer on a test question enters into the model. Here is a slightly modified version of this code intended to estimate the latent ability of each student:

#from pylab import * #Pylab will not install with pip so I just loaded numpy itself
from numpy import *
import numpy
from pymc import *
from pymc.Matplot import plot as mplot

numquestions = 300 # number of test items being simulated
numpeople = 10 # number of participants
numthetas = 1 # number of latent proficiency variables

generating = 0
theta_initial = zeros((numthetas, numpeople))
correctness = np.random.randint(2, size= numquestions * numpeople) == 1 #Produces Error
#correctness = np.random.randint(2, size= numquestions * numpeople) == -1 #all False code runs fine
#correctness = np.random.randint(2, size= numquestions * numpeople) != -1 #all True code throws error message

correctness.shape = (numquestions, numpeople)


# theta (proficiency params) are sampled from a normal distribution
theta = Normal("theta", mu=0, tau=1, value=theta_initial, observed= generating)


# question-parameters (IRT params) are sampled from normal distributions (though others were tried)
a = Normal("a", mu=1, tau=1, value=[[0.0] * numthetas] * numquestions)
# a = Exponential("a", beta=0.01, value=[[0.0] * numthetas] * numquestions)
b = Normal("b", mu=0, tau=1, value=[0.0] * numquestions)

# take vectors theta/a/b, return a vector of probabilities of each person getting each question correct
@deterministic
def sigmoid(theta=theta, a=a, b=b): 
    bs = repeat(reshape(b, (len(b), 1)), numpeople, 1)
    return np.zeros_like(1.0 / (1.0 + exp(bs - dot(a, theta)))) #np.zeros_like fixes error

# take the probabilities coming out of the sigmoid, and flip weighted coins
correct = Bernoulli('correct', p=sigmoid, value=correctness, observed=not generating)

# create a pymc simulation object, including all the above variables
m = MCMC([a,b,theta,sigmoid,correct])

# run an interactive MCMC sampling session
m.isample(iter=20000, burn=15000)


mydict = m.stats()
print(mydict['theta']['mean']) #Get ability parameters for each student

When I run the script I get the error message:

pymc.Node.ZeroProbability: Stochastic correct's value is outside its support,
 or it forbids its parents' current values.`

which traces back to the line:

correct = Bernoulli('correct', p=sigmoid, value=correctness, observed=not generating)

I checked the original script (which toggles between generating results from latent values and calculating latent values from results) and the correctness variable, which I thought of as the X matrix of test results described above, is full of False values. When I set correctness to be full of False values then the script completes. However, this would seem to say every student got every question wrong, which would not make a whole lot of sense. I thought it might be the correct answer for the question, so I set all values in correctness to True but that produced the same error. What am I doing wrong and how can I estimate latent ability using an IRT model from an X matrix of whether or not a student got the correct answer on a test question using pymc?

Answer 1

You have been bit by a sneaky part of Python. The global import of pymc replaced your numpy exp with a different exp. To get the exp you want, you can use np.exp in your sigmoid deterministic. (Where did np. come from, I wonder?)

return np.exp(1.0 / (1.0 + np.exp(bs - dot(a, theta))))

It looks like you still have some debugging to do, but I hope this gets you unstuck. It is a good example of why I favor the pattern:

import numpy as np, pymc as pm