Python 中的贝叶斯 IRT 模型使用 pymc

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

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

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 矩阵使用 IRT 模型估计学生的潜在能力，以确定学生是否使用 pymc 在测试问题上得到正确答案？

Answer 1

你被 Python 的一个偷偷摸摸的部分咬了。 pymc 的全局导入用不同的 exp 替换了您的 numpy 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

Python 中的贝叶斯 IRT 模型使用 pymc

Bayesian IRT Model in Python using pymc

python

bayesian

pymc