CSDN专帖系列之一： 根据某一特殊规律的概率生成随机数

CSDN专帖系列之一： 根据某一特殊规律的概率生成随机数

Author：Jeff   2006-7-7

关键字：C++ 算法 随机数 概率

环境：Window XP Professional + SP2, VC6.0

 

CSDN原帖：

主　　题：求一算法，随机数和概率相关

作　　者：adintr (www.adintr.com)

链        接：http://community.csdn.net/Expert/topic/4753/4753180.xml?temp=.4426691

描        述：
有一整数 base ，其值大概在 10 - 100 左右， 现在需要生成一随机整数 x，其范围在 0 - n ( 0<n<17) 之间，要求这个随机数出现的概率为 ：
base ^ x / (base ^0 + base ^ 1 + base ^ 2 + ... + base ^ n)
也就是说如果假设生成的随机数为 0 的概率为 1 点的话，则生成的随机数为 1 的概率就是 base 点，为 2 的概率就是 base * base 点，为 3 的概率就是 base ^ 3 点。

注意 base 和 n 的取值达到较大时 100^16 是超过整数的表示范围的。
效率当然是越高越好， base 和 n 的值可以预先确定下来

 

在帖子中楼主已经对他的问题作了解答，他的回答算法是这样的：

首先生成 n 个 [0 ，base + 1) 之间的随机数，则每一位上出现 0 的几率是 1/(base + 1), 为非 0 的几率是 base / (base + 1)。

现在这 n 个随机数全为 0， 则让最后的结果 x = 0, 这种情况出现的几率是 1 / (base + 1)^n

如果第一位为非 0， 其余全为0， 则 x = 1, 出现的几率为 base / (base + 1) ^ n

如果前两位为非 0， 其余为 0， 则 x = 2, 出现几率为 base ^ 2 / (base + 1) ^ n

同理，当前 x 位为非 0， 其余为 0，则 x = x, 出现几率就是 base ^ x / (base + 1)^n

这样就保证了每以个数出现的几率都比小一的那个数大 base 点。

 

关于判断位为不为0时，是不是要求末尾的0连续？当时没有想清楚，后来明白了，这是一种排列情况。虽然有很多种组合方式，但是一旦判断条件定下来，情况也就确定了，每次都必须符合这种特定的排列，是作为组合的一种。因此对于末尾的0并不需要连续。可以这样理解：
        前两位为非 0， 其余为 0， 则 x = 2, 出现几率为 base ^ 2 / (base + 1) ^ n
        如果第一位为非0，末尾为非0，其余位为0，出现几率为
            [base/(base + 1)] * [1/(base + 1)] * ··· * [1/(base + 1)] * [base/(base + 1)]  
         = base ^ 2 / (base + 1) ^ n

两个几率是一样的。在程序中只处理一种情况就够了，且只能处理一种情况。

 

另外，当n值越大，程序就会递归得越深(因为连续随机到n个0的情况很少)。因此层次也不确

定，可能很深， 有可能栈崩掉。

 

        后来，经过分析，对我的算法进行修改，并加以实现。呵呵，发现速度不错。*_*

        经过修改后的算法：

        1． 先计算(0, RAND_MAX]内base的幂M，使得M满足K = base ^ M + base ^ (M - 1) + … + base <= RAND_MAX。(因为随机数最大为RAND_MAX)

        2． 计算各段区间的上限和下限。

        3． 如果n > M，直接到4；否则，跳到6。

4． 生成[0，K]范围内的一个随机数S。

5． 如果S ≠ 0, 计算S位于哪块区间，将其记为H, 则生成的随机数为H + n – M – 1。否则

n ß n – M, 跳到3。

6． 计算K2 = base ^ n + base ^ (n - 1) + … + base。

7． 生成[1，K2]范围内的一个随机数S2。

8． 计算S2位于哪块下限区间内，记为H2, 则生成的随机数为H2 – 1。

 

将base的幂作为分割点把整个线形空间(即[0，base ^ n + base ^ (n - 1) + … + base]空间)分成n + 1段，那么区间的分割状况就是：
          [0，0]，

 [1，base]，

 [base + 1，base ^ 2 + base]，

 ……，

  [base ^ (n - 1) + base ^ (n - 2) + … + base + 1，base ^ n + base ^ (n - 1) + … + base]

 则各段区间对应的比例为：
          1 / (base ^ n + base ^ (n - 1) + … + base + 1)

   base / (base ^ n + base ^ (n - 1) + … + base + 1)

   …….

   base ^ n / (base ^ n + base ^ (n - 1) + … + base + 1)

 和产生随机数的概率是一一对应的。因此，只需要随机数能够分布到整个线形空间，落在哪个区间，则生成了对应的随机数。

 但是base ^ n可能很大，不只rand()函数产生不了这么大范围的值，连int32都可能表示不了。

 

 这个时候，应该对区间进行适当的缩放，满足rand()函数能够产生收缩后的区间内的随机数，当然这个区间应该尽可能得大。计算一下区间的缩小比例，

 L = base ^ M + base ^ (M – 1) + … + base + 1应该尽可能的接近RAND_MAX，且比RAND_MAX小或者相等。在允许误差的情况下，整个线形空间[0，base ^ M + … + base + 1 )各个区间对应的产生的随机数为：

[0，1)： 0 < x < n – M

[1, base + 1)： n – M + 1

……

[base ^ (M – 1) + … + base + 1, base ^ M + base ^ (M – 1) + … + base + 1)： n

当随机数为0时，递归处理余下的线形空间，直到没有余下空间为止。用图例说明如下：

 

 

    下面的是代码及注释：

int  GetPower(long base)  // 也就是求上面所说的M。

{

    int    pow = -1; 

    long   sum = 1L , old = sum;

    long   total = 0L ;

    // 这里没有计算最后的+1，因为那是一个开区间。

    while (32767 > total && old <= sum) {

        old = sum;

        sum *= base;

        total += sum;

        ++pow;

    }   

   

    return pow; //返回base的幂。

}

 

long  power(long x, long y) // 计算x ^ y。(几个函数名都意义重了，不管了~~~)

{

    long   sum = 1L ;

    while (y-- > 0L ) {

        sum *= x;

    }

    return sum;

}

 

int main(void) {

    srand((unsigned)time(NULL));

 

// 下面的代码对下面宏的边界情况没有处理

#define  NUM_BASE    10  // 即base（不要大于RAND_MAX）

#define  NUM_POWER  6   // 即n（不要小于0）

 

    unsigned long  count[NUM_POWER];

    memset(count, 0, sizeof(count));

 

    int  pow = GetPower(NUM_BASE); // 获得幂M

    int  smaller = pow > NUM_POWER ? NUM_POWER : pow;     

   

    long  region[NUM_POWER + 1];

    long  lower_limit[NUM_POWER + 1];

    memset(region, 0, sizeof(region)); //区间上限

    memset(lower_limit, 0, sizeof(lower_limit));  // 区间下限

    for (int i = 1; i <= smaller; ++i) {       

        region[i] = region[i - 1] + power(NUM_BASE, i);  // 根据base计算上限

        lower_limit[i] = region[i - 1] + 1L ;       // 下限

    }   

 

    long  k = 100000000L ;

    while (--k >= 0L ) {

        int   n = NUM_POWER;

            while (0 < n - pow) {

            long seed = rand() % (region[smaller] + 1); //生成随机数[0， region[smaller] ]

            int  s = smaller + 1;

            while (--s >= 0 && lower_limit[s] > seed); //查找.

            if (s > 0) {  //找到了，（下限不是0）

                count[n - pow + s - 1]++;

                break;

            }

            n -= pow; //seed为0的情况，循环处理剩下的空间。

        }

        if (0 >= n - pow) {        //处理余下的一小段空间。

            long seed = (long)rand() % region[n] + 1; // [1， region[n] ]

            int  s = n + 1;

            while (--s >= 0 && lower_limit[s] > seed);

            if (s > 0) {

                count[s - 1]++;               

            }

        }

    }

 

    for (k = 0L ; k < NUM_POWER; ++k) {

        printf("%d: %d/n", k, count[k]);

    }

 

        return 0;

}

 

 

第一次的测试结果：（5秒）

0: 847

1: 8344

2: 91541

3: 915915

4: 9153458

5: 89829895

 

第二次的测试结果：（5秒）

0: 843

1: 8316

2: 91336

3: 915133

4: 9155684

5: 89828688

 

  可以看出规律大概是对的，只是0，1和其他数值相比，比例降低得更快。因为0和1合起来才占据了区间的一个单位，现在把0和1对应的值加起来差不多900，比较符合。

  理论上来说，楼主的算法接近真实值一些，但是速度太慢。而这个算法相对来说，误差大，速度快。时间基本上都花在了最外围的循环上。