机器算法验证 - 在 PyMC 中使用经验先验 - 吾爱随笔录

在 PyMC 中使用经验先验

机器算法验证贝叶斯马尔可夫链蒙特卡罗事先的后部 pymc

2022-04-11 12:06:49

我正在使用 PyMC 对后验分布进行采样，并且在使用来自样本而不是模型的先验时遇到了障碍。我的情况如下：

我有一些参数的经验数据 $z$ 我从中计算概率分布 $p(z)$ . 没有已知的模型/参数化分布 $z$ . 我所拥有的只是经验值。众所周知 $z$ 范围在 0 到 30 之间。
我有一些新的观察 $x$ 我计算可能性 $p(x|z)$ , 也是凭经验的。
我希望找到用上面的新数据更新的后验分布。（这个后验成为下一组观察的新先验。冲洗并重复。）

例如，考虑一下：

import numpy as np
import matplotlib.pyplot as plt

# Ignore the fact that I'm using a mixture model. For all practical
# purposes, I do not know how this is generated. 
old_data = np.array([3 * np.random.randn(1000) + 20, 
                     3 * np.random.randn(1000) + 5,
                     4 * np.random.randn(1000) + 10])

new_data = np.array([4 * np.random.randn(50) + 8, 
                     2 * np.random.randn(50) + 17])

plt.hist(filter(lambda x: 0 <= x <= 30, old_data.flatten()), 
         bins=range(0, 30), normed=True, alpha=0.5)
plt.hist(filter(lambda x: 0 <= x <= 30, new_data.flatten()), 
         bins=range(0, 30), normed=True, alpha=0.5)
plt.legend(['p(z)', 'p(x|z)'])

来到 PyMC，我发现的所有现有资源（例如，来自贝叶斯黑客的这一章）对观察值 ( observations = pm.Normal("obs", center_i, tau_i, value=data, observed=True)) 使用正态分布，对精度和均值分别使用均匀和正态先验分布。

我不确定如何在 PyMC 中断言我的先验是这里的分布而不是模型。我也尝试过将@Stochastic装饰器与一起使用observed=True，但我认为我并不完全理解它。此外，我似乎仍然无法找到避免指定模型的方法。

我是否从根本上误解了 MCMC 库的目的？如果是这样，我应该如何进行？

我真正想要的只是用新的观察来更新我先前的信念，但我认为解决方案并不像将两个直方图相乘（和归一化）那么简单。

1个回答

如果你已经有一个 $p(\theta)$ 和一个可能性 $p(x|\theta)$ ，那么你可以很容易地找到后验 $p(\theta|x)$ 通过将这些相乘并归一化：

p (θ | x) = \frac{p (θ) p (x | θ)}{p (x)} \propto p (θ) p (x | θ)

$p(\theta|x)=\frac{p(\theta)p(x|\theta)}{p(x)}\propto p(\theta)p(x|\theta)$

https://en.wikipedia.org/wiki/Posterior_probability

以下代码演示了估计表示为直方图的后验，因此它可以用作下一个先验：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

# using Beta distribution instead on Normal to get finite support
support_size=30
old_data=np.concatenate([np.random.beta(70,20,1000), 
                     np.random.beta(10,40,1000),
                     np.random.beta(80,80,1000)])*support_size
new_data=np.concatenate([np.random.beta(20,10,1000), 
                     np.random.beta(10,20,1000)])*support_size

# convert samples to histograms
support=np.arange(support_size)
old_hist=np.histogram(old_data,bins=support,normed=True)[0]
new_hist=np.histogram(new_data,bins=support,normed=True)[0]

# obtain smooth estimators from samples
soft_old=gaussian_kde(old_data,bw_method=0.1)
soft_new=gaussian_kde(new_data,bw_method=0.1)

# posterior histogram (to be used as a prior for the next batch)
post_hist=old_hist*new_hist
post_hist/=post_hist.sum()

# smooth posterior
def posterior(x):
    return soft_old(x)*soft_new(x)/np.sum(soft_old(x)*soft_new(x))*x.size/support_size

x=np.linspace(0,support_size,100)

plt.bar(support[:-1],old_hist,alpha=0.5,label='p(z)',color='b')
plt.bar(support[:-1],new_hist,alpha=0.5,label='p(x|z)',color='g')
plt.plot(x,soft_old(x),label='p(z) smoothed',lw=2)
plt.plot(x,soft_new(x),label='p(x|z) smoothed',lw=2)
plt.legend(loc='best',fontsize='small')
plt.show()

plt.bar(support[:-1],post_hist,alpha=0.5,label='p(z|x)',color='r')
plt.plot(x,soft_old(x),label='p(z) smoothed',lw=2)
plt.plot(x,soft_new(x),label='p(x|z) smoothed',lw=2)
plt.plot(x,posterior(x),label='p(z|x) smoothed',lw=2)
plt.legend(loc='best',fontsize='small')
plt.show()

但是，如果您想将经验先验与一些 MCMC 模型结合起来，我建议您看一下PyMC 的潜力，它的主要应用之一是“软数据”。如果您需要针对该问题的答案，请更新您的问题。

其它你可能感兴趣的问题

上一篇从单个样本估计几何分布的参数下一篇随机分布偏差背后的机制