小文
公众号 小文的数据之旅在前面的两期中，我们介绍了如何对统计数据进行整理与描述，使得我们对数据的概貌有一个简单的了解，然而简单的描述只能实现对统计数据粗浅的利用，要挖掘出数据中隐含的信息还需要用推断统计的方法，其中涉及到一些概率、概率分布等等。比如说一个班级中平均身高为166.5，标准差为14.3，共有50个人，那么随机抽取一个人，身高在170以上的概率有多大呢？本文包括的知识点如下：概率的相关概念以及运算离散型随机变量及其分布：伯努利分布、二项分布、泊松分布连续型随机变量及其分布：正态分布1、概率1.1 概率的相关概念古典概率：源于赌博，其基本的前提是游戏遵循公平规则，而且所有基本结果发生的概率相同，比如说一个均匀的骰子，每个点数的出现概率是相等的。相对频率：多次重复做一个试验，那么其概率就是事件最终发生的次数比率。比如扔一个硬币，共扔了5次，其中有3次正面，2次反面，那么其正面发生的概率为60%。根据大数法则，相对频率最终会收敛于其发生的真实概率。个别概率：生活中大部分事件不会重演，个别概率是个人对某种结果发生的主观估计。比如说昨天下了一整天大雨，个人认为下雨的概率大于50%时，那么我今天就会带伞出门。条件概率：是指事件A在另外一事件B已经发生条件下的发生概率。比如说昨天下了一整天雨为事件A，在此事件A发生的基础上今天下雨(事件B)的概率，表示为P(B|A)。1.2 概率的运算一般的加法法则：P(E \cup F) = P(E)+P(F)-P(E \cap F)当E与F相互排斥时，特殊的加法法则： P(E \cup F) = P(E)+P(F) 减法法则：
P(E) = 1-P(\bar{E}) 一般的乘法法则： P(E \cap F) = P(E|F)P(F) 当E与F相互独立时，特殊的乘法法则： P(E \cap F) = P(E)P(F)，即 P(E|F) = P(E) 贝叶斯法则： P(E|F) = \frac{P(E)P(F|E)}{P(E)P(F|E)+P(\bar{E})P(F|\bar{E})}=\frac{P(E \cap F)}{P(F)} 2、 一试验两定理伯努利试验：如果一个随机的试验只包含两个可能结果，第一种可能结果的发生概率是p，另一种是1-p，则该试验称为伯努利试验。抛硬币就是最为经典的伯努利试验。大数法则：当样本数量越多，则其平均就越趋近期望值，它描述的是样本均值收敛到总体均值。用数学的语言表达就是： \lim_{n \rightarrow \infty}{P(|\frac{S_{n}}{n}-p|<\epsilon)} =1
其中n为试验次数，p为每次试验成功的概率，Sn为总成功次数， \epsilon 代表一个任意小的正数。中心极限定理：当样本数量足够大的时候，样本分布符合 X\sim N(μ,σ^2) ，即正态分布，它描述的是当样本数量足够大时，样本均值的分布慢慢变成正态分布。3、 随机变量分类及其分布3.1 随机变量分类如果随机变量X的所有取值都可以逐个列举出来，则称X为离散型随机变量，其分布主要有两种：二项分布和泊松分布，常写作P(X=x)。如果随机变量X的所有取值无法逐个列举出来，而是取数轴上某一区间的任一点，则称X为连续型随机变量，其分布主要有：正态分布，常写作P(X<=x)。问：为什么离散型随机变量写作P(X=x)，而连续型随机变量写作P(X<=x)？这个其实很好理解，因为连续随机变量有无穷多个取值，那么单个取值的概率趋向于0，这样的概率对我们来说没有意义，因此对于连续随机变量，往往会计算一个区间内的概率分布，记作P(X<=x)。常用概率密度函数描述连续随机变量的概率分布。3.2 期望E与方差D期望：表示随机变量本身的平均水平或集中程度，计算方式类似于加权算术平均数，即每个随机变量取值乘以其对应的概率, E(X)=\sum_{i=1}^{n}{x_{i}p_{i}}。 方差：表示随机变量取值的离散程度，D(X)=E[X-E(X)]^{2}=E(X^{2})-[E(x)]^{2}，后者比前者计算更为便捷。这里推导一下：因为 \sum_{i=1}^{n}p_{i}=1 ，设 \bar{u}=E(X) ,那么 期望与方差的一些运算（1）线性变换的通用公式，这个很好证明，这里就不证明了，只要记住基础概率没变，只是随机变量进行了线性变化。E(aX+b) = aE(X)+b D(aX+b) = a^{2}D(X) （2）如果X1,X2,...Xn是X的独立观测值，那么每个观测的期望值、方差值一致，则E(X1+X2+X3+...+Xn)=nE(X) D(X1+X2+X3+...+Xn)=nD(X) （3）如果X,Y为独立随机变量，则E(X \pm Y) = E(X) \pm E(Y) D(X \pm Y) = D(X) + D(Y) 无论随机变量做加法还是减法，其方差都是增大的！3.3 几何分布应用条件：进行一系列独立试验，每一次试验或成功或失败，每一次试验的成功概率相同。解决的问题：为了取得第一次成功，需要进行多少次试验。若每次试验中事件A发生的概率为p，试验进行到事件A出现时停止，此时所进行的试验次数为k，那么发生的概率： P(X=k) = p(1-p)^{k-1}期望： E(X) = \frac {1}{p} ，即希望 \frac {1}{p} 次成功。方差： D(X) = \frac {1-p}{p^{2}}从图中可以看到当k等于1时，P(X=k)达到最大值，随着k增大，P(X=k)逐渐下降。结合几何分布的期望，我们也可以得到相同的结果，当p=100%时，仅需尝试一次即可成功！当p趋向于0时，k趋向于正无穷，即基本不会成功。3.4 二项分布应用条件：进行一系列次数有限的独立试验，每一次试验或成功或失败，每一次试验的成功概率相同。解决的问题：在n次试验中能成功多少次。若重复进行n次试验，事件A发生了k次且每次试验中事件A发生的概率为p，那么：发生的概率为： P(X=k) =C_{n}^{k}p^{k}(1-p)^{n-k} 期望： E(X) = np 方差： D(X) = np(1-p) 二项分布的形状根据n与p的大小会发生变化，p越接近0.5，图形越对称。看上面的图是不是觉得有点像正态分布的形状呢？没错，在满足一定条件下，可以用正态分布近似二项分布！当n=1时,也称为伯努利分布。 这里可能就有朋友会问，人家二项分布好好地有自己的计算公式，为什么要用正态分布来近似替代？两个字，简单！因为正态分布的计算过程比二项分布更为简单！3.5 泊松分布应用条件：单个事件在给定区间内随机、独立地发生，已知给定区间内的事件平均发生次数，或者叫发生率，且这个发生次数是有限的。 解决的问题：在给定区间内的事件发生次数。比如说某企业每月发生事故的次数，某仪器每月出现故障的次数、保险公司每天理赔的次数等等，设λ为给定的时间内时间的事件发生的平均次数，k为事件发生的次数。其发生的概率为： P(X=k)=\frac{\lambda^{k}e^{-\lambda}}{k!} 期望： E(X) = \lambda 方差： D(X) = \lambda λ的大小会影响泊松分布的形状，λ小时向右偏斜，随着λ增大，分布逐渐变得对称，即当λ大到一定程度时，可以用正态分布近似替代！另外看上图是不是跟二项分布的形状也有几分类似？没错，当 P\leq0.25,n>20,np\leq5 时，可以用泊松分布（计算更简单）近似二项分布的效果良好。3.6 正态分布连续随机变量中最为重要的概率分布为正态分布，上述说到二项分布与泊松分布时，都提到了正态分布，正态分布就那么牛逼吗？它之所以被称为正态，是因为它的形态看起来合乎理想，是自然界中的 “理想”分布。正态分布具有钟形曲线，曲线对称，均值和中位数均位于中央，具有最大概率密度，越是偏离均值，概率密度减小。记作 X\sim N(\mu,\sigma^{2}) 。正态分布的概率可以通过两种方式就得，分别是计算分布曲线下方的面积和借助概率表进行查找。其概率密度函数： f(x) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}期望： E(X) = u 方差： D(X) = \sigma^{2}通过计算分布曲线下方的面积可以求出概率，但是计算难度很大。不过幸运的是我们还可以借助概率表进行查找。 步骤如下：（1） 确定随机变量的分布与范围（2） 标准化（转为均值为0，标准差为1的标准正态分布）（3） 借助概率表查找概率 这里举个例子说明一下，有个女生高64英寸，穿上5英寸高跟鞋去约会，男生的身高分布X~N(71,20.25)，那么她的约会者比她高的概率是多少？答：（1） 确定分布与范围：女生穿上高跟鞋后身高为69英寸，男生身高分布为X~N(71,20.25)，即求P(X>69)。（2） 标准化： Z = \frac{(x-\mu)}{\sigma}=\frac {69-71}{\sqrt{20.25}} = -0.44(保留两位小数)（3）通过查找概率表，得到P(Z <-0.44)=0.3300，那么 P(X>69) = P(Z>-0.44) = 1-P(Z<-0.44) = 1-0.3300 = 0.67 正态分布的应用在很多场景中，都会通过小概率帮助我们做决策，正态分布在质量管理中的应用比较有名的就是所谓的6西格玛的奇迹了，下图描述了一个标准正态分布。在标准正态分布中，约有99.7%数值分布在距离平均值有3个标准差之内的范围，也就是说仅有0.3%的数值分布在3σ之外，那么分布在6σ之外的概率有多小呢？大概为十亿分之二。这小到什么概念？通常人们认为奇迹发生的概率为百万分之一，也就是说分布在6σ之外的概率比奇迹还奇迹。但是在大量的试验中发现，实际生产过程中均值可能会发生1.5σ的偏移，所以我们通常说的6σ实际上指的是4.5σ，也就是有名的百万分之三点四，发生的概率也与奇迹相当了，在现实的生活中有理由认为概率低于百万分之三点四的事情是不会发生的。正态分布近似替代二项分布、泊松分布在一些特定情况下，可以用正态分布近似替代二项分布。如果X~B(n,p)，且np>5,n(1-p)>5，则可以用X~N(np,np(1-p))近似代替X。但是要注意一个问题，前面我们说到二项分布是离散变量，而正态分布是连续变量，因此用正态分布近似二项分布时，需要做连续性修正。原因如下：因此，一般用P(X<b+0.5) 近似替代 P(X\leq b) P(a-0.5 <X<b+0.5) 近似替代 P(a\leq X \leq b) P( X>a-0.5) 近似替代 P(X\geq a) 随着λ变大，泊松分布图的外形看起来越来越像正态分布。当λ大于15时，就可以通过正态分布近似替代泊松分布，如X~ Po(λ)且λ>15，就可以用X~N(λ，λ)近似替代，同样别忘了连续性修正！— end —小文的数据之旅戳右上角「+关注」获取最新share如果喜欢，请分享or点赞 }

这篇回答节选自我在专栏《机器学习中的数学：概率统计》中的一篇文章，我们来谈一下大数定律。也欢迎关注我的知乎账号 @石溪 ，将持续发布机器学习数学基础及算法应用等方面的精彩内容。1.一个背景话题我们首先从一个大家都非常熟悉的场景切入来展开我们本节内容的讨论。比如说，我们想获得本省15岁男生的平均身高，这时你会怎么做？显然你不会也不可能真的去统计全省所有15岁男生的身高，然后再求平均值，这样做不太现实。因此，你会去找一些样本，也就是找一部分本省15岁的男生，取他们身高的平均值，用这个样本的平均值去近似的估计所有15岁男生的平均身高。没错，一般就是这么干的，那接下来我再问你，找100个样本取得的平均值和1000个样本所取得的平均值，哪一个你认为更有可能接近真实的全省男生的平均身高（也就是期望）呢？你会说应该是1000个吧，毕竟样本数量多，上下偏差相互抵消，应该会更接近一些。你的直觉没有错。在数据分析的应用中，经常会有上述类似的应用场景，我们需要分析一类对象，常常得去获取他的关键参数，就比如上面所提到的全体男生身高的均值，但是现实中我们不可能去穷尽全部的研究对象，而是只能取得一部分的样本，通过计算这部分样本的参数值去近似的估计总体的目标参数，样本数量越大，近似效果越好。这里的理论依据就是我们下面要详细讲解的大数定理，大数定理是一个非常底层的基础性原理，大量的机器学习理论和算法实际上都建立在这个基础之上。我们常常是理所当然的直觉上感受到他的存在，却很少仔细想过背后的原因。那么通过这篇回答，我们会深入透彻的理解大数定理、中心极限定理背后的极限思想。具体内容安排如下：首先：我们介绍大数定理的原理，并用python语言对其进行模拟，给大家一个更直观的感受；接着：我们会介绍中心极限定理以及他在工程实践中的应用价值和意义，同时也会用python语言来对其进行模拟；最后：我们会讲解大数定理在机器学习和数据分析中的一个重要应用：蒙特卡罗方法，并仔细剖析他的应用场景和内涵。2.大数定理2.1.原理介绍好了，下面我们开始正式进入到大数定理的内容中。我们有如下的随机变量：X_1,X_2,...,X_n 他们彼此之间满足独立同分布，因此他们拥有相同的均值 \mu 和方差 \sigma^2 。此时，我们重点来研究这一组随机变量的均值：M_n=\frac{X_1+X_2+...+X_n}{n} ，显然 M_n 也是一个随机变量。那么，M_n 的期望和方差就是我们此时重点关心的问题。首先，我们从期望的定义入手，来观察一下随机变量M_n 的期望E[M_n] ：E[M_n]=E[ \frac{X_1+X_2+...+X_n}{n}]=\frac{1}{n}(E[X_1]+E[X_2]+...+E[X_n])=\frac{1}{n}\cdot n\cdot \mu=\mu=E[X_i]不难发现，一组独立同分布随机变量均值的期望就等于随机变量的期望，这个结论很直观。下面我们再来看看 M_n 的方差var[M_n] ：var[M_n]=var[ \frac{X_1+X_2+...+X_n}{n}]=\frac{1}{n^2}var[X_1+X_2+...+X_n]=\frac{1}{n^2}(var[X_1]+var[X_2]+...+var[X_n])=\frac{1}{n^2}\cdot n\cdot \sigma^2=\frac{\sigma^2}{n}我们从推导中发现，n 个独立同分布随机变量的均值的方差，是单一随机变量方差的\frac{1}{n}。没错，均值的方差变小了，并且随机变量 X的个数 n 越多，方差越小，他们的分布更加紧密的围绕在了期望的周围。特别的，当 n \rightarrow \infty 时，随机变量均值的方差趋近于0： var[M_n]=\frac{\sigma^2}{n} \rightarrow 0结合前前后后的一大段推导和论述，我们可以得出这么一个结论：独立同分布的随机变量 X_1,X_2,...,X_n ，他们的均值 M_n 的分布会更加接近于实际分布的均值 \mu ，随着样本量 n 的增大，他逐渐收敛于 \mu ，当 n \rightarrow \infty 时，也就是说当样本量非常大的时候，通过抽样样本计算所得到的平均值可以说就是 E[X] 了。独立同分布的随机变量序列的样本均值，在大样本的情况下，以很大的概率与随机变量的均值非常接近。这也就是为什么说，当独立同分布的样本数量 n 充分大时，样本均值（频率）是概率 P 的一个非常好的估计。这就回到本节最初提到的那个小问题了，样本数量到底是选100还是选1000，相信大家都会有明确的理论支撑了。这里我们还简单的提两个结论，其实说起来他们都比较直观。2.2.两个重要的不等关系谈到这里，我们来看概率统计中的两个非常重要的不等关系，以及他们背后的直观逻辑。第一个是马尔科夫不等式。用简洁直白的语言描述他就是：对于一个非负的随机变量 X ，如果他的均值很小的话，那么这个随机变量取到一个大值的概率是非常小的。想想确实很直观，因为随机变量的分布大部分都集中在均值附近，越远离均值，概率就越小。描述这个现象，有一个专门的不等式，叫做马尔科夫不等式：P(|M_n-\mu
\ge \epsilon) \le \frac{\sigma^2}{n \epsilon^2}第二个是切比雪夫不等式。我们同样简单的描述一下这个不等式的内涵：如果一个随机变量的方差非常小的话，那么这个随机变量取到远离均值 \mu 的概率也是非常小的，这个说实话也非常直观，同样有一个专门的不等式来描述他，叫切比雪夫不等式：P(|X-\mu|\ge c) \le \frac{\sigma^2}{c^2}2.3.大数定理的模拟下面我们来实际模拟一下大数定理。在第一个模拟的例子中，我们给大家一个感性的认识：我们生成3组各15000个服从参数为(10,0.4)的二项分布随机变量，随机变量的期望为 n*p=4 ，然后观察随着样本数目的增大，样本均值和实际分布期望之间的关系。代码片段：import numpy as np
from scipy.stats import binom
import matplotlib.pyplot as plt
import seaborn
seaborn.set()
n = 10
p = 0.4
sample_size = 15000
expected_value = n*p
N_samples = range(1, sample_size, 10)
for k in range(3):
binom_rv = binom(n=n, p=p)
X = binom_rv.rvs(size=sample_size)
sample_average = [X[:i].mean() for i in N_samples]
plt.plot(N_samples, sample_average,
label='average of sample {}'.format(k))
plt.plot(N_samples, expected_value * np.ones_like(sample_average),
ls='--', label='true expected value:np={}'.format(n*p), c='k')
plt.legend()
plt.show()运行结果：图1.大数定理模拟我们设置了三个相同的试验组，从试验结果中我们可以发现，在每一组试验中，随着样本数量的逐渐增大，样本均值都会越来越收敛于随机变量的期望。接下来我们再来看看第二个用于模拟大数定理的例子，我们从大数定理的定义出发，我们先生成1000000个服从均值为0，标准差为20正态分布的样本。依次进行三种不同的处理，并观察对应的三组分布图像：图像1：原始正态分布的样本分布图像，颜色为蓝色；图像2：从1000000个原始正态分布样本中，每次随机选取5个数，计算他们的均值，重复操作10000次，观察这10000个均值的分布，颜色为红色；图像3：从1000000个原始正态分布样本中，每次随机选取50个数，计算他们的均值，重复操作10000次，观察这10000个均值的分布，颜色为绿色；代码片段：import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import seaborn
seaborn.set()
norm_rvs = norm(loc=0, scale=20).rvs(size=1000000)
plt.hist(norm_rvs, normed=True, alpha=0.3, color='b', bins=100, label='original')
mean_array = []
for i in range(10000):
sample = np.random.choice(norm_rvs, size=5, replace=False)
mean_array.append(np.mean(sample))
plt.hist(mean_array, normed=True, alpha=0.3, color='r', bins=100, label='sample size=5')
for i in range(10000):
sample = np.random.choice(norm_rvs, size=50, replace=False)
mean_array.append(np.mean(sample))
plt.hist(mean_array, normed=True, alpha=0.3, color='g', bins=100, label='sample size=50')
plt.gca().axes.set_xlim(-60, 60)
plt.legend(loc='best')
plt.show()运行结果：图2.大数定理的模拟这个程序采样的规模比较大，可能运行的时间会比较长。从图中我们发现，随着每次选取的样本数量的增多，样本均值分布的图像越来越向期望集中，再一次佐证了大数定理。3.中心极限定理3.1.原理介绍下面接着来看另外一个现象，我们还是获取随机变量序列： X_1,X_2,...,X_n ，这 n 个随机变量满足独立同分布，均值为 \mu ，方差为 \sigma^2 。我们在这组随机变量序列的基础之上得到一个新的随机变量：Z_n=\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma}针对随机变量 Z_n ，我们很容易计算出：E[Z_n]=E[\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma}]=\frac{E[X_1+X_2+...+X_n]-n\mu}{\sqrt{n}\sigma}=0var[Z_n]=var[\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma}]=\frac{var[X_1]+var[X_2]+...+var[X_n]}{n\sigma^2}=1即，随机变量 Z_n 的期望为0，方差为1。关键的一点是，随着样本个数 n 的增大，随机变量Z_n 的分布逐渐趋向于一个标准正态分布，当 n \rightarrow \infty 时，随机变量的分布收敛于一个标准正态分布。更为重要的一点是，这个定理对随机变量 X 的原始分布没有任何要求，非常具有一般性。3.2.中心极限定理的工程意义实际上，中心极限定理中的随机变量Z_n=\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma} 是经过标准化处理的，如果单单只考虑 n 个随机变量的和，我们很容易得到：S_n=X_1+X_2+...+X_n 的分布趋近于一个均值为 n\mu ，方差为 n\sigma ^2 的正态分布。中心极限定理的意义在于，大量样本的独立随机因素的叠加是趋近于一个正态分布的，这一点在很多工程领域很常见也很关键。更重要的一点是，他不需要我们去搞明白随机变量X 的分布列或者概率密度函数，这往往是非常复杂的，甚至根本就无从得知，我们只需要知道他的均值和方差就可以进行后续的处理和分析了。3.3.中心极限定理的模拟下面，我们还是来举个例子，模拟验证一下中心极限定理。我们从一个服从参数p=0.3的几何分布中进行采样，共分三组试验，分别每次采样2个、5个、50个样本，每组试验各重复100000次，然后按照Z_n=\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma} 进行标准化，得到三组试验对应的结果，最后对试验结果进行可视化观察：代码片段：import numpy as np
from scipy.stats import geom
import matplotlib.pyplot as plt
import seaborn
seaborn.set()
fig, ax = plt.subplots(2, 2)
geom_rv = geom(p=0.3)
geom_rvs = geom_rv.rvs(size=1000000)
mean, var, skew, kurt = geom_rv.stats(moments='mvsk')
ax[0, 0].hist(geom_rvs, bins=100, normed=True)
ax[0, 0].set_title('geom distribution:p=0.3')
n_array = [0, 2, 5, 50]
for i in range(1, 4):
Z_array = []
n = n_array[i]
for j in range(100000):
sample = np.random.choice(geom_rvs, n)
Z_array.append((sum(sample) - n * mean) / np.sqrt(n * var))
ax[i//2, i%2].hist(Z_array, bins=100, normed=True)
ax[i//2, i%2].set_title('n={}'.format(n))
ax[i//2, i%2].set_xlim(-3, 3)
plt.show()运行结果：图3.中心极限定理的模拟左上第一幅图是几何分布的原始图像，我们发现，随着单次采样个数的逐渐增加，随机变量Z_n=\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma} 的分布图像越来越趋近于一个标准正态分布。4.大数定理的应用：蒙特卡罗方法4.1.机器学习中的应用背景用大样本数据计算出来的频率去估计概率，这就是大数定理的本质，而大数定理思想的一个非常典型的应用就是蒙特卡罗方法。蒙特卡罗方法，又叫统计模拟方法，名字很洋气，思想很粗暴，真的很管用。它使用随机数来进行场景的模拟或者过程的仿真，其思想核心就是通过模拟出来的大量样本集或者随机过程去近似我们想要研究的实际问题对象，这是一类非常重要的数值计算方法。该方法的名字来源于世界著名的赌城蒙特卡罗。赌博和概率，二者相视一笑、不谋而合。这种方法最初应用于20世纪40年代美国的曼哈顿原子弹计划，如今在数据分析和机器学习领域中到处都有他的身影。以下是蒙特卡罗方法的几类非常典型的应用：1.近似计算不规则面积/体积/积分2.模拟随机过程，预测随机过程可能性结果的区间范围3.利用马尔科夫链-蒙特卡罗方法（MCMC）进行未知参数的统计推断要知道，蒙特卡洛方法不仅仅只是一种方法技巧，更是一种思考问题的方式。在这篇回答中，我们主要介绍第一个应用点，让大家进一步理解大数定理的原理和蒙特卡罗方法的要点。第二个和第三个应用如果大家有兴趣可以关注我的专栏。4.2.利用蒙特卡罗方法计算不规则面积蒙特卡罗方法可以近似计算出不规则图形的面积，特别对于那些难以用解析方法计算的图像，非常有效。这里，我们利用蒙特卡罗方法来近似计算一个圆的面积，然后估计出 \pi 的近似值，选择圆作为例子的原因不是因为圆无法通过解析法进行计算，而是因为大家对他的面积比较熟悉，方便我们进行结果的对比。代码片段：import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from scipy.stats import uniform
import seaborn
seaborn.set()
n = 100000
r = 1.0
o_x, o_y = (0., 0.)
uniform_x = uniform(o_x-r,2*r).rvs(n)
uniform_y = uniform(o_y-r,2*r).rvs(n)
d_array = np.sqrt((uniform_x - o_x) ** 2 + (uniform_y - o_y) ** 2)
res = sum(np.where(d_array < r, 1, 0))
pi = (res / n) /(r**2) * (2*r)**2
fig, ax = plt.subplots(1, 1)
ax.plot(uniform_x, uniform_y, 'ro', markersize=0.3)
plt.axis('equal')
circle = Circle(xy=(o_x, o_y), radius=r, alpha=0.5)
ax.add_patch(circle)
print('pi={}'.format(pi))
plt.show()运行结果：pi=3.14096图4.蒙特卡洛方法近似计算圆面积我们结合这张图，详细的分析一下这段程序。我们近似计算的目标就是这个蓝色的半径为r=1的单位圆面积，而这个单位圆的外接正方形的边长为l=2r=2，因此外接正方形的面积为4。我们生成100000个在外接正方形内均匀分布的点，均匀的撒下去，这里面的核心原理就是：\frac{圆形面积}{正方向面积}\approx \frac{圆内点的个数}{点的总个数}这样就可以估算出单位圆的面积了。而为了估算\pi ，我们回到这个例子中的参数，可以得到这么一个公式：\frac{\pi r^2}{(2r)^2}\approx \frac{圆内点的个数}{点的总个数}\rightarrow \pi \approx 4 * \frac{圆内点的个数}{点的总个数}而且大数定理告诉我们，随着样本数量的增大，我们用这种方式模拟出来的 \pi 值应该是越来越趋近于真实值，样本无穷大的时候收敛于真值。这就证明了应用大数定理的蒙特卡罗方法的合理性和有效性。至于有同学会问了，是怎么判定一个点是否位于圆当中？这里用的是距离法，也就是计算每个点到原点的距离，如果小于等于半径，就说明这个点在圆内。这里是不是有一个bug？就是说这个规则只对圆这种特殊图像有效，如果是一个完全不规则的图像，我们如何来计算这个图像的面积？我们可以借助计算机图像中的像素来进行。比如把不规则图像内部都涂黑，图像外部都留白，我们还是均匀的撒下大样本量的点，如果某个点位于的坐标，他的像素是黑色的，则证明这个点在图像内部，反之就在图像外部，利用这个方法就能统计出位于图像内部的点的个数。5.小结在这一篇回答中，我们从样本均值和总体期望的关系入手，详细的讨论了概率统计中的重要概念：大数定理与中心极限定理，目的是让大家能够建立好概率统计中极限思维的概念基础。在实践环节中，我们介绍了大数定理的重要应用：蒙特卡罗方法，列举了他在机器学习中的几类重要应用，希望能有助于大家强化大数定理的思想方法。此内容节选自我的专栏《机器学习中的数学：概率统计》，前三小节免费试读，欢迎订阅当然还有《机器学习中的数学（全集）》系列专栏，欢迎大家阅读，配合食用，效果更佳~有订阅的问题可咨询微信：zhangyumeng0422}

这篇回答节选自我在专栏《机器学习中的数学：概率统计》中的一篇文章，我们来谈一下大数定律。也欢迎关注我的知乎账号 @石溪 ，将持续发布机器学习数学基础及算法应用等方面的精彩内容。1.一个背景话题我们首先从一个大家都非常熟悉的场景切入来展开我们本节内容的讨论。比如说，我们想获得本省15岁男生的平均身高，这时你会怎么做？显然你不会也不可能真的去统计全省所有15岁男生的身高，然后再求平均值，这样做不太现实。因此，你会去找一些样本，也就是找一部分本省15岁的男生，取他们身高的平均值，用这个样本的平均值去近似的估计所有15岁男生的平均身高。没错，一般就是这么干的，那接下来我再问你，找100个样本取得的平均值和1000个样本所取得的平均值，哪一个你认为更有可能接近真实的全省男生的平均身高（也就是期望）呢？你会说应该是1000个吧，毕竟样本数量多，上下偏差相互抵消，应该会更接近一些。你的直觉没有错。在数据分析的应用中，经常会有上述类似的应用场景，我们需要分析一类对象，常常得去获取他的关键参数，就比如上面所提到的全体男生身高的均值，但是现实中我们不可能去穷尽全部的研究对象，而是只能取得一部分的样本，通过计算这部分样本的参数值去近似的估计总体的目标参数，样本数量越大，近似效果越好。这里的理论依据就是我们下面要详细讲解的大数定理，大数定理是一个非常底层的基础性原理，大量的机器学习理论和算法实际上都建立在这个基础之上。我们常常是理所当然的直觉上感受到他的存在，却很少仔细想过背后的原因。那么通过这篇回答，我们会深入透彻的理解大数定理、中心极限定理背后的极限思想。具体内容安排如下：首先：我们介绍大数定理的原理，并用python语言对其进行模拟，给大家一个更直观的感受；接着：我们会介绍中心极限定理以及他在工程实践中的应用价值和意义，同时也会用python语言来对其进行模拟；最后：我们会讲解大数定理在机器学习和数据分析中的一个重要应用：蒙特卡罗方法，并仔细剖析他的应用场景和内涵。2.大数定理2.1.原理介绍好了，下面我们开始正式进入到大数定理的内容中。我们有如下的随机变量：X_1,X_2,...,X_n 他们彼此之间满足独立同分布，因此他们拥有相同的均值 \mu 和方差 \sigma^2 。此时，我们重点来研究这一组随机变量的均值：M_n=\frac{X_1+X_2+...+X_n}{n} ，显然 M_n 也是一个随机变量。那么，M_n 的期望和方差就是我们此时重点关心的问题。首先，我们从期望的定义入手，来观察一下随机变量M_n 的期望E[M_n] ：E[M_n]=E[ \frac{X_1+X_2+...+X_n}{n}]=\frac{1}{n}(E[X_1]+E[X_2]+...+E[X_n])=\frac{1}{n}\cdot n\cdot \mu=\mu=E[X_i]不难发现，一组独立同分布随机变量均值的期望就等于随机变量的期望，这个结论很直观。下面我们再来看看 M_n 的方差var[M_n] ：var[M_n]=var[ \frac{X_1+X_2+...+X_n}{n}]=\frac{1}{n^2}var[X_1+X_2+...+X_n]=\frac{1}{n^2}(var[X_1]+var[X_2]+...+var[X_n])=\frac{1}{n^2}\cdot n\cdot \sigma^2=\frac{\sigma^2}{n}我们从推导中发现，n 个独立同分布随机变量的均值的方差，是单一随机变量方差的\frac{1}{n}。没错，均值的方差变小了，并且随机变量 X的个数 n 越多，方差越小，他们的分布更加紧密的围绕在了期望的周围。特别的，当 n \rightarrow \infty 时，随机变量均值的方差趋近于0： var[M_n]=\frac{\sigma^2}{n} \rightarrow 0结合前前后后的一大段推导和论述，我们可以得出这么一个结论：独立同分布的随机变量 X_1,X_2,...,X_n ，他们的均值 M_n 的分布会更加接近于实际分布的均值 \mu ，随着样本量 n 的增大，他逐渐收敛于 \mu ，当 n \rightarrow \infty 时，也就是说当样本量非常大的时候，通过抽样样本计算所得到的平均值可以说就是 E[X] 了。独立同分布的随机变量序列的样本均值，在大样本的情况下，以很大的概率与随机变量的均值非常接近。这也就是为什么说，当独立同分布的样本数量 n 充分大时，样本均值（频率）是概率 P 的一个非常好的估计。这就回到本节最初提到的那个小问题了，样本数量到底是选100还是选1000，相信大家都会有明确的理论支撑了。这里我们还简单的提两个结论，其实说起来他们都比较直观。2.2.两个重要的不等关系谈到这里，我们来看概率统计中的两个非常重要的不等关系，以及他们背后的直观逻辑。第一个是马尔科夫不等式。用简洁直白的语言描述他就是：对于一个非负的随机变量 X ，如果他的均值很小的话，那么这个随机变量取到一个大值的概率是非常小的。想想确实很直观，因为随机变量的分布大部分都集中在均值附近，越远离均值，概率就越小。描述这个现象，有一个专门的不等式，叫做马尔科夫不等式：P(|M_n-\mu
\ge \epsilon) \le \frac{\sigma^2}{n \epsilon^2}第二个是切比雪夫不等式。我们同样简单的描述一下这个不等式的内涵：如果一个随机变量的方差非常小的话，那么这个随机变量取到远离均值 \mu 的概率也是非常小的，这个说实话也非常直观，同样有一个专门的不等式来描述他，叫切比雪夫不等式：P(|X-\mu|\ge c) \le \frac{\sigma^2}{c^2}2.3.大数定理的模拟下面我们来实际模拟一下大数定理。在第一个模拟的例子中，我们给大家一个感性的认识：我们生成3组各15000个服从参数为(10,0.4)的二项分布随机变量，随机变量的期望为 n*p=4 ，然后观察随着样本数目的增大，样本均值和实际分布期望之间的关系。代码片段：import numpy as np
from scipy.stats import binom
import matplotlib.pyplot as plt
import seaborn
seaborn.set()
n = 10
p = 0.4
sample_size = 15000
expected_value = n*p
N_samples = range(1, sample_size, 10)
for k in range(3):
binom_rv = binom(n=n, p=p)
X = binom_rv.rvs(size=sample_size)
sample_average = [X[:i].mean() for i in N_samples]
plt.plot(N_samples, sample_average,
label='average of sample {}'.format(k))
plt.plot(N_samples, expected_value * np.ones_like(sample_average),
ls='--', label='true expected value:np={}'.format(n*p), c='k')
plt.legend()
plt.show()运行结果：图1.大数定理模拟我们设置了三个相同的试验组，从试验结果中我们可以发现，在每一组试验中，随着样本数量的逐渐增大，样本均值都会越来越收敛于随机变量的期望。接下来我们再来看看第二个用于模拟大数定理的例子，我们从大数定理的定义出发，我们先生成1000000个服从均值为0，标准差为20正态分布的样本。依次进行三种不同的处理，并观察对应的三组分布图像：图像1：原始正态分布的样本分布图像，颜色为蓝色；图像2：从1000000个原始正态分布样本中，每次随机选取5个数，计算他们的均值，重复操作10000次，观察这10000个均值的分布，颜色为红色；图像3：从1000000个原始正态分布样本中，每次随机选取50个数，计算他们的均值，重复操作10000次，观察这10000个均值的分布，颜色为绿色；代码片段：import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import seaborn
seaborn.set()
norm_rvs = norm(loc=0, scale=20).rvs(size=1000000)
plt.hist(norm_rvs, normed=True, alpha=0.3, color='b', bins=100, label='original')
mean_array = []
for i in range(10000):
sample = np.random.choice(norm_rvs, size=5, replace=False)
mean_array.append(np.mean(sample))
plt.hist(mean_array, normed=True, alpha=0.3, color='r', bins=100, label='sample size=5')
for i in range(10000):
sample = np.random.choice(norm_rvs, size=50, replace=False)
mean_array.append(np.mean(sample))
plt.hist(mean_array, normed=True, alpha=0.3, color='g', bins=100, label='sample size=50')
plt.gca().axes.set_xlim(-60, 60)
plt.legend(loc='best')
plt.show()运行结果：图2.大数定理的模拟这个程序采样的规模比较大，可能运行的时间会比较长。从图中我们发现，随着每次选取的样本数量的增多，样本均值分布的图像越来越向期望集中，再一次佐证了大数定理。3.中心极限定理3.1.原理介绍下面接着来看另外一个现象，我们还是获取随机变量序列： X_1,X_2,...,X_n ，这 n 个随机变量满足独立同分布，均值为 \mu ，方差为 \sigma^2 。我们在这组随机变量序列的基础之上得到一个新的随机变量：Z_n=\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma}针对随机变量 Z_n ，我们很容易计算出：E[Z_n]=E[\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma}]=\frac{E[X_1+X_2+...+X_n]-n\mu}{\sqrt{n}\sigma}=0var[Z_n]=var[\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma}]=\frac{var[X_1]+var[X_2]+...+var[X_n]}{n\sigma^2}=1即，随机变量 Z_n 的期望为0，方差为1。关键的一点是，随着样本个数 n 的增大，随机变量Z_n 的分布逐渐趋向于一个标准正态分布，当 n \rightarrow \infty 时，随机变量的分布收敛于一个标准正态分布。更为重要的一点是，这个定理对随机变量 X 的原始分布没有任何要求，非常具有一般性。3.2.中心极限定理的工程意义实际上，中心极限定理中的随机变量Z_n=\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma} 是经过标准化处理的，如果单单只考虑 n 个随机变量的和，我们很容易得到：S_n=X_1+X_2+...+X_n 的分布趋近于一个均值为 n\mu ，方差为 n\sigma ^2 的正态分布。中心极限定理的意义在于，大量样本的独立随机因素的叠加是趋近于一个正态分布的，这一点在很多工程领域很常见也很关键。更重要的一点是，他不需要我们去搞明白随机变量X 的分布列或者概率密度函数，这往往是非常复杂的，甚至根本就无从得知，我们只需要知道他的均值和方差就可以进行后续的处理和分析了。3.3.中心极限定理的模拟下面，我们还是来举个例子，模拟验证一下中心极限定理。我们从一个服从参数p=0.3的几何分布中进行采样，共分三组试验，分别每次采样2个、5个、50个样本，每组试验各重复100000次，然后按照Z_n=\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma} 进行标准化，得到三组试验对应的结果，最后对试验结果进行可视化观察：代码片段：import numpy as np
from scipy.stats import geom
import matplotlib.pyplot as plt
import seaborn
seaborn.set()
fig, ax = plt.subplots(2, 2)
geom_rv = geom(p=0.3)
geom_rvs = geom_rv.rvs(size=1000000)
mean, var, skew, kurt = geom_rv.stats(moments='mvsk')
ax[0, 0].hist(geom_rvs, bins=100, normed=True)
ax[0, 0].set_title('geom distribution:p=0.3')
n_array = [0, 2, 5, 50]
for i in range(1, 4):
Z_array = []
n = n_array[i]
for j in range(100000):
sample = np.random.choice(geom_rvs, n)
Z_array.append((sum(sample) - n * mean) / np.sqrt(n * var))
ax[i//2, i%2].hist(Z_array, bins=100, normed=True)
ax[i//2, i%2].set_title('n={}'.format(n))
ax[i//2, i%2].set_xlim(-3, 3)
plt.show()运行结果：图3.中心极限定理的模拟左上第一幅图是几何分布的原始图像，我们发现，随着单次采样个数的逐渐增加，随机变量Z_n=\frac{X_1+X_2+...+X_n-n\mu}{\sqrt{n}\sigma} 的分布图像越来越趋近于一个标准正态分布。4.大数定理的应用：蒙特卡罗方法4.1.机器学习中的应用背景用大样本数据计算出来的频率去估计概率，这就是大数定理的本质，而大数定理思想的一个非常典型的应用就是蒙特卡罗方法。蒙特卡罗方法，又叫统计模拟方法，名字很洋气，思想很粗暴，真的很管用。它使用随机数来进行场景的模拟或者过程的仿真，其思想核心就是通过模拟出来的大量样本集或者随机过程去近似我们想要研究的实际问题对象，这是一类非常重要的数值计算方法。该方法的名字来源于世界著名的赌城蒙特卡罗。赌博和概率，二者相视一笑、不谋而合。这种方法最初应用于20世纪40年代美国的曼哈顿原子弹计划，如今在数据分析和机器学习领域中到处都有他的身影。以下是蒙特卡罗方法的几类非常典型的应用：1.近似计算不规则面积/体积/积分2.模拟随机过程，预测随机过程可能性结果的区间范围3.利用马尔科夫链-蒙特卡罗方法（MCMC）进行未知参数的统计推断要知道，蒙特卡洛方法不仅仅只是一种方法技巧，更是一种思考问题的方式。在这篇回答中，我们主要介绍第一个应用点，让大家进一步理解大数定理的原理和蒙特卡罗方法的要点。第二个和第三个应用如果大家有兴趣可以关注我的专栏。4.2.利用蒙特卡罗方法计算不规则面积蒙特卡罗方法可以近似计算出不规则图形的面积，特别对于那些难以用解析方法计算的图像，非常有效。这里，我们利用蒙特卡罗方法来近似计算一个圆的面积，然后估计出 \pi 的近似值，选择圆作为例子的原因不是因为圆无法通过解析法进行计算，而是因为大家对他的面积比较熟悉，方便我们进行结果的对比。代码片段：import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from scipy.stats import uniform
import seaborn
seaborn.set()
n = 100000
r = 1.0
o_x, o_y = (0., 0.)
uniform_x = uniform(o_x-r,2*r).rvs(n)
uniform_y = uniform(o_y-r,2*r).rvs(n)
d_array = np.sqrt((uniform_x - o_x) ** 2 + (uniform_y - o_y) ** 2)
res = sum(np.where(d_array < r, 1, 0))
pi = (res / n) /(r**2) * (2*r)**2
fig, ax = plt.subplots(1, 1)
ax.plot(uniform_x, uniform_y, 'ro', markersize=0.3)
plt.axis('equal')
circle = Circle(xy=(o_x, o_y), radius=r, alpha=0.5)
ax.add_patch(circle)
print('pi={}'.format(pi))
plt.show()运行结果：pi=3.14096图4.蒙特卡洛方法近似计算圆面积我们结合这张图，详细的分析一下这段程序。我们近似计算的目标就是这个蓝色的半径为r=1的单位圆面积，而这个单位圆的外接正方形的边长为l=2r=2，因此外接正方形的面积为4。我们生成100000个在外接正方形内均匀分布的点，均匀的撒下去，这里面的核心原理就是：\frac{圆形面积}{正方向面积}\approx \frac{圆内点的个数}{点的总个数}这样就可以估算出单位圆的面积了。而为了估算\pi ，我们回到这个例子中的参数，可以得到这么一个公式：\frac{\pi r^2}{(2r)^2}\approx \frac{圆内点的个数}{点的总个数}\rightarrow \pi \approx 4 * \frac{圆内点的个数}{点的总个数}而且大数定理告诉我们，随着样本数量的增大，我们用这种方式模拟出来的 \pi 值应该是越来越趋近于真实值，样本无穷大的时候收敛于真值。这就证明了应用大数定理的蒙特卡罗方法的合理性和有效性。至于有同学会问了，是怎么判定一个点是否位于圆当中？这里用的是距离法，也就是计算每个点到原点的距离，如果小于等于半径，就说明这个点在圆内。这里是不是有一个bug？就是说这个规则只对圆这种特殊图像有效，如果是一个完全不规则的图像，我们如何来计算这个图像的面积？我们可以借助计算机图像中的像素来进行。比如把不规则图像内部都涂黑，图像外部都留白，我们还是均匀的撒下大样本量的点，如果某个点位于的坐标，他的像素是黑色的，则证明这个点在图像内部，反之就在图像外部，利用这个方法就能统计出位于图像内部的点的个数。5.小结在这一篇回答中，我们从样本均值和总体期望的关系入手，详细的讨论了概率统计中的重要概念：大数定理与中心极限定理，目的是让大家能够建立好概率统计中极限思维的概念基础。在实践环节中，我们介绍了大数定理的重要应用：蒙特卡罗方法，列举了他在机器学习中的几类重要应用，希望能有助于大家强化大数定理的思想方法。此内容节选自我的专栏《机器学习中的数学：概率统计》，前三小节免费试读，欢迎订阅当然还有《机器学习中的数学（全集）》系列专栏，欢迎大家阅读，配合食用，效果更佳~有订阅的问题可咨询微信：zhangyumeng0422}