蒙特卡洛系列算法需要進行隨機變量的采樣，先理解如何生成（仿真）隨機變量（這里只討論連續(xù)型隨機變量的生成）。

逆變換法（Inverse Transform Method）

• 原理：如果我們有一個隨機變量U符合[0，1]均勻分布，并且,那么X就是一個分布函數(shù)（CDF）為F的隨機變量。也就是說如果我們需要生成符合標分布F的隨機變量X，先用均勻分布生成隨機變量U，將其輸入，得到的函數(shù)輸出結果就是我們需要的X。

一句話總結：嚴格單調遞增（保證有逆函數(shù)）的累積分布函數(shù)服從均勻分布

image.png

• 證明：見參考文獻

• 算法

  
  
  image.png

• 算法實踐：采樣符合logistic分布的隨機變量X

1. logistic分布的CDF是
   求其逆函數(shù)得
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，

對于另外一些隨機變量，比如正態(tài)分布，不可求，需要用其他方法，
主要有構成法、Box-Muller法、接受-拒絕法

Box Muller 方法 -僅針對正態(tài)分布

• 原理：如題意描述，這個聯(lián)合隨機分布是一個二元正態(tài)分布，其邊緣分布 X，Y是是兩個獨立的標準正態(tài)分布
  
  image.png
• 證明：略

• 算法：python中默認的正態(tài)分布生成方法即是這個

  
  
  image.png

其中包含了 --》如何利用逆變換法生成符合指數(shù)分布的隨機變量T：

假設 則有, ，由于, 則，所以算法中對于直接生成

-算法實踐：利用box-muller算法生成標準正態(tài)分布

# 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
• 算法實踐：利用接受拒絕法從指數(shù)分布中生成標準正態(tài)分布

求出c：

1. 的概率密度函數(shù)（PDF）為
2. 假設，目標分布為 (因為指數(shù)分布的定義中,所以先生成標準正態(tài)分布的絕對值)
3. 
4. c =
   其中sup是上確界的意思，sup和max都是定義在實數(shù)集的函數(shù)，不同點在于sup是對這個集合取上確界，max是取最大值，最大值存在時，最大值就是上確界

算法過程:

1. 使用逆變換法生成指數(shù)分布
2. 生成均勻分布
3. 代入得
   if , x = y
4. 得到 之后,令
   則

def acc_rej():
    results = []
    N = 10000
    for i in range(N):
        #先生成符號q分布即指數(shù)分布的隨機變量y
        #u = np.random.uniform(0,1)
        #y = -np.log(1-u) 利用逆變換法生成指數(shù)分布
        y = np.random.exponential(1) #直接利用numpy中的函數(shù)生成指數(shù)分布
        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數(shù)組 取出其中的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函數(shù)
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算法所需要的迭代次數(shù)將會減少，從而導致計算量減小。
  但是，如果q足夠接近p，那么采樣分布q也變得困難（之所以使用接收拒絕法正是因為難以采樣分布p）因此如何平衡容易采樣的分布q和足夠小的c是關鍵。一方面，我們想要c足夠小，接近1以減小計算量，另一方面我們又希望q容易采樣。

逆變換法 vs 接受拒絕法

• 逆變換法是在CDF上操作,需要CDF逆函數(shù)可求
• 接受拒絕法是在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/ 作者總結了生成正態(tài)分布的各種方法，包含相關證明。