機器學習項目中不可忽視的一個密辛 - 大數定理、中心極限定理

時間 2019-11-06

原文原文鏈接

1. 前言

0x1：機率論的起源

大數定律和中心極限定律是伴隨着古典統計，古典機率論，函數分析，極限理論，現代機率論與數理統計這些理論一塊兒發展起來的。科學家在提出和發展機率論理論體系的過程當中，一個最大也是最根本的挑戰就是，如何從理論上強有力地證實機率論是一個理論合理的數學理論，而不只僅是基於大量的實際實驗而獲得的概括假設，誠然，不少科學理論都是從科學實驗裏概括總結獲得的通常性推理，例如亞里士多德的重力理論，可是最終它們都必須被合理地嵌入在某個理論體系中，和其餘理論互相推導成立，才能真正成爲一門科學。html

在機率論發展的初期，人們主要仍是以實驗概括爲主，人們發現，當大量重複某一相同的實驗的時候，其最後的實驗結果可能會穩定在某一數值附近。就像拋硬幣同樣，當咱們不斷地拋，拋個上千次，甚至上萬次，咱們會發現，正面或者反面向上的次數都會接近一半。除了拋硬幣，現實中還有許許多多這樣的例子，像擲骰子，最著名的實驗就是泊松拋針實驗。這些實驗都像咱們傳達了一個共同的信息，那就是大量重複實驗最終的結果都會比較穩定。那穩定性究竟是什麼？怎樣去用數學語言把它表達出來？這其中會不會有某種規律性？是必然的仍是偶然的？ios

17 世紀中葉，人們開始對機會性遊戲的數學規律進行探討。它的發展與數學史上一些偉大的名字相聯繫，如帕斯卡、費馬、惠更斯、詹姆斯、伯努利、棣莫弗、拉普拉斯等。算法

1654 年，費馬與帕斯卡的通訊中關於分賭注問題的討論被公認爲是機率論誕生的標誌。問題是這樣的：併發

「兩個賭徒相約賭若干局，誰先贏 s 局就算贏了，當賭徒 A 贏 a 局 (a < s)，而賭徒 B 贏 b 局(b < s) 時，賭博被迫停止，應該怎樣分配賭注才合理？」app

在三年後，惠根斯亦用本身的方法解決了這一問題，並寫成了《論賭博中的計算》一書，這就是機率論最先的論著，他們三人提出的解法中，都首先涉及了數學指望（mathematical expectation）這一律念，並由此奠基了古典機率的基礎。dom

通過幾百年的發展，大數定律體系已經很完善了，也出現了更多更普遍的大數定律，例如切比雪夫大數定律，辛欽大數定律，泊松大數定律，馬爾科夫大數定律等等。正是這些數學家們的不斷研究，大數定律才得以如此迅速發展，才得以完善。機器學習

在這篇文章中，筆者將嘗試從時間線的維度，介紹幾個核心的大數定理和中心極限定律的發展脈輪，並討論其應用場景。函數

Relevant Link: 工具

https://baike.baidu.com/item/%E5%A4%A7%E6%95%B0%E5%AE%9A%E5%BE%8B

2. 伯努利大數定律（1713年）

天然科學的發展老是從簡入難，從物理、生活中的平常現象發現規律，並進行概括總結逐步發展起來的。post

歷史上首個極限定理是由瑞士數學家雅各布.伯努利（1654年-1705年），在其遺著中發表了機率論中的第一個極限定理——伯努利大數定理，即「在屢次重複獨立的試驗中，事件發生的頻率有愈來愈穩定的趨勢」。

這正是頻率穩定性的定理形式，它的出現意味着機率論由創建走向發展的階段。

0x1：數學公式定義

設 na 是n重伯努利試驗中，事件A發生的次數，p是事件A在每次事件中發生的機率，則對任意的ε>0，有：

或

在n重伯努利試驗中，設：

則 X1，X2，....，Xn 是相互獨立同分布的隨機變量，且Xi~B(1，p)，i=1，2，...，n，其中 p=P(A)。

由伯努利大數定律知，當n充分大時，有以下等價近似公式：

A發生的頻率

這是機率論第一個很是深入的認識，當大量相互獨立重複試驗中，能夠用某個事件A發生的頻率來近似每次試驗中事件A發生的機率。這就是伯努利大數定律的直觀意義。

當n充分大時，頻率與其機率能以任意接近的機率趨向於1。

所以實際中，只要試驗次數足夠多，能夠用頻率做爲機率的估計。

同時伯努利大數定律也解釋了機率存在的客觀意義，即爲何「大數次」重複試驗下，事件的機率是存在的，正是由於頻率的這種穩定性，咱們才意識到機率的存在，纔有了機率論這門學科。

Relevant Link:

https://baike.baidu.com/item/%E4%BC%AF%E5%8A%AA%E5%88%A9%E5%A4%A7%E6%95%B0%E5%AE%9A%E5%BE%8B/7922139

3. 棣莫弗-拉普拉斯中心極限定理（1733年）

1730 年，法國數學家棣莫弗（1677年-1754年）出版的著做《分析雜論》中包含了著名的棣莫弗─拉普拉斯定理。他使用正態分佈取估計n（很大）時拋擲硬幣出現正面次數的分佈，即二項分佈B(n，0.5)。這就是機率論中第二個基本極限定理的雛形。

將近80年後，拉普拉斯（1749年-1827年）在 1812 年出版的《機率的分析理論》中，首先明確地對機率做了古典的定義。他指出當n很大時，二項分佈B(n，p)(0<p<1)均可以用正態分佈逼近。因此後人稱之爲棣莫弗-拉普拉斯中心極限定理。

另外，拉普拉斯他又和數個數學家創建了關於正態分佈及最小二乘法的理論。

0x1：數學公式定義

定理核心公式以下：

設隨機變量序列X1，X2，..... 是一個獨立同分布的隨機變量序列（注意這是二項分佈的前提假設），且Xi ~ B(1，p)，i = 1，2，....，則對任意實數x，有：

仔細看這個公式，會發現幾點：

分子就是二項分佈隨機變量序列的累加和，減去隨機變量的指望np
分佈就是二項分佈的標準差
分子除分母的整個公式，能夠當作是隨機變量指望值標準化的過程
同時括號裏的<=符號，表達說標準化的指望偏離0原點的機率，也即偏差機率，這個偏離偏差機率是符合正態分佈的

德莫弗-拉普拉斯中心極限定理不只在機率論發展的早期起太重要的做用，並且在工程實踐中依然被大量使用，在大數據前提下，基於該定理的近似等價計算方法在複雜工程問題中很是方便。

0x2：由該定理獲得二項分佈的正態分佈近似解

由二項分佈的可加性知道，所以咱們有：

這就是二項分佈的原始機率計算公式，精確的機率值理論上是能夠精確算出的。

可是，實際問題中當n較大時，計算並不方便，這時候就須要藉助極限、近似等價這2個數學工具來幫助咱們簡化計算。

泊松定理曾經告訴咱們，當p ≤ 0.1時，能夠用泊松分佈做近似計算，這個近似等價其實已經很好了，可是這個假設仍是有些侷限，若是p比較大，則沒法繼續用泊松近似。

如今棣莫弗-拉普拉斯中心極限定理告訴咱們，也能夠用正態分佈做近似計算，它的優勢是不受「p≤0.1」的限制，只需n足夠大。

由德莫弗 -拉普拉斯中心極限定理推得，若是隨機變量Y~B(n，p)，那麼，當n較大時，有以下近似等價公式：

顯然，棣莫弗-拉普拉斯中心極限定理是列維-林德伯格中心極限定理的特例，由於服從二項分佈，因此這個定理又稱爲二項分佈的正態近似。

0x3：棣莫弗-拉普拉斯中心極限定理對伯努利大數定理的更細緻證實

前面由伯努利大數定理知，當n充分大時，能夠用做爲p的近似，可是至於近似程度如何，伯努利大數定律沒有給出。

中心極限定律對近似的程度進行了註釋：

中心極限定理對頻率趨近於機率的論證更爲細緻。

0x4：棣莫弗-拉普拉斯定理應用舉例

1. 題目

一本20萬字的長篇小說進行排版。假定每一個字被錯排的機率爲10-5。試求這本小說出版後發現有6個以上錯字的機率，假定各個字是否被錯排是相互獨立的。

2. 解題分析

顯然，這道題符合二項分佈，咱們應該用二項分佈公式來解題，可是注意這裏n有20w，是一個比較大的數字，所以咱們採用棣莫弗-拉普拉斯定義或者泊松定理來計算近似等價結果。

設錯字總數爲隨機變量X ，X~B (200000，10^-5 )，np=2，

因此有正態分佈近似等價結果：

以及泊松分佈近似等價結果（λ=2）：

Relevant Link:

https://en.wikipedia.org/wiki/De_Moivre%E2%80%93Laplace_theorem

4. 列維-林德伯格中心極限定理（1920年）

咱們將棣莫弗-拉普拉斯定理推廣到通常化的相互獨立同部分隨機變量中，即列維-林德伯格中心極限定理。列維（1886年-1971年）是法國數學家，對極限理論和隨機過程理論作出了傑出的貢獻，林德伯格（1876年-1932年）是芬蘭數學家因中心極限定理而聞名於世。

0x1：數學公式定義

設X1，X2，... 是任意一個獨立同分布的隨機變量序列，且均值和方差存在（隨機變量序列收斂），即：

則對任意一個x，-∞ < x < ∞，總有：

其中，Φ(x)是 N(0，1) 的正態分佈函數。

上式也被稱爲列維 -林德伯格(Levy-Lindberg)中心極限定理。

因爲：

所以，定理中的機率其實是的標準化隨機變量的分佈函數值，任意的機率分佈函數均可以進行標準化。

列維-林德伯格中心極限定理告訴咱們，不論X1，X2，... 原來服從什麼分佈，當n足夠大時，總能夠近似地認爲：

或者

在實際問題中，若n較大，能夠利用正態分佈近似求得機率：

這在工程實踐中特別有用。

0x2：公式的合理性分析

1. 微觀少數樣本的隨機構成宏觀統計意義上的有序

咱們在討論正態分佈的時候，曾經闡述過正態分佈存在的合理性。當一個變量，受到大量微小的，獨立的隨機因素（即X1，X2，... Xn）的累計影響時，這個變量通常服從正態分佈，也能夠反過來理解爲，正是由於這些大量微小的，獨立的隨機因素（即X1，X2，... Xn）的累計影響，致使最終從宏觀上呈現出正態分佈的結果。中心極限定理正是這種直觀經驗的嚴格數學表達。

定理的條件要求隨機變量相互獨立而且服從同一分佈。

這裏相互獨立意味着隨機變量之間不相互影響；

同分布是指每一個隨機變量在隨機變量序列的前n項部分和中的地位相同，也即每一個隨機變量對前n項部分和的影響都是微小的。

這就解釋了天然界中一些現象受到許多相互獨立且微小的隨機因素影響，總的影響就能夠看做服從或近似服從正態分佈。

2. 高爾頓釘板實驗

爲了更好地闡明「微觀少數樣本的隨機性」以及「宏觀大數統計意義上的正態有序」這2個概念，咱們來看一個歷史上著名的例子，高爾頓釘板實驗。

弗朗西斯·高爾頓（Francis Galton，1822年2月16日－1911年1月17日），英國科學家和探險家。他曾到西南非洲探險，因樹立功績而知名並被選爲英國皇家地理學會會員，三年後又入選英國皇家學會，晚年受封爲爵士。
他的學術研究興趣普遍，包括人類學、地理、數學、力學、氣象學、心理學、統計學等方面。他是查爾斯·達爾文的表弟，深受其進化論思想的影響，把該思想引入到人類研究。
他着重研究個別差別，從遺傳的角度研究個別差別造成的緣由，開創了優生學。他關於人類官能的研究開闢了個體心理和心理測驗研究的新途徑。
他在統計學方面也有突出的貢獻，高爾頓在1877年發表關於種子的研究結果，指出迴歸到平均值（regression toward the mean）現象的存在，這個概念與現代統計學中的「迴歸」並不相同，可是倒是迴歸一詞的起源。

下圖是高爾頓釘板實驗的示意圖，實驗過程以下：

每一圓圈點表示釘在板上的一顆釘子，它們彼此的距離均相等，上一層的每一顆的水平位置剛好位於下一層的兩顆正中間，總共有n排釘子，至關於作n次伯努利實驗。

從入口處放進一個直徑略小於兩顆釘子之間的距離的小圓玻璃球，當小圓球向降低落過程當中，碰到釘子後皆以1/2的機率向左或向右滾下，因而又碰到下一層釘子。如此繼續下去，直到滾到底板的一個格子內爲止。

把許許多多一樣大小的小球不斷從入口處放下，只要球的數目至關大，它們在底板將堆成近似於正態的密度函數圖形（即：中間高，兩頭低，呈左右對稱的古鐘型）。

來稍微分析一下小球下落的結果，小球堆積的形態取決於小球最終下落在底部隔板的位置的分佈。設隨機變量X爲「小球最終下落在底部隔板中的位置」，X的機率分佈由一系列的隨機變量序列組成，即：

Xi = -1：小球碰到第 i 排釘子後向左下落
Xi = 1：小球碰到第 i 排釘子後向右下落

顯然對單次實驗（投擲小球）來講，小球會落在哪一個方格Xi徹底是隨機的，幾乎看不出什麼規律，這就是所謂的「微觀少數樣本的隨機性」，這也是機率論的精髓思想，表面隨機的事物背後都蘊含着某種肯定的統計規律性。

從機率分佈的角度來看，顯然有，隨機變量X表示最終小球下落的結果，是由全部隨機變量序列累加獲得的，和的分佈計算很複雜。可是高爾頓通過試驗發現，隨着試驗次數的增長，小球的堆積形態呈現出正態分佈的形態。這就是所謂的「宏觀統計意義上的有序性」。

通過疊加，本來取值任意一點的可能性由相同變爲了向中心位置聚攏，這和咱們的直覺是多麼不一樣。

Relevant Link:

https://baike.baidu.com/item/%E9%AB%98%E5%B0%94%E9%A1%BF%E9%92%89%E6%9D%BF/6765470
https://baike.baidu.com/item/%E6%9E%97%E5%BE%B7%E4%BC%AF%E6%A0%BC%E5%88%97%E7%BB%B4%E5%AE%9A%E7%90%86/10388983?fr=aladdin

5. 辛欽大數定律（相互獨立同部分大數定律）

辛欽大數定律繼承並發展了伯努利大數定律。

辛欽大數定律面向的是通常化的相互獨立同分布的隨機變量狀況，並不侷限於伯努利隨機試驗。伯努利大數定律是相互獨立同分布大數定律的特例，這點咱們接下來會討論。

0x1：數學公式定義

則隨機變量序列 X1，X2，....，Xn 相互獨立同分布（任意兩個隨機變量線性無關），若均值存在，即

則對任意ε>0，有：

也能夠表示成：

在許多實際問題中，方差存在不必定知足，蘇聯數學家辛欽（1894年-1959年）證實了相互獨立同分布情形下，僅指望存在、方差不存在時結論仍然成立，所以相互獨立同分布大數定律又稱爲辛欽大數定律。

0x2：算數平均值法則的理論依據

相互獨立同分布（辛欽）大數定律是咱們平常生活中常用的算術平均值法則的理論依據。

爲了精確稱量物體的質量μ，可在相同條件下重複稱n次，結果可記爲x1，x2，....，xn。可看作n個相互獨立同分布的隨機變量（X1，X2，....，Xn）的一次觀測值。

X1，X2，....，Xn服從同一分佈，它們共同的指望記爲物體的真實質量μ。

由相互獨立同分布（辛欽）大數定律可知，當n充分大時有：

。

這意味着逐漸趨向於μ，也即隨機變量的算數平均值具備穩定性。

在物理實驗中咱們就是採用這種方法測得物體質量的，例如，爲獲得一顆鑽石的真實質量，咱們測n次取其算法平均值便可。

算法平均值法則提供了一條切實可操做的途徑來獲得物體的真實值。大數定律從理論上給出了這個結論的嚴格證實，而不是僅僅靠直覺。

Relevant Link:

https://baike.baidu.com/item/%E8%BE%9B%E9%92%A6%E5%A4%A7%E6%95%B0%E5%AE%9A%E5%BE%8B/7922188

6. 切比雪夫大數定律

切比雪夫大數定律是對伯努利大數定律和辛欽大數定律的繼承和發展。是由俄國數學家切比雪夫（1821年-1894年）發表的。

切比雪夫大數定律是面向全部兩兩不相關且方差一致有界的隨機變量，並不要求隨機變量獨立同分布。相互獨立同分布的大數定律是切比雪夫大數定律的特例，由於切比雪夫大數定律中的在相互獨立同分布大數定律的條件下爲μ。

切比雪夫大數定律比較理論化，在討論具體的數學公式以前，咱們須要先介紹幾個數學概念。

0x1：切比雪夫（Chebyshev）不等式

隨機變量X的取值老是圍繞着其指望變更，若X的分佈已知時，能夠計算事件的機率。

切比雪夫不等式給出了任意機率分佈X時，的上限計算公式。

設隨機變量X的數學指望E(X)及方差D(X)存在，則對於任意的ε>0，有：

事件理解爲：「隨機變量X關於其指望發生了較大誤差」，切比雪夫不等式給出了此事件的機率上界，它與方差成正比：

方差越大，此上界就越大
方差越小，X在其指望附近取值的密集程度就越高，那麼遠離指望的區域的機率上界就越小

咱們對切比雪夫公式進行一個變換使其更好地體現出「機率和誤差的反比關係」，

，令，則有，

從上式中，能夠很清楚看到，X偏離E[X]兩個標準差距離的機率小於1/4。

切比雪夫不等式進一步說明了方差的機率意義，即方差是隨機變量取值與其中心位置的偏離程度的一種度量指標。

值得注意的是，顯然利用切比雪夫不等式估計「隨機變量X關於其指望μ發生了較大誤差」的機率是粗糙的，引入切比雪夫不等式的另外一個目的是，它是證實大數定律的工具之一。

0x2：依機率收斂

在前面第二小節討論到伯努利大數定律的時候，曾經提到頻率的穩定性。

設隨機事件A的機率P(A)=p，在n重伯努利試驗中事件A發生的頻率爲fn(A)，當n很大時，fn(A)將與p很是接近。

天然會想到，應該用極限概念來描述這種穩定性，可是不能簡單地使用高等數學中數列的極限，由於fn(A)本質上是一個隨機變量，它隨着不一樣的n次試驗可能取不一樣的值，這就須要對隨機變量序列引進新的收斂性定義。

設X1，X2，... 是一個隨機變量序列，若是存在一個常數c，使得對任意一個ε>0，總有：

那麼，稱隨機變量序列X1，X2，... 依機率收斂於c，記做

依機率收斂性的直觀意義是，當n足夠大時，隨機變量Xn幾乎老是取接近於常數c的值。

利用求對立事件的機率計算公式，依機率收斂性也能夠等價地表示成：

通常地，Xi 不必定相互獨立，也不必定服從0-1分佈，把具備這種形式的依機率收斂性的結論統稱爲大數定律，即通常化的大數定律。

0x3：切比雪夫大數定律數學公式

設 X1，X2，... 是兩兩不相關（注意兩兩不相關並非獨立同分布，只是兩兩不線性相關）的隨機變量序列。若是存在常數c，使得 D(Xi)≤c，i=1，2，...，那麼有：

由切比雪夫不等式推得，對任意一個ε>0，當n → ∞時，有：

要特別注意的是，切比雪夫大數定律的主要條件是「方差有界」。

在特殊條件下，即隨機變量相互獨立同分布，即E(Xi) = μ，i=1，2，... ，則切比雪夫大數定律等價於辛欽大數定律，即：

0x4：三個大數定律之間的條件關係

三個大數定律條件是不一樣的：

切比雪夫大數定律不要求隨機變量序列同分布，甚至不要求相互獨立，只要兩兩不線性相關、方差一致有界便可
辛欽大數定律要求隨機變量相互獨立且同分布，但不要求方差存在，僅指望存在便可
努利大數定律要求隨機變量相互獨立且同分布，但限定於伯努利兩點分佈

7. 大數定律和中心極限定理的異同

0x1：中心極限定理和大數定律的區別

中心極限定理是隨機變量和的分佈收斂到正態分佈的一類定理，而隨機變量的和又和隨機變量的均值有密切的聯繫，而大數定律論證的主要部分就是隨機變量均值的收斂性特色，所以，中心極限定理和大數定律之間有千絲萬縷的聯繫。

大數定律是說，n只要愈來愈大，把這n個獨立同分布的數加起來去除以n獲得的這個樣本均值（也是一個隨機變量）會依機率收斂到真值u，可是樣本均值的分佈是怎樣的咱們不知道
中心極限定理是說，n只要愈來愈大，這n個數的樣本均值會趨近於正態分佈，而且這個正態分佈以u爲均值，sigma^2/n爲方差
綜上所述，這兩個定律都是在說樣本均值性質。隨着n增大，大數定律說樣本均值幾乎必然等於均值。中心極限定律說，它愈來愈趨近於正態分佈，而且這個正態分佈的方差愈來愈小

0x2：中心極限定理在機率論中的核心地位

不一樣的中心極限定理的差別就在於對隨機變量序列作出了不一樣的假設，因爲中心極限定理的有力支撐使正態分佈在機率論與數理統計中佔據了獨特的核心地位，它是20世紀初機率論研究的中心內容。

在大數n的狀況下，全部機率分佈函數都會收斂到正態分佈這個中心極限的形式上。同時各個機率分佈函數之間也存在互相推導的關係，這也是爲何GMM從理論上能夠擬合任意複雜的機率分佈函數的緣由。

每當這時，筆者就很是喜歡祭出一張很是出名且有趣的圖：

8. 大數定律的理論應用

大數定律支撐了不少應用理論的發展，例如：

1）算術平均法則
2）頻率估計機率方法
3）數理統計中參數的點估計思想，這也是不少傳統機器學習模型訓練的理論支撐
4）中心極限定理在數理統計的區間估計與假設檢驗問題中的應用

9. 大數定律存在的一個實驗例證

0x1：泊松隨機變量的收斂 - 在大數N時，大數定理廣泛存在的一個例證

下面是三個不一樣的泊松隨機變量在大數定理上的實例，用三個不一樣的曲線表示。三個泊松分佈的參數𝜆 = 4.5，樣本個數 sample_size=100000 個。

泊松分佈自己也是一個隨機過程，參數 𝜆 同樣，只是代表它們符合同一個機率分佈，可是泊松分佈自己每次的取值依然是一個隨機過程。這裏用三個泊松分佈來進行實驗，只是想說明無論過程如何隨機，可是其指望依然會收斂到一個具體的值上，就像有一隻魔法之手在牽引着每次的取值。

下面代碼計算𝑛 = 1到 sample_size 時樣本的均值：

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt

figsize(12.5, 5)
import pymc as pm

sample_size = 100000
expected_value = lambda_ = 4.5
poi = pm.rpoisson
N_samples = range(1, sample_size, 100)
print N_samples

for k in range(3):

    samples = poi(lambda_, size=sample_size)

    partial_average = [samples[:i].mean() for i in N_samples]   # 從[1:sample_size]不一樣樣本數下，三個泊松分佈的均值

    plt.plot(N_samples, partial_average, lw=1.5, label="average \
of  $n$ samples; seq. %d" % k)


plt.plot(N_samples, expected_value * np.ones_like(partial_average),
         ls="--", label="true expected value", c="k")

plt.ylim(4.35, 4.65)
plt.title("Convergence of the average of \n random variables to its \
expected value")
plt.ylabel("average of $n$ samples")
plt.xlabel("# of samples, $n$")
plt.legend()

plt.show()

從上圖中，咱們能夠清晰地看出，當 n 很小時，樣本均值的變更很大，隨着 n 的增大都在逐漸趨向平緩，最終在 4.5 水平線附近小幅擺動。數學家和統計學家將這種「小幅擺動」稱爲收斂。

咱們所熟知的定理：泊松分佈的均值等於其參數𝜆。這裏所謂的泊松分佈指的是一個真實機率分佈，可是真實機率分佈只有上帝知道咱們並不知道，可是當樣本數足夠大時，根據大數定理，咱們能夠用必定數量的樣本均值來近似代替真實分佈的均值，進而獲得真實分佈的指望。

1. 收斂到指望值的速度

明白了分佈的收斂特性以後，另外一個與此相關的問題是：收斂到指望值的速度有多快？

1）經過實驗方式來觀察收斂速度

咱們選擇一個特定的 n，而後重複實驗數千次，並計算與實際指望值的距離的均值：

以上公式能夠理解爲：前 n 個採樣的樣本均值與真實值的距離（平均意義上的）

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import pymc as pm

figsize(12.5, 4)

sample_size = 100000
expected_value = lambda_ = 4.5
poi = pm.rpoisson
N_samples = range(1, sample_size, 100)

N_Y = 250  # use this many to approximate D(N)
N_array = np.arange(1000, 50000, 2500)  # use this many samples in the approx. to the variance.
D_N_results = np.zeros(len(N_array))

def D_N(n):
    """
    This function approx. D_n, the average variance of using n samples.
    """
    Z = poi(lambda_, size=(n, N_Y))
    average_Z = Z.mean(axis=0)
    return np.sqrt(((average_Z - expected_value) ** 2).mean())


for i, n in enumerate(N_array):
    D_N_results[i] = D_N(n)


plt.xlabel("$N$")
plt.ylabel("expected squared-distance from true value")
plt.plot(N_array, D_N_results, lw=3,
         label="expected distance between\n\
expected value and \naverage of $N$ random variables.")
plt.plot(N_array, np.sqrt(expected_value) / np.sqrt(N_array), lw=2, ls="--",
         label=r"$\frac{\sqrt{\lambda}}{\sqrt{N}}$")
plt.legend()
plt.title("How 'fast' is the sample average converging? ")

plt.show()

如咱們預期，隨着 N 增大，樣本均值與實際值間距的指望逐漸減少。但也注意到，這個收斂速率也在下降。

樣本點從 10000 增長到 20000 時，均值和指望的距離從 0.020 下降到 0.015，減小了 0.005；當樣本點從 20000 增長到 40000 時，樣本均值和指望之間的距離從 0.015 下降到 0.010，僅僅減小了 0.005。說明收斂的速率在逐漸變小。

2）經過公式方式來描述收斂速率

實際上收斂的速率是能夠經過公式衡量的，在上圖的橙色曲線中，是經過下面的公式生成的：

該公式模擬了樣本的均值和真實指望之間的距離。

這個特性是很是有趣的，由於對於一個給定的𝑁，就能夠知道估計出的樣本平均值和真實的指望值之間的偏差大小。也很容易直觀地想象，若是當前方差很小了，說明擺動幅度已經很小了，基本上確定是快要進入收斂狀態了，因此收斂速率天然也減少了。

2. 鏈接機率和頻率之間的橋樑 - 大數定理

頻率統計 <---->(大數定理) <----> 指望 <----> 機率

在這裏，大數定理充當了一個橋樑的做用，將頻率統計和機率鏈接了起來。而鏈接的水泥就是指望。

同時，咱們也要明白，大數定理在 N 很小的時候會失效，在 N 很小的時候，此時獲得的近似結果是不可靠的。

0x2：小數據的無序性 - 當 N 很小時大數定理會失效

只有在 N 足夠大時，大數定理才成立，而後數據量並不是老是足夠大的。若是任意運用這必定理，就會致使很荒謬的錯誤，咱們這章來討論這個問題。

1. 實例1 - 地理數據的聚合

有些數據會以某一屬性進行分組。例若有些數據按照國家，縣區或者城市進行分組。固然人口數量也會根據地理區域的不一樣而不一樣。若是要對某一地理區域中的一些特徵進行統計平均時，必需要注意大數定理適用範圍，並瞭解在人口數較少時它是如何失效的。
假設在數據集中包含有 5000 個縣的數據。除此以外人口數量在整個國家的分佈符合 100 到 1500 的均勻分佈。咱們關心的是每一個縣中人口的平均高度。

咱們假設不論生活在哪裏，每一個人的身高都服從相同的分佈：height~Normal(150,15)

把數據按照縣級別進行分組，因此只能對每一個縣裏的數據作平均。數據以下圖：

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import pymc as pm

figsize(12.5, 4)
std_height = 15
mean_height = 150

n_counties = 5000
pop_generator = pm.rdiscrete_uniform
norm = pm.rnormal

# generate some artificial population numbers
population = pop_generator(100, 1500, size=n_counties)

average_across_county = np.zeros(n_counties)
for i in range(n_counties): # 5000個縣
    # generate some individuals and take the mean
    average_across_county[i] = norm(mean_height, 1. / std_height ** 2,  # 每一個縣都服從相同的分佈: height~Normal(150,15)
                                    size=population[i]).mean()　　# 不一樣的縣人口

# located the counties with the apparently most extreme average heights.
i_min = np.argmin(average_across_county)
i_max = np.argmax(average_across_county)

# plot population size vs. recorded average
plt.scatter(population, average_across_county, alpha=0.5, c="#7A68A6")
plt.scatter([population[i_min], population[i_max]],
            [average_across_county[i_min], average_across_county[i_max]],
            s=60, marker="o", facecolors="none",
            edgecolors="#A60628", linewidths=1.5,
            label="extreme heights")

plt.xlim(100, 1500)
plt.title("Average height vs. County Population")
plt.xlabel("County Population")
plt.ylabel("Average height in county")
plt.plot([100, 1500], [150, 150], color="k", label="true expected \
height", ls="--")
plt.legend(scatterpoints=1)

plt.show()

從上圖中能夠觀察到什麼現象？若是不考慮人口數量大小，推斷出的結果就可能存在錯誤。

例如，若是不考慮人口數大小，咱們可能會說有最高身高人口的縣也有最低身高的人口，反之亦然（不考慮人口數量，圖中的 x 軸就變成一個點了，紅色圓圈的點就會在一條豎線上了，因此就反映出最高人口的縣也會有最低人口，最低人口的縣也會有最高人口)。

經過這個例子，咱們彷佛得出了一個結論，小數據集狀況下大數定理失效了！可是這種推斷結果然的是對的嗎？

從貝葉斯不肯定性理論的角度出發，準確的答案應該是：不肯定！

理論上來講，確實存在可能在說圖中紅圈所在的 x 軸對應的縣確實存在極值的身高，可是在這種小數據狀況下，不肯定性是很是高的，高到咱們幾乎不能將其採納做爲一種推斷結果。換句話說，這種高風險的推斷對咱們來講是沒有意義的。

從大數定理的角度來講，這個縣不必定有最極端高度的人存在。錯誤的緣由是人口數量少時計算的均值並不能反映出人口真實的指望值(真實是 μ = 150)。樣本的數量 N 越小，大數定理就會失效的越厲害。

2. 實例2 - Kaggle的美國人口普查反饋比例預測

下面是 2010 年美國人口普查數據，它不是按照縣對數據分區，而是按照街區進行分區。

本例中的數據來自於 Kaggle 機器學習競賽中使用的數據集。競賽的目的是爲了預測各個街區對人口普查郵件的回覆率，回覆率的取值範圍爲 0 到 100，使用的預測特徵包括:人口平均收入，女性數量，停車場的個數，小孩的平均個數等。

下面畫出人口調查郵件的回覆率和街區中人口數量的關係圖（多維自變量不方便在一張圖上展現）

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import pymc as pm

figsize(12.5, 6.5)
data = np.genfromtxt("./data/census_data.csv", skip_header=1,
                     delimiter=",")
plt.scatter(data[:, 1], data[:, 0], alpha=0.5, c="#7A68A6")
plt.title("Census mail-back rate vs Population")
plt.ylabel("Mail-back rate")
plt.xlabel("population of block-group")
plt.xlim(-100, 15e3)
plt.ylim(-5, 105)

i_min = np.argmin(data[:, 0])
i_max = np.argmax(data[:, 0])

plt.scatter([data[i_min, 1], data[i_max, 1]],
            [data[i_min, 0], data[i_max, 0]],
            s=60, marker="o", facecolors="none",
            edgecolors="#A60628", linewidths=1.5,
            label="most extreme points")

plt.legend(scatterpoints=1)

plt.show()

上圖反映了典型的統計現象。這裏的典型指的是散點圖的形狀。散點圖是一個典型的三角形，隨着樣本數量的增長，這個三角形變得越緊實（由於愈來愈符合大數定理的要求）。

能夠這麼理解，隨着 N 不斷增大，後驗機率愈來愈收斂到某一個固定值，而在 N 很小的時候，後驗機率的擺動是很劇烈的。

3. 實例3 - 如何對Reddit網站上的評論進行排序

上面的兩個例子都旗幟鮮明地表達了一個觀點：大數據！big data！機率統計和機器學習只有在大數據時才能體現出威力，若是你沒有大數據，那就別癡心妄想作機器學習項目了。

可是世上之事，不如意者十之，咱們不可能永遠在項目中能拿到大數據，就算你身在BAT、google這樣的大數據公司，你所擁有的也不必定就是大數據。

參考Andrew Gelman的一句名言：樣本歷來都不是足夠大的。若是 N 太大不足以進行足夠精確的估計，你須要得到更多的數據。但當 N 「足夠大」，你能夠開始經過劃分數據研究更多的問題，例如在民意調查中，當你已經對全國的民意有了較好的估計，你能夠開始分性別、地域、年齡進行更多的統計。N 歷來都沒法作到足夠大，由於當它一旦大了，你老是能夠開始研究下一個問題從而須要更多的數據。

回到咱們這個小節的問題，考慮在線商品的評分問題，這在生活中很常見，各類商城類app裏都會有針對某個商品的歷史購買者評論。

你如何評價只有一我的，兩我的或者三我的對一個商品的五星級評價結果。咱們多少都會明白，這樣少的評分結果沒法真正的反映出商品的真實價值。這種評價信息不充分的狀況給排序帶來了困難。不少人都知道在線搜索圖書，視頻時，若是用它們的評分或者評論對返回結果進行排序，最終用戶獲得的結果並不會很好。一般狀況下搜索結果中排在最前面的視頻或評論每每有很好的評分，而這些好的評分都是由粉絲評出的。還有一些真正高質量的視頻或者評論排在了稍微靠後的位置，它們的評分在 4.8 左右。咱們如何對這種狀況進行修正？

1）咱們如何對評論進行排序

Reddit是一個流行的網站。它主要提供一些小說或者圖片的連接，這個網站最流行的部分是它對每一個連接的評論。Reddit能夠對每一個評論進行投票(被稱爲同意票或否決票)。默認狀態下Reddit將把最好的評論排在前面。

如今假設咱們是Reddit後臺的算法工程師，咱們會如何設計這種排序算法呢？

須要明白的是這是一個很是複雜且牽一髮而動全身的事情，一個簡單的例子是消費者公衆是具備從衆效應的，即大多數時候咱們會優先看大多數的意見，且很容易接受所謂的大多數人的意見。這種現象表如今不少方面，例如百度搜索的首頁的訪問量會遠大於後面的頁面，且這種差距一旦創建，會愈來愈大；一個商品的評論一旦上了淘寶的爆款熱搜，它的購買熱度會愈來愈大，由於公衆會優先購買它們最早看到的所謂爆款。

而評論排序這個問題也是同樣的，一旦一個評論被排到靠前的位置，用戶會優先看到這些評論，而有更大的概率接受該評論的觀點，從而造成從衆效應，最終致使靠前的評論的點贊數會愈來愈多。

基於這樣的複雜性緣由，咱們接下來對評論排序這件事展開多維度的討論。

熱度

認爲點贊數越多的評論越好。可是一個問題在於，當一評論有幾百個贊，可是被反對了上千次時，雖然熱度很高，可是把這條評論當成最佳評論是有爭議的。

可是這裏的問題也很是複雜，一個觀點有上百文贊同，有上千人反對，那麼這個觀點就必定是對，或者必定就是錯的嗎？都不必定，不少時候，這取決於下游的商業和市場策略。

差數

用點贊數減去反對數來打分。這解決了以上熱度指標的問題。可是依然有缺陷，它不能解決評論的時間屬性致使的問題。

由於一個連接在線上會存在不少，不斷有用戶來點贊或反對，更老的評論會隨着時間積累更多的點贊數，但並不意味着他們是最佳評論。這個問題的根源顯然是沒有進行歸一化。

可是反過來想，這個問題也十分複雜，爲何必定要歸一化呢？試想，一個評論若是在一個至關長的時間區間內，持續不斷地被其餘用戶所接受並點贊，這不就意味着這個評論的觀點接受度很高嗎？

時間調節

考慮用差數除以評論的壽命（至關於歸一化）。這樣獲得了一種相似每秒差數或每分鐘差數這樣的速率值。

看起來足夠合理了，可是一樣存在缺陷，好比一個一秒前發佈的評論，只須要一票，就能擊敗一百秒前的有用九十九票的評論。出現這種問題的緣由是沒有把持續時間的跨度考慮進去。

避免這一問題的一種解法是隻針對 t 秒前的評論進行統計（過濾掉剛剛出現的新鮮的評論）。可是 t 怎麼選呢？並且這樣是否是意味着 t 秒內的評論都是很差的呢？

讀者看出來了嗎？咱們是在對不穩定量和穩定量進行比較（新評論 PK 老評論）

好評率

對某一評論的好評數除以對該評論的好評和差評的總和，獲得好評率，做爲排名權重。這解決了 3 中的問題，若是用好評率做爲排名結果，即便一個新的評論也會和老的評論同樣獲得較高的排名得分。

可是依然存在缺陷，在這種度量標準下，若是對一個評論的好評個數只有一個，即好評率爲 1.0，而另外一個評論的好評個數爲 999，差評個數爲 1，此時這個評論的好評率爲 0.999，顯而後一個評論的質量要比前一個高，這「可能」不太客觀。

我說「可能」是有理由的，前者雖然只有一個同意票，但仍是有可能因爲有 999 個同意票的後者，只是說這個可能性有多少罷了，咱們不能武斷地下結論。

2）真實的好評率如何獲得？ - 貝葉斯推斷來幫忙

咱們真正目標就是估計真實的好評率。注意，真實的好評率並非觀測到的好頻率，真實的好評率隱藏起來了，咱們能觀測到的僅僅是每一個評論的好評和差評。

當有 999 個好評，1 個差評時，咱們確實能夠說真實的好評率接近於 1，咱們能這樣說的緣由是由於有大數定理的存在；

反之，當只有一次好評時，咱們沒法斷言真實的好評率就是 1，由於在小數據下，後驗估計時不肯定的，這聽起來像個貝葉斯問題，咱們顯然沒法獲得一個對好評率的準確估計，咱們只能對好評率可能的機率分佈進行一個貝葉斯推斷。

貝葉斯方法須要好評率的先驗知識，得到先驗知識的一個方法是基於好評率的歷史分佈。能夠用好評率的歷史分佈做爲定好評率的先驗分佈。只要對Reddit的評論進行挖掘就能夠獲得好評率的歷史分佈。計算好評率的歷史分佈時還應注意如下問題：

1. 傾斜數據: 大部分的評論都沒有投票數據，所以有許多評論的好評率都很是小 (極端值)，這就使分佈大多數集中在極端值附件。能夠爲投票的個數設個閾值，大於 該閾值的狀況才予以考慮。閾值的選擇還要考慮其對好評率精度的影響。
2. 誤差數據: Reddit有許多不一樣的子頁面，稱爲subreddits。例若有兩個頁面，第一個爲:r/aww 頁面粘貼許多可愛動物的照片，第二個爲:r/politics 討論政治問題。用戶對這兩個頁面的評論行爲是不 相同的。瀏覽第一個頁面的人的性格會更友好，更有愛心，也就會有更多的好評(與第二個頁面相比)，而第二個頁面的評論可能會有更多的爭論或不一樣的觀點出現。所以並非全部的評論都是一致的。

鑑於此，我以爲用均勻分佈做爲評論好評率的先驗分佈是最佳的作法。

假設共有𝑁個投票數，好評率爲𝑝，好評的個數就是二項隨機變量（好評或差評），參數爲 𝑝 和𝑁。構造一個函數對 𝑝 進行貝葉斯推斷。

import pymc as pm


def posterior_upvote_ratio(upvotes, downvotes, samples=20000):
    """
    This function accepts the number of upvotes and downvotes a particular submission received, 
    and the number of posterior samples to return to the user. Assumes a uniform prior.
    """
    N = upvotes + downvotes
    upvote_ratio = pm.Uniform("upvote_ratio", 0, 1)
    observations = pm.Binomial("obs", N, upvote_ratio, value=upvotes, observed=True)
    # do the fitting; first do a MAP as it is cheap and useful.
    map_ = pm.MAP([upvote_ratio, observations]).fit()
    mcmc = pm.MCMC([upvote_ratio, observations])
    mcmc.sample(samples, samples / 4)
    return mcmc.trace("upvote_ratio")[:]

下面是後驗分佈的結果：

figsize(11., 8)
posteriors = []
colours = ["#348ABD", "#A60628", "#7A68A6", "#467821", "#CF4457"]
for i in range(len(submissions)):
    j = submissions[i]
    posteriors.append(posterior_upvote_ratio(votes[j, 0], votes[j, 1]))
    plt.hist(posteriors[i], bins=18, normed=True, alpha=.9,
             histtype="step", color=colours[i % 5], lw=3,
             label='(%d up:%d down)\n%s...' % (votes[j, 0], votes[j, 1], contents[j][:50]))
    plt.hist(posteriors[i], bins=18, normed=True, alpha=.2,
             histtype="stepfilled", color=colours[i], lw=3, )

plt.legend(loc="upper left")
plt.xlim(0, 1)
plt.title("Posterior distributions of upvote ratios on different submissions")

從上圖中能夠看出，某些分佈很窄，其餘分佈表現爲長尾。這體現了咱們對真實好評率的不肯定性。

3）如何針對評論的後驗分佈的進行排序

回到咱們這個小節的主題，咱們的目標是：如何對評論進行從好到壞的排序？

固然，咱們是沒法對分佈進行排序的，排序只能是標量值。有不少種方法可以從分佈中提取出標量值，用指望或均值來表示分佈就是一個方法，可是均值並非一個好辦法，由於它沒有考慮到分佈的不肯定性。

咱們這裏建議使用 95%最小可信值，定義爲真實參數只有 5% 的可能性低於該值（即 95% 置信區間的下界），下圖是根據 95%最小可信值畫的後驗分佈。

N = posteriors[0].shape[0]
lower_limits = []

for i in range(len(submissions)):
    j = submissions[i]
    plt.hist(posteriors[i], bins=20, normed=True, alpha=.9,
             histtype="step", color=colours[i], lw=3,
             label='(%d up:%d down)\n%s...' % (votes[j, 0], votes[j, 1], contents[j][:50]))
    plt.hist(posteriors[i], bins=20, normed=True, alpha=.2,
             histtype="stepfilled", color=colours[i], lw=3, )
    v = np.sort(posteriors[i])[int(0.05 * N)]
    # plt.vlines( v, 0, 15 , color = "k", alpha = 1, linewidths=3 )
    plt.vlines(v, 0, 10, color=colours[i], linestyles="--", linewidths=3)
    lower_limits.append(v)
    plt.legend(loc="upper left")

plt.legend(loc="upper left")
plt.title("Posterior distributions of upvote ratios on different submissions");
order = np.argsort(-np.array(lower_limits))
print(order, lower_limits)

爲什麼基於這個量進行排序是個好主意？

由於基於 95% 最小可信值進行排序，是對最佳評論的一個最爲保守的估計。即使在最差狀況下，也就是說咱們明顯高估了好評率時，也能確保最佳的評論排在頂端。在這一排序思路中，咱們利用瞭如下很天然的特性：

1. 給定兩個好評率相同的評論，咱們會選擇票數最多的做爲更加評論，由於票數更多的評論在貝葉斯後驗機率分佈中更集中，落在 95%置信區間的機率也更大，也即咱們更肯定其好評率更好。
2. 給定兩個票數同樣的評論，咱們會選擇好評數更多的。

10. 大數定律在工程領域的具體應用 - 保險行業的賺錢原理

0x1：保險行業的業務背景

保險是對風險的保障，它提供一份幫助人們分散危險、分攤損失的機制，這就是保險的本質。

本質是投保人經過不一樣數額的保費，將不一樣類型的風險打包轉賣給了保險公司，因此保險費是創建保險基金的來源。

0x2：保費制定的理論支撐 - 大數定律

大數定律是保險業保險費計算的科學理論基礎，當投保人（保單）的數量足夠大時，由切比雪夫大數定律知，被保人繳納的純保費與其能得到的賠款的指望值是相同的。這個結論反過來，也能夠說明保險公司應該如何合理的收取保費。

假設有n個投保人購買了n份相互獨立的保險，每一個人出險的機率爲p，每一個人可能得到的賠償金爲a，應繳納的保費爲b。

當n足夠大時，設購買該險種的投保人出險的人數爲X，由已知得 X~B(n，p)。則：

nb = a * E(X) = a * np

=> b = ap

即保費的制定是每一個人可能得到賠償金的指望。

11. 咱們怎樣更好應用大數定理

儘管大數定理很是酷，可是正如其名，它只適用於足夠大的數據量。咱們已經看到，若是不考慮數據的構造，那麼估計結果會受到很大影響。

1. 經過（低成本地）得到大量後驗樣本，能夠確保大數定理適用於指望的近似估計
2. 在樣本數量很小的貝葉斯推斷中，後驗分佈包括更多的隨機性。從後驗分佈中將能看出這些隨機性的存在。後驗分佈越平坦隨機性越大，越陡峭隨機性越小。所以推斷的結果是須要調校的。
3. 不考慮樣本量會帶來很大的影響，此時排序的依據每每是不穩定的，並會致使很是病態的排序結果。