拒绝采样总结：rand7生成rand10、单位圆内部均匀采样、Box-Muller、MCMC采样法、吉布斯采样Gibbs

几乎所有采样方法都是以某个均匀分布为基础，经过变换+拒绝采样得到新的采样方式。拒绝采样的一大痛点在于拒绝测次数太多效率低，实际中需结合概率密度函数pdf和累积概率密度函数cdf来设计期待包络的缝隙可以尽量小。cdf是一个单调递增的函数，实际上中存在着以下关系，主要是解决部分pdf不好求的问题
$$cd{f}^{-1}(x)=pdf(x)$$
一般来说，如果随机变量u和随机变量x存在某种变换关系，那么概率密度函数p(u)和p(x)存在以下关系，其中J是Jacobian行列式，一会儿在单位圆内部均匀采样的例子中可以看到效果
$$p(u)\ast J=p(x)$$

rand7生成rand10

这个问题的主要点在于被拒绝后怎么优化次数，以下是一个比较好的思路，刚好是Leetcode 470题

class Solution:
    def rand10(self):
        row,col = rand7(),rand7()
        res1 = (row-1)*7+col
        if (res1 <= 40):
            return res1 % 10 + 1
        res2 = (res1-41)*7+rand7() #(大于40的随机数-40-1)*7+rand7()
        if (res2 <= 60):
            return res2 % 10 + 1
        res3 = (res2-61)*7+rand7()
        while (res3 == 21):
            res3 = (res2-61)*7+rand7() #(大于60的随机数-60-1) *7+rand7()
        return res3 % 10 + 1

在单位圆内部均匀采样

错误做法

我们假设有两个从0到1的随机数a和b，
$$\begin{gathered} r=a \\ theta=b\ast 2\ast Pi \end{gathered}$$
上述随机点显然可以布满整个圆，但是并不是均匀分布的。简单证明如下，如果上述算法可以均匀生成随机点，那么我们假设r0为0到1之间的一个固定数值。由于r=a，而theta不影响点到圆心的距离，所以点在以r0为半径的圆内部的概率为P(a < r0) = r0。而从面积角度讲，点在以r0为半径的圆内部的概率应该是该面积与单位圆的面积比例，即r0r0。由于r0在0和1之间，显然r0!=r0r0。所以我们的假设是不成立的，即上述算法不是均匀分布的。

极坐标思路是对的但局部性不好

还是假设有两个从0到1的随机数a和b，注意这下给a开根号了
$$\begin{gathered} r=\sqrt{a} \\ theta=b\ast 2\ast Pi \end{gathered}$$
方法正确性可以从周长和面积两个角度去解释：

• 周长： 考虑单位圆内部的每一个圆环，生成的点落在半径为r的圆环上的概率应当与圆环的周长成正比，这点没有问题。因为$pdf(r)$在定义域(0,1)内积分为1，而${\int }_0^{1}2rdr=1$，所以$pdf(r)=2r$
• 面积：$cdf(r)$刚好是面积，${\int }_0^{r}pdf(r)dr=r^2$，这回和平方成正比，所以$r^2=a$，a是(0,1)之间均匀分布，那么$r=\sqrt{a}$

代码如下

import matplotlib.pyplot as plt
import random
from math import sqrt
from typing import List
import math
import matplotlib.pyplot as plt

def randpoints1(xc,yc,radius) -> List[float]:
    u, theta = random.random(), random.random() * 2 * math.pi
    r = sqrt(u)
    return [xc + r * math.cos(theta) * radius, yc + r * math.sin(theta) * radius]

def plt1():
    x_axis, y_axis = [], []
    for i in range(1000):
        x, y = randpoints1(0, 0, 1)
        x_axis.append(x)
        y_axis.append(y)

    plt.plot(x_axis, y_axis, '+')
    plt.show()


if __name__ == '__main__':
    plt1()

局部性也好的Concentric Map方案

思路源自Peter Shirley很早就有一篇文章A Low Distortion Map Between Disk and Square，左边是原始随机数(a,b)生成的坐标系（我们假设有这样的坐标系），右边是经过了上述变换后的均匀采样，我们发现图片中的F已经完全走样了，基本上看不出来了

但下面变化就避免了这个问题
$$\begin{gathered} r=a \\ theta=Pi\ast b/(4\ast a) \end{gathered}$$
通过上诉变换，我们可以把正方形压缩成圆的形状

刚才F的局部性也好很多

算法正确性的证明用到了最开始提到的Jacobian行列式
我们只需要找到p(a,b)与p(r,theta)之间的关系就好了，而这两者的关系是不难计算的：
p(a,b) = 0.25（因为a和b都是-1到1的均匀分布,p(a,b)=p(a)×p(b)=1/2×1/2=1/4）
p(u,v) = p(a,b) / |J_T|
其中
u = a * cos( Pi * b / ( 4 * a ) )
v = a * sin( Pi * b / ( 4 * a ) )
J_T为Jacobian行列式，通过计算我们得到
$$\left|\begin{array}{cc}\frac{\partial u}{\partial a}& \frac{\partial u}{\partial b}\\ \frac{\partial v}{\partial a}& \frac{\partial v}{\partial b}\end{array}\right|=\left|\begin{array}{cc}\cos (\frac{\pi b}{4a})+\frac{\pi b}{4a}\sin (\frac{\pi b}{4a})& -\frac{\pi }{4}\sin (\frac{\pi b}{4a})\\ \sin (\frac{\pi b}{4a})-\frac{\pi b}{4a}\cos (\frac{\pi b}{4a})& \frac{\pi }{4}\cos (\frac{\pi b}{4a})\end{array}\right|$$
所以p(u,v) = 0.25 / ( Pi / 4 ) = 1 / Pi，即在圆内部每一点被采样到的概率都为1/Pi，所以上述变换为均匀分布。
注意，Concentric Map方案需要转一下角度，也就是说只有上述推导1,3符合要求，实现可以参考http://psgraphics.blogspot.com/2011/01/improved-code-for-concentric-map.html和https://github.com/erich666/jgt-code/blob/master/Volume_02/Number_3/Shirley1997/disk.cc：

void lowd_square_to_disk(float x_in, float y_in, float *x_out, float *y_out) {
    float theta,r;
    float a = 2*x_in - 1;
    float b = 2*y_in - 1;
    if (a == 0 && b == 0) {
        r = theta = 0;
    }
    else if (a*a> b*b) { 
        r = a;
        theta = (M_PI/4)*(b/a);
    } else {
        r = b;
        theta = (M_PI/2) - (M_PI/4)*(a/b);
    }
    *x_out = r*cos(theta);
    *y_out = r*sin(theta);
}

用Python实现：

class Solution:

    def __init__(self, radius: float, x_center: float, y_center: float):
        self.radius = radius
        self.xc = x_center
        self.yc = y_center

    def randPoint(self) -> List[float]:
        a,b = random.random()*2-1,random.random()*2-1
        if (a == 0.0 or b == 0.0):
            return [0.0, 0.0]

        if (a*a < b*b):
            pi_4 = math.pi/4*a/b
            return [self.xc+math.cos(pi_4)*self.radius*b, self.yc+math.sin(pi_4)*self.radius*b]
        else:
            pi_4 = math.pi/2-math.pi/4*b/a
            return [self.xc+math.cos(pi_4)*self.radius*a, self.yc+math.sin(pi_4)*self.radius*a]

笔者自己构思不用转角度的思路

假设a,b是0,1之间均匀分布的随机数，那么
$$\begin{gathered} u=max(a,b) \\ v=\cos (2\pi \ast min(a,b)/max(a,b)) \end{gathered}$$
a>b和b>a完美对称，我们先考虑a>b的情况，那么就是
$$\begin{gathered} u=a\ast \cos (2\pi \ast b/a) \\ v=a\ast \sin (2\pi \ast b/a) \end{gathered}$$
同样利用雅各比行列式：
$$\left|\begin{array}{cc}\frac{\partial u}{\partial a}& \frac{\partial u}{\partial b}\\ \frac{\partial v}{\partial a}& \frac{\partial v}{\partial b}\end{array}\right|=\left|\begin{array}{cc}\cos (\frac{2\pi b}{a})+\frac{2\pi b}{a}\sin (\frac{2\pi b}{a})& -2\pi \sin (\frac{2\pi b}{a})\\ \sin (\frac{2\pi b}{a})-\frac{2\pi b}{a}\cos (\frac{2\pi b}{a})& 2\pi \cos (\frac{2\pi b}{a})\end{array}\right|=2\pi$$
因为p(a,b)=p(a)*p(b)=1(a,b是0,1之间均匀分布的随机数)，所以p(u,v) = 1 / (2Pi)，即在圆内部一半点被采样到的概率都为1 / (2pi)，对称下另外一半也是1/(2pi)

代码实现如下：

import matplotlib.pyplot as plt
import random
from math import sqrt
from typing import List
import math
import matplotlib.pyplot as plt

def randpoints1(xc,yc,radius) -> List[float]:
    u, theta = random.random(), random.random() * 2 * math.pi
    r = sqrt(u)
    return [xc + r * math.cos(theta) * radius, yc + r * math.sin(theta) * radius]

def randpoints2(xc,yc,radius) -> List[float]:
    #r=a theta=Pi∗b/(4∗a)
    a,b = random.random(),random.random()
    max_ab, min_ab = max(a,b),min(a,b)
    random_2pi = (math.pi*2)*(min_ab/max_ab)
    theta = random_2pi
    delta_x = max_ab * math.cos(theta)*radius
    detta_y = max_ab * math.sin(theta)*radius
    return [xc+delta_x,yc+detta_y]


def plt1():
    x_axis, y_axis = [], []
    for i in range(1000):
        x, y = randpoints1(0, 0, 1)
        x_axis.append(x)
        y_axis.append(y)

    plt.plot(x_axis, y_axis, '+')
    plt.show()

def plt2():
    x_axis, y_axis = [], []
    for i in range(30000):
        x, y = randpoints2(0, 0, 1)
        x_axis.append(x)
        y_axis.append(y)

    plt.plot(x_axis, y_axis, '+')
    plt.show()

if __name__ == '__main__':
    plt2()

高斯分布的采样（当然是说截断到某个范围内）

如果范围是负无穷到正无穷，高斯分布的概率密度函数是高中学过的那个。但是如果截断到某个范围内，概率密度函数是$er{f}^{-1}(x)$不是初等函数：
$$erf(x)=\frac{2}{\sqrt{\pi }}{\int }_0^t{e}^{-t^2}dt$$
所以弄俩独立变量x,y再求上面的，转换到极坐标，这种方式叫做Box-Muller方法
但是Box-Muller方法本身要有三角函数，也不是最快的，效率更高的Ziggurat方法

Box-Muller方法

以下转载自 https://eipi10.cn/statistics/2020/07/17/box-muller/

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

def box_muller(SampleSize=1000000):
    U1 = np.random.uniform(0,1,SampleSize)
    U2 = np.random.uniform(0,1,SampleSize)
    X = np.cos(2*np.pi*U1)*np.sqrt(-2*np.log(U2))
    return X

X = box_muller()
plt.hist(X, bins=np.linspace(-4,4,81),facecolor="blue")
plt.title('X')
 
# 检验分布是否符合正态分布
x_test = stats.kstest(X, 'norm')
# 大多数情况下，p值>>0.05，证明不能拒绝X是正态分布。
print(x_test.pvalue)

ziggurat算法

以下转载自 https://www.zhihu.com/question/29971598，下图中提到的超越方程即解lnx = x这种非初等函数的方程。ziggurat算法即划分成等面积的块，然后拒绝采样

马尔科夫蒙特卡洛采样法（MCMC）

构造一个马尔可夫链,使得该马尔可夫链的平稳分布就是目标分布;然后,从任何一个初始状态出发,沿着马尔可夫链进行状态转移,最终得到的状态转移序列会收敛到目标分布,由此可以得到目标分布的一系列样本

非周期的马尔科夫链模型的状态转移矩阵收敛到的稳定概率分布与我们的初始状态概率分布无关

先看什么是周期的马尔科夫链：

非周期的马尔科夫链数学表示如下：

马尔科夫链的细致平稳条件

MCMC采样公式

解释一下上面公式，$x_t$表示第t轮采样得到的样本，根据$Q(x|x_t)$采样出样本$x_{\ast }$，$\pi (x_{\ast })$是根据要采样的那个平稳分布的$\pi (x)$的pdf算出来的概率，而$Q(x_{\ast },x_t)$是根据$Q(x|x_t)$对应的概率密度pdf算出来的概率，这俩乘了一下，就是那个$\alpha (x_t,x_{\ast })$。

Metropolis-Hastings采样

Metropolis-Hastings采样法是MCMC采样法的一个改进，防止MCMC中被拒绝的概率过高。而MH的核心在于需要目前分布的pdf也就是下图里的$\pi (x)$和一个状态转移矩阵Q（这个Q的转移是当下采样值$x_t$转转移到$x_{\ast }$，实际中经常用正态分布来采样）

结合code来看一下：

import random
import math
from scipy.stats import norm
import matplotlib.pyplot as plt
# matplotlib inline
def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y
T = 5000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t<T-1:
    t = t + 1
    # 把Q(x,x_t)定义为均值为pi[t-1]，方差为1的正态分布，从而根据pi[t-1]采样出pi_star
    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)
    # norm_dist_prob是根据采样得到的值和pdf曲线算一个概率
    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))
    # 这个采样的公式就是图里的简化版pi(j)/pi(i)，其中i就是之前轮采样的值，j是根据概率转移Q得到的下一轮的值
    u = random.uniform(0, 1)
    if u < alpha:
        pi[t] = pi_star[0]
    else:
        pi[t] = pi[t - 1]

plt.scatter(pi, norm.pdf(pi, loc=3, scale=2))
num_bins = 50
plt.hist(pi, num_bins, facecolor='red', alpha=0.7)
plt.show()

再来看一下完整的包含Q的code，转载自 https://github.com/abdulfatir/sampling-methods-numpy/blob/master/Metropolis-Hastings.ipynb

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm as tqdm

# 希望采样P这个分布
P = lambda x: 3 * np.exp(-x*x/2) + np.exp(-(x - 4)**2/2)
Z = 10.0261955464

x_vals = np.linspace(-10, 10, 1000)
y_vals = P(x_vals)

f_x = lambda x: x
g_x = lambda x: np.sin(x)
true_expected_fx = 10.02686647165
true_expected_gx = -1.15088010640

# 概率转移定义成正态分布
Q = lambda t1, t0: np.exp(-(t1 - t0)**2/2)/(np.sqrt(2 * np.pi))
x0 = 0
xt = x0
samples = []
for i in range(100000):
    xt_candidate = np.random.normal(xt, 1)
    accept_prob = (P(xt_candidate) * Q(xt, xt_candidate))/(P(xt) * Q(xt_candidate, xt))
    if np.random.uniform(0, 1) < accept_prob:
        xt = xt_candidate
    samples.append(xt)
# burn_in表示前burn_in个样本不稳定
burn_in = 1000
samples = np.array(samples[burn_in:])
expected_f_x = np.mean(f_x(samples))
expected_g_x = np.mean(g_x(samples))
expected_f_x *= Z
expected_g_x *= Z
print('E[f(x)] = %.5f, Error = %.5f' % (expected_f_x, abs(expected_f_x - true_expected_fx)))
print('E[g(x)] = %.5f, Error = %.5f' % (expected_g_x, abs(expected_g_x - true_expected_gx)))

plt.hist(samples, bins=50, histtype='bar', facecolor='g', alpha=0.75, label='bins')
plt.plot(x_vals, y_vals/Z, 'r', label='P(x)')
plt.title('Metropolis Hastings')
plt.legend(loc='upper right', shadow=True)
plt.show()

完整转载

来自https://zhuanlan.zhihu.com/p/30003899

吉布斯采样法

吉布斯采样法是MH采样法的一个特例，特点是只对样本的一个维度进行采样和更新
来自https://zhuanlan.zhihu.com/p/30003899的解释比较通俗：

下面的推导都转载自https://kexue.fm/archives/8084：

转载鸣谢：

1. https://blog.youkuaiyun.com/codeboycjy/article/details/6225886
2. https://leetcode.cn/problems/generate-random-point-in-a-circle/solution/zai-yuan-nei-sui-ji-sheng-cheng-dian-by-qp342/
3. 葫芦书