在最后一章，小弟为大家介绍两个统计方法，说老实话，这两个方法可以说的是上古秘籍。统计的起源最早可追溯到1910年，那时候计算机还没出生呢，当中国还是清朝伪满洲的时候，有这么一群数学家利用纸和笔在条件很艰苦的情况下，想出了一些统计方法，这些方法就是本章要介绍的。我们不能忘记过去数学家们的智慧和对如今大数据时代的贡献。

仿真法/统计模拟方法(simulation based methods)

在18世纪，也就是中国乾隆时期，西方有这么个问题—蒲丰投针问题:
假设我们有一个以平行且等距木纹(黑线)铺成的地板（如图），

随意抛一支长度比木纹之间距离小的针(红线段)，求针和其中一条木纹相交的概率。
经过蒲丰计算，如果条纹间距为a，针的长度为$l$,那么针与其中平行线中任意一条相交的概率为$\frac{2l}{\pi a}$.这里推导也是比较简单的(相关推导网上均有)，但小弟在这还是写几笔:

设X表示针的中点到平行线的距离，X服从均匀分布(0,$\frac{a}{2}$),也就是(0<X<$\frac{a}{2}$)的概率是均匀的(数据分析模型第二章，连续均匀分布)
$$f(X)=\{\begin{array}{ll}\frac{1}{\frac{a}{2}-0}=\frac{2}{a},& 0<X<\frac{a}{2}\\ 0,& 其他\end{array}$$

设Y表示针与平行线的夹角。Y服从均匀分布(0,$\frac{\pi }{2}$)
$$f(Y)=\{\begin{array}{ll}\frac{1}{\frac{\pi }{2}-0}=\frac{2}{\pi },& 0<Y<\frac{\pi }{2}\\ 0,& 其他\end{array}$$

X和Y相互独立,那么X和Y的联合概率为:
$$f(X,Y)=\{\begin{array}{ll}\frac{2}{a}\ast \frac{2}{\pi }=\frac{4}{\pi a},& 0<X<\frac{a}{2},0<Y<\frac{\pi }{2}\\ 0,& 其他\end{array}$$

当$X<\frac{l}{2}sinY$时，针与线相交。那么:
P { X < l 2 s i n Y } = ∫ 0 π 2 ∫ 0 1 2 s i n ( y ) f ( X , Y ) d x d y = ∫ 0 π 2 ∫ 0 1 2 s i n ( y ) 4 π a d x d y = 2 l π a P\{X<\frac{l}{2}sinY\}=\int_{0}^{\frac{\pi}{2}}\int_{0}^{\frac{1}{2}sin (y)}f(X,Y)dxdy=\int_{0}^{\frac{\pi}{2}}\int_{0}^{\frac{1}{2}sin(y)}\frac{4}{\pi a}dxdy=\frac{2l}{\pi a} P{X<2l​sinY}=∫02π​​∫021​sin(y)​f(X,Y)dxdy=∫02π​​∫021​sin(y)​πa4​dxdy=πa2l​
（根据数据分析模型第二章知识点和内容）
因为$X<\frac{l}{2}sinY$,所以${\int }_0^{\frac{1}{2}sin(y)}$,因为$0<Y<\frac{\pi }{2}$所以${\int }_0^{\frac{\pi }{2}}$. 先算${\int }_0^{\frac{1}{2}sin(y)}(\cdot )dx$，后算${\int }_0^{\frac{\pi }{2}}(\cdot )dy$。

有位叫蒙特卡洛(Monte Carlo)同学给了另外一个解决方法，那就是实验，硬干它。经过大量实验，记录多少次针与条纹相交的次数，从而得到那么针与其中平行线中任意一条相交的概率。我们也知道根据弱大数理论，如果实验多次，那么它估计的概率其实也是很准确的。

但当时人家蒙特卡洛同学本意不是为了估计针与条纹相交的概率，人家统计这个概率是为了估计$\pi$值。利用实验得到的概率具体值和蒲丰利用微积分计算的概率$\frac{2l}{\pi a}$,来反向估计出$\pi$值。所以蒲丰投针问题也是用来计算$\pi$值的。

以现在看，当时人家蒙特卡洛同学通过实验的方法来估计某一事件的概率，大家觉得这不就是小儿科么。但这就是早期统计的起源，也是最早的统计模拟方法。

有的时候，对于某件事，我们很难纯用数学来计算它的概率，举几个例子:

1.如果$X_1,X_2～N(0,1)$，那么$P(\sqrt{|X_1|}+\sqrt{|X_1|}>2)$是多少？

2.如果$X_1～Poi({\lambda }_1),X_2～Poi({\lambda }_2)$,那么$P(X_1-X_2^{\frac{3}{2}}>0)$是多少？

这些问题，要是靠纯数学计算推导是很难的，这个时候我们需要模拟这个事件，也就是做实验来计算它的概率。

实验概率/经验概率收敛(Convergence of Empirical Probabilities)
我们可以利用统计模拟方法来解决这些问题:

第一步:
反复m次，每一次从它的分布中产生$X_1,X_2..,X_p$个数据。

第二步:估计概率，利用实验
$$P(X_1,X_2,..,X_p\in Z)\approx \frac{1}{m}\sum _{i=1}^{m}I(X_1,X_2,..,X_p\in Z)$$
$I(\cdot )$为指标方程，即满足条件则为1，不满足条件则为0

举个例子:
$X～N(0,1)$,计算 P { ∣ X ∣ > 1 } P\{\sqrt{|X|}>1\} P{∣X∣ ​>1}的概率

我们先利用纯数学来计算下:
设:$Q=\sqrt{|X|}$,那么根据正态分布$N(0,1)$,得:
$$P(Q=q)=(\frac{2^3}{\pi }{)}^{\frac{1}{2}}qexp(-\frac{q^4}{2})$$
那么:
$$P(Q>1)={\int }_1^{\infty }P(q)dq=1-\frac{2}{\sqrt{\pi }}{\int }_0^{x}exp(-t^2)dt=0.3171$$

我们在通过实验概率/经验概率收敛的方法来算下:

#MATLAB或者R
X=rnorm(m,0,1)
mean(sqrt(abs(X))>1)
>m=100, 0.3300
>m=1000, 0.3310
>m=100000, 0.3176
>m=10^9, 0.317102

那么根据抽取$10^9$次数据(每次实验抽一个数据)，我们利用实验估计 P { ∣ X ∣ > 1 } P\{\sqrt{|X|}>1\} P{∣X∣ ​>1}的概率大约为0.317102.

我们可以发现其实通过实验来计算概率和纯数学计算的结果几乎一样，当我们的实验次数足够多的情况下。

当然了这个方法也有缺点:
那就是实验多少次才算足够，当我们计算出的概率结果比$\frac{1}{m}$还要小，说明我们的计算已经被$\frac{1}{m}$所影响了，那样计算出的概率结果是不准确的。因为如果当m足够大时，我们的$\frac{1}{m}$会趋近于0.

但我们有解决这个缺点的办法:
我们会先叫实验的次数和每次实验抽取多少数据均为m，如果最后计算的概率结果小于$\frac{1}{m}$,那就说明估计不准确，对于较小值的概率结果，我们则需要较大的实验次数，试验次数至少是$\frac{1}{概率结果}$的10倍。

伪随机数(Pseudo Random Numbers)
在做仿真法/统计模拟法时，我们会需要一些随机数，例如:根据分布随机生成数据，随机生成数字。
虽然计算机可以随机生成数字，但如果计算机是根据一些特定规则进而随机生成数字，这类随机数被称为伪随机数。伪随机数可以用计算机大量生成，是模拟研究中重要的一步，可以提高效率。

伪随机数最初是利用均匀分布的思想，但如今小弟就不知道了，因为创造伪随机数的方法多了，例如马特塞特旋转演算，线性同余法，乘同余法，后面的两种方法均是离散数学范畴但还是用到了均匀分布还是比较简单的。据小弟所知，大部分编程语言创造伪随机数还是利用均匀分布的思想。

小弟在此就分享下利用均匀分布来创建伪随机数，若有其他创造伪随机数的方法，那基本上是大同小异。

均匀分布随机变量产生其他分布的随机变量
首先我们要有个均匀分布$U(0,1)$,选取变量∈(0,1)。那么我们可以利用均匀分布变量∈(0,1)来得到其他分布的变量。

第一步:给定其他分布或者概率密度函数$Z=F(Y)$

第二步:计算它的反函数$F^{-1}(Z)=Y$,那么我们都知道这个$Z$∈(0,1)。有的没有反函数就分段。

第三步：生成随机变量X从均匀分布$U(0,1)$

第四步: 带入反函数$F^{-1}(X)=Y$那么我们便利用均匀分布的变量$X$，得到了其他分布的变量$Y$了。

我们既然可以利用均匀分布随机变量得到其他分布的随机变量。那么，如何叫它随机呢？

随机化:
随机化主要是在于随机选取均匀分布的变量上。
$$f(X_0,X_1,...,X_{n-2},X_{n-1})=X_n$$
举个例子:
我们先初始给一个$X_0$的值，该值被称为种子(seed)。$f(\cdot )$被称为转换方程。这里转换方程就很多了，各有各的优点，例如有的方法可以使随机数更均匀，有的影响随机数列的周期长短，有的需要参数间互为质数等等。
当我们有了$X_0$,我们利用$f(\cdot )$得到$X_1$,再利用$X_0,X_1$带入$f(\cdot )$中得到$X_2$。但无论怎样我们的$X_0,X_1,...,X_{n-2},X_{n-1}\in (0,1)$.

这样我们的实现了在均匀分布上随机选取变量，从而随机产生其他分布上的变量。

伪随机数缺点:
好的种子比较关键。另外如果使用相同的种子，那么它们产生的伪随机数是相同的。

自助法(Bootstrap method)

自助法也经常被用在统计模拟方法中，自助法和统计模拟法均为大家族多次采样方法中的一员。我们在前两章其实已经了解了它们家族的另外一员，那就是CV(交叉验证)。

对于交叉验证，我们估计对于未来数据(测验数据)预测的差值。

对于自助法, 它是用来看看我们的估计是否稳定(例如标准差，置信区间等等)

在第三章中，我们经过多次取样，我们会得到样本分布，即如下图:

Population, 则是我们的总体数据，多次取样(samples)，取相同数量的数据，再利用每次取样得来的样本数据估计总体参数。

举个例子:
我们目前有一组很感兴趣的样本，那么我们可以利用上图的自助法分布来估计这组由感兴趣样本所估计出的参数$\bar{\theta }$偏差和方差

假如我们已经反复取样了M次，分别估计它们的参数，偏有了M个${\bar{\theta }}^{(i)}$，那么我们可以估计出该感兴趣样本参数的$\bar{\theta }$的偏差和方差

估计$\bar{\theta }$的偏差:
$$bias=\frac{1}{M}\sum _{i=1}^M(\bar{\theta }-{\bar{\theta }}^{(i)})$$

估计$\bar{\theta }$的方差:
$$Var=(\frac{1}{M-1})\sum _{i=1}^M({\bar{\theta }}^{(i)}-\bar{\theta }{)}^2$$

除此之外，经过多次取样后，我们自然也会得到样本的分布。就像我们之前第三章和第四章内容所讲，经过利用样本的估计参数后(就是一些数学转换)，我们自然也会得到估计参数的分布。

举个例子：
我们抽一次样本得$y=(y_1,y_2,..,y_n)$,该样本有n个数据，并作为我们感兴趣的样本。

我们利用该感兴趣的样本来估计参数$\bar{\theta }$

利用自助法逻辑：
1.抽取M次数样本，每个样本均也需要有n个数据
2.每个样本分别估计参数记作${\bar{\theta }}^{(i)},i=1,2,3..,M$

那么我们便可以做很多事情了，估计我们兴趣样本的偏差，方差。当然了，因为我们有了${\bar{\theta }}^{(i)},i=1,2,3..,M$的分布，我自然也可以计算这个分布的置信区间。

大家不要被小弟上述例子从而局限了想象力，自助法的本质就是多次取样。还记得第八章中的一幅图么，即如下:

我们不看绿线，就看图中橘色虚线和线段，即平均拟合模型和95%平均拟合模型。这里95%其实就是95%${\bar{\theta }}^{(i)},i=1,2,3..,M$的分布置信区间，然后把估计的参数依次带入到模型里去，然后在把你想预测的数据代入进去，这样你的预测就会有个范围即95%。至于橘色线段即平均拟合，则是算的预测值的均值罢了。

在第6章，我们也接触过简单线性回归对于预测的置信区间，还记得么$(Y-{\bar{\beta }}_0-{\bar{\beta }}_1{x}_0)/[\sqrt{\frac{n+1}{n}+\frac{(x_0-\bar{x}{)}^2}{({\sum }_{i=1}^n{x}_i^2-n{\bar{x}}^2)}}\sqrt{\frac{RSS}{n-2}}]～t_{n-2}$，$x_0$则是测试数据。对没错，有了分布我们便可以得到关于简单线性回归的预测95%置信区间。你可以用纯数学的方式来找到关于某个模型的分布，但现实很残酷，有些模型过于复杂，利用纯数学的方式我们是很难做的，旦旦一个简单线性回归的分布都如此复杂，就不要说比简单新型回归更复杂的模型了。还是以上图为例，该图其实是多项式分布模型，绿色线为真多项式分布，并且有5项，要想用数学找到该模型的分布是很困难的。所以才会有蒙特塔罗的思想，就是多做实验，多抽样，用简单粗暴的方式。

总之一句话，多做实验，多动手，多取样，再看分布，有时候你觉得能用数学推，到最后你会发现现实很残酷。

自助法的优点:
简单粗暴，很容易做，也很容易理解。
估计准确度会随着样本的增加而增加
自助法经常用在计算偏差，方差，预测差，置信区间。

自助法缺点:
算起来太慢了。
对于某些问题需要转换，例如拉索

排列检验(Permutation tests)
在最后，小弟介绍一个另一个统计方法，即排列检验，经常用在样本分布上面。其实大家早就知道了，它就是p值。

简单回顾下
$$H_0:\beta =0,H_1\ne 0$$
计算p值，看是否拒绝原假设。

那么如何利用p值用在多次抽样上面呢？
我们用例子来解释:

我们要用排列检验计算p值，来看是否$\bar{\beta }=1.4310$

现在我们从这些数据(总体)里取样。取个m次，分别估计参数${\bar{\beta }}^{(i)}$.
此时我们计算p值:
$$p\approx \frac{1}{m}\sum _{i=1}^{m}I(|{\bar{\beta }}^{(i)}|>|\bar{\beta }|)$$
$I(\cdot )$符合条件为1，不符合条件为0

我们来张图：

经过m=10000次的抽样，我们估计10000个${\bar{\beta }}^{(i)}$即我们的x轴.计算频率，则有了y轴。我们利用所有数据估计的$\bar{\beta }$=1.4310.可以看到$P(|{\bar{\beta }}^{(i)}|>|\bar{\beta }|)$即上图±1.4310的两个外端概率是很小的约为0.0013(p值太小)拒绝原假设。我们其实能够清楚看到，也侧面说明了，$\bar{\beta }\ne 1.4310$。

排列检验与我们之间的假设检验区别在于，排列检验不需要基于原假设，但假设检验需要基于原假设去计算。

在排列检验里，其实它更注重的是多次取样后的数据分布。然后根据条件，计算最小概率即p值来进行观察原假设或者备择假设。

排列检验的优点:
对于非正态分布的模型计算p值更准确
好写代码，很好应用

排列检验的缺点
结果很容易受到采样的次数影响。
可能会运行很慢

三. 结语

有公式推导错误的或理论有谬误的，请大家指出，我好及时更正，感谢。

感谢大家来阅读小弟我分享的一些关于基础统计的知识。如果将来有时间，小弟会陆续分享更多知识，例如高等数据分析，机器学习，计算机视觉，神经网络，优化等几门课的知识分享。也希望小弟可以坚持分享下来吧。