概率統計18——再看大數定律

　　在對不了解概率的人解釋期望時，我總是敷衍地將期望解釋為均值。這種敷衍的說法之所以行得通，正是由於大數定律起了作用。

　　

　　人們在實踐中發現，盡管每個隨機變量的取值不同，但當隨機變量大量出現時，它們的均值卻相對恆定，這個規律就是大數定律。

一個公平的骰子

　　我們有一個公平的骰子，每個點數出現的概率都是1/6，如果只投擲一次，完全無法預測它的點數，但是如果把連續投擲20次看作一次試驗，卻發現每次試驗的點數的均值總是在3.5附近徘徊。每次試驗投擲的次數越多，點數的均值越穩定。

import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(10, 5))
plt.subplots_adjust(hspace=0.5, wspace=0.3) 

m = 100 # 試驗的次數
for i in range(1, 5):
    n = 2 * 10 ** i # 每次試驗投擲n次骰子
    mean_list = [np.random.randint(1, 7, n).mean() for i in range(1, m + 1)] # 每次試驗的均值
    ax = fig.add_subplot(2, 2, i)
    ax.plot(list(range(1, m + 1)), mean_list, label='均值')
    plt.yticks(list(range(0, 7))) # 重置y軸的坐標
    ax.set_xlabel('試驗次數')
    ax.set_ylabel('均值')
    ax.set_title('每次試驗投{}次骰子'.format(n))
    ax.legend(loc='upper right')

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標簽
plt.show()

　　

均值的期望和均值的方差

　　我們用隨機變量X₁表示第一次投骰子的結果，X₂表示第二次……，一個公平的骰子可以得到下面的結論：

　　連續投擲后將形成一個隨機變量序列{ X₁, X₂, ……, X_n}，這個序列的均值是：

　　序列中的每個變量都是隨機的，因此它們的和及均值也是隨機的，也就是說然是隨機的。以n=20為例：

　　而均值的期望則不同：

　　每次擲骰子的結果都是獨立同分布的，它們有相同的期望和方差。用μ和σ²表示隨機變量的期望和方差，於是我們得到了均值的期望：

　　均值的期望等於整體的期望，是一個定值。對於公平的骰子來說，整體的期望與整體的均值相等。

　　均值的方差是：

辛欽大數定律和伯努利大數定律

　　n個樣本均值的方差是總體方差的1/n，隨着n的增大，該方差也越來越小，逐漸趨近於0。方差為0表示隨機變量相對於的波動0，即不含隨機性。換句話說，隨着n的增大，樣本的均值逐漸收斂於μ，這就是所謂的大數定律。這回有意思了，隨着n的增大，作為均值的隨機變量居然逐漸變得不含隨機性。

　　這里的大數定律也稱為辛欽大數定律或弱大數定律。抄一下教科書上的定義：

　　設隨機變量X₁, X₂, …, X_n,…相互獨立，服從同一分布且具有數學期望E(X_k)=μ(k=1,2,…)，則序列 依概率收斂於μ，即。

　　

　　樣本均值的標准差是，標准差表示隨機變量與總體之間的差異度。想要把差異度縮減10倍，那么n的數量要增加10²倍。

　　

　　伯努利大數定律是辛欽大數定律的一個重要推論。

　　設f_A是n次獨立重復試驗中事件A發生的次數，p是事件A在每次試驗中發生的概率，則對於任意正數ε > 0，有：

　　別被定義唬住了，它想解釋的問題很簡單：當n足夠大時，f_A/n與P的差值趨近於0，也就是頻率趨近於概率，這也是我們用比例估計作為點估計量的基礎。在實際問題中，當試驗次數很大時，便可以用事件的頻率來代替事件的概率。

　　其實辛欽大數定律也可以用①的方式表示：

　　對於概率的符號，有時候有大括號，有時候用小括號，這個還是別太糾結，實際上沒什么區別，用大括號只是為了強調事件是一個集合。

不存在的期望

　　大數定律並不是對所有問題都生效，比如在不存在期望的情況下。

　　現在我們把骰子換成一枚公平的硬幣，參加一個關於硬幣第幾次正面朝上的游戲。把隨機變量X=x看作在第x次投硬幣時會得到第1次正面朝上的結果：

　　如果第1次就得到正面朝上的結果（X=1），得2分

　　如果第2次才得到正面朝上的結果（X=2），得4分

　　如果第3次才得到正面朝上的結果（X=3），得8分

　　……

　　如果第n次才得到正面朝上的結果（X=n），得2ⁿ分

　　現在我們需要計算得分Y=2^X的期望。顯然X符合幾何分布，每個X又和Y一一對應，因此Y出現的概率等同於X出現的概率：

　　

　　這個期望是發散的，即沒法得出具體的數值，也沒法求極限，E[Y]=∞表示不存在期望。

蒙特卡羅積分法

　　假設我們遇到了一個難以計算的積分：

　　包括換元法、三角替換、分部積分在內的各種方法都無濟於事，此時可以求助於暴力的大數定律。

　　我們知道積分的幾何意義是曲線與x軸圍成的面積：

　　曲線f(x)把矩形分割成兩部分，U是曲線上方的面積，L是陰影部分的面積，R是整個矩形的面積，是已知量。現在我從們這個矩形內隨機選擇一點r_i=(x, y)，那么根據幾何概型，這個點落入L中的概率是L/R。既然積分的是在計算面積L，我們就可以借助積分對落在L中的某一點的概率進行精確的表達：

　　上式同時附帶得出積分的結果：

　　R是已知量，只要求得p，就能得到積分的結果。於是積分問題變成了概率問題。

　　現在讓每個隨機點r_i都對應一個隨機變量X_i，X_i滿足下面的特性：

　　每個隨機變量的期望是：

　　將大數定律應用於X_i上，當n→∞時，X_i的均值依概率收斂於p，即：

　　n是落在矩形中的隨機點的數量，如果隨機點落在L中，X_i=1，否則X_i=0，這個結果實際是告訴我們，當n足夠大時，落在L中的點的概率趨近於落在L中點的數量與落在矩形中的總點數的比值。

　　將這個結果代入積分式：

　　於是我們可以借助物理實驗和計算機的模擬求解這個難以對付的積分。

　　

　　我們用程序模擬蒙特卡羅積分法，計算下面四分之一圓的面積：

　　圓的曲線是 x² + y² = 1，如果一個隨機點(x_i, y_i)滿足x² + y² < 1，則表示該點落在了圓內。

 

import matplotlib.pyplot as plt
import numpy as np

def create_ri(n):
    ''' 成 n 個隨機點 '''
    xs, ys = np.random.random_sample(n), np.random.random_sample(n)  # 隨機點的橫縱坐標
    return xs, ys

def cnt_circle(xs, ys):
    ''' 圓中隨機點的數量 '''
    ps = xs ** 2 + ys ** 2
    return ps[ps < 1].size

fig = plt.figure()
plt.subplots_adjust(hspace=0.5)
ax = fig.add_subplot(2, 1, 1)
x = y = np.arange(0, 1, 0.001)
x, y = np.meshgrid(x,y)
ax.contour(x, y, x**2 + y**2, [1])
xs, ys = create_ri(100)
ax.scatter(xs, ys, s=10)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.axis('scaled')

xl, yl = [], []
n_max = 10000
for n in range(100, n_max, 10):
    xs, ys = create_ri(n)
    cnt = cnt_circle(xs, ys)
    Area = cnt / n
    xl.append(n)
    yl.append(Area)
    n *= 10
ax = fig.add_subplot(2, 1, 2)
ax.hlines(np.pi / 4, xmin=-1000, xmax=n_max + 1000, label='$\pi/4$')
ax.plot(xl, yl, label='count($r_i$)/n')
ax.legend(loc='upper right')
ax.set_xlabel('n')
ax.set_ylabel('Area')
plt.show()

 

 

 

　　隨着n的增加，藍色曲線的波動越來越小，計算的結果越來越接近真實值。

　　正像前面提到的，如果想讓計算結果的精度提升10倍，需要增加100倍的試驗次數，因此蒙特卡羅法的效率並不高。

　　盡管大數定律下的蒙特卡羅法很好用，但仍然需要小心謹慎，關鍵之處在於大數定律沒有為我們明確指出到底什么才是“大數”，到底需要多少次試驗才能充分接近我們所期待的極限。無論n有多大，我們仍然不能否認存在這樣的情況：所拋出的硬幣全部正面朝上，盡管這種情況發生的概率很小。確定這個概率的方式是中心極限定理的內容。

---

　　出處：微信公眾號 "我是8位的"

　　本文以學習、研究和分享為主，如需轉載，請聯系本人，標明作者和出處，非商業用途！ 

　　掃描二維碼關注作者公眾號“我是8位的”

ps = xs ** 2 + ys ** 2
return ps[ps < 1].size