1

我正在尝试在我的研究中遇到的各种数据上练习使用 pymc3，但是当每个人给我多个数据点并且每个人都来自不同的组时，我无法思考如何拟合模型（所以尝试分层模型）。

这是我正在使用的练习场景：假设我们有 2 组人，每组 N = 30。所有 60 人都经过 10 个问题的调查，每个人可以回答（“1”）或不回答（“0”）每个问题。因此，对于每个人，我都有一个长度为 10 的数组，其中包含 1 和 0。

为了对这些数据建模，我假设每个人都有一些潜在的特征“theta”，每个项目都有一个“歧视”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.Data 做一些事情来解决这个问题，但是文档对我来说并不完全清楚。

基本上，我习惯于在 JAGS 中构建模型，在其中循环遍历每个数据点，但 pymc3 似乎并没有那样工作。我对如何在模型中构建/索引我的随机变量感到困惑，以确保概率改变我希望它们从试验到试验的方式，并确保我估计的参数对应于正确的组中的正确的人。

提前感谢您的帮助。我对 pymc3 很陌生，并试图掌握它，并想尝试与 JAGS 不同的东西。

编辑：我能够通过首先通过循环试验构建我需要的数组来解决这个问题，然后使用以下方法转换数组：

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

然后我把这个新变量放在伯努利实例的“p”参数中，它起作用了！这是更新后的完整模型：（下面，我将 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
4

0 回答 0