統計-stats

時間 2019-11-11

標籤統計 stats 简体版

原文原文鏈接

統計-stats

SciPy的stats模塊包含了多種機率分佈的隨機變量[1]，隨機變量分爲連續和離散兩種。全部的連續隨機變量都是rv_continuous的派生類的對象，而全部的離散隨機變量都是rv_discrete的派生類的對象。html

Footnotes數組

[1]	本節中的隨機變量是指機率論中的概念，不是Python中的變量

連續和離散機率分佈

可使用下面的語句得到stats模塊中全部的連續隨機變量：dom

>>> from scipy import stats
>>> [k for k,v in stats.__dict__.items() if isinstance(v, stats.rv_continuous)]
['genhalflogistic','triang','rayleigh','betaprime', ...]

連續隨機變量對象都有以下方法：函數

rvs：對隨機變量進行隨機取值，能夠經過size參數指定輸出的數組的大小。
pdf：隨機變量的機率密度函數。
cdf：隨機變量的累積分佈函數，它是機率密度函數的積分。
sf：隨機變量的生存函數，它的值是1-cdf(t)。
ppf：累積分佈函數的反函數。
stat：計算隨機變量的指望值和方差。
fit：對一組隨機取樣進行擬合，找出最適合取樣數據的機率密度函數的係數。

03-scipy/scipy_stats.pyspa

機率密度函數、直方圖統計和累積分佈函數code

下面以正規分佈爲例，簡單地介紹隨機變量的用法。下面的語句得到缺省正規分佈的隨機變量的指望值和方差，咱們看到缺省狀況下它是一個均值爲0，方差爲1的隨機變量：orm

>>> stats.norm.stats()
(array(0.0), array(1.0))

norm能夠像函數同樣調用，經過loc和scale參數能夠指定隨機變量的偏移和縮放參數。對於正態分佈的隨機變量來講，這兩個參數至關於指定其指望值和標準差[2]：htm

>>> X = stats.norm(loc=1.0, scale=2.0)
>>> x.stats()
(array(1.0), array(4.0))

下面調用隨機變量X的rvs()方法，獲得包含一萬次隨機取樣值的數組x，而後調用NumPy的mean()和var()計算此數組的均值和方差，其結果符合隨機變量X的特性：對象

>>> x = X.rvs(size=10000) # 對隨機變量取10000個值
>>> np.mean(x) # 指望值
1.0181259658732724
>>> np.var(x)  # 方差
4.00188661640059

咱們也可使用fit()方法對隨機取樣序列x進行擬合，它返回的是與隨機取樣值最吻合的隨機變量的參數：排序

>>> stats.norm.fit(x) # 獲得隨機序列指望值和標準差
array([ 1.01810091,  2.00046946])

接下來比較隨機變量X的機率密度函數和對數組x進行直方圖統計的結果：

>>> t = np.arange(-10, 10, 0.01)
>>> pl.plot(t, X.pdf(t)) # 繪製機率密度函數的理論值
>>> p, t2 = np.histogram(x, bins=100, normed=True)
>>> t2 = (t2[:-1] + t2[1:])/2
>>> pl.plot(t2, p) # 繪製統計所獲得的機率密度

其中histogram()對數組x進行直方圖統計，它將數組x的取值範圍分爲100個區間，並統計x中的每一個值落入各個區間中的次數。histogram()返回兩個數組p和t2，其中p表示各個區間的取樣值出現的頻數，因爲normed參數爲True，所以p的值是正規化以後的結果。t2表示區間，因爲其中包括區間起點和終點，所以t2的長度爲101。【圖:正規分佈的機率密度函數(左)和累積分佈函數(右)】(左)顯示了機率密度函數和直方圖統計的結果，能夠看出兩者是一致的。

下面的程序繪製隨機變量X的累積分佈函數和數組p的累加結果，其結果如【圖:正規分佈的機率密度函數(左)和累積分佈函數(右)】(右)所示。

>>> pl.plot(t, X.cdf(t))
>>> pl.plot(t2, np.add.accumulate(p)*(t2[1]-t2[0]))

正規分佈的機率密度函數(左)和累積分佈函數(右)

有些隨機分佈除了loc和scale參數以外，還須要額外的形狀參數。例如伽瑪分佈可用於描述等待個獨立的隨機事件發生所需的時間，就是伽瑪分佈的形狀參數。下面計算形狀參數爲1和2時的伽瑪分佈的指望值和方差：

>>> stats.gamma.stats(1.0)
(array(1.0), array(1.0))
>>> stats.gamma.stats(2.0)
(array(2.0), array(2.0))

伽瑪分佈的尺度參數 $\theta$ 和隨機事件發生的頻率相關，由scale參數指定：

>>> stats.gamma.stats(2.0, scale=2)
(array(4.0), array(8.0))

根據伽瑪分佈的數學定義可知其指望值爲 $k \theta$ ，方差爲 $k \theta^2$ 。上面的程序驗證了這兩個公式。

當隨機分佈有額外的形狀參數時，它所對應的rvs()、pdf()等方法都會增長額外的參數接收形狀參數。例以下面的程序調用rvs()對 $k=2,\theta=2$ 的伽瑪分佈取4個隨機值：

>>> x = stats.gamma.rvs(2, scale=2, size=4)
>>> x
array([ 2.20814048,  3.56652153,  4.30088176,  0.68262888])

接下來調用pdf()查看上面4個隨機值所對應的機率密度：

>>> stats.gamma.pdf(x, 2, scale=2)
array([ 0.18301012,  0.1498734 ,  0.12519094,  0.12130919])

也能夠先建立將形狀參數和尺度參數固定的隨機變量，而後再調用其pdf()計算機率密度：

>>> X = stats.gamma(2, scale=2)
>>> X.pdf(x)
array([ 0.18301012,  0.1498734 ,  0.12519094,  0.12130919])

當分佈函數的值域爲離散時咱們稱之爲離散機率分佈。例如投擲有六個面的骰子時，只能得到1到6的整數，所以所獲得的機率分佈爲離散的。對於離散隨機分佈，一般使用機率質量函數(PMF)描述其分佈狀況。

在stats庫中全部描述離散分佈的隨機變量都從rv_discrete類繼承。也能夠直接用rv_discrete類自定義離散機率分佈。例如假設有一個不均勻的骰子，它的各點出現的機率不相等。咱們能夠用下面的數組x保存骰子的全部可能值，數組p保存每一個值出現的機率：

>>> x = range(1,7)
>>> p = (0.4, 0.2, 0.1, 0.1, 0.1, 0.1)

因而能夠用下面的語句定義表示這個特殊骰子的隨機變量，並調用其rvs()方法投擲此骰子20次，得到符合機率p的隨機數：

>>> dice = stats.rv_discrete(values=(x,p))
>>> dice.rvs(size=20)
array([2, 5, 1, 2, 1, 1, 2, 4, 1, 3, 1, 1, 4, 3, 1, 1, 1, 2, 6, 4])

Footnotes

[2]	標準差是方差的算術平方根，所以標準差爲2.0時，方差爲4.0

二項、泊松、伽瑪分佈

本節用幾個實例程序對機率論中的二項分佈、泊松分佈以及伽瑪分佈進行一些實驗和討論。

二項分佈是最重要的離散機率分佈之一。假設有一種只有兩個結果的試驗，其成功機率爲，那麼二項分佈描述了進行次這樣的獨立試驗，成功次的機率。二項分佈的機率質量函數公式以下：

$f(k;n,p) = \frac{n!}{k!(n-k)!} p^k(1-p)^{n-k}$

例如能夠經過二項分佈的機率質量公式計算投擲5次骰子出現3次6點的機率。投擲一次骰子，點數爲6的機率(即試驗成功的機率)爲，試驗次數爲。使用二項分佈的機率質量函數pmf()能夠很容易計算出現次6點的機率。和機率密度函數pdf()相似，pmf()的第一個參數爲隨機變量的取值，後面的參數爲描述隨機分佈所需的參數。對於二項分佈來講，參數分別爲和，而取值範圍則爲0到之間的整數。下面的程序計算爲0到6所對應的機率：

>>> stats.binom.pmf(range(6), 5, 1/6.0)
array([0.401878, 0.401878, 0.160751, 0.032150, 0.003215, 0.000129])

由結果可知：出現0或1次6點的機率爲40.2%，而出現3次6點的機率爲3.215%。

在二項分佈中，若是試驗次數很大，而每次試驗成功的機率很小，其乘積比較適中時，那麼試驗成功的次數的機率能夠用泊松分佈近似描述。

在泊松分佈中使用 $\lambda$ 描述單位時間(或單位面積)中隨機事件的平均發生率。若是咱們將二項分佈中的試驗次數看做單位時間中所作的試驗次數，那麼它和事件出現的機率的乘積就是事件的平均發生率，即 $\lambda=n p$ 。泊松分佈的機率質量函數公式以下：

$f(k;\lambda)=\frac{e^{-\lambda}\lambda^k}{k!}$

下面的程序分別計算二項分佈和泊松分佈的機率質量函數，其結果如【圖:當n足夠大時二項分佈和泊松分佈近似相等】所示，能夠看出當足夠大時，兩者是十分接近的。

03-scipy/scipy_binom_poisson.py

比較二項分佈和泊松分佈的機率質量函數

程序中事件平均發生率 $\lambda$ 恆等於10。根據二項分佈的試驗次數，計算每次事件出現的機率 $p=\lambda/n$ 。程序中運算部分大體以下：

>>> _lambda = 10.0
>>> k = np.arange(20)
>>> possion = stats.poisson.pmf(k, _lambda)  # 泊松分佈
>>> binom100 = stats.binom.pmf(k, 100, _lambda/100)    # 二項分佈 n=100
>>> binom1000 = stats.binom.pmf(k, 1000, _lambda/1000) # 二項分佈 n=1000
>>> np.max(np.abs(binom100-possion)) # 計算最大偏差
0.006755311103353312
>>> np.max(np.abs(binom1000-possion)) # n爲1000時，偏差較小
0.00063017540509099912

當n足夠大時二項分佈和泊松分佈近似相等

泊松分佈適合描述單位時間內隨機事件發生的次數的分佈狀況。例如某一設施在必定時間內的使用次數，機器出現故障的次數，天然災害發生的次數等等。

爲了加深讀者對泊松分佈概念的理解，下面咱們使用隨機數模擬泊松分佈，並與其機率質量函數進行比較，其結果如【圖:模擬泊松分佈】所示。圖中，每秒內事件的平均發生次數爲10，即 $\lambda=10$ 。其中左圖的觀察時間爲1000秒，而右圖的觀察時間爲50000秒。能夠看出觀察時間越長，每秒內事件發生的次數越符合泊松分佈。

03-scipy/scipy_poisson_sim.py

模擬泊松分佈

因爲上面的程序中包含了許多繪圖方面的代碼，下面咱們直接在IPython中介紹泊松分佈的模擬過程。首先定義事件發生率 $\lambda$ 和觀察時間:

>>> _lambda = 10
>>> time = 10000

能夠用NumPy的隨機數生成函數rand()，產平生均分佈於0到time之間的_lambda*time個事件所發生的時刻。因爲rand()產生的是0到1之間的平均分佈的隨機數，所以須要對其結果擴大time倍：

>>> t = np.random.rand(_lambda*time)*time

用histogram()能夠統計數組t中每秒以內的事件發生的次數count：

>>> count, time_edges = np.histogram(t, bins=time, range=(0,time))
>>> count
array([10,  9,  8, ..., 11, 10, 18])

根據泊松分佈的定義，count數組中的數值的分佈狀況應該符合泊松分佈。下面統計事件次數在0到20區間內的機率分佈。當histogram()的normed參數爲True而且每一個統計區間的長度爲1時，其結果和機率質量函數相等。

>>> dist, count_edges = np.histogram(count, bins=20, range=(0,20), normed=True)
>>> poisson = stats.poisson.pmf(x, _lambda)
>>> np.max(np.abs(dist-possion)) # 最大偏差很小，符合泊松分佈
0.0088356241037075706

還能夠換一個角度看隨機事件的分佈問題。咱們能夠觀察相鄰兩個事件之間的時間間隔的分佈狀況，或者隔k個事件的時間間隔的分佈狀況。根據機率論，事件之間的時間間隔應符合伽瑪分佈，因爲時間間隔能夠是任意數值，所以伽瑪分佈是一種連續機率分佈。伽瑪分佈的機率密度函數公式以下，它描述第個事件發生所需的等待時間的機率分佈。 $\Gamma(k)$ 是伽瑪函數，當爲整數時，它的值和的階乘相等。

$f(X;k, \lambda)=\frac{X^{(k-1)}\lambda^{k} e^{(-\lambda X)}}{\Gamma(k)}$

下面的程序模擬了事件的時間間隔的伽瑪分佈，其結果如【圖:模擬伽瑪分佈】所示。圖中的觀察時間爲1000秒，平均每秒產生10個事件。左圖中「k=1」，它表示相鄰兩個事件間的間隔的分佈，而「k=2」則表示相隔一個事件的兩個事件間的間隔的分佈，能夠看出它們都符合伽瑪分佈。

03-scipy/scipy_gamma_sim.py

模擬伽瑪分佈

下面咱們直接在IPython中模擬伽瑪分佈。首先在10000秒以內產生100000個隨機事件發生的時刻。所以事件的平均發生次數爲每秒10次：

>>> _lambda = 10
>>> time = 10000
>>> t = np.random.rand(_lambda*time)*time

爲了計算事件先後的時間間隔，須要先對隨機時刻進行排序：

>>> t.sort()

而後分別計算「k=1」和「k=2」時的時間間隔：

>>> s1 = t[1:] - t[:-1]    #相鄰兩事件的時間間隔
>>> s2 = t[2:] - t[:-2]    #相隔一個事件的兩個事件的時間間隔

對s1和s2分別調用histogram()進行機率統計，設置normed爲True能夠直接統計機率密度：

>>> dist1, x1 = np.histogram(s1, bins=100, normed=True)
>>> dist2, x2 = np.histogram(s2, bins=100, normed=True)

histogram()返回的第二個值爲統計區間的邊界，下面gamma.pdf()計算伽瑪分佈的機率密度時，使用各個區間的中值進行計算。pdf()的第二個參數爲k值，scale參數爲 $1/\lambda$ ：

>>> gamma1 = stats.gamma.pdf((x1[:-1]+x1[1:])/2, 1, scale=1.0/_lambda)
>>> gamma2 = stats.gamma.pdf((x2[:-1]+x2[1:])/2, 2, scale=1.0/_lambda)
>>> np.max(np.abs(gamma1 - dist1))
0.13557317865888141
>>> np.max(np.abs(gamma2 - dist2))
0.087375030861794656

因爲機率密度函數的值自己比較大，所以上面的偏差已經很小了：

>>> np.max(gamma1), np.max(gamma2)
(9.3483221580498537, 3.6767953241013656)

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。