Sampling 随即取样问题归纳

《编程珠玑》第12章读书笔记：
这里把常见的随机取样问题归纳了下，以下所有算法和程序都经过实际测试过，如果有问题，请留言给我。转载请注明出处。

大纲：

问题1  如何随即产生1-N之间一个随机数K.

问题2  从N中随即取K个，N已知。

问题3  从N中随即取K个，N未知。

问题4  从N中随即取1个，N未知。

问题1,   如何随即产生1-N之间一个随机数K. rand[1,N]

          简单做法是,通过调用库函数

rand()% N + 1

但是这样并不严格，因为rand()产生的随机数在[0, RAND_MAX]之间，而RAND_MAX值一般不大，比如在我机器上RAND_MAX为32767，因此，如果需要产生一个大随机数的话则需要：

rand() * RAND_MAX + rand()

 

问题1扩展： 如果给随即函数rand[1,N] 如何产生rand[1,M]

      1.1 如果M< N, 则对rand[1,N] 直接舍去大于M的数即可。

      1.2 如果M> N, 例如已知rand[1,5], 求rand[1,7].

           解法一：

int rand7() {     while (1)    {      int num = 5*(rand5() -1) + rand5（）- 1；      if (num < 21) return num % 7;    } } 

           解法很多，这里给个我的版本：

     step1 [1, 5] ---> [0, 1] (大于3为1，小于3为0，等于3重做）

     step2 [0, 1] ---> [1, 7] (根据[0,1] 分别产生3位二进制的每一位，如果全为0，重做）

#include <stdio.h> #include <stdlib.h> /** * rand[1, 5] **/ int rand5() {     return rand() % 5 + 1; // 0-4 + 1 = 1- 5 }  /** *rand[0, 1] **/ int rand1() {     int temp;     while( (temp = rand5())  == 3)     {         temp = rand5();     }    if(temp > 3) return 1;    else return 0; }  /** *rand[1, 7] **/ int rand7() {     int t1 = rand1() ;     int t2 = rand1();     int t3 = rand1();     while( t1 == t2 && t2 == t3 && t3 == 0)     {              t1  =  rand1();              t2  =  rand1();              t3  =  rand1();     }     return (t1<<2) + (t2<<1) + t3; }  int main() {    int i;      float n1 ,  n2 , n3 ,n4 , n5, n6 , n7;     n1 = n2 = n3 = n4 = n5 = n6 = n7 =0;      int temp;     for( i = 0; i < 1000000; i++)     {         temp = rand7();         if(temp == 1)       n1++;         else if(temp== 2)  n2++;         else if(temp== 3)  n3++;         else if(temp== 4)  n4++;         else if(temp== 5)  n5++;         else if(temp == 6) n6++;         else if(temp == 7) n7++;         else printf("error/n");     }     int num = n1 + n2 + n3 + n4 + n5 + n6 + n7 ;     printf("1/7 = %0.3f/n", 1.0/7);     printf("%d 次rand[1,7]概率为/n1: %0.3f  2: %0.3f  3: %0.3f 4: %0.3f  5: %0.3f  6: %0.3f  7:%0.3f/n", num, n1/num, n2/num, n3/num, n4/num, n5/num, n6/num, n7/num);     return 0; }

运行结果为：

1/7 = 0.143

1000000 次rand[1,7]概率为

1: 0.144 2: 0.143 3: 0.143 4: 0.143 5: 0.143 6: 0.142 7:0.143

 

问题2       从N中随即取K个N已知

 

方法一：TAOCP 第二卷3.4.2中的S算法, 例如从N=5中取K=3个。选择1个概率为3/5，当选择了1后，选择2的概率为2/4；如果1没有被选中，选择2个概率为3/4。

int randNK(int N, int K) {     int selected = 0;      int i;     for( i = 0; i < N; i++)     {         if( (K - selected) >= (rand() % (N-i) + 1 ) )         {             printf("%d ", i+1);             selected++;         }     }     return 0; }

 

方法二：类似于从一个装有N个球的盒子里抽取K个球，首先产生[1, N]的随机数m，如果m没有被选中过，则选取M,否则重新选取。

void randNK(int N, int K) {     set S;      while(S.size() < K)         S.insert( rand() % N + 1 );      set::iterator i;     for( i = S.begin(); i != S.end(); i++)     {         cout<< * i << " ";     }     cout<<endl; }

 

当要选取的K很接近N的时候，方法二的效率会受到影响。比如从N=100个取K=100个，当已经选取了99个时候，产生第100个满足条件的概率为(100-99)/100=1/100,也就是平均需要100次才能得到满足条件的第100个随机数。

Bob Floyd发现了这个问题，并提出了方法二的改进方法：

方法三: Folyd采用递归的思想，为了从N中选取K个，则先从N-1中选取K-1个，再从N个中选取第K个。

递归下去，从N-2选K-2个，…，从N-K+2中选取2个，从N-K+1中选取1个。

很容易看出每个元素被选中的概率为1/(N-K+1)。

非递归的写法：

"for( j = N - K + 1; j <= N; j++) 
{ 
S.insert(rand()%(j+1) + 1); 
}" 
这个 for 循环只能执行 K 次，但是那个插入语句好像不能保证每次都插入不同的数，也就是说最后得到的随机数的个数可能小于 K。根据 Floyd 的想法，我认为应该是这样的： 
for( j = N - K + 1; j <= N; j++ ) { 
int t = rand() % (j - 1) + 1; 
if( S.find(t) == S.end() ) { 
S.insert(t); 
} else { 
S.insert(j); 
} 
} 

上面的 t 是落在 1 到 j - 1 范围内的，肯定不会等于 j。因此，上面的方法能保证 for 循环的每次执行都会插入一个不同的数，从而保证最后得到的是 K 个数。

方法四：洗牌算法

洗牌算法就是把一个数组全部打乱.

void swap(int * a, int * b) {         int temp = *a;         *a = *b;         *b = temp; } void shuffle(int * array, int n) {         for( int i = 0; i < n; i++)         {                 swap(array+i, array + rand()%(n-i));         } }

如果洗牌算法中只打乱前面K个元素，则成了从N个随即选取K个

void rand_shuffle(int N, int K) {         int * array = new int[N];         for(int i = 0; i < N; i++)         {                 array[i] = i+1;         }         for(int i = 0; i < K; i++)         {                 swap(array+i, array + rand()%(N-i));         }          //print         for(int i = 0; i < K; i++)         {                 cout<<array[i] <<" ";         }         cout<<endl ; }

 

问题3       从N中随即取K个N未知

 

这个问题我以前写在javaeye的,这个解法同样适用于问题2:

reservoir sampling
//Init : a reservoir with the size： K
for(int i = K; i < N; ++i){
   int M = rand() % i;
   if( M < K)
      SWAP the Mth value and ith value
}

解决方案就是蓄水库抽样（reservoid sampling）。主要思想就是保持一个集合（这个集合中的每个数字出现），作为蓄水池，依次遍历所有数据的时候以一定概率替换这个蓄水池中的数字。

（1）初始情况。出现在水库中的k个元素的出现概率都是一致的，都是1。这个很显然。
（2）第一步。第一步就是指，处理第k+1个元素的情况。

分两种情况：1,元素全部都没有被替换；2，其中某个元素被第k+1个元素替换掉。
对于情况2：第k+1个元素被选中的概率是k/(k+1)（根据公式k/i），所以这个新元素在水库中出现的概率就一定是k/(k+1)（不管它替换掉哪个元素，反正肯定它是以这个概率出现在水库中）。下面来看水库中剩余的元素出现的概率，也就是1-P(这个元素被替换掉的概率)。水库中任意一个元素被替换掉的概率是：(k/k+1)*(1/k)=1/(k+1)，意即首先要第k+1个元素被选中，然后自己在集合的k个元素中被选中。那它出现的概率就是1-1/(k+1)=k/(k+1)。可以看出来，旧元素和新元素出现的概率是相等的。
情况1：当元素全部都没有替换掉的时候，每个元素的出现概率肯定是一样的，这很显
然。但具体是多少呢？就是1-P(第k+1个元素被选中)=1-k/(k+1)=1/(k+1)。

（3）归纳法：重复上面的过程，只要证明第i步到第i+1步，所有元素出现的概率是相
等的即可。

当把K取为1的时候，就得到了问题3

问题4  从N中随即取1个，其中N未知

把问题三中K取1，即第K个以1/K的概率选取。

 

参考文献:

《编程珠玑》 第12章，《编程珠玑2》 第13章，《计算机程序设计艺术》第二卷 第三章

come from http://wansishuang.appspot.com/?p=59003#Q1