探索常用的數(shù)據(jù)分析統(tǒng)計(jì)分布
本文用Python統(tǒng)計(jì)模擬的方法，介紹四種常用的統(tǒng)計(jì)分布，包括離散分布：二項(xiàng)分布和泊松分布，以及連續(xù)分布，指數(shù)分布和正態(tài)分布，最后查看人群的身高和體重?cái)?shù)據(jù)所符合的分布。
# 導(dǎo)入相關(guān)模塊import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport seaborn as sns%matplotlib inline%config InlineBackend.figure_format = 'retina'
01 隨機(jī)數(shù)
計(jì)算機(jī)發(fā)明后，便產(chǎn)生了一種全新的解決問題的方式：使用計(jì)算機(jī)對(duì)現(xiàn)實(shí)世界進(jìn)行統(tǒng)計(jì)模擬。該方法又稱為“蒙特卡洛方法（Monte Carlo method）”，起源于二戰(zhàn)時(shí)美國(guó)研制原子彈的曼哈頓計(jì)劃，它的發(fā)明人中就有大名鼎鼎的馮·諾依曼。蒙特卡洛方法的名字來源也頗為有趣，相傳另一位發(fā)明者烏拉姆的叔叔經(jīng)常在摩洛哥的蒙特卡洛賭場(chǎng)輸錢，賭博是一場(chǎng)概率的游戲，故而以概率為基礎(chǔ)的統(tǒng)計(jì)模擬方法就以這一賭城命名了。
使用統(tǒng)計(jì)模擬，首先要產(chǎn)生隨機(jī)數(shù)，在Python中，numpy.random 模塊提供了豐富的隨機(jī)數(shù)生成函數(shù)。比如生成0到1之間的任意隨機(jī)數(shù)：
np.random.random(size=5) # size表示生成隨機(jī)數(shù)的個(gè)數(shù)array([ 0.32392203, 0.3373342 , 0.51677112, 0.28451491, 0.07627541])又比如生成一定范圍內(nèi)的隨機(jī)整數(shù)：
np.random.randint(1, 10, size=5) # 生成5個(gè)1到9之間的隨機(jī)整數(shù)array([5, 6, 9, 1, 7])計(jì)算機(jī)生成的隨機(jī)數(shù)其實(shí)是偽隨機(jī)數(shù)，是由一定的方法計(jì)算出來的，因此我們可以按下面方法指定隨機(jī)數(shù)生成的種子，這樣的好處是以后重復(fù)計(jì)算時(shí)，能保證得到相同的模擬結(jié)果。
np.random.seed(123)在NumPy中，不僅可以生成上述簡(jiǎn)單的隨機(jī)數(shù)，還可以按照一定的統(tǒng)計(jì)分布生成相應(yīng)的隨機(jī)數(shù)。這里列舉了二項(xiàng)分布、泊松分布、指數(shù)分布和正態(tài)分布各自對(duì)應(yīng)的隨機(jī)數(shù)生成函數(shù)，接下來我們分別研究這四種類型的統(tǒng)計(jì)分布。
np.random.binomial()
np.random.poisson()
np.random.exponential()
np.random.normal()
二項(xiàng)分布是n個(gè)獨(dú)立的是/非試驗(yàn)中成功的次數(shù)的概率分布，其中每次試驗(yàn)的成功概率為p。這是一個(gè)離散分布，所以使用概率質(zhì)量函數(shù)（PMF）來表示k次成功的概率：
最常見的二項(xiàng)分布就是投硬幣問題了，投n次硬幣，正面朝上次數(shù)就滿足該分布。下面我們使用計(jì)算機(jī)模擬的方法，產(chǎn)生10000個(gè)符合（n，p）的二項(xiàng)分布隨機(jī)數(shù)，相當(dāng)于進(jìn)行10000次實(shí)驗(yàn)，每次實(shí)驗(yàn)投擲了n枚硬幣，正面朝上的硬幣數(shù)就是所產(chǎn)生的隨機(jī)數(shù)。同時(shí)使用直方圖函數(shù)繪制出二項(xiàng)分布的PMF圖。
def plot_binomial(n,p):'''繪制二項(xiàng)分布的概率質(zhì)量函數(shù)'''sample = np.random.binomial(n,p,size=10000) # 產(chǎn)生10000個(gè)符合二項(xiàng)分布的隨機(jī)數(shù)bins = np.arange(n+2)plt.hist(sample, bins=bins, align='left', normed=True, rwidth=0.1) # 繪制直方圖#設(shè)置標(biāo)題和坐標(biāo)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， 那么出現(xiàn)正面次數(shù)的分布符合上圖所示的二項(xiàng)分布。該分布左右對(duì)稱，最有可能的情況是正面出現(xiàn)5次。
但如果這是一枚作假的硬幣呢？比如正面朝上的概率p=0.2，或者是p=0.8，又會(huì)怎樣呢？我們依然可以做出該情況下的PMF圖。
fig = plt.figure(figsize=(12,4.5)) #設(shè)置畫布大小p1 = fig.add_subplot(121) # 添加第一個(gè)子圖plot_binomial(10, 0.2)p2 = fig.add_subplot(122) # 添加第二個(gè)子圖plot_binomial(10, 0.8)
這時(shí)的分布不再對(duì)稱了，正如我們所料，當(dāng)概率p=0.2時(shí)，正面最有可能出現(xiàn)2次；而當(dāng)p=0.8時(shí)，正面最有可能出現(xiàn)8次。
03 泊松分布
泊松分布用于描述單位時(shí)間內(nèi)隨機(jī)事件發(fā)生次數(shù)的概率分布，它也是離散分布，其概率質(zhì)量函數(shù)為：
比如你在等公交車，假設(shè)這些公交車的到來是獨(dú)立且隨機(jī)的（當(dāng)然這不是現(xiàn)實(shí)），前后車之間沒有關(guān)系，那么在1小時(shí)中到來的公交車數(shù)量就符合泊松分布。同樣使用統(tǒng)計(jì)模擬的方法繪制該泊松分布，這里假設(shè)每小時(shí)平均來6輛車（即上述公式中l(wèi)ambda=6）。
lamb = 6sample = np.random.poisson(lamb, size=10000) # 生成10000個(gè)符合泊松分布的隨機(jī)數(shù)bins = np.arange(20)plt.hist(sample, bins=bins, align='left', rwidth=0.1, normed=True) # 繪制直方圖# 設(shè)置標(biāo)題和坐標(biāo)軸plt.title('Poisson PMF (lambda=6)')plt.xlabel('number of arrivals')plt.ylabel('probability')plt.show()
04 指數(shù)分布
比如上面等公交車的例子，兩輛車到來的時(shí)間間隔，就符合指數(shù)分布。假設(shè)平均間隔為10分鐘（即1/lambda=10)，那么從上次發(fā)車開始，你等車的時(shí)間就滿足下圖所示的指數(shù)分布。
tau = 10sample = np.random.exponential(tau, size=10000) # 產(chǎn)生10000個(gè)滿足指數(shù)分布的隨機(jī)數(shù)plt.hist(sample, bins=80, alpha=0.7, normed=True) #繪制直方圖plt.margins(0.02)# 根據(jù)公式繪制指數(shù)分布的概率密度函數(shù)lam = 1 / taux = np.arange(0,80,0.1)y = lam * np.exp(- lam * x)plt.plot(x,y,color='orange', lw=3)#設(shè)置標(biāo)題和坐標(biāo)軸plt.title('Exponential distribution, 1/lambda=10')plt.xlabel('time')plt.ylabel('PDF')plt.show()
05 正態(tài)分布
正態(tài)分布是一種很常用的統(tǒng)計(jì)分布，可以描述現(xiàn)實(shí)世界的諸多事物，具備非常漂亮的性質(zhì)，我們?cè)谙乱恢v參數(shù)估計(jì)之中心極限定理時(shí)會(huì)詳細(xì)介紹。其概率密度函數(shù)為：
以下繪制了均值為0，標(biāo)準(zhǔn)差為1的正態(tài)分布的概率密度曲線，其形狀好似一口倒扣的鐘，因此也稱鐘形曲線。
def norm_pdf(x,mu,sigma):'''正態(tài)分布概率密度函數(shù)'''pdf = np.exp(-((x - mu)**2) / (2* sigma**2)) / (sigma * np.sqrt(2*np.pi))return pdfmu = 0 # 均值為0sigma = 1 # 標(biāo)準(zhǔn)差為1# 用統(tǒng)計(jì)模擬繪制正態(tài)分布的直方圖sample = np.random.normal(mu, sigma, size=10000)plt. hist(sample, bins=100, alpha=0.7, normed=True)# 根據(jù)正態(tài)分布的公式繪制PDF曲線x = np.arange(-5, 5, 0.01)y?=?norm_pdf(x,?mu,?sigma)plt.plot(x,y, color='orange', lw=3)plt.show()
06 身高、體重的分布
以上從計(jì)算機(jī)模擬的角度出發(fā)，介紹了四種分布，現(xiàn)在讓我們看一下現(xiàn)實(shí)中的數(shù)據(jù)分布。繼續(xù)上一講數(shù)據(jù)探索之描述性統(tǒng)計(jì)中使用的BRFSS數(shù)據(jù)集，我們查看其中的身高和體重?cái)?shù)據(jù)，看看他們是不是滿足正態(tài)分布。
首先導(dǎo)入數(shù)據(jù)，并編寫繪制PDF和CDF圖的函數(shù) plot_pdf_cdf()，便于重復(fù)使用。
# 導(dǎo)入BRFSS數(shù)據(jù)import brfssdf = brfss.ReadBrfss()height = df.height.dropna()weight = df.weight.dropna()
def plot_pdf_cdf(data, xbins, xrange, xlabel):'''繪制概率密度函數(shù)PDF和累積分布函數(shù)CDF'''fig = plt.figure(figsize=(16,5)) # 設(shè)置畫布尺寸p1 = fig.add_subplot(121) # 添加第一個(gè)子圖# 繪制正態(tài)分布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')# 繪制數(shù)據(jù)的直方圖plt.hist(data, bins=xbins, range=xrange, rwidth=0.9,alpha=0.5, normed=True, label='observables')????#?圖片設(shè)置plt.legend()plt.xlabel(xlabel)plt.title(xlabel +' PDF')p2 = fig.add_subplot(122) #添加第二個(gè)子圖# 繪制正態(tài)分布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')????#?繪制數(shù)據(jù)的CDF曲線plt.hist(data, cumulative=True, bins=1000, range=xrange,normed=True, histtype='step', lw=2, label='observables')????#圖片設(shè)置plt.legend(loc='upper left')plt.xlabel(xlabel)plt.title( xlabel + ' CDF')plt.show()
人群的身高分布比較符合正態(tài)分布。
plot_pdf_cdf(data=height, xbins=21, xrange=(1.2, 2.2), xlabel='height')
但是體重分布明顯右偏，與對(duì)稱的正態(tài)分布存在一定的差異。
plot_pdf_cdf(data=weight, xbins=60, xrange=(0,300), xlabel='weight')
將體重?cái)?shù)據(jù)取對(duì)數(shù)值后，其分布就與正態(tài)分布非常吻合。
log_weight = np.log(weight)plot_pdf_cdf(data=log_weight, xbins=53, xrange=(3,6), xlabel='log weight')
參考資料：
維基百科：蒙特卡羅方法
《Think Stats 2》
《統(tǒng)計(jì)學(xué)》，William Mendenhall著
猜你喜歡
評(píng)論
圖片
表情