PYMC3:NUTS很难从一个分层的零膨胀gamma mod采样

Question

我试图复制a paper from Richard McElreath的数据分析，在这个分析中，他用一个分层的零膨胀Gamma模型来拟合数据。这些数据是关于大约15000个猎人在20年中的狩猎收益。由于许多狩猎旅行都是零回报，因此该模型假设每一次狩猎都有pi的零回报概率，以及{}的正收益概率，其服从参数为alpha和{}。在

预测变量是年龄，模型使用年龄多项式（最多3阶）来建模pi和{}。由于15000旅行属于150个体猎人，因此每个猎人都有自己的系数，所有的系数都遵循一个共同的多元正态分布。有关模型的详细信息，请参阅以下代码。模型规格看起来不错，但NUTS在开始取样时遇到了麻烦：大约20分钟后，它只给出了大约10个样品，取样器就停在那里，告诉我要花几百个小时才能完成取样。我想知道是什么引起了这些问题。在

通常的进口货

import pymc3 as pm
import numpy as np
from pymc3.distributions import Continuous, Gamma
import theano.tensor as tt

数据可以从github获得

^{pr2}$

零膨胀伽马的对数概率密度函数。在

class ZeroInflatedGamma(Continuous):
    def __init__(self, alpha, beta, pi, *args, **kwargs):
        super(ZeroInflatedGamma, self).__init__(*args, **kwargs)
        self.alpha = alpha
        self.beta = beta
        self.pi = pi = tt.as_tensor_variable(pi)
        self.gamma = Gamma.dist(alpha, beta)

    def logp(self, value):
        return tt.switch(value > 0,
                         tt.log(1 - self.pi) + self.gamma.logp(value),
                         tt.log(self.pi))

这是一个索引9X9矩阵之前的相关矩阵，pymc3中的LKJ先验是一维向量

dim = 9
n_elem = dim * (dim - 1) / 2
tri_index = np.zeros([dim, dim], dtype=int)
tri_index[np.triu_indices(dim, k=1)] = np.arange(n_elem)
tri_index[np.triu_indices(dim, k=1)[::-1]] = np.arange(n_elem)

这是模型

with pm.Model() as Vary9_model:

    # hyper-priors
    mu_a = pm.Normal('mu_a', mu=0, sd=100, shape=9)
    sigma_a = pm.HalfCauchy('sigma_a', 5, shape=9)

    # build the covariance matrix
    C_triu = pm.LKJCorr('C_triu', n=2, p=9)    
    C = tt.fill_diagonal(C_triu[tri_index], 1)
    sigma_diag = tt.nlinalg.diag(sigma_a)
    cov = tt.nlinalg.matrix_dot(sigma_diag, C, sigma_diag)

    # priors for each hunter and all the linear components, 9 dimensional Gaussian  
    a = pm.MvNormal('a', mu=mu_a, cov=cov, shape=(n_hunter, 9))

    # linear function  
    mupi = a[:,0][idx_hunter] + a[:,1][idx_hunter] * age + a[:,2][idx_hunter] * age2 + a[:,3][idx_hunter] * age3
    mualpha = a[:,4][idx_hunter] + a[:,5][idx_hunter] * age + a[:,6][idx_hunter] * age2 + a[:,7][idx_hunter] * age3

    pi = pm.Deterministic('pi', pm.math.sigmoid(mupi))
    alpha = pm.Deterministic('alpha', pm.math.exp(mualpha))
    beta = pm.Deterministic('beta', pm.math.exp(a[:,8][idx_hunter]))

    y_obs = ZeroInflatedGamma('y_obs', alpha, beta, pi, observed=y)

    Vary9_trace = pm.sample(6000, njobs=2)

这是模型的状态：

Auto-assigning NUTS sampler...
Initializing NUTS using advi...
Average ELBO = -28,366: 100%|██████████| 200000/200000 [15:36<00:00, 213.57it/s]
Finished [100%]: Average ELBO = -28,365
  0%|          | 22/6000 [15:51<63:49:25, 38.44s/it]

我对这个问题有一些想法，但不确定是什么原因。在

• 九维高斯太难采样了吗？我以前只将mualpha和mupi的截获建模为二元高斯模型，虽然速度慢，但有效（模型拟合大约需要20分钟）

• 是概率密度造成了问题吗？我自己写的密度函数，不确定它是否运行良好。我认为密度函数在零处是不可微的，这会给坚果采样器带来麻烦吗？

• 是因为预测变量高度相关吗？在这个模型中，线性模型的组成部分是年龄的三次多项式，自然地，预测因子是高度相关的。

或者是因为别的原因？在

作为补充说明，我尝试使用Metropolis取样器，我的计算机内存不足，链仍然没有聚合。在

Answer 1

零膨胀伽马看起来很好。密度函数对于π，α和β是可微的。这就是观察到的变量所需要的全部。如果你试图估计这些值，你只需要关于这个值的导数。在

在实施LKJCorr时存在一个问题： https://github.com/pymc-devs/pymc3/pull/1863 你可以在主人身上再试一次。遗憾的是，pymc3不支持在cholesky分解的参数化中使用MVNormal和LKJCorr。这也可能有帮助。github上有一个正在处理的请求，请求对此进行： https://github.com/pymc-devs/pymc3/pull/1875

为了提高收敛性，可以尝试对a进行非中心参数化。一些类似于

a_raw = pm.Normal('a_raw', shape=(9, n_hunter))
a = mu_a[None, :] + tt.dot(tt.slinalg.cholesky(cov), a_raw)

当然，如果我们有那个叫乔尔斯基·伊克科尔的人。。。在