前面兩小節對隨機變量做了一個概述:
數理統計與概率論及Python實現(1)——概率論中基本概念
數理統計與概率論及Python實現(2)——隨機變量
數理統計與概率論及Python實現(3)——隨機變量概述
這一節主要記錄一維離散型隨機變量以及關於它們的一些性質。對於概率論與數理統計方面的計算及可視化,主要的Python包有scipy, numpy和matplotlib等。
scipy是Python中使用最為廣泛的科學計算工具包,再加上numpy和matplotlib,基本上可以處理大部分的計算和作圖任務。下面是對scipy的介紹:
SciPy是一個開源的Python算法庫和數學工具包。SciPy包含的模塊有最優化、線性代數、積分、插值、特殊函數、快速傅立葉變換、信號處理和圖像處理、常微分方程求解和其他科學與工程中常用的計算。與其功能相類似的軟體還有MATLAB、GNU Octave和Scilab。SciPy目前在BSD許可證下發布。它的開發由Enthought資助。
上面的介紹中沒有提到stats模塊,這個模塊中包含了概率論及統計相關的函數。
下面是調用一個分布函數的常用方法(以最常見的正態分布為例):
初始化一個分布函數(也叫作凍結的分布);
調用該分布函數的方法或計算其數值特徵;
import numpy as npfrom scipy import statsimport matplotlib.pyplot as plt
norm_dis = stats.norm(5, 3) # 利用相應的分布函數及參數,創建一個凍結的正態分布(frozen distribution)x = np.linspace(-5, 15, 101) # 在區間[-5, 15]上均勻的取101個點
# 計算該分布在x中個點的概率密度分布函數值(PDF)pdf = norm_dis.pdf(x)
# 計算該分布在x中個點的累計分布函數值(CDF)cdf = norm_dis.cdf(x)
# 下面是利用matplotlib畫圖plt.figure(1)# plot pdfplt.subplot(211) # 兩行一列,第一個子圖plt.plot(x, pdf, 'b-', label='pdf')plt.ylabel('Probability')plt.title(r'PDF/CDF of normal distribution')plt.text(-5.0, .12, r'$\mu=5,\ \sigma=3$') # 3是標準差,不是方差plt.legend(loc='best', frameon=False)# plot cdfplt.subplot(212)plt.plot(x, cdf, 'r-', label='cdf')plt.ylabel('Probability')plt.legend(loc='best', frameon=False)
plt.show()
正態分布N(5,3**2)的概率密度函數和累計分布函數
2. 伯努利分布
伯努利分布應該是所有分布裡面最簡單的分布,也是二項分布的基本單元。其樣本空間中只有兩個點,一般取為{0、1};不同的伯努利分布只是取到這兩個值的概率不同。如果將拋一次硬幣看作是一次伯努利實驗,且將正面朝上記做1,反面朝上記做0。那麼伯努利分布中的參數p就表示硬幣正面朝上的概率,其期望同參數p,也是關注硬幣正面朝上的概率。
每種分布都是一種模型,都有其適用的實例。伯努利分布適合於試驗結果只有兩種可能的單次試驗。例如拋一次硬幣,其結果只有正面或反面兩種可能;一次產品質量檢測,其結果只有合格或不合格兩種可能。
1.1 定義
伯努利分布(英語:Bernoulli distribution,又名兩點分布或者0-1分布,是一個離散型概率分布,為紀念瑞士科學家雅各布·伯努利而命名。)若伯努利試驗成功,則伯努利隨機變量取值為1。若伯努利試驗失敗,則伯努利隨機變量取值為0。記其成功概率為
p(0≤p≤1),失敗概率為 q=1−p。則其概率質量函數(PMF)為:
伯努利分布只有一個參數p,伯努利分布只有一個參數p,記做X∼Bernoulli(p),或X∼B(1,p),讀作X服從參數為p的伯努利分布。
1.2 Python的實現
def binom_dis(n=1, p=0.1): """ 二項分布,模擬拋硬幣試驗 :param n: 實驗總次數 :param p: 單次實驗成功的概率 :return: 試驗成功的次數 """ binom_dis = stats.binom(n, p) simulation_result = binom_dis.rvs(size=5) # 取5個符合該分布的隨機變量 print(simulation_result) # [ 7 11 13 8 13], 每次結果會不一樣 prob_10 = binom_dis.pmf(10) print(prob_10) # 0.117
binom_dis(n=20, p=0.6)
柱狀圖表示的伯努利分布
B(1,0.3)的PMF為了得到比較準確的某個服從伯努利分布的隨機變量的期望,需要大量重複伯努利試驗,例如重複n次,然後利用"正面朝上的次數/n"來估計p。
3. 二項分布
如果把一個伯努利分布獨立的重複n次,就得到了一個二項分布。二項分布是最重要的離散型概率分布之一。隨機變量
X要滿足這個分布有兩個重要條件:各次試驗的條件是穩定的;
各次試驗之間是相互獨立的。
還是利用拋硬幣的例子來比較伯努利分布和二項分布:如果將拋一次硬幣看作是一次伯努利實驗,且將正面朝上記做1,反面朝上記做0。那麼拋n次硬幣,記錄正面朝上的次數Y,Y就服從二項分布。假如硬幣是均勻的,Y的取值應該大部分都集中在n/2附近,而非常大或非常小的值都很少。由此可見,二項分布關注的是計數,伯努利分布關注的是比值(正面朝上的計數/n)。
現實生活中有許多現象程度不同地符合這些條件,例如經常用來舉例子的拋硬幣,擲骰子等。如果每次試驗條件都相同,那麼硬幣正面朝上的次數以及某一個點數出現的次數都是非常典型的符合二項分布的隨機變量。均勻硬幣拋1000次,則正面朝上的次數X∼Binomial(1000,0.5);有六個面的骰子,擲100次,則6點出現的次數X∼Binomial(100,1/6)
2.1 定義
二項分布有兩個參數——試驗次數n和每次試驗成功的概率p. 其概率質量函數為:
2.2 Python的實現
下面的代碼用來模擬拋一枚不均勻的硬幣20次,其中正面朝上的概率為0.6
def binom_dis(n=1, p=0.1): """ 二項分布,模擬拋硬幣試驗 :param n: 實驗總次數 :param p: 單次實驗成功的概率 :return: 試驗成功的次數 """ binom_dis = stats.binom(n, p) simulation_result = binom_dis.rvs(size=5) # 取5個符合該分布的隨機變量 print(simulation_result) # [ 7 11 13 8 13], 每次結果會不一樣 prob_10 = binom_dis.pmf(10) print(prob_10) # 0.117
binom_dis(n=20, p=0.6)上面定義了一個n=20,p=0.6的二項分布,意思是說每次試驗拋硬幣(該硬幣正面朝上的概率大於背面朝上的概率)20次並記錄正面朝上的次數。
第9行"size=5"表示這樣的試驗重複了5次;第10行是試驗結果(第一次試驗,正面朝上出現了7次;第二次試驗,正面朝上出現了11次...);第11行表示計算正面朝上的次數為10的概率,由於每次試驗拋硬幣20次,因此試驗結果從0到20都有可能,只是概率不同而已。下面是該分布的概率質量分布函數圖:
二項分布
B(20,0.6)的PMF從圖2-1中可以明顯看到該分布的概率質量分布函數圖明顯向右邊偏移,在x=12處取到最大概率。這是因為這個硬幣正面朝上的概率大於反面朝上的概率。
為了比較準確的得到某個服從二項分布的隨機變量的期望,需要大量重複二項分布試驗,例如有m個人進行試驗(每人拋n次),然後利用"所有人得到的正面次數之和/m"來估計np。總共相當於做了nm次伯努利實驗。
PMF的實現:
def binom_pmf(n=1, p=0.1): """ 二項分布有兩個參數 https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binom.html#scipy.stats.binom :param n:試驗次數 :param p:單次實驗成功的概率 :return: """ binom_dis = stats.binom(n, p) x = np.arange(binom_dis.ppf(0.0001), binom_dis.ppf(0.9999)) print(x) # [ 0. 1. 2. 3. 4.] fig, ax = plt.subplots(1, 1) ax.plot(x, binom_dis.pmf(x), 'bo', label='binom pmf') ax.vlines(x, 0, binom_dis.pmf(x), colors='b', lw=5, alpha=0.5) ax.legend(loc='best', frameon=False) plt.ylabel('Probability') plt.title('PMF of binomial distribution(n={}, p={})'.format(n, p)) plt.show()
binom_pmf(n=20, p=0.6)
4. 泊松分布日常生活中,大量事件的發生是有固定頻率的。例如某醫院平均每小時出生3個嬰兒,某網站平均每分鐘有2次訪問等。它們的特點就是,我們可以預估這些事件在某個時間段內發生的總次數,但是沒法知道具體的發生時間。已知平均每小時出生3個嬰兒,請問下一個小時,會出生幾個?
有可能一下子出生6個,也有可能一個都不出生。這是我們沒法知道的。
如果某事件以固定強度
λλ,隨機且獨立地出現,該事件在單位時間內出現的次數(個數)可以看成是服從泊松分布。泊松分布適合於描述單位時間內隨機事件發生的次數的概率分布。如某一服務設施在一定時間內受到的服務請求的次數,電話交換機接到呼叫的次數、汽車站臺的候客人數、機器出現的故障數、自然災害發生的次數、DNA序列的變異數、放射性原子核的衰變數、雷射的光子數分布等等。
3.1 定義
泊松分布有一個參數λ(有的地方表示為μ),表示單位時間(或單位面積)內隨機事件的平均發生次數,其PMF表示為:
下面是參數μ=8時的泊松分布的概率質量分布圖(在scipy中將泊松分布的參數表示為μ):
泊松分布
P(8)的PMF
def poisson_pmf(mu=3): """ 泊松分布 https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.poisson.html#scipy.stats.poisson :param mu: 單位時間(或單位面積)內隨機事件的平均發生率 :return: """ poisson_dis = stats.poisson(mu) x = np.arange(poisson_dis.ppf(0.001), poisson_dis.ppf(0.999)) print(x) fig, ax = plt.subplots(1, 1) ax.plot(x, poisson_dis.pmf(x), 'bo', ms=8, label='poisson pmf') ax.vlines(x, 0, poisson_dis.pmf(x), colors='b', lw=5, alpha=0.5) ax.legend(loc='best', frameon=False) plt.ylabel('Probability') plt.title('PMF of poisson distribution(mu={})'.format(mu)) plt.show()
poisson_pmf(mu=8)
5. 泊松分布與二項分布的關係
如果僅僅是看二項分布與泊松分布的概率質量分布圖,也可以發現它們的相似度非常高。事實上這兩個分布內在聯繫十分緊密。泊松分布可以作為二項分布的極限得到。一般來說,若X∼B(n,p),其中nn很大,pp很小,而np=λ不太大時,則X的分布接近於泊松分布P(λ).
從下圖中可以非常直觀的看到兩者的關係:
更多內容,歡迎交流學習