使用 PYMC3 对 RV 求和

Question

我正在尝试从图像中实现模型。我是 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？

Answer 1

C 的正确形状是 (W,D)，由于这一切都是基于张量对象的 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 函数允许任意操作），这可以进一步简化，但我找不到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}的要求。