从任意分布中构造任意分布

2020-07-23

前言

设想我们希望在计算机中编写程序生成满足给定分布的随机数,例如说,生成一些近似服从正态分布或者近似服从均匀分布的随机数,可是,我们仅能够获得产生自不确定分布的随机数,比如说能够接触的随机源仅包括/dev/urandom/dev/random,和/dev/arandom,在这样的条件下,如何能获得满足特定分布的随机数呢?注意到/dev/urandom中的 “u” 并不是 “Uniform distribution” 中的 “u”,而是 “Unlimited” 中的 “u”,表示“能够无限地产生随机数”,即想要多少个就能够产生多少个随机数.而且,我们并不能够确定/dev/urandom在任何时候都是近似服从均匀分布的(近似服从表示相差不太远,表示在$\alpha=0.05$的置信水平下不能拒绝$\chi^2$检验的原假设$\mathsf{H}_0: 分布是均匀的$).

两种思路

两种思路大概可总结概括为:1)修正法,和:2)转换法.

先说修正法,首先我们知道:通过收集一些数据,我们能够对这些数据做一些非参数(Non-parametric)检验(使用非参数统计的方法进行假设检验),这些非参数检验有助于我们得出结论——这些数据是不是产生自某某分布,那么具体来说,我们定义一个全局变量,存储历史生成的数据,或者至少是历史生成的随机数的每个数字的数量统计,然后,每生成一个新的随机数,我们将这个新的随机数与历史数据一并,做某种非参数检验,如果非参数检验没有拒绝原假设,那么我们就接受这个随机数,并且更新历史数据,如果说原假设被拒绝了,就说明如果接受这个随机数的话势必将使得历史分布和理论分布偏离过大换句话来说就是得到的不是我们想要的分布了,所以就拒绝这个随机数并重新生成另一个再进行检验直到通过为止,并且得到能通过假设检验的随机数之后也更新历史数据,这样一来,很容易理解,从始至终,对历史数据进行假设检验都表明服从某种分布,亦即服从我们想要的分布,从而我们的目的就达到了.这就是使用「修正法」来产生任意给定分布的基本思路,就好像是一艘略微失控的宇宙飞船,修正法就好像某种意志,监控该宇宙飞船的飞行状态,在必要时刻对其进行干预,从而确保它是指是在预定航线上飞行而不偏离预定航线太远.或者,你也可以将该方法理解为统计学上的「重抽样」.

再说转换法,我们都听说过中心极限定理,中心极限定理是说,对于一个未知的分布$D$(假设它均值为$\mu$,方差为$\sigma^2$且方差$\sigma^2 > 0$),我对它进行多次抽样,第$i$次从$D$中抽取$m$个随机数($m$相当大),并且对于第$i$次抽样,我们记$Y_i$为

$$ Y_i = \frac{X_1 + X_2 + \cdots + X_m}{m} $$

假设我们前后总共进行了$n$次这样的抽样(同样地,$n$也相当大),我们就这样收集到一组抽样数据$Y_1, Y_2, \cdots, Y_n$,我们从这$n$个$Y_i$数据画出的直方图会看到,$Y_i$似乎是产生至某组参数下的正态分布,事实上,我们假设这个$Y_i$它所来自的分布是${D}^\star$,并且设$Y$是产生自${D}^\star$的随机变量,那么,事实上,必然有

$$ \mathsf{E}(Y) = \mu \\ \mathop{var}(Y) = \frac{\sigma^2}{n} $$

并且${D}^\star$是正态分布.

也就是说这种通过“每次抽样并且取算术平均”构造得来的抽样分布正是正态分布,而不用管原来的分布是什么分布!鉴于此,我们应该知道怎样对未知随机数产生器产生的随机数进行处理并最终得到满足给定条件的正态分布了.

高尔顿板,或者说,当试验次数足够多时,二项分布的试验结果看起来也近似是正态的,试验结果的直方图画出来近似为正态分布的直方图.

非参数检验

在非参数统计(Nonparametric statistics)中,我们一般不知道得到的数据是服从什么分布,举例来说,在与之相反的参数统计(Parametric statistics)中,我们可以对一组数据进行$t$检验,以此判断这组数据是否参数自一个均值为$\mu$的正态分布,然而,这要求数据来自正态分布,因为只有假定数据是来自正态分布之后,才能说该组数据的算术平均是服从$t$分布.换句话说,在非参数统计中,很多情况下,对于数据是服从什么分布,我们根本不知道,一般也不随意做假设,因为假设往往与实际偏差过大,也就是说,非参数统计尝试对来自我们不熟悉的情形的随机变量进行建模,而参数统计所讨论的情形一般都是我们已然熟悉的.

设想,一个骰子(Dice),这个骰子来自一个陌生的赌场,我们不知道该赌场的经营者的道德人品具体如何,我们也不清楚该骰子的具体制作过程和具体质地,换用非参数统计的话来说:我们对这个骰子它所代表的分布知之甚少,然而,如要判断该骰子是否是公平的,即,该骰子掷出六个面当中各个面出现的概率是否均等,我们可以收集数据(掷很多次这个骰子),并且使用一种或者若干种检验方法来检验之,我们先介绍比较简单易懂的卡方检验($\chi^2$-test)法:卡方检验方法对比的是实验数据和理论数据,计算它们的差别,当这两者的差别太大的时候,原假设将被拒绝,原假设拒绝就意味著理论和实验相差过大,就意味著经验分布和理论分布有显著的不一致.

首先我们来看一枚质地均匀的骰子,我们这样来构造它:

fairDice = EmpiricalDistribution[
    {1/6, 1/6, 1/6, 1/6, 1/6, 1/6} -> {1, 2, 3, 4, 5, 6}
]

fairDiceData = RandomVariate[fairDice, 10000];

Counts[fairDiceData]

用这些Mathematica代码我们在计算机中模拟了一枚质地均匀的骰子,并且掷了这枚骰子一万次,并将实验结果打印出来,如图:

质地均匀的骰子投掷实验统计数据

质地均匀的骰子投掷实验统计数据

输出结果的意思是:2出现了1624次,6出现了1719次,4出现了1651次,1出现了1689次,5出现了1654次,3出现了1663次.直观地看来,这枚骰子应该是公平的,因为这些数量并没有相差太多.下面我们开始做卡方检验:

首先写出原假设和备则假设:

$$ \begin{cases} H_0: 这枚骰子出现各个面的概率全都是1/6 \\ H_1: 这枚骰子出现各个面的概率不全是1/6 \end{cases} $$  

为了检验这一组假设,我们首先计算原假设所对应的理论数值,即,假如说原假设成立的话,那么每个面应该出现多少次,其实也很好算,我们总共做了$10000$次试验,假如说原假设成立,那么每个面都应该出现$10000 \times 1/6$次也就大约是$1666$次.现在为了描述实验结果和理论上每个面出现次数的差距,我们引入卡方统计量

$$ X^2 = \sum_{i = 1}^{k} \frac{(o(X=i)-t(X=i))^2}{o(X=i)} $$

求和范围$i=1$到$i=k$对应骰子可能出现的值是从1到6,所以这里$k$等于$6$,括号里的$t(X=i)$表示骰子掷出$i$的次数,在这里就是$t(X=1) = 1689, t(X=2) = 1624$,而$o(X=i)$表示$t(X=i)$的理论值,在这里就是$o(X=i)=10000/6$.总的来说就是“实际减期望求平方再比上期望再求和”,很容易看出当期望频数和实际频数相差较大时,这个值也大,反之同理,因而它可以用来度量实验数据和理论数据的偏差程度.卡方统计量都有自由度,对于$k$个面的骰子来说它对于的卡方统计量的自由度就是$k-1$,那么在这里这个卡方统计量的自由度就是$5$,这个数字可以用来查分位数(也就是临界值).按照这一个公式,我们现在就可以来算一下这个卡方统计量:

$$ \begin{align} X^2 &= \frac{(1689 - 10000/6)^2}{10000/6} + \frac{(1624 - 10000/6)^2}{10000/6} \\ &+ \frac{(1663 - 10000/6)^2}{10000/6} + \frac{(1651 - 10000/6)^2}{10000/6} \\ &+ \frac{(1654 - 10000/6)^2}{10000/6} + \frac{(1719 - 10000/6)^2}{10000/6} \\ &\approx 3.2864 \end{align} $$  

意思就是说,根据这一万次试验收集到的数据,我们同时还得到一个自由度为5的服从卡方分布的统计量的观测值为$3.2864$,如何理解这个数字呢?$X^2$它服从自由度为$5$的卡方分布$\chi^{2}(5)$,它取什么值,是有一个概率的,实际上,我们会得知,这个$X^2$它取值大于或等于$3.2864$的概率是大概是等于

$$ \mathsf{Pr}(X^2 \geq 3.2864) \approx 0.655924 $$

原假设$H_0$同时还意味著:实际频数和理论频数之间的误差是随机产生的,现在,有了这个概率,我们知道,这个误差(现在由这个卡方统计量$X^2$描述)在$[3.2864, +\infty)$的概率是$0.655924$,假如说我们要拒绝原假设,说“这个使得$X^2=3.2864$的误差的产生并非偶然,而是骰子根本就不均匀”,那么,我们将面临$0.655924$这么大的犯Ⅰ类错误的可能性,显然,我们证据不足,不足以斩钉截铁地拒绝原假设,按照无罪推定原则,没有足够有力的证据,万不得已的情况下我们不会轻易地拒绝原假设.现在,我们还不能说骰子是不公平的,至少,仅仅根据现有的数据,看不出来如此.

现在,我们再来看一枚质地不均匀的骰子:

unFairWeight = {4, 6, 1, 7, 10, 9}
probs = Normalize[unFairWeight, Total]
unFairDice = EmpiricalDistribution[probs -> {1, 2, 3, 4, 5, 6}]
unFairDiceData = RandomVariate[unFairDice, 10000];
Counts[unFairDiceData]

输出结果如下

用一枚不公平的骰子做的实验

用一枚不公平的骰子做的实验

会发现,为了故意构造这样一枚不公平的骰子,我们故意地将出现6个面的概率分别设置为$4/37$,$6/37$,$1/37$,$10/37$,和$9/37$,实验结果也不出我们所料,在一万次试验中,该骰子掷出5的次数为$2669$次,掷出6的次数为$2430$次,掷出1的次数为$1080$次,掷出4的次数为$1896$次,掷出3的次数为$277$次,掷出2的概率为$1648$次,很显然,其中很多数字都和一万次试验下平均每个面出现$1666$次这个理论频数相差甚远,我们再设立与上边相似的一组假设(原假设说骰子公平,备择假设说骰子不公平),并且按照类似的方法计算卡方统计量以对假设进行检验:

$$ \begin{align} X^2 &= \frac{(2669-10000/6)^2}{10000/6} + \frac{(2430-10000/6)^2}{10000/6} \\ &+ \frac{(1080-10000/6)^2}{10000/6} + \frac{(1896-10000/6)^2}{10000/6} \\ &+ \frac{(277-10000/6)^2}{10000/6} + \frac{(1648-10000/6)^2}{10000/6} \\ &\approx 2349.39 \end{align} $$  

可现在,如果是这个使得卡方$\chi^{2}(5)$统计量$X^2$取值为$2349.39$的实际——理论误差是随机产生的,是产生于偶然因素的,而不是因为骰子本来就不公平,那么,随机产生这么大的误差的概率有多大呢?

$$ \begin{align} p &= \mathsf{Pr}(X^2 \geq 2349.39) \\ &\leq \mathsf{Pr}(X^2 \geq 234.939) \\ &\leq \mathsf{Pr}(X^2 \geq 23.4939) \\ &\approx 0.000271529 \end{align} $$  

不会超过万分之$2.72$.也就是说,一枚均匀的骰子要能够掷出一组和理论偏差这样大的试验结果,概率顶多不会超过万分之$2.72$,要说这枚骰子是均匀的,那几乎是不可能的!因此,我们已经有足够的证据(试验数据和计算结果)来拒绝原假设,我们可以说:我们可以在犯Ⅰ类错误的概率不超过$0.000271529$的基础上拒绝主张骰子公平的原假设,并且接受主张骰子不公平的备择假设,从而判定该枚骰子是不公平的.

通过上面的两个例子:一个正例和一个反例(构造出来的公平的骰子和不公平的骰子),我们讲完了如何收集数据,并且如何通过卡方检验的方法来判断一枚质地未知的骰子是否是公平的,可看做是非参数统计学讨论过的的其中一种方法的简单介绍,如果你觉得这个卡方检验的例子略有复杂,下面,我还将讲解一个更加简单的例子:如何来判断一枚质地未知的硬币是否是公平的.

我们用到的方法被称为二项式检验(Binomial test),已知某种随机过程$X$能够产生两种可能的输出结果:$X=A$,和$X=B$,并且已知产生$X=A$的概率是$Pr(X=A)=p$,相应地,产生$X=B$的概率就只能是$Pr(X=B)=1-p$,现在,我们进行$n$次这种试验,即,设法使这种随机过程发生$n$次,并且统计$n$次试验中,出现$X=A$的次数并将其记为$Y$,那么这个$Y$,它实际上就是服从二项分布(Binomial distribution),记做$Y \sim B(n, p)$,二项式系数符号${n \choose k}$表示:从$n$个独立个体组成的组合中不放回地抽取$k$个个体,得到的所有可能的由$k$个个体组成的组合的个数,当然,它同时还表示$(1+x)^2$的完全平方展开式中$x$的次数为$k$的项的系数(正因此被称为二项式系数),那么,我们容易证明,$Y=k, 0 \leq k \leq n, k \in \mathbb{Z}$的概率是

$$ \mathsf{Pr}(Y=k) = {n \choose k} p^k (1-p)^{n-k}. $$

回到硬币的话题上来,为检验一枚硬币是否是公平的(公平意味著该硬币掉落在平面上时出现正面和反面的概率都是$0.5$),我们按照惯例,首先建立一组待检验的假设,原假设是我们没有足够证据就不能够也不应当随意拒绝的选项(遵从「无罪推定原则」):

$$ \begin{cases} H_0: 待检验的硬币是公平的,出现正面的概率是0.5 \\ H_1: 待检验的硬币不是公平的,出现正面的概率不是0.5 \end{cases} $$  

作为演示,我们构造两个例子同样是一正一反,分别包含一枚公平的硬币,和一枚不公平的硬币:

从输出结果中我们可以看到,对于那一枚公平的硬币,其10000次抛掷试验中有4963次是出现正面(我们约定数字1代表正面,数字0代表反面),5037次是出现反面,相差并不太大,而对于另外那一枚不公平的硬币,其10000次抛掷试验中有6955次是出现正面,3045次是出现反面,差距是非常明显的.我们知道,这种差距一眼就能看出来,可是怎样套用公式来得出结论呢?

硬币的抛掷结果无法是正反两面,而出现正面的次数和出现反面的次数定然很难严格相等,因此,出现正面的次数和出现反面的次数这两个数当中定然就有比较小的那一个,这两数中较「小」的那一个有多「小」就正好度量了硬币的实际频数与理论频数的偏差程度,正因此,例如当正面次数和反面次数都非常接近时,最小的那个数也不会太小,例如案例中那枚公平的硬币:4963次正面和5037次反面,最小的是4963,但是距离0仍然相差很远,而另外那一枚不公平的硬币:6955次正面和3045次反面,这其中的3045次反面距离6955可就比4963距离5037大得多,距离0也相比4963距离0更近,也就能说明这枚硬币的实际频数与理论频数相差更大,也就更有理由怀疑这枚硬币是不公平的.

所以,我们的检测方法就是,取正面频数和反面频数两数中较小的那个(如果想等那就任取),例如,对于第一枚硬币的$(4963, 5037)$我们就取$4963$,然后我们可以计算“出现的正面次数小于或等于4963次的概率”:

$$ \begin{align} \mathsf{Pr}(Y \leq 4963) &= \sum_{i=0}^{4963} {10000 \choose i} p^{i}(1-p)^{10000-i} \\ &\approx 0.232696 \end{align} $$  

由于要处理的数字过于庞大,我们是使用Mathematica数学软件进行计算,计算的代码是

CDF[BinomialDistribution[10000, 0.5], 4963]

另外,鉴于$n$非常大,还可用正态分布$N(np, np(1-p))$来近似二项分布$B(n, p)$,这里,我们就用$N(5000, 2500)$来近似$B(10000, 0.5)$,其中$N(5000, 2500)$代表均值为$5000$方差为$2500$的正态分布:

CDF[NormalDistribution[5000, Sqrt[10000*0.5*0.5]], 4963]

结果是$0.22965$,和直接用二项分布累积概率函数算出的$0.232696$只相差一个百分位.

那么,现在我们已经知道了,对于一枚公平的硬币来说,对它抛掷10000次,正面(或者其中一面)出现的次数小于或者等于$4963$次的概率为$0.232696$,这个结果或许还不足以说明什么,但是,我们考虑减少这个数字,例如,我们考虑$\mathsf{Pr}(Y \leq k), k \in \mathbb{Z}, 0 \leq k \leq n$,即其中一面出现的次数小于$k$的概率,经过我们的计算:

Table[
    {
        4963 - d, 
        CDF[BinomialDistribution[10000, 0.5], 4963 - d]
    }, 
    {d, 0, 200, 50}
]

输出结果为:

{
    {4963, 0.232696}, 
    {4913, 0.0418126}, 
    {4863, 0.00316531}, 
    {4813, 0.0000955874}, 
    {4763, 1.11795*10^-6}
}

我们发现,当$k$值逐渐减小时,概率$\mathsf{Pr}(Y \leq k)$是快速地减小,我们看到,其中一面出现的次数小于等于$4913$的概率约是$0.0418126$,其中一面出现的次数小于等于$4863$的概率约是$0.00316531$,其中一面出现的次数小于等于$4813$的概率约是$0.0000955874$,而其中一面出现的次数小于等于$4763$的概率仅约为$1.11795 \times 10^{-6}$,到这里,我们再去看后面那一枚不公平的硬币,其中出现次数较小的是反面,出现了$3045$次,而我们知道不论反面或者正面,对于一个公平的硬币来说(假设它是公平的),在10000次试验下,出现的次数小于等于$4763$的概率是不到$1.11795 \times 10^{-6}$的,换句话来说,如果这枚硬币是公平的,那么抛掷它10000下,正面或者反面仅仅出现不到$4763$次是绝无可能的,更不要说仅仅出现$3045$次反面,因此,我们断定,这枚硬币铁定是有问题的而绝非巧合,如果是巧合,这个概率实在也太小了!就这样,我们就可以拒绝主张“硬币是公平的”的原假设而接受主张“硬币不是公平的”的备择假设.

像上面出现的这一些$0.232696$,$0.0418126$,……,$1.11796 \times 10^{-6}$这样的数字被称为$p$-值,可以简单地理解为原假设成立的概率(尽管这样有失严谨),一般来说,当这个$p$-值过小时,我们就会拒绝原假设并且接受备择假设从而破案,$p$-值多「小」才算是「小」呢?一般会给出一个固定的数字被称为「显著性水平」,比如说给定显著性水平$\alpha=0.05$,如果$p$-值小于$0.05$,那么我们就是备择假设的主张的$\alpha=0.05$显著的,当然了,我们还可以更加谨慎一些,把$\alpha$设得更低,例如$\alpha=0.01$,这样如若我们最终拒绝了原假设,那必然是因为我们收集到了及其有理的证据从而有了充足的理由和依据,从而将犯Ⅰ类错误(错怪原假设的错误)的可能性降到最低,从而最终做出既可信又可靠且有理有据的决策或者说判断.

总结

一开始我们说了,假如说我们得到的随机数生成器所能生成的随机数的分布是不确定的,比如说我们不能确定/dev/urandom是严格的均匀分布,或者说虽说/dev/urandom是均匀分布但是我们想要其他分布,我们可以一边继续用/dev/urandom生成随机数,一边随时进行假设检验,如果通过则接受之,不通过则拒绝之并且再生成一个随机数指定通过为止,按照这种方法,我们可以得到任意想要的分布,不管正态分布也好,二项分布也好,泊松分布也好,两点分布也好,都可以办到;我们还讲了分布和分布之间存在着一些互相转化的关系,例如中心极限定理指出,对任意分布进行抽样再求平均就可以得到正态分布,高尔顿板的案例也体现了如何由二项分布近似得到正态分布.而在后面,我们则详细介绍了假设检验的过程,具体是怎样检验给定的数据是否服从某种分布,我们先后分别是从两枚骰子(一枚公平骰子和一枚非公平骰子)和两枚硬币(一枚公平硬币和一枚非公平硬币)入手,公平的骰子可以看做是只产生整数的离散型的均匀分布(取各个值的概率相等),而公平的硬币可以看做是命中率为$0.5$的两点分布,根据骰子或者硬币的抛掷数据检验骰子或硬币是否公平,就相当于根据历史数据检验未知总体是否服从某种分布,其基本的原理的一样的.

概率论和数理统计方法描述生活和生产活动中的方方面面,是人类对大自然和广阔宇宙进行探究的有力工具,并且能够对随机过程进行准确的建模;随机数生成与生成的随机数的随机性在信息安全领域的工程实践中占用重要的地位,在该领域的工程实践中,密码学安全的随机数生成器和密码学不安全的随机数生成器是严格区分开的并且应用领域也有严格界定,这正体现了高质量的随机数生成器对于构建一个安全的系统来说是多么不可或缺.

这篇文章仅是作为我个人的学习笔记被撰写,其存在的意义在于概念的介绍而非严格且详实的讲解,出现错误或许在所难免,请读者辩证地以及宽容地看待.

参考资料

[1] Central limit theorem - Wikipedia

[2] Types of generators - The Rust Rand Book

仅供参考probability

小小的改进给探索子博客带来了大大的进步