2
我试图从a paper from Richard McElreath复制数据分析，其中他将数据与分层零膨胀Gamma模型拟合。这些数据是关于狩猎返回约15000狩猎之旅从约150二十多年的猎人。由于许多狩猎旅行的回报率为零，因此该模型假定每次旅行具有零回报率的概率为pi，而1 - pi为正数回报率的概率为alpha和beta之间的Gamma分布。PYMC3：NUTS难以从分层零充气伽马模型采样
预测变量是年龄，该模型使用年龄多项式（最高为3）来模拟pi和alpha。而且由于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
n_trip = len(d)
n_hunter = len(d['hunter.id'].unique())
idx_hunter = d['hunter.id'].values
y = d['kg.meat'].values
age = d['age.s'].values
age2 = (d['age.s'].values)**2
age3 = (d['age.s'].values)**3
日志概率密度函数被obtainted为零膨胀的伽玛。
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矩阵的相关矩阵，所述LKJ事先在pymc3被给定为一个维向量
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采样器，我的电脑已经耗尽内存和链还没有收敛。