采樣之蒙特卡羅方法

   隨機采樣方法

    蒙特卡洛（Monte Carlo）方法是二十世紀四十年代中期由於科學技術的發展和電子計算機的發明，而被提出的一種以概率統計理論為基礎的數值計算方法。它的核心思想就是使用隨機數（或更常見的偽隨機數）來解決一些復雜的計算問題。

   模擬方法:是一種基於“隨機數”的計算方法，基於數值采樣的近似推斷方法,也被稱為蒙特卡羅( MonteCarlo )方法、隨機模擬方法。

    通常均勻分布    的樣本，即我們熟悉的類rand()函數，可以由線性同余發生器生成；而其他的隨機分布都可以在均勻分布的基礎上，通過某種函數變換得到，比如，正態分布可以通過Box-Muller變換得到。然而，這種變換依賴於計算目標分布的積分的反函數，當目標分布的形式很復雜，或者是高維分布時，很難簡單變換得到。

當一個問題無法用分析的方法來求精確解，此時通常只能去推斷該問題的近似解，而隨機模擬(MCMC)就是求解近似解的一種強有力的方法。

     隨機模擬的核心就是對一個分布進行抽樣(Sampling)。

   Monte Carlo 方法有這些

    產生獨立樣本
        基本方法:概率積分變換(第一部分已講)
        接受—拒絕(舍選)采樣
        重要性采樣
   產生相關樣本:

      Markov Chain Monte Carlo
      Metropolis-Hastings算法
      Gibbs Sampler 

   直接采樣

    通過對均勻分布采樣，實現對任意分布采樣。

    一般而言均勻分布 Uniform(0,1)的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成[0,1]之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布 Uniform(0,1) 的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。

而我們常見的概率分布，無論是連續的還是離散的帆布，都可以基於Uniform(0,1)的樣本生成。

    直接采樣步驟：

1. 從Uniform(0,1)隨機產生一個樣本z
2. 另z=h(y)，其中h(y)為y的累積概率分布CDF
3. 計算
4. 結果y為對p(y)的采樣

    Note：需要知道累積概率分布的解析表達式，且累積概率分布函數存在反函數。但是如果h(x)不能確定或者沒有無法解析求逆則直接抽樣不再合適。對於復雜的現實模型其實是不常用的。

   接受-拒絕采樣（Acceptance-Rejection Sampling）

    簡稱拒絕采樣，基本思想：假設我們需要對一個分布f(x)進行采樣，但是卻很難直接進行采樣，所以我們想通過另外一個容易采樣的分布g(x)的樣本，用某種機制去除掉一些樣本，從而使得剩下的樣本就是來自與所求分布f(x)的樣本。

   采樣過程

 

   

    下面這張圖很好地闡釋了拒絕采樣的基本思想。假設我們想對 PDF 為 p(x)p(x) 的函數進行采樣，但是由於種種原因（例如這個函數很復雜），對其進行采樣是相對困難的。但是另外一個 PDF 為 q(x)q(x) 的函數則相對容易采樣，例如采用 Inverse CDF 方法可以很容易對對它進行采樣，甚至 q(x)q(x) 就是一個均勻分布（別忘了計算機可以直接進行采樣的分布就只有均勻分布）。那么，當我們將 q(x)q(x) 與一個常數 MM 相乘之后，可以實現下圖所示之關系，即 M⋅q(x)M⋅q(x) 將 p(x)p(x) 完全“罩住”。

    然后重復如下步驟，直到獲得 mm 個被接受的采樣點：

• 從 q(x) 中獲得一個隨機采樣點 x(i)
  
• 對於 xi 計算接受概率（acceptance probability）
• 從 Uniform(0,1) 中隨機生成一個值，用 u表示
• 如果 α≥u，則接受 xi 作為一個來自 p(x)的采樣值，否則就拒絕 xi 並回到第一步

    接受率問題

    現在有一個問題，m怎么求？

    其實也很簡單，看圖有 mq(x) >= p(x)，那么m>= p(x) / q(x)，我們求p(x) / q(x)的最大值，即為m。

   隨機采樣算法的實現

    對截斷正態分布[Truncated normal distribution]采樣

    

# -*- coding: utf-8 -*-
from scipy import optimize
from scipy.stats import norm, uniform
import numpy as np
import matplotlib.pyplot as plt

trunc = [0, 4]  # 實際分布截斷坐標點
p = [1, 1]  # 實際分布參數（均值，標准差）
q = [0, 2]  # 建議分布參數（均值，標准差）

def k(p, q, trunc):
    '''
    求k = max_x{p(x)/q(x)}
    '''
    #  loc代表了均值, scale代表標准差
    #  cdf(x, loc=0, scale=1)；輸入x，返回概率，既密度函數的面積
    #  pdf(x, loc=0, scale=1)；輸入x，返回概率密度函數

    trunc_factor = (norm.cdf(trunc[1], p[0], p[1]) - norm.cdf(trunc[0], p[0], p[1])) / p[1]
    print(trunc_factor)
    exp_k = lambda x: norm.pdf(x, p[0], p[1]) / norm.pdf(x, q[0], q[1]) / trunc_factor
    # exp_k = lambda x: norm.pdf(x, p[0], p[1]) / norm.pdf(x, q[0], q[1])#未截斷的
    max_x = optimize.minimize_scalar(lambda x: -norm.pdf(x, p[0], p[1]) / norm.pdf(x, q[0], q[1]), bounds=trunc)['x']
    k = exp_k(max_x)
    print("max_x = {}\nk = {}\n".format(max_x, k))
    return k, max_x

def show_k(k, max_x, p, q, trunc):
    '''
    求出k后繪制建議分布概率密度和實際分布概率密度圖，看p(x)和k*q(x)是否相切
    '''
    # x = np.linspace(norm.ppf(0.01, loc=p[0], scale=p[1]), norm.ppf(0.99, loc=p[0], scale=p[1]), N)
    x = np.linspace(trunc[0], trunc[1], 100)
    q = k * norm.pdf(x, loc=q[0], scale=q[1])  # 建議分布概率密度
    p = norm.pdf(x, loc=p[0], scale=p[1]) / (
        norm.cdf(trunc[1], p[0], p[1]) - norm.cdf(trunc[0], p[0], p[1]))  # 實際分布概率密度
    plt.plot(x, q, 'r')
    plt.plot(x, p, 'g')
    plt.axvline(max_x, color='b', label=max_x)  # 相切點
    plt.text(max_x, 0, str(round(max_x, 2)))
    plt.show()
def acc_rej_sample(k, p, q, trunc, N):
    '''
    接受拒絕采樣
    :param N: 采樣數
    '''
    # rvs(loc=0, scale=1, size=1, random_state=None)：產生隨機數
    z = norm.rvs(loc=q[0], scale=q[1], size=N)  # 從建議分布采樣
    mu = uniform.rvs(size=N)  # 從均勻分布采樣
    z = z[(mu <= norm.pdf(z, p[0], p[1]) / (k * norm.pdf(z, q[0], q[1])))]  # 接受-拒絕采樣
    z = z[z >= trunc[0]]
    z = z[z <= trunc[1]]
    # print("sampled z = \n{}\n".format(z))
    return z


def show_z(z, p, trunc):
    '''
    采樣得到采樣樣本z后看是否采樣得到實際正態分布的近似
    '''
    # 采樣分布概率密度圖
    cnts, bins = np.histogram(z, bins=500, normed=True)
    bins = (bins[:-1] + bins[1:]) / 2
    plt.plot(bins, cnts, label='sampling dist')
    # plt.hist(z, bins=500, normed=True)

    # 實際分布概率密度圖（截斷后的）
    x = np.linspace(trunc[0], trunc[1], 100)
    trunc_factor = (norm.cdf(trunc[1], p[0], p[1]) - norm.cdf(trunc[0], p[0], p[1])) / p[1]
    plt.plot(x, norm.pdf(x, loc=p[0], scale=p[1]) / trunc_factor, 'r', label='real dist')
    plt.legend()
    plt.show()


if __name__ == '__main__':
    k, max_x = k(p, q, trunc)
    # show_k(k, max_x, p, q, trunc)
    z = acc_rej_sample(k, p, q, trunc, N=10000000)
    show_z(z, p, trunc)

    

截斷面積trunc_factor：0.839994848037

k最大時的x點max_x = 1.333333333395553

求解出的k值：k = 2.812780139369929

 

采樣結果同真實結果的比較：

 

   參考文獻：

           https://blog.csdn.net/u010159842/article/details/78959177

           https://www.jianshu.com/p/3d30070932a8

           https://blog.csdn.net/baimafujinji/article/details/51407703