MH算法拟合标准柯西分布

1. 题目

利用M-H算法从标准柯西分布中产生随机数，丢弃链的前1000个值，比较生成链观测值的十分位数和柯西分布理论10分位数的拟合情况，并画出QQ图和链的直方图。本题中提议分布取为$N(x_t,{\sigma }^2)$。

2. 代码展示

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

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

按照给定目标分布f，提议分布g，按照MH算法我们编写代码：

def Metropolis(m, f, g, sigm):
    k = 0
    x = np.zeros(m)
    u = np.random.rand(m)
    x[0] = g.rvs(0, scale=sigm)

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

    return  x

我们设定马氏链的长度为10000；预烧期为1000；目标分布为标准柯西分布；提议分布为均值为0、方差为2的正态分布：

name = '6_3'
m = 100000
burn = 1000
sigm = [0.05, 0.5, 1, 2, 16]
f = st.cauchy
g = st.norm
x = Metropolis(m, f, g, sigm[2])

运行代码，我们得到如下结果：

reject probablity: 0.22952

拒绝概率落入$10\%-40\%$的区间内，符合马氏链的要求。

3. 模型评价

在此部分，我们可视化了算法过程的Trace Plot图、QQ图和采样分布和目标分布的分布密度图来评价我们的模型。

3.1 轨迹图

我们以迭代次数为横轴，以马尔科夫链的值为纵轴，编写绘制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("x", fontsize=12)

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

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

从中我们看出，马尔科夫链在迭代100000轮后无明显变化趋势，已基本收敛。

3.2 分布直方图

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

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.3 分布密度图

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

def DrawDistribution(x, burn, f):
    x = x[burn:]
    tick = np.arange(x.min(), x.max(), 0.001)
    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.4 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()

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

10 quantitile of sample distribution:-2.5831631465936775
10 quantitile of theoretical distribution:-3.1107905666662963

可视化结果如下：

我们可以观察到，在靠近中心的$(0,0)$附近，两分布拟合状况较好，分布的拟合状况像两侧逐渐偏移。整体来看，量分布拟合状况良好。