拒绝采样(reject sampling)原理详解

文章转自：https://blog.csdn.net/jteng/article/details/54344766
版权归原作者！

---

  蒙特·卡罗方法(Monte Carlo method)也称统计模拟方法，通过重复随机采样模拟对象的概率与统计的问题，在物理、化学、经济学和信息技术领域均具有广泛应用。拒绝采样(reject sampling)就是针对复杂问题的一种随机采样方法。
  首先举一个简单的例子介绍Monte Carlo方法的思想。假设要估计圆周率π
的值，选取一个边长为1的正方形，在正方形内作一个内切圆，那么我们可以计算得出，圆的面积与正方形面积之比为π/4
。现在在正方形内随机生成大量的点，如图1所示，落在圆形区域内的点标记为红色，在圆形区域之外的点标记为蓝色，那么圆形区域内的点的个数与所有点的个数之比，可以认为近似等于π/4
。因此，Monte Carlo方法是通过随机采样的方式，以频率估计概率。

import numpy as np
import matplotlib.pyplot as plt

size = 1e+07
z = np.random.normal(loc = 1.4,scale = 1.2, size = size)
sigma = 1.2
qz = 1/(np.sqrt(2*np.pi)*sigma**2)*np.exp(-0.5*(z-1.4)**2/sigma**2)
k = 3
u = np.random.uniform(low = 0, high = k*qz, size = size)

pz =  0.3*np.exp(-(z-0.3)**2) + 0.7* np.exp(-(z-2.)**2/0.3)
sample = z[pz >= u]
plt.hist(sample,bins=150,normed = True)

import numpy as np
import matplotlib.pyplot as plt

def f1(x):
    return 0.3*np.exp(-(x-0.3)**2) + 0.7* np.exp(-(x-2.)**2/0.3)
def f2(x):
    sigma =1.2
    return 3/(np.sqrt(2*np.pi)*sigma**2)*np.exp(-0.5*(x-1.4)**2/sigma**2)
x = np.arange(-4.,6.,0.01)

plt.plot(x,f1(x),color = "red")
plt.plot(x,f2(x),color = "blue")
plt.xticks([])
plt.yticks([])
plt.ylim(0,0.9)
plt.xlim(-4,6)
plt.plot([0.3,0.3],[0,0.54601532],color = "black")
plt.plot(0.3,0.54601532,'b.')
plt.fill_between(x,f1(x),f2(x),color = (0.7,0.7,0.7))
plt.annotate('$z_0$',xy=(0.,0),xytext=(0.2,-0.04),fontsize=15)
plt.annotate('$u_0$',xy=(0.,0.),xytext=(0.35,0.15),fontsize=15)
plt.annotate('$kq(z_0)$',xy=(0.,0.),xytext=(-0.8,0.55),fontsize=15)
plt.annotate('$p(z)$',xy=(0.,0.),xytext=(2,0.15),fontsize=15)
plt.annotate('$kq(z)$',xy=(0.,0.),xytext=(2.7,0.5),fontsize=15)

另外，参考：https://zhuanlan.zhihu.com/p/34071776