Python中基于pymc的贝叶斯IRT模型

Question

我想用Python来估计一个IRT模型。更具体地说，以学生参加考试的典型IRT为例。对于每个学生，我们观察他们是否对考试中回答的问题给出了正确的答案。这给出了一个观测结果矩阵X，对于每一个问题，我们想要估计(1)一个困难参数α和(2)一个判别参数β，这样我们也可以作为每个测试问题得到正确答案的函数来估计每个学生的潜在能力Y，即α+βX。我能找到的最好的例子是如何用α来估计这类IRT贝叶斯模型。我从这个例子中不明白的是，一个学生是否得到正确答案的X矩阵进入了模型。下面是这个代码的一个稍微修改的版本，旨在估计每个学生的潜在能力：

#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

当我运行脚本时，我会收到错误消息：

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

可以追溯到这条线上：

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

我检查了原始脚本(它在从潜在值生成结果和从结果计算潜在值之间切换)和correctness变量(我认为它是上面描述的测试结果的X矩阵)充满了False值。当我将correctness设置为满是False值时，脚本就完成了。然而，这似乎意味着每个学生都把每个问题都搞错了，这是没有道理的。我认为这可能是这个问题的正确答案，所以我将correctness中的所有值设置为True，但这会产生相同的错误。我做错了什么?如何从X矩阵中估计潜在能力，用pymc判断学生是否在试题上得到正确答案？

Answer 1

你被Python的一个狡猾的部分咬了一口。pymc的全球导入将您的numpy exp替换为不同的exp。要获得所需的exp，可以在sigmoid确定性中使用np.exp。(我想知道np.是从哪里来的？)

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

看起来你还有一些调试要做，但我希望这能让你摆脱困境。这是我为什么喜欢这种模式的一个很好的例子：

import numpy as np, pymc as pm