蒙特卡洛理解与python实现1-如何生成随机变量

蒙特卡洛系列算法需要进行随机变量的采样，先理解如何生成（仿真）随机变量（这里只讨论连续型随机变量的生成）。

逆变换法（Inverse Transform Method）

• 原理：如果我们有一个随机变量U符合[0，1]均匀分布，并且,那么X就是一个分布函数（CDF）为F的随机变量。也就是说如果我们需要生成符合标分布F的随机变量X，先用均匀分布生成随机变量U，将其输入，得到的函数输出结果就是我们需要的X。

一句话总结：严格单调递增（保证有逆函数）的累积分布函数服从均匀分布

image.png

• 证明：见参考文献

• 算法

  
  image.png

• 算法实践：采样符合logistic分布的随机变量X

1. logistic分布的CDF是
   求其逆函数得
2. 给定,则

#下面所有代码均引用这两个包
import numpy as np
import matplotlib.pyplot as plt  

def inverse_transform():
    results = []
    N = 1000
    for i in range(N):
        u = np.random.uniform(0,1)
        x = np.log(u / (1-u))
        results.append(x)
    return results

results_inv = inverse_transform()
plt.hist(results_inv,bins=100,color='green')
plt.show()

image.png

3. 理论验证

其中是指u的均匀分布的CDF，

对于另外一些随机变量，比如正态分布，不可求，需要用其他方法，
主要有构成法、Box-Muller法、接受-拒绝法

Box Muller 方法 -仅针对正态分布

• 原理：如题意描述，这个联合随机分布是一个二元正态分布，其边缘分布 X，Y是是两个独立的标准正态分布
  image.png
• 证明：略

• 算法：python中默认的正态分布生成方法即是这个

  
  image.png

其中包含了 --》如何利用逆变换法生成符合指数分布的随机变量T：

假设 则有, ，由于, 则，所以算法中对于直接生成

-算法实践：利用box-muller算法生成标准正态分布

# Box-Muller  
# Box-Muller  
def box_muller():
    results = []
    N = 1000
    for i in range(N):
        u = np.random.uniform(0, 1, 2) #同时生成2个符合【0，1】均匀分布的随机变量【u1，u2】
        x = np.cos(2*np.pi*u[0])*np.sqrt(-2*np.log(u[1]))
        y = np.sin(2*np.pi*u[0])*np.sqrt(-2*np.log(u[1]))
        results.append(x) #采样x、y中任意一个即可
    return results

results_boxmuller = box_muller()
plt.hist(results_boxmuller,bins=100,color='green')
plt.show()

image.png

接受拒绝法（Acceptance-Rejection Algorithm）

• 原理：如果直接生成目标p分布不容易生成，可以先生成容易得到的q分布，再利用接受拒绝法得到p。
• 证明：见参考文献

• 算法：

  
  image.png
• 算法实践：利用接受拒绝法从指数分布中生成标准正态分布

求出c：

1. 的概率密度函数（PDF）为
2. 假设，目标分布为 (因为指数分布的定义中,所以先生成标准正态分布的绝对值)
3. 
4. c =
   其中sup是上确界的意思，sup和max都是定义在实数集的函数，不同点在于sup是对这个集合取上确界，max是取最大值，最大值存在时，最大值就是上确界

算法过程:

1. 使用逆变换法生成指数分布
2. 生成均匀分布
3. 代入得
   if , x = y
4. 得到 之后,令
   则

def acc_rej():
    results = []
    N = 10000
    for i in range(N):
        #先生成符号q分布即指数分布的随机变量y
        #u = np.random.uniform(0,1)
        #y = -np.log(1-u) 利用逆变换法生成指数分布
        y = np.random.exponential(1) #直接利用numpy中的函数生成指数分布
        u = np.random.uniform(0,1)
        if u <= np.exp(-0.5*(y-1)**2):
            x = y
            z = np.random.choice([x,-x],size = 1,p=[0.5,0.5])
            results.append(z[0]) #z为numpy数组 取出其中的value
    return results
     
results_accrej = acc_rej()
plt.hist(results_accrej, bins=100, color='green')
plt.show()

• 算法实践2:利用接受拒绝法从q:中生成p:

1. Beta分布的PDF：
   其中的B为beta函数
2. 代入得, 均匀分布，则取c=1.5
3. 将c代入，得

#用均匀分布生成beta(2,2)分布x
def acc_rej():
    results = []
    N = 10000
    for i in range(N):
        y = np.random.uniform(0,1) #q分布
        u = np.random.uniform(0,1)
        if u <= (-4)*y**2+4*y:
            x = y
            results.append(x)
    return results
     
results_accrej = acc_rej()
plt.hist(results_accrej, bins=100, color='green')
plt.show()

• 接受拒绝法中隐含的tradeoff
  如果我们选择的分布q越接近p，则c越接近于1，因此生成的随机变量y越容易被接收，那么acc-rej算法所需要的迭代次数将会减少，从而导致计算量减小。
  但是，如果q足够接近p，那么采样分布q也变得困难（之所以使用接收拒绝法正是因为难以采样分布p）因此如何平衡容易采样的分布q和足够小的c是关键。一方面，我们想要c足够小，接近1以减小计算量，另一方面我们又希望q容易采样。

逆变换法 vs 接受拒绝法

• 逆变换法是在CDF上操作,需要CDF逆函数可求
• 接受拒绝法是在PDF上操作
  • 使用该方法，对于要生成的目标p分布不需要进行归一化，算法本身是自我归一化的（self-normalizing）
  • 

参考文献

1. http://www.columbia.edu/~mh2078/MonteCarlo.html 中的第一讲Generating Random Variables and StochasticProcesses。里面的课件和notes给了很多详细的证明过程，以及其他的内容如何生成多维的随机变量、离散的随机变量等。
2. https://cosx.org/2015/06/generating-normal-distr-variates/ 作者总结了生成正态分布的各种方法，包含相关证明。