Python – 伯樂在線
我們提供的服務(wù)有：網(wǎng)站設(shè)計制作、成都網(wǎng)站建設(shè)、微信公眾號開發(fā)、網(wǎng)站優(yōu)化、網(wǎng)站認證、橋西ssl等。為1000多家企事業(yè)單位解決了網(wǎng)站和推廣的問題。提供周到的售前咨詢和貼心的售后服務(wù)，是有科學管理、有技術(shù)的橋西網(wǎng)站制作公司
首頁所有文章觀點與動態(tài)基礎(chǔ)知識系列教程實踐項目工具與框架工具資源Python小組伯樂在線 Python - 伯樂在線 所有文章 實踐項目 如何在Python中實現(xiàn)這五類強大的概率分布如何在Python中實現(xiàn)這五類強大的概率分布
2015/04/25 · 實踐項目 · 概率分布
分享到： 12
本文由 伯樂在線 - feigao.me 翻譯，Daetalus 校稿。未經(jīng)許可，禁止轉(zhuǎn)載！
英文出處：。歡迎加入翻譯組。
R編程語言已經(jīng)成為統(tǒng)計分析中的事實標準。但在這篇文章中，我將告訴你在Python中實現(xiàn)統(tǒng)計學概念會是如此容易。我要使用Python實現(xiàn)一些離散和連續(xù)的概率分布。雖然我不會討論這些分布的數(shù)學細節(jié)，但我會以鏈接的方式給你一些學習這些統(tǒng)計學概念的好資料。在討論這些概率分布之前，我想簡單說說什么是隨機變量（random variable）。隨機變量是對一次試驗結(jié)果的量化。
舉個例子，一個表示拋硬幣結(jié)果的隨機變量可以表示成Python
X = {1 如果正面朝上,
2 如果反面朝上}
12X = {1 如果正面朝上,
2 如果反面朝上}
隨機變量是一個變量，它取值于一組可能的值（離散或連續(xù)的），并服從某種隨機性。隨機變量的每個可能取值的都與一個概率相關(guān)聯(lián)。隨機變量的所有可能取值和與之相關(guān)聯(lián)的概率就被稱為概率分布（probability distributrion）。
我鼓勵大家仔細研究一下scipy.stats模塊。
概率分布有兩種類型：離散（discrete）概率分布和連續(xù)（continuous）概率分布。
離散概率分布也稱為概率質(zhì)量函數(shù)（probability mass function）。離散概率分布的例子有伯努利分布（Bernoulli distribution）、二項分布（binomial distribution）、泊松分布（Poisson distribution）和幾何分布（geometric distribution）等。
連續(xù)概率分布也稱為概率密度函數(shù)（probability density function），它們是具有連續(xù)取值（例如一條實線上的值）的函數(shù)。正態(tài)分布（normal distribution）、指數(shù)分布（exponential distribution）和β分布（beta distribution）等都屬于連續(xù)概率分布。
若想了解更多關(guān)于離散和連續(xù)隨機變量的知識，你可以觀看可汗學院關(guān)于概率分布的視頻。
二項分布（Binomial Distribution）
服從二項分布的隨機變量X表示在n個獨立的是/非試驗中成功的次數(shù)，其中每次試驗的成功概率為p。
E(X) = np, Var(X) = np(1?p)
如果你想知道每個函數(shù)的原理，你可以在IPython筆記本中使用help file命令。 E(X)表示分布的期望或平均值。
鍵入stats.binom?了解二項分布函數(shù)binom的更多信息。
二項分布的例子：拋擲10次硬幣，恰好兩次正面朝上的概率是多少？
假設(shè)在該試驗中正面朝上的概率為0.3，這意味著平均來說，我們可以期待有3次是硬幣正面朝上的。我定義擲硬幣的所有可能結(jié)果為k = np.arange(0,11)：你可能觀測到0次正面朝上、1次正面朝上，一直到10次正面朝上。我使用stats.binom.pmf計算每次觀測的概率質(zhì)量函數(shù)。它返回一個含有11個元素的列表（list），這些元素表示與每個觀測相關(guān)聯(lián)的概率值。
您可以使用.rvs函數(shù)模擬一個二項隨機變量，其中參數(shù)size指定你要進行模擬的次數(shù)。我讓Python返回10000個參數(shù)為n和p的二項式隨機變量。我將輸出這些隨機變量的平均值和標準差，然后畫出所有的隨機變量的直方圖。
泊松分布（Poisson Distribution）
一個服從泊松分布的隨機變量X，表示在具有比率參數(shù)（rate parameter）λ的一段固定時間間隔內(nèi)，事件發(fā)生的次數(shù)。參數(shù)λ告訴你該事件發(fā)生的比率。隨機變量X的平均值和方差都是λ。
E(X) = λ, Var(X) = λ
泊松分布的例子：已知某路口發(fā)生事故的比率是每天2次，那么在此處一天內(nèi)發(fā)生4次事故的概率是多少？
讓我們考慮這個平均每天發(fā)生2起事故的例子。泊松分布的實現(xiàn)和二項分布有些類似，在泊松分布中我們需要指定比率參數(shù)。泊松分布的輸出是一個數(shù)列，包含了發(fā)生0次、1次、2次，直到10次事故的概率。我用結(jié)果生成了以下圖片。
你可以看到，事故次數(shù)的峰值在均值附近。平均來說，你可以預計事件發(fā)生的次數(shù)為λ。嘗試不同的λ和n的值，然后看看分布的形狀是怎么變化的。
現(xiàn)在我來模擬1000個服從泊松分布的隨機變量。
正態(tài)分布（Normal Distribution）
正態(tài)分布是一種連續(xù)分布，其函數(shù)可以在實線上的任何地方取值。正態(tài)分布由兩個參數(shù)描述：分布的平均值μ和方差σ2 。
E(X) = μ, Var(X) = σ2
正態(tài)分布的取值可以從負無窮到正無窮。你可以注意到，我用stats.norm.pdf得到正態(tài)分布的概率密度函數(shù)。
β分布（Beta Distribution）
β分布是一個取值在 [0, 1] 之間的連續(xù)分布，它由兩個形態(tài)參數(shù)α和β的取值所刻畫。
β分布的形狀取決于α和β的值。貝葉斯分析中大量使用了β分布。
當你將參數(shù)α和β都設(shè)置為1時，該分布又被稱為均勻分布（uniform distribution）。嘗試不同的α和β取值，看看分布的形狀是如何變化的。
指數(shù)分布（Exponential Distribution）
指數(shù)分布是一種連續(xù)概率分布，用于表示獨立隨機事件發(fā)生的時間間隔。比如旅客進入機場的時間間隔、打進客服中心電話的時間間隔、中文維基百科新條目出現(xiàn)的時間間隔等等。
我將參數(shù)λ設(shè)置為0.5，并將x的取值范圍設(shè)置為 $[0, 15]$ 。
接著，我在指數(shù)分布下模擬1000個隨機變量。scale參數(shù)表示λ的倒數(shù)。函數(shù)np.std中，參數(shù)ddof等于標準偏差除以 $n-1$ 的值。
結(jié)語（Conclusion）
概率分布就像蓋房子的藍圖，而隨機變量是對試驗事件的總結(jié)。我建議你去看看哈佛大學數(shù)據(jù)科學課程的講座，Joe Blitzstein教授給了一份摘要，包含了你所需要了解的關(guān)于統(tǒng)計模型和分布的全部。
如果隨機變量X的所有取值都可以逐個列舉出來，則稱X為離散型隨機變量。相應的概率分布有二項分布，泊松分布。
如果隨機變量X的所有取值無法逐個列舉出來，而是取數(shù)軸上某一區(qū)間內(nèi)的任一點，則稱X為連續(xù)型隨機變量。相應的概率分布有正態(tài)分布，均勻分布，指數(shù)分布，伽馬分布，偏態(tài)分布，卡方分布，beta分布等。(真多分布，好恐怖~~)
在離散型隨機變量X的一切可能值中，各可能值與其對應概率的乘積之和稱為該隨機變量X的期望值，記作E(X) 。比如有隨機變量，取值依次為：2，2，2，4，5。求其平均值：(2+2+2+4+5)/5 = 3。
期望值也就是該隨機變量總體的均值。 推導過程如下：
= (2+2+2+4+5)/5
= 1/5 2 3 + 4/5 + 5/5
= 3/5 2 + 1/5 4 + 1/5 5
= 0.6 2 + 0.2 4 + 0.2 5
= 60% 2 + 20% 4 + 20%*5
= 1.2 + 0.8 + 1
= 3
倒數(shù)第三步可以解釋為值為2的數(shù)字出現(xiàn)的概率為60%，4的概率為20%，5的概率為20%。 所以E(X) = 60% 2 + 20% 4 + 20%*5 = μ = 3。
0-1分布（兩點分布），它的隨機變量的取值為1或0。即離散型隨機變量X的概率分布為：P{X=0} = 1-p, P{X=1} = p，即：
則稱隨機變量X服從參數(shù)為p的0-1分布，記作X~B（1，p)。
在生活中有很多例子服從兩點分布，比如投資是否中標，新生嬰兒是男孩還是女孩，檢查產(chǎn)品是否合格等等。
大家非常熟悉的拋硬幣試驗對應的分布就是二項分布。拋硬幣試驗要么出現(xiàn)正面，要么就是反面，只包含這兩個結(jié)果。出現(xiàn)正面的次數(shù)是一個隨機變量，這種隨機變量所服從的概率分布通常稱為 二項分布 。
像拋硬幣這類試驗所具有的共同性質(zhì)總結(jié)如下：（以拋硬幣為例）
通常稱具有上述特征的n次重復獨立試驗為n重伯努利試驗。簡稱伯努利試驗或伯努利試驗概型。特別地，當試驗次數(shù)為1時，二項分布服從0-1分布(兩點分布)。
舉個栗子：拋3次均勻的硬幣，求結(jié)果出現(xiàn)有2個正面的概率 。
已知p = 0.5 (出現(xiàn)正面的概率) ，n = 3 ，k = 2
所以拋3次均勻的硬幣，求結(jié)果出現(xiàn)有2個正面的概率為3/8。
二項分布的期望值和方差 分別為：
泊松分布是用來描述在一 指定時間范圍內(nèi)或在指定的面積或體積之內(nèi)某一事件出現(xiàn)的次數(shù)的分布 。生活中服從泊松分布的例子比如有每天房產(chǎn)中介接待的客戶數(shù)，某微博每月出現(xiàn)服務(wù)器癱瘓的次數(shù)等等。 泊松分布的公式為 ：
其中 λ 為給定的時間間隔內(nèi)事件的平均數(shù)，λ = np。e為一個數(shù)學常數(shù)，一個無限不循環(huán)小數(shù)，其值約為2.71828。
泊松分布的期望值和方差 分別為：
使用Python繪制泊松分布的概率分布圖：
因為連續(xù)型隨機變量可以取某一區(qū)間或整個實數(shù)軸上的任意一個值，所以通常用一個函數(shù)f(x)來表示連續(xù)型隨機變量，而f(x)就稱為 概率密度函數(shù) 。
概率密度函數(shù)f(x)具有如下性質(zhì) ：
需要注意的是，f(x)不是一個概率，即f(x) ≠ P(X = x) 。在連續(xù)分布的情況下，隨機變量X在a與b之間的概率可以寫成：
正態(tài)分布（或高斯分布）是連續(xù)型隨機變量的最重要也是最常見的分布，比如學生的考試成績就呈現(xiàn)出正態(tài)分布的特征，大部分成績集中在某個范圍（比如60-80分），很小一部分往兩端傾斜（比如50分以下和90多分以上）。還有人的身高等等。
正態(tài)分布的定義 ：
如果隨機變量X的概率密度為( -∞x+∞)：
則稱X服從正態(tài)分布，記作X~N(μ,σ2)。其中-∞μ+∞，σ0， μ為隨機變量X的均值，σ為隨機變量X的標準差。 正態(tài)分布的分布函數(shù)
正態(tài)分布的圖形特點 ：
使用Python繪制正態(tài)分布的概率分布圖：
正態(tài)分布有一個3σ準則，即數(shù)值分布在(μ-σ,μ+σ)中的概率為0.6827，分布在（μ-2σ,μ+2σ)中的概率為0.9545，分布在(μ-3σ,μ+3σ)中的概率為0.9973，也就是說大部分數(shù)值是分布在(μ-3σ,μ+3σ)區(qū)間內(nèi)，超出這個范圍的可能性很小很小，僅占不到0.3%，屬于極個別的小概率事件，所以3σ準則可以用來檢測異常值。
當μ=0，σ=1時，有
此時的正態(tài)分布N(0,1) 稱為標準正態(tài)分布。因為μ，σ都是確定的取值，所以其對應的概率密度曲線是一條 形態(tài)固定 的曲線。
對標準正態(tài)分布，通常用φ(x)表示概率密度函數(shù)，用Φ(x)表示分布函數(shù)：
假設(shè)有一次物理考試特別難，滿分100分，全班只有大概20個人及格。與此同時語文考試很簡單，全班絕大部分都考了90分以上。小明的物理和語文分別考了60分和80分，他回家后告訴家長，這時家長能僅僅從兩科科目的分值直接判斷出這次小明的語文成績要比物理好很多嗎？如果不能，應該如何判斷呢？此時Z-score就派上用場了。 Z-Score的計算定義 ：
即 將隨機變量X先減去總體樣本均值，再除以總體樣本標準差就得到標準分數(shù)啦。如果X低于平均值，則Z為負數(shù)，反之為正數(shù) 。通過計算標準分數(shù)，可以將任何一個一般的正態(tài)分布轉(zhuǎn)化為標準正態(tài)分布。
小明家長從老師那得知物理的全班平均成績?yōu)?0分，標準差為10，而語文的平均成績?yōu)?2分，標準差為4。分別計算兩科成績的標準分數(shù)：
物理：標準分數(shù) = (60-40)/10 = 2
語文：標準分數(shù) = (85-95)/4 = -2.5
從計算結(jié)果來看，說明這次考試小明的物理成績在全部同學中算是考得很不錯的，而語文考得很差。
指數(shù)分布可能容易和前面的泊松分布混淆，泊松分布強調(diào)的是某段時間內(nèi)隨機事件發(fā)生的次數(shù)的概率分布，而指數(shù)分布說的是 隨機事件發(fā)生的時間間隔 的概率分布。比如一班地鐵進站的間隔時間。如果隨機變量X的概率密度為：
則稱X服從指數(shù)分布，其中的參數(shù)λ0。 對應的分布函數(shù) 為：
均勻分布的期望值和方差 分別為：
使用Python繪制指數(shù)分布的概率分布圖：
均勻分布有兩種，分為 離散型均勻分布和連續(xù)型均勻分布 。其中離散型均勻分布最常見的例子就是拋擲骰子啦。拋擲骰子出現(xiàn)的點數(shù)就是一個離散型隨機變量，點數(shù)可能有1，2，3，4，5，6。每個數(shù)出現(xiàn)的概率都是1/6。
設(shè)連續(xù)型隨機變量X具有概率密度函數(shù)：
則稱X服從區(qū)間(a,b)上的均勻分布。X在等長度的子區(qū)間內(nèi)取值的概率相同。對應的分布函數(shù)為：
f(x)和F(x)的圖形分別如下圖所示：
均勻分布的期望值和方差 分別為：
如何根據(jù)概率密度函數(shù)生成隨機數(shù)
我這里并不是要講“偽隨機”、“真隨機”這樣的問題，而是關(guān)于如何生成服從某個概率分布的隨機數(shù)（或者說 sample）的問題。比如，你想要從一個服從正態(tài)分布的隨機變量得到 100 個樣本，那么肯定抽到接近其均值的樣本的概率要大許多，從而導致抽到的樣本很多是集中在那附近的。當然，要解決這個問題，我們通常都假設(shè)我們已經(jīng)有了一個 生成 0 到 1 之間均勻分布的隨機數(shù)的工具，就好像 random.org 給我們的結(jié)果那樣，事實上許多時候我們也并不太關(guān)心它們是真隨機數(shù)還是偽隨機數(shù)，看起來差不多就行了。 :p
現(xiàn)在再回到我們的問題，看起來似乎是很簡單的，按照概率分布的話，只要在概率密度大的地方多抽一些樣本不就行了嗎？可是具體要怎么做呢？要真動起手 來，似乎有不是那么直觀了。實際上，這個問題曾經(jīng)也是困擾了我很久，最近又被人問起，那我們不妨在這里一起來總結(jié)一下。為了避免一下子就陷入抽象的公式推 導，那就還是從一個簡單的具體例子出發(fā)好了，假設(shè)我們要抽樣的概率分布其概率密度函數(shù)為 p(x) = \frac{1}{9}x^2 ，并且被限制在區(qū)間 [0, 3] 上，如右上圖所示。
好了，假設(shè)現(xiàn)在我們要抽 100 個服從這個分布的隨機數(shù)，直觀上來講，抽出來的接近 3 的數(shù)字肯定要比接近 0 的數(shù)字要多。那究竟要怎樣抽才能得到這樣的結(jié)果呢？由于我們實際上是不能控制最原始的隨機數(shù)生成過程的，我們只能得到一組均勻分布的隨機數(shù)，而這組隨機數(shù) 的生成過程對于我們完全是透明的，所以，我們能做的只有把這組均勻分布的隨機數(shù)做一些變換讓他符合我們的需求。找到下手的點了，可是究竟要怎樣變換呢？有 一個變換相信大家都是很熟悉的，假設(shè)我們有一組 [0,1] 之間的均勻分布的隨機數(shù) X_0 ，那么令 X_1=3X_0 的話，X_1 就是一組在 [0,3] 之間均勻分布的隨機數(shù)了，不難想象，X_1 等于某個數(shù) x^* 的概率就是 X_0 等于 x^*/3 的概率（“等于某個數(shù)的概率”這種說法對于連續(xù)型隨機變量來說其實是不合適的，不過大概可以理解所表達的意思啦）。似乎有一種可以“逆轉(zhuǎn)回去”的感覺了。
于是讓我們來考慮更一般的變換。首先，我們知道 X_1 的概率密度函數(shù)是 f(x) = 1/3, x\in[0,3] ，假設(shè)現(xiàn)在我們令 Y = \phi (X_1) ，不妨先假定 \phi(\cdot) 是嚴格單調(diào)遞增的函數(shù)，這樣我們可以求其逆函數(shù) \phi^{-1}(\cdot) （也是嚴格單調(diào)遞增的）。現(xiàn)在來看變換后的隨機變量 Y 會服從一個什么樣的分布呢？
這里需要小心，因為這里都是連續(xù)型的隨機變量，并不像離散型隨機變量那樣可以說成“等于某個值的概率”，因此我們需要轉(zhuǎn)換為概率分布函數(shù)來處理，也就是求一個積分啦：
\displaystyle F(x) = P(X \leq x) = \int_{-\infty}^x f(t)dt
那么 X_1 的概率分布函數(shù)為 F(x) = \frac{1}{3}x 。很顯然 Y 小于或等于某個特定的值 y^* 這件事情是等價于 X_1=\phi^{-1}(Y)\leq\phi^{-1}(y^*) 這件事情的。換句話說，P(Y\leq y^*) 等于 P(X_1 \leq \phi^{-1}(y^*)) 。于是，Y 的概率分布函數(shù)就可以得到了：
\displaystyle G(y) = P(Y \leq y) = P(X_1 \leq \phi^{-1}(y)) = F(\phi^{-1}(y))
再求導我們就能得到 Y 的概率密度函數(shù)：
\displaystyle g(y) = \frac{dG(y)}{dy} = f(\phi^{-1}(y))\fracsqu6kqw{dy}\phi^{-1}(y)
這樣一來，我們就得到了對于一個隨機變量進行一個映射 \phi(\cdot) 之后得到的隨即變量的分布，那么，回到我們剛才的問題，我們想讓這個結(jié)果分布就是我們所求的，然后再反推得 \phi(\cdot) 即可：
\displaystyle \frac{1}{9}y^2 = g(y) = f(\phi^{-1}(y))\fracsqu6kqw{dy}\phi^{-1}(y) = \frac{1}{3}\fracsqu6kqw{dy}\phi^{-1}(y)
經(jīng)過簡單的化簡就可以得到 \phi^{-1}(y) = \frac{1}{9} y^3 ，亦即 \phi(x) = (9x)^{1/3} 。也就是說，把得到的隨機數(shù) X_1 帶入到到函數(shù) \phi(\cdot) 中所得到的結(jié)果，就是符合我們預期要求的隨機數(shù)啦！ :D 讓我們來驗證一下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
#!/usr/bin/python import numpy as np import matplotlib.pyplot as plot N = 10000 X0 = np.random.rand(N) X1 = 3*X0 Y = np.power(9*X1, 1.0/3) t = np.arange(0.0, 3.0, 0.01) y = t*t/9 plot.plot(t, y, 'r-', linewidth=1) plot.hist(Y, bins=50, normed=1, facecolor='green', alpha=0.75)plot.show()
這就沒錯啦，目的達成啦！讓我們來總結(jié)一下。問題是這樣的，我們有一個服從均勻分布的隨機變量 X ，它的概率密度函數(shù)為一個常數(shù) f(x)=C ，如果是 [0,1] 上的分布，那么常數(shù) C 就直接等于 1 了?，F(xiàn)在我們要得到一個隨機變量 Y 使其概率密度函數(shù)為 g(y) ，做法就是構(gòu)造出一個函數(shù) \phi(\cdot) 滿足（在這里加上了絕對值符號，這是因為 \phi(\cdot) 如果不是遞增而是遞減的話，推導的過程中有一處就需要反過來）
\displaystyle g(y) = f(\phi^{-1}(y))\left|\fracsqu6kqw{dy}\phi^{-1}(y)\right| = C\left|\fracsqu6kqw{dy}\phi^{-1}(y)\right|
反推過來就是，對目標 y 的概率密度函數(shù)求一個積分（其實就是得到它的概率分布函數(shù) CDF ，如果一開始就拿到的是 CDF 當然更好），然后求其反函數(shù)就可以得到需要的變換 \phi(\cdot) 了。實際上，這種方法有一個聽起來稍微專業(yè)一點的名字：Inverse Transform Sampling Method 。不過，雖然看起來很簡單，但是實際操作起來卻比較困難，因為對于許多函數(shù)來說，求逆是比較困難的，求積分就更困難了，如果寫不出解析解，不得已只能用數(shù) 值方法來逼近的話，計算效率就很讓人擔心了?？墒聦嵣弦彩侨绱?，就連我們最常見的一維標準正態(tài)分布，也很難用這樣的方法來抽樣，因為它的概率密度函數(shù)
\displaystyle g(y) = \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}y^2}
的不定積分沒有一個解析形式。這可真是一點也不好玩，費了這么大勁，結(jié)果好像什么都干不了?？磥磉@個看似簡單的問題似乎還是比較復雜的，不過也不要灰心，至少對于高斯分布來說，我們還有一個叫做 Box Muller 的方法可以專門來做這個事情。因為高斯分布比較奇怪，雖然一維的時候概率分布函數(shù)無法寫出解析式，但是二維的情況卻可以通過一些技巧得出一個解析式來。
首先我們來考慮一個二維的且兩個維度相互獨立的高斯分布，它的概率密度函數(shù)為
\displaystyle f(x,y) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\cdot\frac{1}{\sqrt{2\pi}}e^{-\frac{y^2}{2}} = \frac{1}{2\pi}e^{-\frac{x^2+y^2}{2}}
這個分布是關(guān)于原點對稱的，如果考慮使用極坐標 (\theta,r) （其中 \theta\in[0,2\pi), r\in[0,\infty) ）的話，我們有 x = r\cos\theta，y=r\sin\theta 這樣的變換。這樣，概率密度函數(shù)是寫成：
\displaystyle f(\theta,r) = \frac{1}{2\pi}e^{-\frac{r^2}{2}}
注意到在給定 r 的情況下其概率密度是不依賴于 \theta 的，也就是說對于 \theta 來說是一個均勻分布，這和我們所了解的標準正態(tài)分布也是符合的：在一個圓上的點的概率是相等的。確定了 \theta 的分布，讓我們再來看 r，用類似于前面的方法：
\displaystyle \begin{aligned} P(rR) = \int_0^{2\pi}\int_0^R\frac{1}{2\pi}e^{\frac{r^2}{2}}rdrd\theta \ = \int_0^Re^{-\frac{r^2}{2}}rdr \ = 1-e^{-\frac{R^2}{2}} \end{aligned}
根據(jù)前面得出的結(jié)論，我現(xiàn)在得到了 r 的概率分布函數(shù)，是不是只要求一下逆就可以得到一個 \phi(\cdot) 了？亦即 \phi(t) = \sqrt{-2\log (1-t)} 。
現(xiàn)在只要把這一些線索串起來，假設(shè)我們有兩個相互獨立的平均分布在 [0,1] 上的隨機變量 T_1 和 T_2 ，那么 2\pi T_1 就可以得到 \theta 了，而 \phi(T_2) = \sqrt{-2\log(1-T_2)} 就得到 r 了（實際上，由于 T_2 和 1-T_2 實際上是相同的分布，所以通常直接寫為 \sqrt{-2\log T_2}）。再把極坐標換回笛卡爾坐標：
\displaystyle \begin{aligned} x = r\cos\theta = \sqrt{-2\log T_2}\cdot \cos(2\pi T_1) \ y = r\sin\theta = \sqrt{-2\log T_2}\cdot \sin(2\pi T_1) \end{aligned}
這樣我們就能得到一個二維的正態(tài)分布的抽樣了?？梢灾庇^地驗證一下，二維不太好畫，就畫成 heatmap 了，看著比較熱的區(qū)域就是概率比較大的，程序如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
#!/usr/bin/python import numpy as np import matplotlib.pyplot as plot N = 50000 T1 = np.random.rand(N) T2 = np.random.rand(N) r = np.sqrt(-2*np.log(T2)) theta = 2*np.pi*T1 X = r*np.cos(theta) Y = r*np.sin(theta) heatmap, xedges, yedges = np.histogram2d(X, Y, bins=80) extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]] plot.imshow(heatmap, extent=extent) plot.show()
畫出來的圖像這個樣子：
不太好看，但是大概的形狀是可以看出來的。其實有了二維的高斯分布，再注意到兩個維度在我們這里是相互獨立的，那么直接取其中任意一個維度，就是一個一維高斯分布了。如下：
如果 X\sim N(0,1) 即服從標準正態(tài)分布的話，則有 \sigma X+\mu \sim N(\mu, \sigma^2) ，也就是說，有了標準正態(tài)分布，其他所有的正態(tài)分布的抽樣也都可以完成了。這下總算有點心滿意足了。不過別急，還有最后一個問題：多元高斯分布。一般最常 用不就是二元嗎？二元不是我們一開始就推出來了嗎？推出來了確實沒錯，不過我們考慮的是最簡單的情形，當然同樣可以通過 \sigma X+\mu 這樣的方式來處理每一個維度，不過高維的情形還有一個需要考慮的就是各個維度之間的相關(guān)性——我們之前處理的都是兩個維度相互獨立的情況。對于一般的多維正態(tài)分布 X\sim N(\mathbf{\mu}, \Sigma) ，如果各個維度之間是相互獨立的，就對應于協(xié)方差矩陣 \Sigma 是一個對角陣，但是如果 \Sigma 在非對角線的地方存在非零元素的話，就說明對應的兩個維度之間存在相關(guān)性。
這個問題還是比較好解決的，高斯分布有這樣的性質(zhì)：類似于一維的情況，對于多維正態(tài)分布 X\sim N(\mathbf{\mu}, \Sigma)，那么新的隨機變量 X_1=\mathbf{\mu}_1 + LX 將會滿足
\displaystyle X_1 \sim N(\mathbf{\mu}_1+L\mu, L\Sigma L^T)
所以，對于一個給定的高斯分布 N(\mathbf{\mu}, \Sigma) 來說，只要先生成一個對應維度的標準正態(tài)分布 X\sim N(0, I) ，然后令 X_1 = \mu+LX 即可，其中 L 是對 \Sigma 進行 Cholesky Decomposition 的結(jié)果，即 \Sigma = LL^T 。
結(jié)束之前讓我們來看看 matlab 畫個 3D 圖來改善一下心情：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
N = 50000; T1 = rand(1, N); T2 = rand(1, N); r = sqrt(-2*log(T2)); theta = 2*pi*T1; X =[r.*cos(theta); r.*sin(theta)]; mu = [1; 2]; Sigma = [5 2; 2 1]; L = chol(Sigma); X1 = repmat(mu,1, N) + L*X; nbin = 30; hist3(X1', [nbin nbin]); set(gcf, 'renderer', 'opengl'); set(get(gca,'child'), 'FaceColor', 'interp', 'CDataMode', 'auto'); [z c] = hist3(X1', [nbin nbin]); [x y] =meshgrid(c{1}, c{2}); figure; surfc(x,y,-z);
下面兩幅圖，哪幅好看一些（注意坐標比例不一樣，所以看不出形狀和旋轉(zhuǎn)了）？似乎都不太好看，不過感覺還是比前面的 heatmap 要好一點啦！
然后，到這里為止，我們算是把高斯分布弄清楚了，不過這只是給一個介紹性的東西，里面的數(shù)學推導也并不嚴格，而 Box Muller 也并不是最高效的高斯采樣的算法，不過，就算我們不打算再深入討論高斯采樣，采樣這個問題本身也還有許多不盡人意的地方，我們推導出來的結(jié)論可以說只能用 于一小部分簡單的分布，連高斯分布都要通過 trick 來解決，另一些本身連概率密度函數(shù)都寫不出來或者有各種奇怪數(shù)學特性的分布就更難處理了。所以本文的標題里也說了，這是上篇，如果什么時候有機會抽出時間 來寫下篇的話，我將會介紹一些更加通用和強大的方法，諸如 Rejection Sampling 、Gibbs Sampling 以及 Markov Chain Monte Carlo (MCMC) 等方法。如果你比較感興趣，可以先自行 Google 一下解饞！ :D
今天小編就為大家分享一篇python高斯分布概率密度函數(shù)的使用詳解，具有很好的參考價值，希望對大家有所幫助。一起跟隨小編過來看看吧
Python正態(tài)分布概率計算方法，喜歡算法的伙伴們可以參考學習下。需要用到math模塊。先了解一下這個模塊方法，再來寫代碼會更好上手。
def st_norm(u):
'''標準正態(tài)分布'''
import math
x=abs(u)/math.sqrt(2)
T=(0.0705230784,0.0422820123,0.0092705272,
0.0001520143,0.0002765672,0.0000430638)
E=1-pow((1+sum([a*pow(x,(i+1))
for i,a in enumerate(T)])),-16)
p=0.5-0.5*E if u0 else 0.5+0.5*E
return(p)
def norm(a,sigma,x):
'''一般正態(tài)分布'''
u=(x-a)/sigma
return(st_norm(u))
while 1:
'''輸入一個數(shù)時默認為標準正態(tài)分布
輸入三個數(shù)(空格隔開)時分別為期望、方差、x
輸入 stop 停止'''
S=input('please input the parameters:\n')
if S=='stop':break
try:
L=[float(s) for s in S.split()]
except:
print('Input error!')
continue
if len(L)==1:
print('f(x)=%.5f'%st_norm(L[0]))
elif len(L)==3:
print('f(x)=%.5f'%norm(L[0],L[1],L[2]))
else:
print('Input error!')
重慶分公司
重慶分公司
