等概率数字筛选问题

如何等概率的从N个元素中选取出K个元素？

从1....n中随机输出m个不重复的数。(迅雷2011.10.21笔试题)

算法1

[cpp] view plain copy

1. knuth(int n, int m)  
2. {  
3.        for (int i=0; i<n; i++)  
4.        {  
5.                  if ( rand_n()<m ) // rand_n()为生成[0,n）这个区间内的随机数  
6.                        // rand_n()<m的概率就是m/n  
7.                   {  
8.                       cout<<i<<endl;  
9.                    }  
10.           }  
11. }  

算法2

[cpp] view plain copy

1. knuth(int n, int m)  
2. {  
3.        srand((unsigned int)time(0));  
4.        for (int i=0; i<n; i++)  
5.        {  
6.     //假设还要从后面n-i个数中选m个元素  
7.     //那么每个元素被选中的概率为rand()%(n-i)<m  
8.     //因此每个元素的选中概率是等价的  
9.                  if (  rand()%(n-i)<m )  
10.                   {  
11.                       cout<<i<<endl;  
12.                       m--;  
13.                    }  
14.           }  
15. }  

如何等概率的从数据流S中选取出K个元素？

由于数据量大小未知，因此无法使用上面的算法。我们能做的仅仅是保证所有已经输入的数据（假设已经有N个了）已k/N的等概率被选出，从而满足需求。算法如下：

[cpp] view plain copy

1. S为元素数据流  
2. 初始化: 初始化一个大小为K的元素空间buf  
3.     for i =1 to k  
4.         buf[i] = S[i]  
5.     for i= k+1 to N  
6.     {  
7.         //第i个元素以 k/i 的概率留下  
8.             M=random(1, i);//生成一个1到i的随机数  
9.             if( M <= k)  
10.                buf[M]= S[i]  
11.     }  

原理：

假设当前是i+1, 那么i+1这个元素被选中的概率是k/i+1。考虑前i个元素，如果前i个元素被选中的概率都是k/i+1，则算法满足要求。 

证明：

a. N<=k 时 元素全保留，留下来的概率 1.

归纳证明：k < i <=N 

1. 当i=k+1时

   第k+1个元素被选中的概率为k/(k+1)。
   从K个位置中随机选出一个位置放入新数据。有1/k的概率被选出。
   需要该操作的概率是k/(k+1),因此有1/k *  k/(k+1) = 1/(k+1)的概率被选出。
   故k/(k+1)的概率保留在缓存中。即就是前面i个元素和第i+1个元素等概率的存在于缓存中，出现的概率为k/i。

2. 假设当 j=i 的时候结论成立，此时以 k/i 的概率来选择第i个元素，前i-1个元素出现在蓄水池的概率都为k/i。
3. 当j=i+1时：

   元素j留下的概率为k/j =  k/i+1;

   从K个缓存中为元素j找一个位置放入K个元素有1/k的概率被替换

   因此 本次操作 K个元素被替换的概率为（k/i+1） *（1/k） = 1/(i+1),

   本次操作K个元素能保留下来的概率为 1 - 1/(i+1) = i/(i+1

   )此前，我们已经从i个元素中等概率k/i的选了k个元素，

   经过这次操作这K个元素有i/(i+1)保留

    因此有k/i * i/(i+1) = k/（i+1）

综上所述，证明成立。

扰乱一个递增序列。 还未想明白

[cpp] view plain copy

1. for i =[0,N)  
2.     swap(x[i],x[rand(i,n-1)];  

有人证明，只要扰乱前m个就可以。

[cpp] view plain copy

1. void sample_shuf(const int N,const int m)  
2. {  
3.     int i, j;  
4.     int *x = new int[N];  
5.     for(i = 0 ; i <N ; i++)  x[i]=i+1;  

[cpp] view plain copy

1. for(i = 0 ; i < m ; i ++)  
2. {  
3.     j = rand(i,N-1);  
4.     swap(x[i],x[j]);  
5. }  
6. sort(x,x+m);  
7. Print(x,m);  
8. delete []x;  

[cpp] view plain copy

1. x= NULL;  

关于采样的几个问题：

1、生成[1,5]整数的发生器，如何生成区间在[1,7]的发生器?
解答：利用拒绝采样定理
       首先，将（1，5）之间的随机发生器使用两次，按照五进制进行使用，拼成一个[0,24]的随即发生器既：（[gen -1][gen -1]）5,每一[]为一个5进制上的位，换算为十
进制为：x=(gen-1)*5+(gen-1) +1,在十进制上的范围为：1-25; 然后将（1,25）平均分配到7中情况上面，考虑21是7的倍数，因此可以每三个做一个映射（当然，也可以不管，直接截断7后面的数字，但是范围太小，效率不高），1-3--》1,4-6--》2,19-21--》7，此时就是等概率的，如果产生了22-25之间的数字，可以有两种方法决定结果：
     （1）拒绝采样，重新再运算
     （2）如果得到了22-25之间的数字，则此次的随即发生器结果，直接使用上一次得到的结果。这个方法有人证明过，是等概率的，算法Metropolis Algorithm。

2、Generate a random permutation for a deck of cards
解答：
        从后往前，第k步的时候，随机产生一个1 到 k，之间的数字j，然后交换j和k处的数字，可以很容易的最后这个排列就是一个等概率得到的排列。

[cpp] view plain copy

1. for k=N:1  
2.     j = rand(1,k)  
3.     swap(j,k)  
4. end  

      同样的，也可以从前往后进行这个过程，不过产生的范围就是变成k-N之间了。

[cpp] view plain copy

1. for k = 1:N  
2.      j = rand(k,N)  
3.      swap(j,k)  
4. end  

参考资料

1. http://blog.sina.com.cn/s/blog_6344728a0100h6pz.html
2. http://blog.csdn.net/v_july_v/article/details/6880698
   
3. http://blog.csdn.net/hackbuteer1/article/details/7971328