PyMC3 用时变捕获概率估计封闭种群的大小
PyMC3 estimating size of closed population with time-varying capture probabilities
我正在尝试将 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
吗?
以下是我如何将此模型转换为 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)。
我正在尝试将 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
吗?
以下是我如何将此模型转换为 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)。