Байесовская модель IRT в Python с использованием pymc

Question

Я хотел бы оценить модель модели ответа на вопрос (IRT) в Python. В частности, возьмите канонический пример IRT для студентов, сдающих экзамен. Для каждого ученика мы наблюдаем, правильно ли ответили на вопросы, которые они ответили на экзамене. Это дает нам наблюдаемую матрицу результатов X, из которой по каждому вопросу мы хотели бы оценить (1) параметр трудности α и (2) параметр дискриминации β, чтобы мы могли также оценить каждую скрытую способность Y студентов как функцию того, получил правильный ответ или нет по каждому тестовому вопросу, т.е. α + βX. Лучший пример, который я мог найти, как оценить этот тип байесовской модели IRT с использованием MCMC в Python, - это 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?

Answer 1

Вы были немного скрыты частью Python. Глобальный импорт pymc заменил ваш numpyexp с другим exp. Чтобы получить нужный вам exp, вы можете использовать np.exp в своем sigmoid детерминистическом. (Откуда np. пришли, интересно?)

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

Похоже, у вас все еще есть некоторые отладки, чтобы сделать, но я надеюсь, что это получает вас отключился. Это хороший пример того, почему я предпочитаю шаблон:

import numpy as np, pymc as pm