5.1.2 采样

什么是采样？

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

采出来的样本对应的是概率分布的x轴
采出来很多样本之后，画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）

算法：

从均匀分布 $[0,1]$ 采样得到 $u$
从高斯分布的累积分布函数的反函数 $x = F^{-1}(u)$ , $x$ 就是需要获得的样本

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

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

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

采样前的准备：

先找一个提议分布（proposed distribution） $G$ ，概率密度函数 $g(x)$ 已知，可以是任意分布（例如高斯分布）
再找一个常数 $C$ ，使得对于任意的 $x$ ，都有 $C \cdot g(x) \geq p(x)$

算法：

从提议分布 $G$ 中进行采样，获取采样样本 $x$
从 $[0, 1]$ 的均匀分布中进行采样一个数值 $u$
如果 $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采样

算法：

从任意一个 $x_0$ 开始
根据当前的样本 $x_t$ 构建一个提议分布 $g(x | x_t)$ （例如高斯分布），然后采样得到 $x^*$
然后计算接受率 $\alpha = \min \left[ 1, \frac{p(x^*) \cdot g(x_t \mid x^*)}{p(x_t) \cdot g(x^* \mid x_t)} \right]$
从 $[0,1]$ 采样一个值 $u$ ,
- 如果 $u< \alpha$ , 接受该样本 $x_{t+1} = x^*$
- 否则，继续维持原来的样本 $x_{t+1} = x_t$
重复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)')

最后更新于1天前