使用 PYMC3 对 RV 求和
Summing RVs using PYMC3
我正在尝试从图像中实现模型。我是 PyMC3 的新手,我不确定如何正确构建模型。我的尝试如下:
# sample data
logprem = np.array([8.66768002, 8.49862181, 8.60410456, 8.54966038, 8.55910259,
8.56216656, 8.51559191, 8.60630237, 8.56140145, 8.50956416])
with Model() as model:
logelr = Normal('logelr', -0.4, np.sqrt(10), shape=10)
α0 = 0
β9 = 0
α = Normal('α', 0, sigma=np.sqrt(10), shape=9)
β = Normal('β', 0, sigma=np.sqrt(10), shape=9)
a = Uniform('a', 0, 1, shape=10)
σ0 = a[0] + a[1] + a[2] + a[3] + a[4] + a[5] + a[6] + a[7] + a[8] + a[9]
σ1 = a[1] + a[2] + a[3] + a[4] + a[5] + a[6] + a[7] + a[8] + a[9]
σ2 = a[2] + a[3] + a[4] + a[5] + a[6] + a[7] + a[8] + a[9]
σ3 = a[3] + a[4] + a[5] + a[6] + a[7] + a[8] + a[9]
σ4 = a[4] + a[5] + a[6] + a[7] + a[8] + a[9]
σ5 = a[5] + a[6] + a[7] + a[8] + a[9]
σ6 = a[6] + a[7] + a[8] + a[9]
σ7 = a[7] + a[8] + a[9]
σ8 = a[8] + a[9]
σ9 = a[9]
σ = [σ0, σ1, σ2, σ3, σ4, σ5, σ6, σ7, σ8, σ9]
for w in range(10):
for d in range(10):
if w == 0:
if d == 9:
μ = logprem[w] + logelr + α0 + β9
else:
μ = logprem[w] + logelr + α0 + β[d]
else:
if d == 9:
μ = logprem[w] + logelr + α[w-1] + β9
else:
μ = logprem[w] + logelr + α[w-1] + β[d]
C = Lognormal('C', μ, σ[d])
运行 这会导致 TypeError
TypeError: For compute_test_value, one input test value does not have the requested type.
The error when converting the test value to that variable type:
Wrong number of dimensions: expected 0, got 1 with shape (10,).
如何用正确的形状定义 C?
C
的正确形状是 (W,D)
并且由于这一切都是基于 Tensor 对象的 Theano 计算图,因此最好避免循环并将自己限制在 theano.tensor
operations。这是一个类似这样的实现:
import numpy as np
import pymc3 as pm
import theano.tensor as tt
D = 10
W = 10
# begin with (W,1) vector, then broadcast to (W,D)
logprem = tt._shared(
np.array(
[[8.66768002, 8.49862181, 8.60410456, 8.54966038, 8.55910259,
8.56216656, 8.51559191, 8.60630237, 8.56140145, 8.50956416]]) \
.T \
.repeat(D, axis=1))
with pm.Model() as model:
logelr = pm.Normal('logelr', -0.4, np.sqrt(10))
# col vector
alpha = pm.Normal("alpha", 0, sigma=np.sqrt(10), shape=(W-1,1))
# row vector
beta = pm.Normal("beta", 0, sigma=np.sqrt(10), shape=(1,D-1))
# prepend zero and broadcast to (W,D)
alpha_aug = tt.concatenate([tt.zeros((1,1)), alpha], axis=0).repeat(D, axis=1)
# append zero and broadcast to (W,D)
beta_aug = tt.concatenate([beta, tt.zeros((1,1))], axis=1).repeat(W, axis=0)
# technically the indices will be reversed
# e.g., a[0], a[9] here correspond to a_10, a_1 in the paper, resp.
a = pm.Uniform('a', 0, 1, shape=D)
# Note: the [::-1] sorts it in the order specified
# such that (sigma_0 > sigma_1 > ... )
sigma = pm.Deterministic('sigma', tt.extra_ops.cumsum(a)[::-1].reshape((1,D)))
# now everything here has shape (W,D) or is scalar (logelr)
mu = logprem + logelr + alpha_aug + beta_aug
# sigma will be broadcast automatically
C = pm.Lognormal('C', mu=mu, sigma=sigma, shape=(W,D))
关键技巧是
- 将零添加到
alpha
和 beta
上,让所有内容都保持在张量形式中
- 使用
tt.extra_ops.cumsum
方法简明表达步骤5;
- 使步骤 6 中的所有项具有形状 (W,D)
如果可以使用加法运算符在 alpha
和 beta
向量之间执行外积(例如,在 R 中,outer
函数允许任意ops) 但我在 theano.tensor
.
下找不到这样的方法
使用 NUTS 并不能很好地采样,但是一旦您实际观察到要插入的 C
的值,它可能会更好。
with model:
# using lower target_accept and tuning throws divergences
trace = pm.sample(tune=5000, draws=2000, target_accept=0.99)
pm.summary(trace, var_names=['alpha', 'beta', 'a', 'sigma'])
因为这只是先验采样,唯一真正有趣的是转换变量的分布 sigma
:
pm.plot_forest(trace, var_names=['sigma'])
哪个可以看到符合sigma_{d} > sigma_{d+1}
.
的要求
我正在尝试从图像中实现模型。我是 PyMC3 的新手,我不确定如何正确构建模型。我的尝试如下:
# sample data
logprem = np.array([8.66768002, 8.49862181, 8.60410456, 8.54966038, 8.55910259,
8.56216656, 8.51559191, 8.60630237, 8.56140145, 8.50956416])
with Model() as model:
logelr = Normal('logelr', -0.4, np.sqrt(10), shape=10)
α0 = 0
β9 = 0
α = Normal('α', 0, sigma=np.sqrt(10), shape=9)
β = Normal('β', 0, sigma=np.sqrt(10), shape=9)
a = Uniform('a', 0, 1, shape=10)
σ0 = a[0] + a[1] + a[2] + a[3] + a[4] + a[5] + a[6] + a[7] + a[8] + a[9]
σ1 = a[1] + a[2] + a[3] + a[4] + a[5] + a[6] + a[7] + a[8] + a[9]
σ2 = a[2] + a[3] + a[4] + a[5] + a[6] + a[7] + a[8] + a[9]
σ3 = a[3] + a[4] + a[5] + a[6] + a[7] + a[8] + a[9]
σ4 = a[4] + a[5] + a[6] + a[7] + a[8] + a[9]
σ5 = a[5] + a[6] + a[7] + a[8] + a[9]
σ6 = a[6] + a[7] + a[8] + a[9]
σ7 = a[7] + a[8] + a[9]
σ8 = a[8] + a[9]
σ9 = a[9]
σ = [σ0, σ1, σ2, σ3, σ4, σ5, σ6, σ7, σ8, σ9]
for w in range(10):
for d in range(10):
if w == 0:
if d == 9:
μ = logprem[w] + logelr + α0 + β9
else:
μ = logprem[w] + logelr + α0 + β[d]
else:
if d == 9:
μ = logprem[w] + logelr + α[w-1] + β9
else:
μ = logprem[w] + logelr + α[w-1] + β[d]
C = Lognormal('C', μ, σ[d])
运行 这会导致 TypeError
TypeError: For compute_test_value, one input test value does not have the requested type.
The error when converting the test value to that variable type:
Wrong number of dimensions: expected 0, got 1 with shape (10,).
如何用正确的形状定义 C?
C
的正确形状是 (W,D)
并且由于这一切都是基于 Tensor 对象的 Theano 计算图,因此最好避免循环并将自己限制在 theano.tensor
operations。这是一个类似这样的实现:
import numpy as np
import pymc3 as pm
import theano.tensor as tt
D = 10
W = 10
# begin with (W,1) vector, then broadcast to (W,D)
logprem = tt._shared(
np.array(
[[8.66768002, 8.49862181, 8.60410456, 8.54966038, 8.55910259,
8.56216656, 8.51559191, 8.60630237, 8.56140145, 8.50956416]]) \
.T \
.repeat(D, axis=1))
with pm.Model() as model:
logelr = pm.Normal('logelr', -0.4, np.sqrt(10))
# col vector
alpha = pm.Normal("alpha", 0, sigma=np.sqrt(10), shape=(W-1,1))
# row vector
beta = pm.Normal("beta", 0, sigma=np.sqrt(10), shape=(1,D-1))
# prepend zero and broadcast to (W,D)
alpha_aug = tt.concatenate([tt.zeros((1,1)), alpha], axis=0).repeat(D, axis=1)
# append zero and broadcast to (W,D)
beta_aug = tt.concatenate([beta, tt.zeros((1,1))], axis=1).repeat(W, axis=0)
# technically the indices will be reversed
# e.g., a[0], a[9] here correspond to a_10, a_1 in the paper, resp.
a = pm.Uniform('a', 0, 1, shape=D)
# Note: the [::-1] sorts it in the order specified
# such that (sigma_0 > sigma_1 > ... )
sigma = pm.Deterministic('sigma', tt.extra_ops.cumsum(a)[::-1].reshape((1,D)))
# now everything here has shape (W,D) or is scalar (logelr)
mu = logprem + logelr + alpha_aug + beta_aug
# sigma will be broadcast automatically
C = pm.Lognormal('C', mu=mu, sigma=sigma, shape=(W,D))
关键技巧是
- 将零添加到
alpha
和beta
上,让所有内容都保持在张量形式中 - 使用
tt.extra_ops.cumsum
方法简明表达步骤5; - 使步骤 6 中的所有项具有形状 (W,D)
如果可以使用加法运算符在 alpha
和 beta
向量之间执行外积(例如,在 R 中,outer
函数允许任意ops) 但我在 theano.tensor
.
使用 NUTS 并不能很好地采样,但是一旦您实际观察到要插入的 C
的值,它可能会更好。
with model:
# using lower target_accept and tuning throws divergences
trace = pm.sample(tune=5000, draws=2000, target_accept=0.99)
pm.summary(trace, var_names=['alpha', 'beta', 'a', 'sigma'])
因为这只是先验采样,唯一真正有趣的是转换变量的分布 sigma
:
pm.plot_forest(trace, var_names=['sigma'])
哪个可以看到符合sigma_{d} > sigma_{d+1}
.