PyMC：逐步采样？

Question

我想知道为什么采样器在逐步采样时速度非常慢。 例如，如果我运行：

mcmc = MCMC(model)
mcmc.sample(1000)

采样速度很快。但是，如果我运行：

mcmc = MCMC(model)
for i in arange(1000):
    mcmc.sample(1)

采样速度较慢（采样越多，速度越慢）。

如果您想知道我为什么要问这个.. 好吧，我需要逐步采样，因为我想在采样器的每一步之后对变量的值执行一些操作。

有没有办法加快速度？

提前谢谢你！

------------------ 编辑 ---------------------------- ----------------------------------

这里我更详细地介绍了具体问题：

我有两个模型在竞争，它们是一个更大模型的一部分，该模型有一个分类变量作为两者之间的“开关”。

在这个玩具示例中，我有观察到的向量“Y”，可以用泊松或几何分布来解释。分类变量“switch_model”在 = 0 时选择几何模型，在 =1 时选择泊松模型。

在每个样本之后，如果 switch_model 选择几何模型，我希望泊松模型的变量不要更新，因为它们不会影响可能性，因此它们只是逐渐消失。如果 switch_model 选择 Poisson 模型，则相反。

基本上，我在每个步骤中所做的就是通过手动将其后退一步来“更改”未选择模型的值。

我希望我的解释和注释代码足够清楚。如果您需要更多详细信息，请告诉我。

import numpy as np
import pymc as pm
import pandas as pd
import matplotlib.pyplot as plt

# OBSERVED VALUES
Y = np.array([0, 1, 2, 3, 8])

# PRIOR ON THE MODELS
pi = (0.5, 0.5)

switch_model = pm.Categorical("switch_model", p = pi) 
# switch_model = 0 for Geometric, switch_model = 1 for Poisson

p = pm.Uniform('p', lower = 0, upper = 1) # Prior of the parameter of the geometric distribution
mu = pm.Uniform('mu', lower = 0, upper = 10) # Prior of the parameter of the Poisson distribution


# LIKELIHOOD
@pm.observed
def Ylike(value = Y, mu = mu, p = p, M = switch_model):
    if M == 0:
        out = pm.geometric_like(value+1, p)
    elif M == 1:
        out = pm.poisson_like(value, mu)
    return out

model = pm.Model([Ylike, p, mu, switch_model])

mcmc = pm.MCMC(model)

n_samples = 5000

traces = {}
for var in mcmc.stochastics:
    traces[str(var)] = np.zeros(n_samples)


bar = pm.progressbar.progress_bar(n_samples)
bar.update(0)

mcmc.sample(1, progress_bar=False)
for var in mcmc.stochastics:
    traces[str(var)][0] = mcmc.trace(var)[-1]


for i in np.arange(1,n_samples):
    mcmc.sample(1, progress_bar=False)
    bar.update(i)
    for var in mcmc.stochastics:
        traces[str(var)][i] = mcmc.trace(var)[-1]
    if mcmc.trace('switch_model')[-1] == 0: # Gemetric wins
        traces['mu'][i] = traces['mu'][i-1] # One step back for the sampler of the Poisson parameter
        mu.value = traces['mu'][i-1]

    elif mcmc.trace('switch_model')[-1] == 1: # Poisson wins
        traces['p'][i] = traces['p'][i-1] # One step back for the sampler of the Geometric parameter
        p.value = traces['p'][i-1]

print '\n\n'

traces=pd.DataFrame(traces)

traces['mu'][traces['switch_model'] == 0] = np.nan
traces['p'][traces['switch_model'] == 1] = np.nan

print traces.describe()

traces.plot()
plt.show()

Answer 1

之所以这么慢，是因为 Python 的 for 循环非常慢，尤其是与 FORTRAN 循环相比时（PyMC 基本上是用它编写的。）如果您可以显示更详细的代码，它可能是更容易看到您正在尝试做什么并提供更快的替代算法。

Answer 2

实际上，我找到了一个“疯狂”的解决方案，我怀疑它为什么会起作用。我仍然想就我的把戏获得专家意见。

基本上，如果我按以下方式修改 for 循环，每 1000 个循环添加一个“重置 mcmc”，采样就会再次启动：

for i in np.arange(1,n_samples):
    mcmc.sample(1, progress_bar=False)
    bar.update(i)
    for var in mcmc.stochastics:
        traces[str(var)][i] = mcmc.trace(var)[-1]
    if mcmc.trace('switch_model')[-1] == 0: # Gemetric wins
        traces['mu'][i] = traces['mu'][i-1] # One step back for the sampler of the Poisson parameter
        mu.value = traces['mu'][i-1]

    elif mcmc.trace('switch_model')[-1] == 1: # Poisson wins
        traces['p'][i] = traces['p'][i-1] # One step back for the sampler of the Geometric parameter
        p.value = traces['p'][i-1]

    if i%1000 == 0:
        mcmc = pm.MCMC(model)

在实践中，此技巧每 1000 步擦除采样器的轨迹和数据库。看起来采样器不喜欢长数据库，虽然我不太明白为什么。 （当然1000步是任意的，太短会增加太多开销，太长会导致trace和数据库太长）。

我觉得这个 hack 有点疯狂，而且绝对不优雅。有任何专家或开发人员对此发表评论吗？谢谢！