1. 大數定理是什麼？
大數定理是一個非常底層的基礎性原理，大量的機器學習理論和演算法實際上都建立在其之上。只是因為其太底層，也太合理的，合理到我們感覺這就是一個毫無爭議的真理而不需要去單獨討論。我們本文和讀者來一起從貝葉斯推斷的角度來嘗試對大數定理展開一些討論。
0x1：大數定理基本數學定義
設?? 為服從特定概率分佈的?個獨立取樣樣本，按照大數定理理論，只要?的期望?[?]是有限的，則就有如下關係：
用文字來描述：即來自同一分佈的一組隨機變數，其均值收斂域該分佈的期望。
以上等式只有在極限條件下成立，但是隨著參加平均計算的樣本越來越多，均值將越來越接近?的期望值?[?]。除了一些極少數的特殊情況外，大數定理對於大多數分佈都是正確的。
0x2：泊松隨機變數的收斂 - 大數定理普遍存在的一個例證
下面是三個不同的泊松隨機變數在大數定理上的例項，用三個不同的曲線表示。 三個泊松分佈的引數? = 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）通過公式方式來描述收斂速率
實際上收斂的速率是可以通過公式衡量的，在上圖的橙色曲線中，是通過下面的公式生成的：
該公式模擬了樣本的均值和真實期望之間的距離。
這個特性是非常有趣的，因為對於一個給定的?，就可以知道估計出的樣本平均值和真實的期望值之間的誤差大小。也很容易直觀地想象，如果當前方差很小了，說明擺動幅度已經很小了，基本上肯定是快要進入收斂狀態了，所以收斂速率自然也減小了。
0x3：連線概率和頻率之間的橋樑 - 大數定理
在開始這個小節討論之前，我們先丟擲一個公式：
頻率統計 <---->(大數定理) <----> 期望 <----> 概率
在期望值和估計的概率之間存在一種隱性關係。
首先定義如下函式：
，這種0-1分佈有非常好的特性，可以很清晰地表達出期望和概率的關係。
根據大數定理，如果有樣本?? ，可以估計出事件 A 的概率，用?(?)表示
由上面只是函式的定義可知，只有事件發生時1?(?)才取 1 其它情況都取 0， 所以只要用試驗中事件發生的總次數除以試驗的次數就是概率（回想一下平時我們是怎 麼樣用頻率近似概率的）
在這裡，大數定理充當了一個橋樑的作用，將頻率統計和概率連線了起來。而連線的水泥就是期望。
同時，我們也要明白，大數定理在 N 很小的時候會失效，在 N 很小的時候，此時得到的近似結果是不可靠的。
2. 小資料的無序性 - 當 N 很小時大數定理的表現如何？
只有在 N 足夠大時，大數定理才成立，然後資料量並非總是足夠大的。如果任意運用這一定理，就會導致很荒謬的錯誤，我們這章來討論這個問題。
0x1：例項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 軸就變成一個點了，紅色圓圈的點就會在一條豎線上了，所以就反映出最高人口的縣也會有最低人口，最低人口的縣也會有最高人口)。
1. 在小資料集情況下大數定理失效了嗎？
但是這種推斷結果真的是對的嗎？
從貝葉斯不確定性理論的角度出發，準確的答案應該是：不確定！
理論上來說，確實存在可能在說圖中紅圈所在的 x 軸對應的縣確實存在極值的身高，但是在這種小資料情況下，不確定性是非常高的，高到我們幾乎不能將其採納作為一種推斷結果。換句話說，這種高風險的推斷對我們來說是沒有意義的。
從大數定理的角度來說，這個縣不一定有最極端高度的人存在。錯誤的原因是人口數量少時計算的均值並不能反映出人口真實的期望值(真實是 μ = 150)。樣本的數量 N 越小，大數定理就會失效的越厲害。
0x2：例項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 很小的時候，後驗概率的擺動是很劇烈的。
0x3：例項3 - 如何對Reddit網站上的評論進行排序
上面的兩個例子都旗幟鮮明地表達了一個觀點：大資料！big data！概率統計和機器學習只有在大資料時才能體現出威力，如果你沒有大資料，那就別痴心妄想做機器學習專案了。
但是世上之事，不如意者十之，我們不可能永遠在專案中能拿到大資料，就算你身在BAT、google這樣的大資料公司，你所擁有的也不一定就是大資料。
參考Andrew Gelman的一句名言：樣本從來都不是足夠大的。如果 N 太大不足以進行足夠精確的估計，你需要獲得更多的資料。但當 N “足夠大”，你可以開始通過劃分資料研究更多的問題，例如在民意調查中，當你已經對全國的民意有了較好的估計，你可以開始分性別、地域、年齡進行更多的統計。N 從來都無法做到足夠大，因為當它一旦大了，你總是可以開始研究下一個問題從而需要更多的資料。
回到我們這個小節的問題，考慮線上商品的評分問題，這在生活中很常見，各種商城類app裡都會有針對某個商品的歷史購買者評論。
你如何評價只有一個人，兩個人或者三個人對一個商品的五星級評價結果。我們多少都會明白，這樣少的評分結果無法真正的反映出商品的真實價值。這種評價資訊不充分的情況給排序帶來了困難。很多人都知道線上搜尋圖書， 視訊時，如果用它們的評分或者評論對返回結果進行排序，終端使用者得到的結果並不會很好。通常情況下搜尋結果中排在最前面的視訊或評論往往有很好的評分，而這些好的評分都是由粉絲評出的。還有一些真正高質量的視訊或者評論排在了稍微靠後的位置， 它們的評分在 4.8 左右。我們如何對這種情況進行修正？
1. 我們如何對評論進行排序
Reddit是一個流行的網站。它主要提供一些小說或者圖片的連結，這個網站最流行的部分是它對每個連結的評論。Reddit可以對每個評論進行投票(被稱為贊成票或否決票)。預設狀態下Reddit將把最好的評論排在前面。
現在假設我們是Reddit後臺的演算法工程師，我們會如何設計這種排序演算法呢？
需要明白的是這是一個非常複雜且牽一髮而動全身的事情，一個簡單的例子是消費者公眾是具有從眾效應的，即大多數時候我們會優先看大多數的意見，且很容易接受所謂的大多數人的意見。這種現象表現在很多方面，例如百度搜尋的首頁的訪問量會遠大於後面的頁面，且這種差距一旦建立，會越來越大；一個商品的評論一旦上了淘寶的爆款熱搜，它的購買熱度會越來越大，因為公眾會優先購買它們最先看到的所謂爆款。
而評論排序這個問題也是一樣的，一旦一個評論被排到靠前的位置，使用者會優先看到這些評論，而有更大的機率接受該評論的觀點，從而形成從眾效應，最終導致靠前的評論的點贊數會越來越多。
基於這樣的複雜性原因，我們接下來對評論排序這件事展開多維度的討論。
1）熱度
認為點贊數越多的評論越好。但是一個問題在於，當一評論有幾百個贊，但是被反對了上千次時，雖然熱度很高，但是把這條評論當成最佳評論是有爭議的。
但是這裡的問題也非常複雜，一個觀點有上百文贊同，有上千人反對，那麼這個觀點就一定是對，或者一定就是錯的嗎？都不一定，很多時候，這取決於下游的商業和市場策略。
2）差數
用點贊數減去反對數來打分。這解決了以上熱度指標的問題。但是依然有缺陷，它不能解決評論的時間屬性導致的問題。
因為一個連結線上上會存在很多，不斷有使用者來點贊或反對，更老的評論會隨著時間積累更多的點贊數，但並不意味著他們是最佳評論。這個問題的根源顯然是沒有進行歸一化。
但是反過來想，這個問題也十分複雜，為什麼一定要歸一化呢？試想，一個評論如果在一個相當長的時間區間內，持續不斷地被其他使用者所接受並點贊，這不就意味著這個評論的觀點接受度很高嗎？
3）時間調節
考慮用差數除以評論的壽命（相當於歸一化）。這樣得到了一種類似每秒差數或每分鐘差數這樣的速率值。
看起來足夠合理了，但是同樣存在缺陷，比如一個一秒前釋出的評論，只需要一票，就能擊敗一百秒前的有用九十九票的評論。出現這種問題的原因是沒有把持續時間的跨度考慮進去。
避免這一問題的一種解法是隻針對 t 秒前的評論進行統計（過濾掉剛剛出現的新鮮的評論）。但是 t 怎麼選呢？而且這樣是不是意味著 t 秒內的評論都是不好的呢？
讀者看出來了嗎？我們是在對不穩定量和穩定量進行比較（新評論 PK 老評論）
4）好評率
對某一評論的好評數除以對該評論的好評和差評的總和，得到好評率， 作為排名權重。這解決了 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. 給定兩個票數一樣的評論，我們會選擇好評數更多的。
3. 大數定理的一些更多思考
儘管大數定理非常酷，但是正如其名，它只適用於足夠大的資料量。我們已經看到，如果不考慮資料的構造，那麼估計結果會受到很大影響。
1. 通過（低成本地）獲得大量後驗樣本，可以確保大數定理適用於期望的近似估計 2. 在樣本數量很小的貝葉斯推斷中，後驗分佈包括更多的隨機性。從後驗分佈中將能看出這些隨機性的存在。後驗分佈越平坦隨機性越大，越陡峭隨機性越小。因此推斷的結果是需要調校的。 3. 不考慮樣本量會帶來很大的影響，此時排序的依據往往是不穩定的，並會導致非常病態的排序結果。