几率统计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位的”