[PBRT-V3]从MCMC的角度分析MLT

关于MLT，之前有过一篇博文：Q139：PBRT-V3，Metropolis Light Transport (MLT)（16.4章节）

为了更好地理解MLT，小编决定从MCMC的角度对该渲染算法分析一下。

关于MCMC，主要参考如下（感谢Eureka和刘建平Pinard写了这些文章）。
Eureka的博文：
马尔可夫链蒙特卡罗算法（MCMC）
刘建平Pinard的博文：
MCMC(一)蒙特卡罗方法
MCMC(二)马尔科夫链
MCMC(三)MCMC采样和M-H采样

关于MCMC的解析，请参考如上链接。
后文主要包含两部分：小编的MCMC阅读笔记、从MCMC的角度分析MLT

MCMC

Metropolis采样的结果真的满足“重要性采样”吗？？？

（实例内容是基于：https://zhuanlan.zhihu.com/p/37121528）
实例：假设目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵 Q(i,j) 的条件转移概率是以 i 为均值，标准差1的正态分布在位置 j 的值。

python实现代码：

#-*- coding: UTF-8 -*-
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt

def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y

T = 400
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t < T-1:
    t = t + 1
    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)   #状态转移进行随机抽样
    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))   #alpha值

    u = np.random.uniform(0, 1)
    if u < alpha:
        pi[t] = pi_star[0]
    else:
        pi[t] = pi[t - 1]


plt.xlim(-5, 11)
plt.ylim(-0.01, 0.3)

plt.scatter(pi, norm.pdf(pi, loc=3, scale=2),label='Target Distribution')
num_bins = 50
plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7,label='Samples Distribution')

plt.legend()
plt.show()

采样结果如下：

为了更好地体现出采样点的分布情况，引入“直方图”

注意：
1，不要将直方图和条形图搞混。可以参考一下：【Python-matplotlib】画直方图（hist）
2，不要将直方图和“细分，累加求积分”混淆。
每一个矩形的高度表示矩形宽度区间内采样点出现的频数或者频率
从上图可以看到：[1,2]、[3,5]区间内的采样点分布相对密集（矩形相对较高）。
采样点的分布情况貌似和目标分布的函数值并不是严格匹配（重要性采样的“重要性”体现得不够好），原因是：采样点个数太少。

下面是采样点个数分别是1600，6400和25600的情况。

从这些图可以看出：随着采样点数目的增加，采样点的分布情况和目标分布越来越接近（也就是，越来越体现“重要性采样”）。
再次提醒：不要将直方图和“细分，累加求积分”混淆。
每一个矩形的高度表示矩形宽度区间内采样点出现的频数或者频率

考虑到正态分布原本是没有必要使用Metropolis方法进行采样的，而是可以直接进行重要性采样。python中有现成的采样函数。
这里咱对比看一下“Metropolis采样”和“直接重要性采样”的效果。
直接重要性采样的python代码如下。

#-*- coding: UTF-8 -*-
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt

pi_star = norm.rvs(loc=3, scale=2, size=25600, random_state=None)

plt.xlim(-5, 11)
plt.ylim(-0.01, 0.3)

plt.scatter(pi_star, norm.pdf(pi_star, loc=3, scale=2),label='Target Distribution')
num_bins = 50
plt.hist(pi_star, num_bins, normed=1, facecolor='red', alpha=0.7,label='Samples Distribution')

plt.legend()
plt.show()

后面是采样点个数分别是400，1600，6400和25600时“直接重要性采样”的情况。




对比Metropolis采样和直接重要性采样的结果，可以看出：
当采样点个数不是足够大时，直接重要性采样的结果的结果要比Metropolis采样的结果要好。
所以，如果能够进行直接重要性采样，就直接采样吧。
这里顺便提一下，通过“反函数法”进行直接重要性采样的步骤如下：

MLT具体是怎么实现路径空间的Metropolis采样的呢？