2. 直接法——利用rand7做计算
由于舍去法每次调用rand7的次数未知,所以希望能够找到一种直接的方法。当然最直接的方法就是用线性同余这样的随机数生成法直接写一个,但是这就用不上rand7了,与题意不符合。有人希望用rand7的各种组合计算来完成,比如(rand7+rand7+rand7+...)%10+1或者(rand7*rand7*rand7*...)%10+1,有些计算确实能够达到很好的均匀度,但却是近似均匀。
从统计学的角度看, 一个rand7就相当于符合均匀分布的随机变量X1,当n个rand7做运算的时候,相当于X1...Xn是符合均匀分布的边缘随机变量,而F(X1, X2,...,Xn)是他们的联合分布,这个分布并非是均匀的,甚至很复杂,如果不用穷举法,几乎没有简便的方法计算每个数字出现的概率。现在要将多随机变量的F(X1,X2,...,Xn)映射到1到10的均匀分布,显然是有一定难度的, 而在本题来说,是不可能的,这是因为:
1) 对rand7的一元运算只能有7种结果,不可能产生10个随机数
2) 现在有二元运算(X)可以是加减乘除或者任何函数任何映射关系,rand7(X)rand7的可能运算方式是7*7种,,n次二元(X)运算后的可能是运算方式是7^(n+1)种,现在要用7^(n+1)种运算过程得到均匀的10种结果,这是不可能的,(因为7^(n+1)不能被10整除),所以只能是近似均匀。
下面来看怎样获得近似的均匀,
ju136提醒我注意到,其中的一个比较好的方法是:
(rand7()+rand7()+rand7()+rand7()+rand7()+rand7()+rand7()+rand7()+rand7()+rand7())%10+1
他获得的均匀度非常好,1出现的最多,5出现的最少,但是概率上仅仅相差0.00002,人类的感觉已经分辨不出来了。但是这个方法需要调用10次rand7,效率上差一些。
有人希望用这样的方法:调用两次rand7从而生成一个7进制的数,然后转换成0-49,刚好是50个数的均匀分布,再取模10。这个方法貌似可行,可是很遗憾的是,这样生成的7进制0-66对应到10进制是0-48,而不是49,少了一个数。
下面这个方法也比较好,(rand7+(rand7+7) +(rand7+14)+...+(rand7+42))%10,这个表达式生成7个随机数,分别均匀分布在1-7,8-14,...7个区间,相加之后再做模10运算,映射到0-9这10个数字。这7^7种运算,统计每个数字可以得到次数,其中5最高,0最低,但是他们几率的差仅为0.00041,人的感觉几乎分辨不出来了。这个方法需要调用7次rand7,效率比上面的代码高一些,因此有:
代码3
int rand10_7plus()
{
return (rand7()+rand7()+rand7()+rand7()+rand7()+rand7()+rand7()+147)%10+1;
}
这种直接方法,无论怎样在理论上都做不到均匀分布,所以,我们要想一些办法来提高。
3. 直接法——利用组合方法进一步提高
考虑数组
{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
{6, 5, 4, 3, 2, 1, 10, 9, 8, 7 }
用rand7_plus在第一行中选择的时候6的机会最大,1的机会最小,但是在第二行选择的时候1的机会最大,6的机会最小,这两行有一定的互补性,所以如果轮流选择这两行,会得到更均匀的分布,根据上面的计数结果,可以得知,最小几率与最大几率的差仅为万分之一。推广一下,如果在数组
int seed10[10][10] = {
{1, 2, 3, 4, 5, 6, 7, 8, 9, 10},
{10, 1, 2, 3, 4, 5, 6, 7, 8, 9},
{9, 10, 1, 2, 3, 4, 5, 6, 7, 8 },
{8, 9, 10, 1, 2, 3, 4, 5, 6, 7 },
{7, 8, 9, 10, 1, 2, 3, 4, 5, 6 },
{6, 7, 8, 9, 10, 1, 2, 3, 4, 5 },
{5, 6, 7, 8, 9, 10, 1, 2, 3, 4 },
{4, 5, 6, 7, 8, 9, 10, 1, 2, 3 },
{3, 4, 5, 6, 7, 8, 9, 10, 1, 2 },
{2, 3, 4, 5, 6, 7, 8, 9, 10, 1 },
};
之中,依次在每一行用rand10_7plus选择,由于每个数字出现的几率均等,并且在每行和每列上出现的几率均等,在多次调用之后就可以得到均匀分布。
代码4
unsigned int gi = 0;
int rand10_matrix()
{
//用rand10_7plus生成一个数,选择列
int n = rand10_7plus() - 1;
//轮流选择seed10的行
return seed10[gi++ % 10][n];
}
如果利用模运算的循环性质,就不需要seed10这个矩阵,可以验证下面的代码与rand10_matrix是等效的:
代码5
unsigned int gi = 0;
int rand10_mod()
{
return ((rand10_7plus() - 1) + gi++) % 10 + 1;
}
这个算法每次调用7次rand7,做一次模运算,不需要额外的内存以及循环,从统计意义上说,已经是个比较好的伪随机数生成器。
4. 随机数测试
最初我在些这个文章的时候,在测试方面出了问题,所以在这里强调两点:
1). 随机度测试,需要每隔10次调用计数一次,以验证每个数字在各个位置上出现的次数是均等的。
代码要与如下类似:
for (k = 0; k < 100000; k++) {
n = 0;
for (j = 0; j < 10; j++)
n = rand10_7plus10();
result[n - 1]++;
}
注意result[n - 1]++是在j=0~10这个循环之外的。
2). 概率计算
在编写代码之前,一定要先证明1-10的均匀分布。例如前面10个rand连加的算法,在人的感觉之中已经分辨不出来概率的差别了,因此需要仔细统计一下10个数字,写一个简单的10层循环计数就可以了。很多朋友忽视了这点,并且,我们上面在第2小节证明了这是不可能做到均匀的。如果想不清楚或者很难证明你的方法,于是这个最简单的代码就非常有用:
for (i1 = 1; i1 <= 7; i1++)
for (i2 = 1; i2 <= 7; i2++)
for (i10 = 1; i10 <= 7; i10++)
result[(i1 + i2 + + i10) % 10]++;
5. 直接法——分布变换与连续随机变量的分布——更实际的应用
除了这道题的一些技巧,题目本身在实战中没有任何应用,比较实际的问题是,假设一个班的学生成绩符合正态分布,如何模拟生成考试成绩。
我们先看如果rand7是1-7的连续均匀分布,如何获得1-10的均匀分布。答案很简单,从几何的角度上看,我们可以把[a,b]线段上的点按照一对一映射到另一个线段[c,d]上去,只需要做一个线性变换y=(x-a)/(b-a)*(d-c)+c. 那么,若x=rand()~U(a,b),则y=~U(c,d),也就是如果rand()是a到b上的均匀分布,则y=(d-c)(x-a)/(b-a)+c是c到d上的均匀分布。对于本例rand10=(rand()-1)/6*9+1. 下面是证明,更一般的情况同理可证:
另外有一个重要的定理来表明变换之后的分布。这可以处理如Y=X^2, Y=e^X等多种变换。定理如下:
这个定理还可以更强一些,f(x)是分段还是也可以,甚至只是一个覆盖(包括)就可以了。从符合一种分布的随机数生成另外一种分布的随机数是统计模拟的课题,其中有非常有趣的变换方法,例如,如果X是(0,1)上的均匀分布,则Y=-a*log(X)是指数分布。
现在来回答如何按照正态分布模拟生成一个班的学生成绩。这个方法被称为Box-Muller算法,如果U1和U2服从 (0,1)区间的均匀分布,做变换R=sqrt(-2*logU1), alpha=2*pi*U2,则X=Rcos(alpha)和Y=Rsin(alpha)是一对独立的 标准正态分布n(0,1)。证明从略。按照这个方法,C语言的rand(),可以模拟生成近似的(0,1)均匀分布,只要rand()/MAX_RAND就可以了,做Box-Muller变换到n(0,1),就可以做出符合学生成绩分布了,概率统计课的基础内容就有。
6. 其他方法
舍去法也是非常重要的一类随机,用来生成各种分布的随机数,另外的方法:比如Metropolis算法,比较著名的还有Markov Chain Monte Carlo (MCMC)算法,这类方法可以看成是一个黑盒子,要求在算法内部通过几次运算很快收敛到一种概率分布,然后返回一个随机数。
7. 参考文献
[2]
Kunth第2卷Seminumerical Algorithms, Random Numbers.