随机游动Metropolis算法拟合标准拉普拉斯分布

1. 题目分析

使用随机游动的Metropolis抽样方法产生标准拉普拉斯分布的随机数，其中使用正态分布产生增量，比较由提议分布不同方差生成的链的差异，并比较每个链的接收概率。

2. 代码展示

我们首先先导入所需要的库：

import numpy as np
import scipy.stats as st
import seaborn as sns
import matplotlib.pyplot as plt

按照给定目标分布f，提议分布g，按照随机游动的Metropolis算法我们编写代码：

def RandomWalkMetropolis(m, f, g, x0):
    k = 0
    np.random.rand(0)
    x = np.zeros(m)
    u = np.random.random(m)
    x[0] = x0

    for i in range(1, m):
        xt = x[i-1]
        y = g.rvs()
        num = f.pdf(y)
        den = f.pdf(xt)
        if u[i] <= num / den:
            x[i] = y
            k += 1
        else:
            x[i] = xt
    print('acceptance rate: ' + str(k / m))

    return  x

我们设定马氏链的长度为100000；预烧期为10000；目标分布为标准拉普拉斯分布；提议分布为均值为0、标准差为$\sigma$的正态分布，其中$\sigma \in [0.05,0.5,1,4,16,32]$，初始化点$x_0=0$：

name = '6_6_sigm0.05'
m = 100000
burn = 10000
sigm = [0.05, 0.5, 1, 4, 16, 32]
f = st.laplace
g = st.norm(scale=sigm[0])
x = RandomWalkMetropolis(m, f, g, 0)

3. 模型评价

在此部分，我们系统的比较了提议分布不同标准差$\sigma$条件下接受概率、Trace Plot图、QQ图和采样分布-目标分布分布密度图来评价我们的算法模型。

3.1 接受概率

我们采用不同的标准差重复运行2中Metropolis算法得到以下结果：

sigma=0.05, acceptance rate: 0.98319
sigma=0.5, acceptance rate: 0.84878
sigma=1, acceptance rate: 0.72404
sigma=4, acceptance rate: 0.34838
sigma=16, acceptance rate: 0.09858
sigma=32, acceptance rate: 0.04895

我们可以看到，随着标准差的增大，接受概率逐渐下降，采样效率明显降低。

3.2 轨迹图

我们以迭代次数为横轴，以马尔科夫链的值为纵轴，编写绘制TracePlot的函数

def DrawTracePlot(x, name):
    dataX = np.arange(1, x.shape[0]+1)
    dataY = x

    plt.figure(figsize=(12, 5))
    plt.plot(dataX, dataY, linewidth=0.4, color="red")
    plt.title("Trace Plot", fontsize=20)
    plt.xlabel("Iterations", fontsize=12)
    plt.ylabel("Sample Values", fontsize=12)

    plt.savefig(name + "TracePlot.png", dpi=200)
    plt.show()

我们调用上述代码，得到TracePlot图。

从中我们看出，随着提议分布标准差的增大，采样的效率大大降低。

3.3 分布直方图

接下来，我们绘制抽样的频率直方图，观察分布的情况。我们编写如下代码：

def DrawHistogram(x, burn, f, name):
    x = x[burn:]
    plt.figure(figsize=(10,5))
    plt.title("Histogram")
    plt.xlabel('Observed Values')
    plt.ylabel('Distribution Density')
    plt.hist(x, bins=100, color="blue", density=True, edgecolor="black", alpha=0.5)
    y = f.rvs(size=x.shape[0])

    plt.savefig(name + " Histogram.png", dpi=200)
    plt.show()

我们传入参数，运行上述代码，得到如下结果：






我们通过肉眼观察，认为抽样分布基本符合标准拉普拉斯分布。

3.4 分布密度图

更进一步的，我们画出理论分布与抽样分布的分布密度图，进行比较。我们编写以下代码：

def DrawDistribution(x, burn, f):
    x = x[burn:]
    tick = np.arange(x.min(), x.max(), 0.001)
    plt.figure(figsize=(10,5))
    plt.plot(tick, f.pdf(tick), color="red")
    p, y = np.histogram(x, bins=100, normed=True)
    y = (y[:-1] + y[1:])/2
    plt.plot(y, p, '--b')
    plt.title('Distribution')
    plt.legend(['target distribution', 'fitted distribution'])
    plt.xlabel('Observed Values')
    plt.ylabel('Distribution Density')
    plt.savefig(name + "Distribution.png", dpi=200)
    plt.show()

我们传入参数，运行上述代码，得到如下结果：

通过观察，我们可以发现随着方差的逐渐增大样本分布密度与理论分布密度呈现逐渐接近的趋势，这一点反映了提议分布的支撑集必须包含目标分布的支撑集这一前提。

3.5 QQ图

最后为了系统性的量化样本分布与理论分布的差异，我们画出了QQ图，我们编写以下代码：

def DrawQQPlot(x, burn, scale):
    x = x[burn:]
    x.sort()
    y = st.cauchy.rvs(size=x.shape[0])
    y.sort()
    quantileX = [np.percentile(x, i) for i in range(101)]
    quantileY = [np.percentile(y, i) for i in range(101)]

    print('10 quantitile of sample:{}'.format(quantileX[10]))
    print('10 quantitile of theoretical distribution:{}'.format(quantileY[10]))

    plt.scatter(quantileY, quantileX, c='blue')
    plt.xlim(scale)
    plt.ylim(scale)
    plt.xlabel('Theoretical Quantitiles')
    plt.ylabel('Sample Quantitiles')
    plt.plot(scale, scale, linewidth=2, color='red')
    plt.savefig(name + "QQPlot.png", dpi=200)
    plt.show()

我们传入参数，运行上述代码，得到如下结果：






我们可以观察到，随着提议分布标准差的增大，抽样分布和原分布越来越接近。

至此，我们的任务就完成了。