Блог

Я пытаюсь попрактиковаться в использовании pymc3 для тех типов данных, с которыми я сталкиваюсь в своих исследованиях, но мне трудно продумать, как подогнать модель, когда каждый человек дает мне несколько точек данных, и каждый человек происходит из другой группы (поэтому я пытаюсь использовать иерархическую модель).

Вот практический сценарий, который я использую: предположим, у нас есть 2 группы людей, N = 30 в каждой группе. Все 60 человек проходят опрос из 10 вопросов, где каждый человек может ответить («1») или не отвечать («0») на каждый вопрос. Итак, для каждого пользователя у меня есть массив длиной 10 с 1 и 0.

Для моделирования этих данных я предполагаю, что у каждого человека есть какая-то скрытая черта «тета», и каждый элемент имеет «различение» a и «сложность» b (это всего лишь базовая модель ответа на элемент), а вероятность ответа («1») задается: (1 exp(-a(theta — b)))^(-1). (Логистика применяется к a (theta — b) .)

Вот как я попытался подогнать его, используя pymc3:

 traces = {}

for grp in range(2):

group = prac_data["Group ID"] == grp
data = prac_data[group]["Response"]

with pm.Model() as irt:
    # Priors
    a_tmp = pm.Normal('a_tmp',mu=0, sd = 1, shape = 10)
    a = pm.Deterministic('a', np.exp(a_tmp))
    # We do this transformation since we must have a >= 0 
    
    b = pm.Normal('b', mu = 0, sd = 1, shape = 10)
    
    # Now for the hyperpriors on the groups:
    theta_mu = pm.Normal('theta_mu', mu = 0, sd = 1)
    theta_sigma = pm.Uniform('theta_sigma', upper = 2, lower = 0)
    
   
    
    theta = pm.Normal('theta', mu = theta_mu,
                             sd = theta_sigma, shape = N)
    
    p = getProbs(Disc, Diff, theta, N)
    
    y = pm.Bernoulli('y', p = p, observed = data)
        
    traces[grp] = pm.sample(1000)
  

Предполагается, что функция «getProbs» дает мне массив вероятностей для случайной величины Бернулли, поскольку вероятность ответа 1 изменяется в зависимости от испытаний / вопросов опроса для каждого человека. Но этот метод выдает мне ошибку, потому что в нем говорится «указать один из p или logit_p», но я думал, что я сделал с функцией?

Вот код для «getProbs» на случай, если это полезно:

 def getProbs(Disc, Diff, THETA, Nprt):
# Get a large array of probabilities for the bernoulli random variable
n = len(Disc)
m = Nprt
probs = np.array([])
for th in range(m):
    for t in range(n):
        p = item(Disc[t], Diff[t], THETA[th])
        probs = np.append(probs, p)
return probs
  

Я добавил параметр Nprt, потому что, если бы я попытался получить длину THETA, это выдало бы мне ошибку, поскольку это объект FreeRV. Я знаю, что могу попытаться векторизовать функцию «item», которая является просто логистической функцией, которую я указал выше, вместо того, чтобы делать это таким образом, но это также привело к ошибке, когда я попытался ее запустить.

Я думаю, что могу что-то сделать с pm.Данные, чтобы исправить это, но документация мне не совсем понятна.

В принципе, я привык строить модели в JAGS, где вы перебираете каждую точку данных, но pymc3, похоже, так не работает. Я не понимаю, как создавать / индексировать мои случайные переменные в модели, чтобы убедиться, что вероятности изменяются так, как я бы хотел, от испытания к испытанию, и чтобы убедиться, что параметры, которые я оцениваю, соответствуют нужному человеку в нужной группе.

Заранее спасибо за любую помощь. Я довольно новичок в pymc3 и пытаюсь освоиться с ним, и хотел попробовать что-то отличное от JAGS.

РЕДАКТИРОВАТЬ: я смог решить эту проблему, сначала создав необходимый мне массив, выполнив цикл испытаний, а затем преобразовав массив с помощью:

 p = theano.tensor.stack(p, axis = 0)
  

Затем я поместил эту новую переменную в аргумент «p» экземпляра Bernoulli, и это сработало! Вот обновленная полная модель: (ниже я импортировал theano.tensor как T)

 group = group.astype('int') 
data = prac_data["Response"]


with pm.Model() as irt:
# Priors

# Item parameters:
a = pm.Gamma('a', alpha = 1, beta = 1, shape = 10) # Discrimination

b = pm.Normal('b', mu = 0, sd = 1, shape = 10) # Difficulty

# Now for the hyperpriors on the groups: shape = 2 as there are 2 groups
theta_mu = pm.Normal('theta_mu', mu = 0, sd = 1, shape = 2)
theta_sigma = pm.Uniform('theta_sigma', upper = 2, lower = 0, shape = 2)


# Individual-level person parameters:
# group is a 2*N array that lets the model know which
# theta_mu to use for each theta to estimate
theta = pm.Normal('theta', mu = theta_mu[group],
                         sd = theta_sigma[group], shape = 2*N)


# Here, we're building an array of the probabilities we need for 
# each trial:
p = np.array([])
for n in range(2*N):
    for t in range(10):
        x = -a[t]*(theta[n] - b[t])
        p = np.append(p, x)
  
# Here, we turn p into a tensor object to put as an argument to the
# Bernoulli random variable
p = T.stack(p, axis = 0)

y = pm.Bernoulli('y', logit_p = p, observed = data)


# On my computer, this took about 5 minutes to run. 
traces = pm.sample(1000, cores = 1)
    
print(az.summary(traces)) # Summary of parameter distributions