0

我正在尝试在 PyMC3 中构建一个多层次、多维的贝叶斯模型。对于这个问题，我将使用具有以下图形结构的较小玩具模型： 成绩单捕获效率的图形模型

其中G代表基因、K细胞类型和C_k细胞类型的细胞k。总体而言，该模型表示从不同细胞类型的细胞集合中采样的基因转录物，其中存在一些由细胞类型平均表达水平参数化的二项式分布mu_gk，以及细胞特异性捕获效率，p_kc。

当我用 NUTS 对这个玩具模型进行采样时，它表现良好并恢复了合理的后验分布：

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


# Generative model for data simulation
def sample_data(G=1, K=2, C_k=100):
    mu_gk = np.random.randint(1, 1000, size=(G, K))
    p_kc = np.random.beta(5, 95, (K, C_k))
    N_gkc = np.random.binomial(mu_gk[:, :, np.newaxis], p_kc[np.newaxis, :, :])

    return N_gkc


G = 10    # genes
K = 5     # cell types
C_k = 20  # cells per type

data = sample_data(G, K, C_k)

with pm.Model() as capture_efficiency:

    # Genes expression levels per cell type
    mu_gk = pm.Uniform('mu_gk', 1, 1000, shape=(G, K, 1))

    # Cell capture efficiencies
    p_kc = pm.Beta('p_kc', shape=(1, K, C_k), alpha=5, beta=95)

    # Captured transcripts
    N_gkc = pm.Binomial('N_gkc', shape=(G, K, C_k),
                        n=tt.tensordot(mu_gk, np.ones((C_k, 1)), [[2], [1]]),
                        p=tt.tensordot(np.ones((G, 1)), p_kc, [[1], [0]]),
                        observed=data)

    trace = pm.sample(5000, tune=10000, target_accept=0.99)

但是，当我尝试使用 Metropolis 进行采样时，例如，

    trace = pm.sample(5000, tune=10000, step=pm.Metropolis())

我收到以下堆栈跟踪和错误消息：

Traceback (most recent call last):
  File "/Applications/PyCharm.app/Contents/helpers/pydev/pydev_run_in_console.py", line 52, in run_file
    pydev_imports.execfile(file, globals, locals)  # execute the script
  File "/Applications/PyCharm.app/Contents/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "/Users/mfansler/Projects/pymc3/intro/capture-efficiency-celltypes.py", line 46, in <module>
    trace = pm.sample(5000, tune=10000,  step=pm.Metropolis(),
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/step_methods/arraystep.py", line 65, in __new__
    step.__init__([var], *args, **kwargs)
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py", line 136, in __init__
    self.delta_logp = delta_logp(model.logpt, vars, shared)
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py", line 624, in delta_logp
    [logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared)
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/theanof.py", line 245, in join_nonshared_inputs
    xs_special = [theano.clone(x, replace, strict=False) for x in xs]
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/theanof.py", line 245, in <listcomp>
    xs_special = [theano.clone(x, replace, strict=False) for x in xs]
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/scan_module/scan_utils.py", line 247, in clone
    share_inputs)
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 232, in rebuild_collect_shared
    cloned_v = clone_v_get_shared_updates(outputs, copy_inputs_over)
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 93, in clone_v_get_shared_updates
    clone_v_get_shared_updates(i, copy_inputs_over)
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 93, in clone_v_get_shared_updates
    clone_v_get_shared_updates(i, copy_inputs_over)
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 93, in clone_v_get_shared_updates
    clone_v_get_shared_updates(i, copy_inputs_over)
  [Previous line repeated 9 more times]
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 96, in clone_v_get_shared_updates
    [clone_d[i] for i in owner.inputs], strict=rebuild_strict)
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/gof/graph.py", line 246, in clone_with_new_inputs
    new_node = self.op.make_node(*new_inputs)
  File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/tensor/elemwise.py", line 230, in make_node
    % (self.input_broadcastable, ib)))
TypeError: The broadcastable pattern of the input is incorrect for this op. Expected (False, False, True), got (False, False, False).

我确实发现了一个 GitHub issue 提交了一些类似的内容，但我不清楚有人为他们的特定模型提出的“解决方法”在我的案例中将如何转化。

我怀疑这个模型与遇到的错误最相关的方面是在实例化二项式随机变量时手动广播参数：

n=tt.tensordot(mu_gk, np.ones((C_k, 1)), [[2], [1]]),
p=tt.tensordot(np.ones((G, 1)), p_kc, [[1], [0]])

它将 2D 张量“挤压”成与所需输出形状匹配的 3D 张量。

这个模型应该如何实现才能避免在运行 Metropolis 时出现错误？

4

1 回答 1

0

对于这个玩具模型，使所有参数采样步骤平坦似乎就足够了：

with pm.Model() as capture_efficiency:

    # Genes expression levels per cell type
    mu_gk_flat = pm.Uniform('mu_gk', 1, 1000, shape=G*K)
    mu_gk = mu_gk_flat.reshape((G, K, 1))

    # Cell capture efficiencies
    p_kc_flat = pm.Beta('p_kc', shape=K*C_k, alpha=5, beta=95)
    p_kc = p_kc_flat.reshape((1, K, C_k))

    # Captured transcripts
    N_gkc = pm.Binomial('N_gkc', shape=(G, K, C_k),
                        n=tt.tensordot(mu_gk, np.ones((C_k, 1)), [[2], [1]]),
                        p=tt.tensordot(np.ones((G, 1)), p_kc, [[1], [0]]),
                        observed=data)

    trace = pm.sample(5000, tune=10000, step=pm.Metropolis())

N_gkc张量形式似乎很好，也许是因为它只涉及似然计算而不是实际的采样步骤。

仅从这个玩具示例来看，尚不清楚是否还需要展平额外的层级。

于 2018-07-17T02:16:48.757 回答