zidingxiangyu 2018-07-11
马尔可夫链蒙特卡罗(MCMC)是贝叶斯统计中广泛流行的技术。它用于后验分布采样,因为分析形式通常是不可跟踪的。但是,在这篇文章中,我们将使用它从静态图像/徽标生成动画。顺便提一下,它可以作为MCMC和拒绝采样的介绍。这个想法基于一个伟大的开源软件包imcmc,它基于PyMC3构建。
MCMC适用于任何维度,我们将重点关注2D发行版。为什么?看看下面的图表。
我们从RGB图像(最好是一个简单的标志)开始,将其转换为灰度,并根据强度为每个像素分配一个概率。在我们的例子中,像素越黑,概率就越高。这导致了离散的二维分布。
一旦我们有了这个分布,我们就可以从中抽取样本(像素)。一般来说,我们寻找以下两个属性:
第二个性质总是在旁观者的眼中,第一个性质可以通过适当的采样算法得到。具体地说,让我们来仔细看看3种非常适合当前问题的抽样算法:Rejection, Gibbs和 Metropolis-Hastings抽样。
第一种方法属于IID(独立且同分布)抽样方法的一类。这本质上意味着当前样本对下一个样本没有影响。该算法假设我们能够从所谓的提案分配中抽取样本。可以使用许多不同的提议分布,但最常见的是Uniform(有限支持)和Normal(无限支持)分布。
为了尽量减少拒绝样本的可能性,我们必须选择一个与我们的目标相似的提案分配。同样,我们希望标度常数M越低越好。见下图(1 d)计划。
在我们的设置中,我们可以采用二维均匀分布作为提案
def rejection_sampling(image, approx_samples):
image_pdf = image / image.sum()
pdf_max = image_pdf.max()
height, width = image_pdf.shape
p_success = 1 / (height * width * pdf_max)
actual_samples = min(int(approx_samples / p_success), int(1e8))
samples_height = np.random.randint(0, high=height, size=actual_samples)
samples_width = np.random.randint(0, high=width, size=actual_samples)
samples_uniform = np.random.uniform(0, 1, size=actual_samples)
result = [(h, w) for (h, w, u) in zip(samples_height, samples_width, samples_uniform) if
(image_pdf[h, w] >= pdf_max * u)]
return result
在较低的维度,拒绝采样表现非常好。然而,随着维数的增加,它遭受了臭名昭着的维度诅咒。对我们来说幸运的是,2D并没有被诅咒,因此拒绝采样是一个很好的选择。
吉布斯抽样属于通过构造马尔可夫链生成样本的第二种抽样类型。因此,这些样本不是独立的。事实上,直到链达到其平稳分布,它们才会均匀分布。出于这个原因,通常的做法是丢弃前x个样本以确保链“忘记”初始化(burn-in)。
算法本身假设我们能够沿着目标的每个维度从条件分布中提取样本。每当我们从这些单变量分布中进行采样时,我们就会考虑剩余分量的最新值。
def gibbs_sampling(image, w_start, samples):
image_pdf = image / image.sum()
height, width = image_pdf.shape
result = []
w_current = w_start
for _ in range(samples):
# sample height
h_given_w = image_pdf[:, w_current] / image_pdf[:, w_current].sum()
h_current = np.random.choice(np.array(range(height)), size=1, p=h_given_w)[0]
# sample width
w_given_h = image_pdf[h_current, :] / image_pdf[h_current, :].sum()
w_current = np.random.choice(np.array(range(width)), size=1, p=w_given_h)[0]
result.append((h_current, w_current))
return result
与吉布斯类似,Metropolis-Hastings的抽样也创造了马尔可夫链。但是,它更通用(吉布斯是一个特例)并且灵活。
考虑到目前的样本,提案的分发给了我们一个新的样本的建议。然后,我们通过检查它比当前样本多多少(少多少)的可能性来评估它的资格,并考虑到提案分布对这个样本可能存在的偏差。所有东西结合起来,我们计算接受新样本的概率然后让随机性来决定。
毫无疑问,提案分配的作用是至关重要的,会影响性能和收敛性。更常见的选择之一是使用以当前状态为中心的Normal。
让我们看看一些结果。对于每个徽标,运行拒绝和Metropolis-Hastings(使用PyMC3)并创建可视化(gif)。最重要的超参数:
对于多模式分布,Metropolis-Hastings可能会陷入某些区域。这是一个非常常见的问题,可以通过更改提案或对更多/更长的链进行抽样来解决。
如果你想深入挖掘并查看源代码,可以访问(https://github.com/jankrepl/creating-animations-with-MCMC/blob/master/main.ipynb)