使用机器学习算法MCMC创建动画

zidingxiangyu 2018-07-11

马尔可夫链蒙特卡罗(MCMC)是贝叶斯统计中广泛流行的技术。它用于后验分布采样,因为分析形式通常是不可跟踪的。但是,在这篇文章中,我们将使用它从静态图像/徽标生成动画。顺便提一下,它可以作为MCMC和拒绝采样的介绍。这个想法基于一个伟大的开源软件包imcmc,它基于PyMC3构建。

初步研究

MCMC适用于任何维度,我们将重点关注2D发行版。为什么?看看下面的图表。

使用机器学习算法MCMC创建动画

我们从RGB图像(最好是一个简单的标志)开始,将其转换为灰度,并根据强度为每个像素分配一个概率。在我们的例子中,像素越黑,概率就越高。这导致了离散的二维分布。

一旦我们有了这个分布,我们就可以从中抽取样本(像素)。一般来说,我们寻找以下两个属性:

  • 样本来自目标分布
  • 样品的连续和它们的可视化在美学上是令人愉快的

第二个性质总是在旁观者的眼中,第一个性质可以通过适当的采样算法得到。具体地说,让我们来仔细看看3种非常适合当前问题的抽样算法:Rejection, Gibbs和 Metropolis-Hastings抽样。

Rejection sampling(拒绝采样)

第一种方法属于IID(独立且同分布)抽样方法的一类。这本质上意味着当前样本对下一个样本没有影响。该算法假设我们能够从所谓的提案分配中抽取样本。可以使用许多不同的提议分布,但最常见的是Uniform(有限支持)和Normal(无限支持)分布。

使用机器学习算法MCMC创建动画

为了尽量减少拒绝样本的可能性,我们必须选择一个与我们的目标相似的提案分配。同样,我们希望标度常数M越低越好。见下图(1 d)计划。

使用机器学习算法MCMC创建动画

在我们的设置中,我们可以采用二维均匀分布作为提案

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并没有被诅咒,因此拒绝采样是一个很好的选择。

吉布斯采样Gibbs Sampling

吉布斯抽样属于通过构造马尔可夫链生成样本的第二种抽样类型。因此,这些样本不是独立的。事实上,直到链达到其平稳分布,它们才会均匀分布。出于这个原因,通常的做法是丢弃前x个样本以确保链“忘记”初始化(burn-in)。

使用机器学习算法MCMC创建动画

算法本身假设我们能够沿着目标的每个维度从条件分布中提取样本。每当我们从这些单变量分布中进行采样时,我们就会考虑剩余分量的最新值。

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采样

与吉布斯类似,Metropolis-Hastings的抽样也创造了马尔可夫链。但是,它更通用(吉布斯是一个特例)并且灵活。

使用机器学习算法MCMC创建动画

考虑到目前的样本,提案的分发给了我们一个新的样本的建议。然后,我们通过检查它比当前样本多多少(少多少)的可能性来评估它的资格,并考虑到提案分布对这个样本可能存在的偏差。所有东西结合起来,我们计算接受新样本的概率然后让随机性来决定。

毫无疑问,提案分配的作用是至关重要的,会影响性能和收敛性。更常见的选择之一是使用以当前状态为中心的Normal。

动画

让我们看看一些结果。对于每个徽标,运行拒绝和Metropolis-Hastings(使用PyMC3)并创建可视化(gif)。最重要的超参数:

  • burn-in samples: 500
  • samples: 10000

使用机器学习算法MCMC创建动画

使用机器学习算法MCMC创建动画

使用机器学习算法MCMC创建动画

使用机器学习算法MCMC创建动画

使用机器学习算法MCMC创建动画

使用机器学习算法MCMC创建动画

使用机器学习算法MCMC创建动画

对于多模式分布,Metropolis-Hastings可能会陷入某些区域。这是一个非常常见的问题,可以通过更改提案或对更多/更长的链进行抽样来解决。

如果你想深入挖掘并查看源代码,可以访问(https://github.com/jankrepl/creating-animations-with-MCMC/blob/master/main.ipynb)

相关推荐