python实现大数定律与中心极限定理

一. 大数定律

设 { X k , k = 1 , 2 , . . . } {\{X_k,k = 1, 2,...\}} {Xk​,k=1,2,...}是独立同分布的随机变量序列，且随机变量$X$的数学期望为$\mu$，方差为${\delta }^2$。则对于任意给定的$ϵ>0$有
$$\underset{n\to \infty }{\lim }P(|\overline{X_i}-u|\ge ϵ)=\underset{n\to \infty }{\lim }P(|\frac{1}{n}(X_1+X_2+...+X_n)-\mu |\ge ϵ)=0$$

大数定律表明随机变量序列的算术平均值依概率收敛于随机变量的期望。

假如随机变量$X$服从二项分布，随机取出一个长度为$n$的随机变量序列$X_k$，随着$n$的增大，该随机变量序列$X_k$的均值的极限值就是随机变量$X$的数学期望。下面用代码来演示这个收敛过程。

import matplotlib.pyplot as plt
from scipy.stats import binom
import numpy as np

plt.rcParams["font.family"] = "SimHei"  # 设置字体
plt.rcParams["axes.unicode_minus"] = False  # 正常显示负号

if __name__ == '__main__':
    fig, ax = plt.subplots(figsize=(10,8))
    # 二项分布的参数
    n, p = 10, 0.4
    # 二项分布的期望
    expected_mean = n * p
    # 样本大小从10到100000,每次递增50
    sample_size_list = range(10, 100000, 50)
    mean_list = []
    for sample_size in sample_size_list:
        sample = binom(n = n, p = p).rvs(size = sample_size)
        mean = np.mean(sample)
        mean_list.append(mean)
    ax.plot(sample_size_list, mean_list, label = 'sample_mean')
    expected_mean_list = [expected_mean for i in sample_size_list]
    ax.plot(sample_size_list, expected_mean_list, label = 'expected_mean')

    ax.legend()
    plt.show()

运行结果：

从图中可以看出随着样本量的增大，样本的均值逐渐收敛于总体期望。

二. 中心极限定理

设 { X k , k = 1 , 2 , . . . } \{X_k, k=1,2,...\} {Xk​,k=1,2,...}是独立同分布的随机变量序列，且随机变量$X$的数学期望为$\mu$，方差为${\delta }^2\ne 0$，则当$n\to \infty$时，随机变量$X^{\ast }=\frac{\overline{X_k}-\mu }{\delta /\sqrt{n}}$分布服从标准正太分布，$n$为随机变量序列的长度，数学表达式如下：

$$\underset{n\to \infty }{\lim }P(\frac{\overline{X_k}-\mu }{\delta /\sqrt{n}}<\alpha )=\Phi (\alpha )=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\alpha }{e}^{-t^2/2}dt$$

假如随机变量$X$服从二项分布，则随机取出一个长度为$n$的随机变量序列$X_k$，随着抽样次数的增大，随机变量序列$X^{\ast }=\frac{\overline{X_k}-\mu }{\delta /\sqrt{n}}$的分布会越来越趋近于标准正太分布。下面用代码来演示这个过程。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom
from scipy.stats import norm
import math

plt.rcParams["font.family"] = "SimHei"  # 设置字体
plt.rcParams["axes.unicode_minus"] = False  # 正常显示负号


if __name__ == '__main__':
    fig, ax = plt.subplots(4, 1, figsize=(10, 8))

    # 二项分布的参数
    n, p = 10, 0.4
    # 二项分布的期望
    expected_mean = n * p
    # 二项分布的标准差
    expected_std = math.sqrt(n * p * (1 - p))
    # 单次抽样的大小
    sample_size = 1000
    # 抽样次数
    sample_times_list = [500, 10000, 100000]

    i = 0
    while i < len(sample_times_list):
        sample_mean_list = []
        for k in range(sample_times_list[i]):
            sample = binom(n=n, p=p).rvs(size = sample_size)
            # 计算抽样样本的均值
            sample_mean = np.mean(sample)
            sample_mean_list.append(sample_mean)

        # 对均值进行标准化处理
        sample_size_list_normalized = [ (mean - expected_mean) / (expected_std / math.sqrt(sample_size)) for mean in sample_mean_list]
        # 画出抽样均值的分布
        ax[i].hist(sample_size_list_normalized, bins=100, density=True, color='blue', label='抽样次数:{}, 单次抽样大小:{}'.format(sample_times_list[i], sample_size))
        ax[i].legend()
        i += 1

    # 画出正太分布
    x = np.linspace(-5, 5, 1000)
    y = norm(loc = 0, scale= 1).pdf(x)
    ax[i].plot(x,y, color = 'red', label = '标准正太分布')
    ax[i].legend()

    plt.show()

运行结果：

从图中可以看出随着抽样次数的增加，随机变量$X^{\ast }=\frac{\overline{X_k}-\mu }{\delta /\sqrt{n}}$的分布越来越跟标准正太分布吻合。