解决等概率随机抽样问题

 本篇讨论蓄水库抽样算法 的 原理及其实现

最近在优快云技术贴区看到一个帖子讨论这个问题：

问题：100个苹果完全随机分给4人，每人可能得0~100个。设计一个随机分配算法。要求：在结果随机（不可预知）基础上每种分配概率均等。如（25，25，25，25），（0，0，0，100）都是分配结果，机率一样。

看到下面有大量的讨论，我感觉都不是很对，所以写了这篇那博文，来讨论下这个问题。

下面我们先看看随机抽样问题，理解了这个算法，上述题目也可以解决。

要求从N个元素中随机的抽取k个元素，其中N无法确定。

这种应用的场景一般是数据流的情况下，由于数据只能被读取一次，而且数据量很大，并不能全部保存，因此数据量N是无法在抽样开始时确定的；但又要保持随机性，于是有了这个问题。所以搜索网站有时候会问这样的问题。这里的核心问题就是“随机”，怎么才能是随机的抽取元素呢？我们设想，买彩票的时候，由于所有彩票的中奖概率都是一样的，所以我们才是“随机的”买彩票。那么要使抽取数据也随机，必须使每一个数据被抽样出来的概率都一样。

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

[cpp] view plain copy

1. Init : a reservoir with the size： k  
2.   
3. for   i= k+1 to N  
4.        M=random(1, i);  
5.        if( M < k)  
6.                SWAP the Mth value and ith value  
7.  end for  

上面是算法的伪代码实现。 解释一下：程序的开始就是把前k个元素都放到水库中，然后对之后的第i个元素，以 k/i 的概率替换掉这个水库中的某一个元素。

下面来具体证明一下：每个水库中的元素出现概率都是相等的。

【证明】

（1）初始情况。出现在水库中的k个元素的出现概率都是一致的，都是1。这个很显然。

（2）第一步。第一步就是指，处理第k+1个元素的情况。分两种情况：元素全部都没有被替换；其中某个元素被第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步，所有元素出现的概率是相等的即可。

上面是我从网上找的一段比较好的证明过程。其实简单说就是可能被替换的概率1/k，可能能选到的替换体概率 1/(i+1)*k.然后两个相乘，结果是 1/i+1。所有的都是等概率的。

有了上面这个算法的支持，我们就可以很容易相处上面题目的解法，从上1-100个数中使用上面算法等概率的取出三个数，这三个数会把1-100 分成四分，然后这四分就是我们的分法。

下面是我实现的代码：

[cpp] view plain copy

1. #include<iostream>  
2. #include<time.h>  
3. #include<stdio.h>  
4. using namespace std;  
5.   
6. const int  N=100;  
7. const int  pool=3;  
8.   
9. void div(void);  
10. void swap(int *a,int *b);  
11. int  random(int min, int max);  
12.   
13. int main(void)  
14. {  
15.     div();  
16.     system("pause");  
17.     return 0;  
18. }  
19. void div()  
20. {  
21.       int  I_rus[N];  
22.       int  I_cun[pool];  
23.       //初始化原始资源  
24.       for(int i=0;i<N;i++)  
25.           I_rus[i]=1+i;  
26.       //初始化缓冲池  
27.       for(int j=0;j<pool;j++)  
28.           I_cun[j]=1+j;  
29.       //算法开始  
30.       for(int k=pool+1;k<N;k++)  
31.       {  
32.           int tem=random(1,k);  
33.           if(tem<pool)  
34.               swap(I_cun[tem],I_rus[k]);  
35.       }  
36.       cout<<I_cun[0]<<endl;  
37.       cout<<I_cun[1]<<endl;  
38.       cout<<I_cun[2]<<endl;  
39.    
40.   
41. }  
42. void swap(int *a,int *b)  
43. {  
44.     int tem;  
45.     tem =*a;  
46.     *a  =*b;  
47.     *b  =tem;  
48. }  
49. //产生随机数1-i  
50. int  random(int min, int max)  
51. {  
52.     srand( (unsigned)time( NULL ) );  
53.     return (min+rand() % (max-min+1))-1;  
54. }  

下面是一次运行结果，使用这三个数把1-100 分成四块就可以产生想要的结果。（我没有做排序）