# 導入相關模塊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之間的任意隨機數：array([ 0.32392203, 0.3373342 , 0.51677112, 0.28451491, 0.07627541])np.random.randint(1, 10, size=5) 計算機生成的隨機數其實是偽隨機數，是由一定的方法計算出來的，因此我們可以按下面方法指定隨機數生成的種子，這樣的好處是以後重複計算時，能保證得到相同的模擬結果。在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) 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) 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) plt.hist(sample, bins=80, alpha=0.7, normed=True) plt.margins(0.02)
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 sigma = 1
sample = np.random.normal(mu, sigma, size=10000)plt. hist(sample, bins=100, alpha=0.7, normed=True)
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()，便於重複使用。import brfssdf = 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) 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) 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') 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')