使用 PyMC 拟合 3 个法线：简单数据的错误收敛

Posted 2023-03-12

【中文标题】使用 PyMC 拟合 3 个法线：简单数据的错误收敛

【英文标题】：Fitting 3 Normals using PyMC: wrong convergence on simple data

【发布时间】：2013-10-07 12:43:57

【问题描述】：

我编写了一个 PyMC 模型，用于将 3 个法线拟合到数据中（类似于 this question 中的那个）。

import numpy as np
import pymc as mc
import matplotlib.pyplot as plt

n = 3
ndata = 500

# simulated data
v = np.random.randint( 0, n, ndata)
data = (v==0)*(10+ 1*np.random.randn(ndata)) \
   + (v==1)*(-10 + 2*np.random.randn(ndata)) \
   + (v==2)*3*np.random.randn(ndata)

# the model
dd = mc.Dirichlet('dd', theta=(1,)*n)
category = mc.Categorical('category', p=dd, size=ndata)
precs = mc.Gamma('precs', alpha=0.1, beta=0.1, size=n)
means = mc.Normal('means', 0, 0.001, size=n)

@mc.deterministic
def mean(category=category, means=means):
    return means[category]

@mc.deterministic
def prec(category=category, precs=precs):
    return precs[category]

obs = mc.Normal('obs', mean, prec, value=data, observed = True)

model = mc.Model('dd': dd,
              'category': category,
              'precs': precs,
              'means': means,
              'obs': obs)

M = mc.MAP(model)
M.fit()
# mcmc sampling
mcmc = mc.MCMC(model)
mcmc.use_step_method(mc.AdaptiveMetropolis, model.means)
mcmc.use_step_method(mc.AdaptiveMetropolis, model.precs)
mcmc.sample(100000,burn=0,thin=10)

tmeans = mcmc.trace('means').gettrace()
tsd = mcmc.trace('precs').gettrace()**-.5
plt.plot(tmeans)
#plt.errorbar(range(len(tmeans)), tmeans, yerr=tsd)
plt.show()

我从中采样数据的分布明显重叠，但有 3 个明显不同的峰（见下图）。将 3 个法线拟合到此类数据应该是微不足道的，我希望它能够在 99% 的 MCMC 运行中产生我从 (-10, 0, 10) 采样的均值。 我期望的结果示例。这发生在十分之二的案例中。 10 个案例中有 6 个发生意外结果的示例。这很奇怪，因为在 -5 上，数据中没有峰值，所以我不能真正找到采样可能陷入的严重局部最小值（从 (-5,-5) 到 (-6,-4)应该提高合身性，等等）。

在大多数情况下（自适应 Metropolis）MCMC 采样卡住的原因可能是什么？有哪些可能的方法来改进它没有的采样程序？

所以运行确实收敛，但并没有真正探索正确的范围。

---

更新： Using different priors，我在 5/10 得到了正确的收敛（大约第一张图片），而在另一个 5/10 得到了错误的收敛（大约第二张图片）。基本上，更改的行是下面的行，并删除了 AdaptiveMetropolis 步骤方法：

precs = mc.Gamma('precs', alpha=2.5, beta=1, size=n)
means = mc.Normal('means', [-5, 0, 5], 0.0001, size=n)

【问题讨论】：

能否也提供型号代码？我想我知道为什么会发生这种情况，但我也想要完整的上下文。

型号代码与链接完全相同。

啊，网址被隐藏了 cc: @***（我希望这有效）

您是否可以将原始代码放在问题中或另一个馅饼中？我认为您的链接包含更新的代码，我正在尝试根据 cmets 和答案重建您最初拥有的内容（也许我应该编写一个 PyMC 模型来帮助我）。

完成。发布的代码显示了与解释完全相同的行为，最后是具有不同先验的更新代码的链接（在 50% 的情况下仍然只产生正确的结果）。

【参考方案1】：

您想使用AdaptiveMetropolis 有什么特别的原因吗？我想香草MCMC 不起作用，你得到了这样的东西：

是的，那不好。我可以做几个cmets。下面我用的是原版 MCMC。

您的means 先前方差0.001 太大。这对应于大约 31 ( = 1/sqrt(0.001) ) 的标准偏差，这太小了。你真的强迫你的手段接近0。你想要一个更大的标准。偏差以帮助探索该地区。我将值减小到 0.00001 并得到了这个：

完美。当然，我先验地知道真正的均值是 50,0 和 -50。通常我们不知道这一点，因此将该值设置得非常小总是一个好主意。

2。你真的认为所有的法线都在 0 处排列，就像你之前的 mean 建议的那样吗？ （你将它们的平均值设置为 0）这个练习的重点是发现它们是不同的，所以你的先验应该反映这一点。类似的东西：

means = mc.Normal('means', [-5,0,5], 0.00001, size=n)

更准确地反映了您的真实信念。这实际上也有助于通过向 MCMC 建议方法应该在哪里来帮助收敛。当然，您必须使用您的最佳估计来得出这些数字（我在这里天真地选择了 -5,0,5）。

【讨论】：

即使更改了之前，它也会爆炸 10 次中的 9 次（注意：我犯了一个错误，数据均值是 [-10,0,10]；但 50 年代的结果是相同的）。我现在使用的确切代码链接在这里：pastebin.com/qauDYXM8

【参考方案2】：

该问题是由category 变量的低接受率引起的。请参阅answer 我给了一个类似的问题。