使用PYMC3对RV求和 [英] Summing RVs using PYMC3

问题描述

我正在尝试从图像中实现模型.我是PyMC3的新手，我不确定如何正确构建模型.我的尝试如下:

I am attempting to implement the model from the image. I am new to PyMC3 and I am not sure how to structure the model correctly. My attempt is below:

# 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

Running this leads to a 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?

How to I define C with the correct shape?

推荐答案

C的正确形状是(W,D)，由于这一切都是基于Tensor对象的Theano计算图，因此最好避免循环并限制自己 theano.tensor操作.这是一种类似的实现方式:

The correct shape for C is (W,D) and since underlying this all is a Theano computational graph on Tensor objects, it is best to avoid loops and restrict oneself to theano.tensor operations. Here is an implementation along such lines:

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)

• prepending and appending the zeros onto alpha and beta, allow everything to be kept in a tensor form
• using the tt.extra_ops.cumsum method to concisely express step 5;
• getting all the terms in Step 6 to have the shape (W,D)

如果可以使用加法运算符在alpha和beta向量之间执行外积，则可以进一步简化此操作(例如，在R中，outer函数允许任意操作)，但我找不到theano.tensor下的这种方法.

This could be simplified even further if one could perform an outer product between the alpha and beta vectors using an addition operator (e.g., in R the outer function allows arbitrary ops) but I couldn't find such a method under theano.tensor.

使用NUTS并不能很好地进行采样，但是一旦您实际观察到要插入的C值，也许会更好.

This doesn't really sample well using NUTS, but perhaps it might be better once you have actually observed values of C to plug in.

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的分布:

Since this is just the prior sampling, the only thing actually interesting are the distributions of the transformed variable sigma:

pm.plot_forest(trace, var_names=['sigma'])

与sigma_{d} > sigma_{d+1}的要求一致.

这篇关于使用PYMC3对RV求和的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！