pym在Python中的贝叶斯IRT模型

Question

我想评估一下Python中的项目响应理论（IRT）模型。更具体地说，以学生参加考试的典型IRT为例。对于每个学生，我们观察他们在考试中回答的问题是否给出了正确的答案。这给了我们一个观察结果矩阵X，从中我们可以估计出（1）一个难度参数α和（2）一个判别参数β，这样我们就可以估计每个学生的潜在能力Y，作为他们在每个测试问题上是否得到正确答案的函数，即α+βX。我能做的最好的例子如何使用Python中的MCMC来估计这种IRT贝叶斯模型是example。从这个例子中我不明白的是，一个学生是否得到正确答案的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

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

^{pr2}$

可以追溯到这条线：

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

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

Answer 1

你被Python的一个狡猾的部分咬了一口。pymc的全局导入将您的numpyexp替换为另一个{}。为了得到你想要的exp，你可以在你的sigmoid中使用np.exp。（我想知道np.是从哪里来的？）在

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

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

^{pr2}$