利用python進行蒙特卡羅模擬

本文用Python統計模擬的方法，介紹四種常用的統計分布，包括離散分布：二項分布和泊松分布，以及連續分布：指數分布和正態分布，最後查看人群的身高和體重數據所符合的分布。

# 導入相關模塊import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport seaborn as sns %matplotlib inline %config InlineBackend.figure_format = 'retina'

隨機數

計算機發明後，便產生了一種全新的解決問題的方式：使用計算機對現實世界進行統計模擬。該方法又稱為「蒙特卡洛方法（Monte Carlo method）」，起源於二戰時美國研製原子彈的曼哈頓計劃，它的發明人中就有大名鼎鼎的馮·諾依曼。蒙特卡洛方法的名字來源也頗為有趣，相傳另一位發明者烏拉姆的叔叔經常在摩洛哥的蒙特卡洛賭場輸錢，賭博是一場概率的遊戲，故而以概率為基礎的統計模擬方法就以這一賭城命名了。

使用統計模擬，首先要產生隨機數，在Python中，numpy.random 模塊提供了豐富的隨機數生成函數。比如生成0到1之間的任意隨機數：

np.random.random(size=5)  # size表示生成隨機數的個數

array([ 0.32392203,  0.3373342 ,  0.51677112,  0.28451491,  0.07627541])

又比如生成一定範圍內的隨機整數：

np.random.randint(1, 10, size=5)  # 生成5個1到9之間的隨機整數

array([5, 6, 9, 1, 7])

計算機生成的隨機數其實是偽隨機數，是由一定的方法計算出來的，因此我們可以按下面方法指定隨機數生成的種子，這樣的好處是以後重複計算時，能保證得到相同的模擬結果。

np.random.seed(123)

在NumPy中，不僅可以生成上述簡單的隨機數，還可以按照一定的統計分布生成相應的隨機數。這裡列舉了二項分布、泊松分布、指數分布和正態分布各自對應的隨機數生成函數，接下來我們分別研究這四種類型的統計分布。

np.random.binomial()

np.random.poisson()

np.random.exponential()

np.random.normal()

二項分布

二項分布是n個獨立的是/非試驗中成功的次數的概率分布，其中每次試驗的成功概率為p。這是一個離散分布，所以使用概率質量函數（PMF）來表示k次成功的概率：

最常見的二項分布就是投硬幣問題了，投n次硬幣，正面朝上次數就滿足該分布。下面我們使用計算機模擬的方法，產生10000個符合（n，p）的二項分布隨機數，相當於進行10000次實驗，每次實驗投擲了n枚硬幣，正面朝上的硬幣數就是所產生的隨機數。同時使用直方圖函數繪製出二項分布的PMF圖。

def plot_binomial(n,p):    '''繪製二項分布的概率質量函數'''    sample = np.random.binomial(n,p,size=10000)  # 產生10000個符合二項分布的隨機數    bins = np.arange(n+2)    plt.hist(sample, bins=bins, align='left', normed=True, rwidth=0.1)  # 繪製直方圖    #設置標題和坐標    plt.title('Binomial PMF with n={}, p={}'.format(n,p))      plt.xlabel('number of successes')    plt.ylabel('probability') plot_binomial(10, 0.5)

投10枚硬幣，如果正面或反面朝上的概率相同，即p=0.5， 那麼出現正面次數的分布符合上圖所示的二項分布。該分布左右對稱，最有可能的情況是正面出現5次。

但如果這是一枚作假的硬幣呢？比如正面朝上的概率p=0.2，或者是p=0.8，又會怎樣呢？我們依然可以做出該情況下的PMF圖。

fig = plt.figure(figsize=(12,4.5)) #設置畫布大小p1 = fig.add_subplot(121)  # 添加第一個子圖plot_binomial(10, 0.2) p2 = fig.add_subplot(122)  # 添加第二個子圖plot_binomial(10, 0.8)

這時的分布不再對稱了，正如我們所料，當概率p=0.2時，正面最有可能出現2次；而當p=0.8時，正面最有可能出現8次。

泊松分布

泊松分布用於描述單位時間內隨機事件發生次數的概率分布，它也是離散分布，其概率質量函數為：

比如你在等公交車，假設這些公交車的到來是獨立且隨機的（當然這不是現實），前後車之間沒有關係，那麼在1小時中到來的公交車數量就符合泊松分布。同樣使用統計模擬的方法繪製該泊松分布，這裡假設每小時平均來6輛車（即上述公式中lambda=6）。

lamb = 6sample = np.random.poisson(lamb, size=10000)  # 生成10000個符合泊松分布的隨機數bins = np.arange(20) plt.hist(sample, bins=bins, align='left', rwidth=0.1, normed=True) # 繪製直方圖# 設置標題和坐標軸plt.title('Poisson PMF (lambda=6)') plt.xlabel('number of arrivals') plt.ylabel('probability') plt.show()

指數分布

指數分布用以描述獨立隨機事件發生的時間間隔，這是一個連續分布，所以用質量密度函數表示：

比如上面等公交車的例子，兩輛車到來的時間間隔，就符合指數分布。假設平均間隔為10分鐘（即1/lambda=10)，那麼從上次發車開始，你等車的時間就滿足下圖所示的指數分布。

tau = 10sample = np.random.exponential(tau, size=10000)  # 產生10000個滿足指數分布的隨機數plt.hist(sample, bins=80, alpha=0.7, normed=True) #繪製直方圖plt.margins(0.02) # 根據公式繪製指數分布的概率密度函數lam = 1 / tau x = np.arange(0,80,0.1) y = lam * np.exp(- lam * x) plt.plot(x,y,color='orange', lw=3)#設置標題和坐標軸plt.title('Exponential distribution, 1/lambda=10') plt.xlabel('time') plt.ylabel('PDF') plt.show()

正態分布

正態分布是一種很常用的統計分布，可以描述現實世界的諸多事物，具備非常漂亮的性質，我們在下一講參數估計之中心極限定理時會詳細介紹。其概率密度函數為：

以下繪製了均值為0，標準差為1的正態分布的概率密度曲線，其形狀好似一口倒扣的鐘，因此也稱鐘形曲線。

def norm_pdf(x,mu,sigma):    '''正態分布概率密度函數'''    pdf = np.exp(-((x - mu)**2) / (2* sigma**2)) / (sigma * np.sqrt(2*np.pi))    return pdf mu = 0    # 均值為0sigma = 1 # 標準差為1# 用統計模擬繪製正態分布的直方圖sample = np.random.normal(mu, sigma, size=10000) plt. hist(sample, bins=100, alpha=0.7, normed=True)# 根據正態分布的公式繪製PDF曲線x = np.arange(-5, 5, 0.01) y = norm_pdf(x, mu, sigma) plt.plot(x,y, color='orange', lw=3) plt.show()

身高、體重的分布

以上從計算機模擬的角度出發，介紹了四種分布，現在讓我們看一下現實中的數據分布。繼續上一講數據探索之描述性統計中使用的BRFSS數據集，我們查看其中的身高和體重數據，看看他們是不是滿足正態分布。

首先導入數據，並編寫繪製PDF和CDF圖的函數 plot_pdf_cdf()，便於重複使用。

# 導入BRFSS數據import brfss df = brfss.ReadBrfss() height = df.height.dropna() weight = df.weight.dropna()

def plot_pdf_cdf(data, xbins, xrange, xlabel):    '''繪製概率密度函數PDF和累積分布函數CDF'''    fig = plt.figure(figsize=(16,5)) # 設置畫布尺寸    p1 = fig.add_subplot(121)  # 添加第一個子圖    # 繪製正態分布PDF曲線    std = data.std()    mean = data.mean()    x = np.arange(xrange[0], xrange[1], (xrange[1]-xrange[0])/100)    y = norm_pdf(x, mean, std)    plt.plot(x,y, label='normal distribution')    # 繪製數據的直方圖    plt.hist(data, bins=xbins, range=xrange, rwidth=0.9,             alpha=0.5, normed=True, label='observables')    # 圖片設置    plt.legend()    plt.xlabel(xlabel)    plt.title(xlabel +' PDF')    p2 = fig.add_subplot(122)  #添加第二個子圖    # 繪製正態分布CDF曲線    sample = np.random.normal(mean, std, size=10000)    plt.hist(sample, cumulative=True, bins=1000, range=xrange,             normed=True, histtype='step', lw=2, label='normal distribution')    # 繪製數據的CDF曲線    plt.hist(data, cumulative=True, bins=1000, range=xrange,             normed=True, histtype='step', lw=2, label='observables')    #圖片設置    plt.legend(loc='upper left')    plt.xlabel(xlabel)    plt.title( xlabel + ' CDF')    plt.show()

人群的身高分布比較符合正態分布。

plot_pdf_cdf(data=height, xbins=21, xrange=(1.2, 2.2), xlabel='height')

但是體重分布明顯右偏，與對稱的正態分布存在一定的差異。

plot_pdf_cdf(data=weight, xbins=60, xrange=(0,300), xlabel='weight')

將體重數據取對數值後，其分布就與正態分布非常吻合。

log_weight = np.log(weight) plot_pdf_cdf(data=log_weight, xbins=53, xrange=(3,6), xlabel='log weight')

參考資料：