机器算法验证 - 在 pymc 中实现有序概率模型 - 吾爱随笔录

在 pymc 中实现有序概率模型

机器算法验证马尔可夫链蒙特卡罗 Python 序数数据罗吉特 pymc

2022-03-29 13:00:23

我正在尝试在 pymc 中实现有序概率模型，但我被卡住了。该模型类似于Welinder 的“群体的多维智慧”，有编码器（由 i 索引）和文档（由 j 索引）。编码人员将代码分配给文档，但编码过程很嘈杂。

我们希望估计两件事。第一的， $z_j$ ，每个文档在某个潜在维度上的真实、潜在价值。第二， $\beta_i$ ，偏差项揭示了每个编码员的评估与组平均值的准确程度。

如果代码是连续的，这将非常容易，但我拥有的数据是有序的。代码在 1,2,...,5 范围内。

形式上，我们可以这样表示：

$x_{ij} \sim Normal( \alpha_i + \beta_iz_j, 1)$

$code_{ij} = cut(x_{ij}, w)$

在哪里 $cut$ 是一个截止函数，它根据最高截止点的索引分配值， $w$ , 超过 $x_{ij}$ .

到目前为止，一切都很好。该模型的问题在于，切点函数是确定性的，并且可以观察到代码。但是在 pymc（和其他 MCMC 程序，例如 JAGS）中，也无法观察到确定性节点。所以这个模型不能直接在pymc中构建。

好像有办法治疗 $x$ 作为确定性的，并且 $codes$ 作为一个随机函数 $x$ . 这可能会使 $code$ 涉及分类节点。但我不确定如何指定概率函数，我有点担心我的整个方法可能会失败。谁能让我直截了当？

另外——长镜头——如果有人在 pymc 中编码，很高兴看到它的源代码。Ordinal probit/logit 是一个非常标准的模型，但我在网上的任何地方都找不到 pymc 示例。

1个回答

这需要相当多的工作，但我最终得到了它。请注意，我使用了来自 github 的开发版本（pymc 2.2grad），而不是 pypi 上可用的旧版本。

此外，这运行得非常缓慢，因为它没有很好地利用 numpy 的数组操作。对于合理大小的数据集，一些智能预处理和对 cutpoints 节点函数的重写可能会解决这个问题。

这是完整的源代码：

from scipy.stats import norm, pearsonr
from pymc import Normal, Lambda, Uniform, Exponential, stochastic, deterministic, observed, MCMC, Matplot
from numpy import mean, std, log
import numpy as np

#Set array dimensions
(I, J, K, M, N) = (5, 20, 3, 4, 1000)

#Set simulation parameters
alpha_star = np.random.normal(0, 1, size=(I,K))
beta_star = np.random.normal(1, 1, size=(I,K))
z_star = np.random.normal(0, 1, size=(J,1))
w_star = np.array([0,1,3])

#Generate data
coder = np.random.randint(I, size=(N))
doc = np.random.randint(J, size=(N))
item = np.random.randint(K, size=(N))

code = np.zeros(shape=(N))
for n in range(N):
    i, j, k = coder[n], doc[n], item[n]
    m = alpha_star[i,k] + beta_star[i,k] * z_star[j] + np.random.normal(0, 1)
    code[n] = 1+sum(m > w_star)
#    print "\t".join([str(x) for x in [i, j, k, m, code[n]]])

#Set GLM parameters
alpha = Normal('alpha', mu=0.0, tau=0.01, value=np.zeros(I*K))
beta = Normal('beta', mu=1.0, tau=0.01, value=np.ones(I*K))

z = Normal('z', mu=0.0, tau=0.01, value=np.random.normal(0,1,J))

w = Exponential('w', .1, value=np.ones(M-3))

#Link functions
mu = Lambda('mu', lambda alpha=alpha, beta=beta, z=z, i=coder, j=doc, k=item: alpha[i+I*k]+beta[i+I*k]*z[j])

@deterministic(plot=False)
def cutpoints(w=w):
    w2 = [-np.inf, 0.0, 1.0]
    v = 1
    for i in w:
        v += i
        w2.append(v)
    w2.append(np.inf)

    cp = np.array( w2 )
    return cp


@stochastic(dtype=int, observed=True)
def y(value=code, mu=mu, cp=cutpoints):

    def logp(value, mu, cp):
        d = norm.cdf(cp[value]-mu)-norm.cdf(cp[value-1]-mu)
        lp = sum(log(d))
        return lp

#Run chain
M = MCMC([alpha, beta, z, mu, w, cutpoints, y])
M.isample(10000, 5000, thin=5, verbose=0)

#Summarize results
Matplot.summary_plot([alpha], name="alpha", path="./graphs/")
Matplot.summary_plot([beta], name="beta", path="./graphs/")
Matplot.summary_plot([z], name="z", path="./graphs/")
Matplot.summary_plot([w], name="w", path="./graphs/")

print pearsonr( alpha_star.transpose().reshape((I*K,)), alpha.stats()['mean'])
print pearsonr( beta_star.transpose().reshape((I*K,)), beta.stats()['mean'])
print pearsonr( z_star.reshape((J,)), z.stats()['mean'])

其它你可能感兴趣的问题

上一篇R：如何“控制”线性混合效应回归模型中的另一个变量？下一篇传达回归模型结果