多观测值的pymc3分层模型，在MCMC期间不计算可能性？

Question

我试图用pymc3来拟合一个层次模型，结果遇到了麻烦，我的前驱被抽样了，但我的可能性没有。我尝试过复制描述的here方法，它适用于多个观察。在

我使用的是R中未标记包中的mallard数据集，作为pandas数据帧加载。你可以从我的GitHub页面here获取我正在使用的CSV。简言之，每一行代表沿一个横断面收集的最多三个观测值（y1、y2和y3列）。我尝试过删除y列或协变量（海拔、森林覆盖率、调查长度）中缺少数据（NaN）的行，这似乎没有改变我遇到的问题。在

mallard = pd.read_csv('mallard.csv')
elev, length, forest = mallard.elev, mallard.length, mallard.forest

# transpose ys to have same shape as series of each covariate (elev, length...)
observations = mallard[['y1','y2','y3']].T
num_transects = len(mallard)

with pm.Model() as mallard_model:
    # priors
    # detection probability considered constant across all transects and plots
    prob = pm.Beta('prob', alpha=1, beta=1)
    # priors for deterministic linear model of local abundance at a transect 
    # that applies to all plots within a transect
    b0 = pm.Normal('b0', mu=0, sd=10) # intercept
    b1 = pm.Normal('b1', mu=0, sd=10) # elevation
    b2 = pm.Normal('b2', mu=0, sd=10) # forest cover
    b3 = pm.Normal('b3', mu=0, sd=10) # survey length

    # linear model of local abundance using log link
    lam = pm.Deterministic('lam', pm.math.exp(b0 + b1*elev + b2*forest + b3*length))

    # likelihood of observations
    # Ni is abundance at a transect, with binomial distribution across 
    # plots within a transect

    Ni = pm.Poisson('Ni', mu=lam, shape=num_transects)
    Y_obs = pm.Binomial('Y_obs', n=Ni, p=prob, observed=observations)

    # inference, use default step functions for each parameter
    trace = pm.sample(draws=5000, init='ADVI', n_init=10000)#, step=pm.Metropolis())

    plt.figure(figsize=(7, 7))
    pm.traceplot(trace[100:]) # leave out first 100 draws as burn-in
    plt.tight_layout()

结果是：
将螺母分配给prob斨u logodds 分配给b0的螺母
分配给b1的螺母
分配给b2的螺母
分配给b3的螺母
分配给Ni的大都市
分配给你的大都会失踪
100%| |∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮∮5000/5000[00:07<；00:00，675.88it/s]

我注意到没有使用ADVI报告初始化输出，这似乎令人不安。在

我还发现，如果我强制模型使用特定的step函数（例如，我在上面对pm.sample的调用中注释掉的pm.Metropolis()），我确实可以对一些参数（但不是所有参数）进行采样。仍然没有初始化报告。在

有什么办法让它工作吗？在

Answer 1

我也有类似的问题，只是我的观察数据不是number of observationsxnumber of variables的格式。结果只使用了第一次观察。在

尝试observations = mallard[['y1','y2','y3']]

如果你想在以后的使用中取样sample_ppc

Answer 2

我认为你的模型有几个问题：

• 离散参数不适用于advi和nuts，我怀疑metropolis采样器也能处理所有这些离散参数。在许多情况下，您可以将它们边缘化，但在这里，您可能希望使用一个连续变量来表示总体规模。也许是这样的（这也考虑到人口规模低于观察到的数量是不可能的）

  sd = pm.HalfCauchy('Ni_sd', beta=2.5)
  trafo = pm.distributions.transforms.lowerbound(observations.max(axis=0))
  Ni = pm.Gamma('Ni', mu=lam, sd=sd, shape=num_transects,
                transform=trafo, testval=observations.max(axis=0) + 1)
  

• 我认为你的模型是无法确定的：你不能总是增加所有的人口规模和减少道具？我看不出这个模型怎么能学到什么有用的道具。这些信息从何而来？

• 为什么人口规模取决于调查长度？这不应该是影响力的支柱吗？还是一个区域？