機率統計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):
    ''' 圓中隨機點的數量 '''
    return len(list(filter(lambda x: x < 1, xs**2 + ys**2)))

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位的」