Metropolis-specific TypeError: The broadcastable pattern of the input is incorrect for this op

Question

我正在尝试在 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).

我确实找到了 a GitHub issue filed for something along these lines，但我不清楚 "workaround" 有人为他们的特定模型提出的建议在我的情况下会如何翻译。

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

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

将 2D 张量 "extrudes" 转换为匹配所需输出形状的 3D 张量。

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

Answer 1

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

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

仅从这个玩具示例来看，尚不清楚是否还需要展平其他层次结构级别。

Metropolis-specific TypeError: The broadcastable pattern of the input is incorrect for this op

Metropolis-specific TypeError: The broadcastable pattern of the input is incorrect for this op

python

theano

pymc3

tensor

hierarchical-bayesian