我想用 Python 估算项目 react 理论 (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 但产生了相同的错误。我做错了什么？如何使用 IRT 模型从 X 矩阵估计学生是否使用 pymc 在测试问题上得到正确答案的潜在能力？
请您参考如下方法:
您被 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
