代码之家 › 专栏 › 技术社区 › merv

metropolis-specific-typeerror:输入的可广播模式对此操作不正确

hierarchical-bayesian pymc3 pymc tensor theano

0

merv · 技术社区 · 7 年前

我试图在pymc3中建立一个多层次、多维的贝叶斯模型。对于这个问题,我将使用具有以下图形结构的较小玩具模型:

在哪里? G 代表基因, K 单元类型,以及 C_k 细胞类型的细胞 k 是的。总的来说,该模型代表从不同细胞类型的细胞集合中取样的基因转录本,其中存在一些由细胞类型平均表达水平参数化的二项式分布, mu_gk ,以及特定于细胞的捕获效率, p_kc 是的。

当我用坚果对这个玩具模型进行抽样时,它做得很好,并恢复了合理的后验分布:

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)

然而,当我尝试与大都会,例如,

    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).

我确实发现 a GitHub issue filed for something along these lines ,但我不清楚有人为他们的特定模型提出的“变通方法”在我的案例中会如何翻译。

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

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

它将二维张量“挤出”成与所需输出形状匹配的三维张量。

如何实施这一模式,以避免大都市运行中的错误?

1 回复 | 直到 7 年前

1

0

merv 7 年前

使所有参数采样步骤变平似乎足以用于此玩具模型:

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 张量形式似乎很好,可能是因为它只涉及似然计算,而不是实际的采样步骤。

仅从这个玩具例子来看,还不清楚是否还需要平展额外的层次结构级别。