5.1.2 采样

什么是采样？

采样（sampling）指的是从一个概率分布中抽取样本的过程。

1. 采出来的样本对应的是概率分布的x轴

2. 采出来很多样本之后，画histogram的图，应该近似概率分布

采样方法

一般来说，我们在python或matlab等语言中可以很轻松的借助内置函数来从一些常见分布中采样：

from scipy.stats import norm,gamma
a = norm.rvs() # 从高斯分布中采样
print(a)
b = gamma.rvs(1) # 从伽马分布中采样
print(b)

但关键问题是：假如不允许使用这些函数，我们又该如何采样呢？

接下来我们以高斯采样的例子介绍几种常见的采样方法。

问题：只给定均匀分布$[0,1]$的样本，如何利用该均匀分布从一个高斯分布中采样?

方法1：逆变化采样（Inverse Transform Sampling）

算法：

1. 从均匀分布$[0,1]$采样得到$u$

2. 从高斯分布的累积分布函数的反函数$x = F^{-1}(u)$, $x$就是需要获得的样本

几乎所有的已知分布都可以用这种方法采样（matlab和python内部就是这么做的）。

但是，如果一个复杂分布，其累积分布函数的反函数未知，如何采样?

方法2：接受-拒绝采样（Accept-Reject Sampling）

采样前的准备：

1. 先找一个提议分布（proposed distribution）$G$，概率密度函数$g(x)$已知，可以是任意分布（例如高斯分布）

2. 再找一个常数$C$，使得对于任意的$x$，都有$C \cdot g(x) \geq p(x)$

算法：

1. 从提议分布$G$中进行采样，获取采样样本$x$

2. 从$[0, 1]$的均匀分布中进行采样一个数值$u$

3. 如果$u \leq \frac{p(x)}{C \cdot g(x)}$ ,那么我们就接受这个样本$x$，保存下来；如果不满足，就拒绝并丢弃这个采样值，重来

以下是代码示例：

首先，设置目标采样分布的概率密度函数：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm,uniform
# 目标采样分布的概率密度函数(已知但是很复杂)
def p(x):
    return ((0.3 * np.exp(-(x-0.3)**2))+0.7*np.exp(-(x-2.)**2/0.3))/1.2113

选定$C$值和提议分布函数$g(x)$，并画出分布。

# 设定C值
C = 2.5

# g(x)为均值=1.4, 标准差=1.2的高斯分布

# 画出两个分布
x = np.arange(-4., 6., 0.01)
plt.plot(x, p(x), color='r', label='p(x)')
plt.plot(x, C*norm.pdf(x, loc=1.4, scale=1.2), color='b', label='C*g(x)')
plt.legend()
plt.xlabel('x')
plt.ylabel('pdf')

现在我们来采样:

sample = []
    
nSample = 10000
for i in range(nSample):
    x = np.random.normal(loc=1.4, scale=1.2) 
    u = np.random.rand()
    if u <= p(x)/(C*norm.pdf(x, loc=1.4, scale=1.2)):
        sample.append(x)

完成之后画出分布：

x = np.arange(-3., 5., 0.01)
plt.plot(x, p(x), color='r', label='p(x)')
plt.plot(x, C*norm.pdf(x, loc=1.4, scale=1.2), color='b', label='C*g(x)')
plt.hist(sample, color='k', bins=150, density=True)
plt.legend()
plt.ylabel('p(x)')

以上采样过程还可以写成矩阵形式

nSample = 10000
# norm和uniform函数都是产生一个正态或均匀分布的对象，可以搭配rvs函数等使用
uniform_rv = uniform(loc=0, scale=1) 
norm_rv = norm(loc=1.4, scale=1.2)
yy = norm_rv.rvs(size=nSample)
uu = uniform_rv.rvs(size=nSample)
mask = (uu <= p(yy)/(C*norm_rv.pdf(yy)))
sample = yy[mask] 

接受拒绝采样的直观理解：

• 在每一个$x$的点上, 都有$p(y) \leq C \cdot g(y)$，即$0 < \frac{f(y)}{C \cdot g(y)} \leq 1$

• 从图上看，就是蓝线和红线的比值。

但是，接受拒绝采样需要找到一个准确的常数$C$，如果这个常数$C$无法直接得到怎么办？

方法3：Metropolis-Hasting采样

算法：

1. 从任意一个$x_0$开始

2. 根据当前的样本$x_t$构建一个提议分布$g(x | x_t)$ （例如高斯分布），然后采样得到$x^*$

3. 然后计算接受率 $\alpha = \min \left[ 1, \frac{p(x^*) \cdot g(x_t \mid x^*)}{p(x_t) \cdot g(x^* \mid x_t)} \right]$

4. 从$[0,1]$采样一个值$u$,

   - 如果 $u< \alpha$, 接受该样本 $x_{t+1} = x^*$

   - 否则，继续维持原来的样本 $x_{t+1} = x_t$

5. 重复2-4步直到达到设定的样本总数

在MCMC的证明中会详细介绍该算法并解释它为什么成立，在这里，我们先来看看代码如何实现。

首先，设置目标分布的概率密度函数：

from scipy.stats import norm
import numpy as np 

# 目标采样分布的概率密度函数
def p(x):
    return ((0.3 * np.exp(-(x-0.3)**2))+0.7*np.exp(-(x-2.)**2/0.3))/1.2113
    
nSample = 100000
sample2 = []

y = 1 # 设置初始值

#进行MH采样
for i in range(nSample):
    
    y2 = norm(loc=y).rvs() # 从提议分布采样一个
    
    a = np.min([1, p(y2)/p(y)]) # 因为提议分布高斯是对称的, 所以可以上下约掉

    u = np.random.rand() # 从[0,1]均匀分布获取一个值
    
    if u <= a:
        sample2.append(y2) # 接受并跳转到新的sample
        y = y2
    else:
        sample2.append(y) # 保留原来那个sample不变
        
# 绘制采样结果
x = np.arange(-3., 5., 0.01)
plt.plot(x, p(x), color='r', label='p(x)')
a,_,_=plt.hist(sample2, color='k', bins=150, density=True)
plt.legend()
plt.xlabel('X')
plt.ylabel('p(x)')