我正试图使用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)），但结果相似，收敛仍然很慢。
看起来您试图建模的是过度分散的数据。事实上，负二项分布只是一个泊松，其平均值根据伽马分布（指数的一般形式）分布。因此，绕过定义1000个变量的一种方法是直接使用负二项式。不过，请记住，尽管名义上有1000个变量，但变量的有效数量在1到1000之间，这取决于均值分布的约束程度。你在这里基本上是在定义随机效应。