0

我正在尝试使用 PyMC 2.3 来获得对复合模型参数的估计。“复合”是指遵循分布的随机变量，其参数是另一个随机变量。（“嵌套”或“分层”有时用于指代这种情况，但我认为它们不太具体，并且在这种情况下会产生更多混乱）。

让我们举一个具体的例子。“真实”数据是从复合分布生成的，该分布是泊松，其参数分布为指数。在普通 scipy 中，数据生成如下：

import numpy as np
from scipy.stats import distributions
np.random.seed(3)  # for repeatability

nsamples = 1000
tau_true = 50
orig_lambda_sample = distributions.expon(scale=tau_true).rvs(nsamples)
data = distributions.poisson(orig_lambda_sample).rvs(nsamples)

我想获得模型参数的估计值tau_true。到目前为止，我在 PyMC 中对这个问题进行建模的方法如下：

tau = pm.Uniform('tau', 0, 100)
count_rates = pm.Exponential('count_rates', beta=1/tau, size=nsamples)
counts = pm.Poisson('counts', mu=count_rates, value=data, observed=True)

请注意，我size=nsamples曾经为每个样本设置一个新的随机变量。

最后，我将 MCMC 运行为：

model = pm.Model([count_rates, counts, tau])
mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000)

该模型收敛（虽然缓慢，> 10^5 次迭代）到以 50 ( tau_true) 为中心的分布。然而，定义 1000 个随机变量来拟合具有单个参数的单个分布似乎有点过头了。

有没有更好的方法来描述 PyMC 中的复合模型？

PS我也尝试了更丰富的先验（tau = pm.Normal('tau', mu=51, tau=1/2**2)），但结果相似，收敛仍然很慢。

4

1 回答 1

1

看起来您要建模的是过度分散的数据。事实上，负二项分布只是一个泊松，其均值根据伽马分布（指数的一般形式）分布。因此，绕过定义 1000 个变量的一种方法是直接使用负二项式。请记住，尽管名义上有 1000 个变量，但变量的有效数量在 1 到 1000 之间，具体取决于均值分布的约束程度。您本质上是在此处定义随机效应。

于 2014-10-14T19:44:14.880 回答