编程珠玑第一章生成随机数随机取样的问题

Posted 2020-10-21 小僧尤里

一、利用随机数函数生成随机数

问题1（《编程珠玑》习题12.1后半段）：

给定一个rand()，可以产生从0到RAND_MAX的随机数，其中RAND_MAX很大（常见值:16位int能表示的最大整数32767），写出利用rand()生成[a,b]中任意整数的函数，其中a>=0, b<=RAND_MAX，且b-a<<RAND_MAX。

分析：

这是在编程工作最常见的随机函数的应用，在这里做一个起点再合适不过。把随机数区间的起点从0变为a，同时把一共RAND_MAX+1个数的区间缩小至只含有b-a+1个数的区间，写为 a + rand()%(b-a+1)，此时显然最大值是a+(b-a)= b。进一步地，这个b-a<<RAND_MAX的条件虽然看上去不起眼，其实很重要。

附加思考：

如果b-a和RAND_MAX很接近会发生什么情况？读者不妨先做思考，问题2的分析会做出解答。这个rand()，其实相当于《编程珠玑》提到的bigrand()。

 

问题2（笔试题变形）：

给定一个随机数函数rand7()，它能以等概率生成1~7这7个整数。请根据rand7()写出类似的rand5()。

分析：

如果直接像问题1中一样，把1+rand7()%5作为rand5()会有什么情况发生？这时确实能产生1~5的随机数没错，可是各个数的概率相等吗？

对于随机数2，既有可能来自于1+1%5，也有可能来自于1+6%5，显然其概率是2/7，而不是1/5，不满足rand5()等概率产生各随机数这一隐含要求。不同于问题1，问题1中一个很大的区间收缩成较小区间时，各个元素映射后的新元素概率虽然概率可能不完全一样，但却是近似相同的。

为了满足等概率的要求，可以这么做：

int rand5() {
    int res;
    do {
        res = rand7();
    } while(t>=6);
    return res;
}

虽然保证了1~5的概率都变成了1/5，但是有一个无法避免的缺点是，每当产生了6或者7都要抛弃，相当于这一次运行是“空转”，浪费了时间。如果对1/5这个概率不明白，可以有两种理解：每次产生6或7就被抛弃，剩下数的概率相等，必然为1/5；或者用更严密的推理：产生1~5的随机数，最终得到某一个的概率为：1/7+(1/7)*(2/7)+(1/7)*(2/7)2+...，无限项求和，结果是1/5。

问题3（笔试题原题）

给定一个随机数函数rand7()，它能以等概率生成1~7这7个整数。请根据rand7()写出类似的rand10()。

分析：

有了问题2的概率基础，把rand7()变成rand10()仅仅需要一点点思考了。

int rand10() {
    int t1,t2,res;
    do {
        t1 = rand7();      //t1以等概率1/5成为{1,2,3,4,5}中的一个数字
    } while(t1>=6);
   
    do {
        t2 = rand7();    //t2以1/2概率成为奇数或偶数，因为{1,2,3,4,5,6]中3奇3偶
    } while(t2==7); 
    res = t1+5*(t2%2);       //res是1~10中的任意一个数的概率都是1/10，注意到%和*具有相同的优先级，这里去掉括号结果是错的    
return res;
}

解释：

t1的概率是{1~5}中每个元素概率为1/5，t2是1/2概率为奇或者偶，所以(t2%2)是以1/2概率选择0或者1，因此5*(t2%2)是以1/2概率为0或者5。这样，选择1的概率是1/5 * 1/2=1/10，...，选择6的概率是1/5 * 1/2为1/10....由此可见。

问题2和问题3是对《编程珠玑》上使用范围不大的randint(i,j)生成其他范围随机数方法的解答。在掌握了问题2和问题3的解法后，你已经学会随机数区间的收缩和扩张，类似的问题迎刃而解。

问题4（《编程珠玑》习题12.1前半段）

C库函数rand()常返回15个随机位，写出bigrand()，能够返回30个随机位。

分析：

其实问题4和问题3有点像，但是不同之处在于，这次我们的视角是从位出发的，把rand()看做了将15个位每一个位以1/2概率设为0或1，从而生成0~RAND_MAX。从某种意义上来说，按这种理解方式来解这个问题更轻松一些，但是仅局限于2的幂减1这样的数值的区间，比如从0到11...11。在此基础之上把两部分的位拼接起来即可。

//《编程珠玑》的答案
//int bigrand() {return RAND_MAX *rand() +rand();}
//最大值不是30个1，在另一个笔记中有叙述，怀疑有错我的答案，把先生成的部分左移15位作为高15位
long bigrand(){return ((long)RAND_MAX+1)*rand() +rand();}

 扩展：拒绝采样

这是一个简单粗暴有效的方法，使用这个方法你可以不用考虑复杂的区间扩展和收缩的问题。以rand7()产生rand10()为例，构造一个二维数组并把它填成：

然后使用两次rand7()分别生成行号和列号，如果对应元素是0，抛弃重来；否则即是结果。其实这种方法依然用到了区间扩展：把7扩展成7*7，并把不符合的部分抛弃。这样横向生成一个1~7的随机数，纵坐标生成1~7的随机数，这样两次随机数便能定位一个数组元素，若该元素为1~10之间的数，则被选中，否则重新选择。这能保证1~10之间的每个数被选中的概率都是1/10。从这里可以看出这种方法其实还是有缺陷的：如果用rand7()生成rand50()，那这50个数可是填满这个表格后还有填不下的。

 

二、利用随机数函数产生随机事件

用随机事件表示随机事件？这个问题具体化为：使用rand()表示以M/N的概率发生的随机事件，其中M<=N，并可用作：if(事件A发生)  ，其中P(A) = M/N，那么表示为：

if(rand()%N< M)

    ...

通过这种方式，我们就可以做出让程序“以M/N的概率执行某个命令”这样的设计了。

 

三、取样问题:从n个元素中随机选取m个

 

1.从概率角度出发

考虑整数0,1,2,...,n-1，可以用上节的方法以概率m/n选取0（推导方式略）。但是对于1，必须考虑之前0是否被选取而以(m-1)/(n-1)或m/(n-1)的概率选取1，后续就更加麻烦。好在迭代是计算机的长项，只需要把这个是否选择当前数的随机事件稍作修改即可，使之变成从r个其余的数中选择s个：

int gen(int m,int n){
    int i, select = m,remaining = n;
    for(i=0;i<n;i++) {
        if(rand() % remaining <select) {
            printf("%d\\n",i);
            select--;
        }
        remaining--;
    }
    return 0;
}

其概率证明可见于Knuth的The Art of Computer Programming第2卷。进一步地可以优化为

int genknuth(int m,int n){
    int i;
    for(i=0;i<n;i++)
        if(rand()%(n -i) < m) {
            printf("%d\\n",i);
            m--;
        }
    return 0;
}

《编程珠玑》提示，这种算法是“所有n的所有m元子集被选中的概率相等”，这个条件强于“所有元素被选中的概率相同”。下面是习题12.2中提到的“以等概率选择所有元素，但是有些m元子集被选中的概率比其他子集大”的算法：直接选择1个数，则这个m元集合为它本身即后续的一共m个数，可能包括回绕。

对于这种方法，总要产生n次随机数。进一步可以写为for(i=0;i<n&&m>0;i++)，但程序的平均运行时间是否变得更快需要权衡。对应地，习题12.7提供了一种稍微减少随机数产生次数的递归函数：

int randselect(int m,int n) {
    int r;
    //assert(m<=n && m>=0);
    if(m>0) {
        r = rand()%n;
        if(r < m) {
            printf("%d\\n",n-1);
            randselect(m-1,n-1);
        }
        else
            randselect(m,n-1);
    }
    return 0;
}

 

2.从集合插入出发

由于集合元素不重复，如果按等概率选择一个随机数，不在集合中就把它插入，反之直接抛弃，直到集合元素个数达到m个，同样可以满足要求，并且用C++的STL很容易实现：

void gensets(int m,int n) {
    set<int> S;
    while(S.size() < m)
        S.insert(rand()%n);
    set<int>::iterator i;
    for(i = S.begin();i!=S.end();++i)
        cout<<*i<<"\\n";
}

这个算法的主要问题是，如果抛弃已存在的元素的次数过多，相当于多次产生随机数并进行集合操作，性能将明显下降。比如当n=100而m=99，取第99个元素时，算法“闭着眼睛乱猜整数，直到偶然碰上正确的那个为止”（《编程珠玑（续）》，13.1节）。虽然这种情况会在“从一般到特殊”提供解决方案，但下面的Floyd算法明显规避了产生随机数超过m次的问题。

习题12.9提供了一种基于STL集合的随机数取样方法，可以在最坏情况下也只产生m个随机数：限定当前从中取值的区间的大小，每当产生重复的随机数，就把这一次迭代时不会产生的第一个随机数拿来替换。

int genfloyd(int m,int n){
    set<int> S;
    set<int>::iterator i;
    for(int j = n-m; j<n;j++) {
        int t = rand()%(j+1);
        if(S.find(t) == S.end())
            S.insert(t);
        else
            S.insert(j);
    }
    for(i=S.begin();i!=S.end();++i)
        cout<<*i<<"\\n";
}

不必基于集合而实现，我自己写了一个原理类似的算法，思想是把“用过”的元素“扔”到下一次迭代的随机数取样区间之外：

int gen(int m,int n) {
    int *array;
    int i,j;
    array = (int *)malloc(sizeof(int) * n);
    for(i=0;i<n;i++)
        array[i] = i;
    while(m>=1) {
        j = rand()%n;
        printf("%d\\n",array[j]);
        if(j<n-1)
            array[j] = array[n-1];
        m--;
        n--;
    }
    return 0;
} 

3、从“打乱顺序”出发

这是个来源于实际的想法：将所有n个元素打乱，取出前m个。更快的做法是，打乱前m个即可。对应的C++代码如下：

int genshuf(int m,int n){
    int i,j;
    int *x = new int[n];
    for(i = 0;i<n;i++)
        x[i] = i;
    for(i = 0;i<m;i++) {
        j = randint(i,n-1);  //randint产生i到n-1之间的随机数
        int t = x[i];x[i] = x[j];x[j] = t;
    }
    //sort(x,x+m);            //sort是为了按序输出 
    for(i=0;i<m;i++)
        cout<<x[i]<<"\\n";
}       

把sort注释掉的这段代码，可以作为随机不重复序列产生器。类似的还有Floyd的算法P。（《编程珠玑（续）》，13.3节）

4、从一般到特殊

以上讨论的几种方式都不限定m和n的取值，只需m<=n即可。对于特殊的取值，有特殊的解决方案，以下是编程珠玑上的两例：

　　1.n = 106而m=n-10，这时m很大较为接近n，可以先生成10个元素，然后输出其余的元素。进行这种处理的额外代码可以提高算法的平均速度。

　　2.n=231而m = 107，m<<n，这时可以先生成1.1*107个元素，排序后去掉重复的（由于n很大，m中出现重复元素的概率很低），得到 107个元素的样本。

附：

（《编程珠玑（续）》习题13.5）Doug McIlory的求N个元素中取M种的第G种情况的组合的算法，原书这个算法我没有理解，也没有看到比较满意的解释，可能用到了某些我不熟悉的组合数性质。原书中没有对其进行进一步解释，并且似乎只是作者的题外话而已。如果看着有困难，这部分代码完全可以跳过。介绍这个算法的原因是用它能在获得随机数G后，直接获得第G种N个元素取M个的取法，相当于只产生了一次随机数。

int comb(int n,int m,int g,int* array){
    int d = 1,t;
    while(m>0) {
        t = combination(n-d,m-1);
        //combination(n,m)是计算n中取m个的组合数
        if(g-t<0) {
            m--;
            array[m] = d;
            printf("%d\\n",d);
        }
        else
            g-=t;
        d++;
    }
    return 0;
}

三、取样问题:从未知总数的元素中选择一个

从事先未知总数为n的对象中随机选择一个的方法。有两种常见的具体问题：

　　1.读取一个未知行数的文件，随机输出其中的的一行，同时最多只能缓冲一行的内容（《编程珠玑（续）》习题15.3利用了这种形式）；

　　2.对链表进行一趟遍历，随机输出其中的一个结点的元素，只能使用一个临时指针。

解法是，以概率1选择第一个元素存入缓存，以概率1/2用第二个元素替换掉缓存，...，直至遍历完所有元素，最后输出缓存的内容。可以分析出此时所有元素留在缓存的概率均为1/n，比如1，是1*(1/2)*(2/3)*...*((n-1)/n)。

int random_select(void){
    int i,num=1;
    for(i=1;i<n;i++)
       if(rand()%i ==0)     //等于0或者任何一个0~i-1之间的数，表示概率1/i。
           num = i;         //以1/i的概率替换num的缓存值。
    return num;
}

对于这个问题的常见应用：马尔科夫文本生成器里选择哈希表中某项对应链表的任意一个结点。