PyMC3 用时变捕获概率估计封闭种群的大小

Question

我正在尝试将 PyMC2 中的潜在多项式模型 this example 移植到 PyMC3。

模型描述

目的是根据T次观察的时间序列，估计封闭种群中个体的实际数量N .

每个人有 2^T 个可能的序列 - 所以如果 T=2，一个人要么永远不会被观察到 (0, 0) [我们称这个序列为 1]，只有第一次 (1, 0 ) [sequence 2]，仅第二次 (0, 1) [sequence 3] 或两次 (1, 1) [sequence 4].

我们知道一个人是否产生了序列 2 到 4，所以我们有这些计数。

但是，该模型取决于我们根据 N.

的采样值确定存在但从未被观察到的个体数量（序列 1）

在最初的 PyMC2 模型中，这是通过简单地将 N - sum(f) 附加到 f 的开头来完成的，其中 f 是观察到的序列计数的 np.array 2 到 2^T。在 PyMC3 中，这必须使用 Theano 来完成，并且其他一些变量也定义为 Theano 共享变量。

PyMC2 中的原始示例

原始 PyMC2 代码如下所示：

from pymc import Lambda, Beta, DiscreteUniform, stochastic, multinomial_like
from numpy import reshape, ma, prod, array, append
import pdb

# Constants
T=2
Cells=4
# Observations seen by both observers, just the first observer, or just the second
f=(22,363,174)
# Inidicators for observer sightings
omegas=reshape((0,0, 1,0, 0,1, 1,1), (Cells, T))

# Capture probabilities
p = Beta('p', alpha=1, beta=1, size=T)

# Cell probabilities
c = Lambda('c', lambda p=p: p**omegas * (1-p)**(1.-omegas))

# Joint probabilities
pi = Lambda('pi', lambda c=c: prod(c, 1))

# Bound total N by the number of observed individuals
N = DiscreteUniform('N', lower=sum(f), upper=10000)

@stochastic(observed=True, dtype=int)
def F(value=f, n=N, p=pi):
    """Complete multinomial for observed and unobserved quatities"""
    return multinomial_like(append(n-sum(value), value), n, p)

已尝试 PyMC3 版本

我第一次尝试 PyMC3 版本（在 Python 2.7 中）如下所示：

from __future__ import division, absolute_import, print_function

import numpy as np
import pymc3 as pm
import theano.tensor as tt
from theano import shared

# Constants

T = 2
CELLS = 4

with pm.Model() as model_mt:

    f = shared(np.array((22, 363, 174)))

    omegas = shared(np.reshape((0,0, 1,0, 0,1, 1,1), (CELLS, T)))

    p = pm.Beta('p', alpha=1, beta=1, shape=T)

    # Treating c and pi as separate steps, as per his example:
    # c = p ** omegas
    # pi = tt.prod(c, axis=1)

    # Shorter line of code:
    pi = tt.prod(p ** omegas, axis=1)

    N = pm.DiscreteUniform('N', lower=f.sum(), upper=10000)

    observed_values = tt.concatenate(([N-tt.sum(f)], f))

    y_obs = pm.Multinomial('y_obs', n=N, p=pi, observed=observed_values)

with model_mt:
    trace = pm.sample(5000)

但结果轨迹看起来难以置信。 p_0 的平均值为 0.27，p_1 的平均值为 0.98，据我所知，这是在 t=1 和 t=2 时看到一个人的值（分别），而 N 的平均值为 930。但是 p_1 不应该接近 (363+174)/N 吗？

Answer 1

以下是我如何将此模型转换为 PyMC3（我一直想这样做，谢谢）：

with pm.Model() as Mt2:

    p = pm.Uniform('p', 0, 1, shape=T)

    c = p**omega * (1 - p)**(1-omega)

    π = pm.Deterministic('π', tt.prod(c, 1))

    π_obs = π[1:] / (1 - π[0])

    f_obs = pm.Potential('f_obs', pm.Multinomial.dist(n=n, p=π_obs).logp(f[1:]))

    N = pm.Uniform('N', 0, 10000)

    n_obs = pm.Potential('n_obs', pm.Binomial.dist(n=tt.floor(N), p=1-π[0]).logp(n))

    trace = pm.sample(1000, tune=1000)

将 N 保持为连续变量并使用 partially-observed 可观察量的因子势允许 PyMC3 使用 NUTS 而不是回退到 Metropolis。

我得到的 N 的估计值刚好超过 611：

我不记得这与 PyMC2（或 Link 和 Barker）输出相比如何，但它比您的估计要小得多。

以下是我对捕获概率的估计：

p[0] 的估计值为 0.325 (0.290, 0.367)，p[1] 的估计值为 0.894 (0.853, 0.936)。

PyMC3 用时变捕获概率估计封闭种群的大小

PyMC3 estimating size of closed population with time-varying capture probabilities

python

numpy

theano

bayesian

pymc3

模型描述

PyMC2 中的原始示例

已尝试 PyMC3 版本