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

初步研究

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

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

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

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

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

Rejection sampling（拒绝采样）

第一种方法属于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并没有被诅咒，因此拒绝采样是一个很好的选择。

吉布斯采样Gibbs Sampling

吉布斯抽样属于通过构造马尔可夫链生成样本的第二种抽样类型。因此，这些样本不是独立的。事实上，直到链达到其平稳分布，它们才会均匀分布。出于这个原因，通常的做法是丢弃前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采样

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

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

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

动画

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

• burn-in samples: 500
• samples: 10000

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

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