等概率随机抽样问题 || 蓄水池抽样算法

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

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

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

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

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

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

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

                       Init : a reservoir with the size： k
 
                       for   i= k+1 to N
                              M=random(1, i);
                              if( M < k)
                                      SWAP the Mth value and ith value
                        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 分成四分，然后这四分就是我们的分法。

下面是我实现的代码：

#include<iostream>
#include<time.h>
#include<stdio.h>
using namespace std;

const int  N=100;
const int  pool=3;

void div(void);
void swap(int *a,int *b);
int  random(int min, int max);

int main(void)
{
	div();
	system("pause");
	return 0;
}
void div()
{
	  int  I_rus[N];
	  int  I_cun[pool];
	  //初始化原始资源
	  for(int i=0;i<N;i++)
		  I_rus[i]=1+i;
	  //初始化缓冲池
	  for(int j=0;j<pool;j++)
		  I_cun[j]=1+j;
	  //算法开始
	  for(int k=pool+1;k<N;k++)
	  {
		  int tem=random(1,k);
		  if(tem<pool)
			  swap(I_cun[tem],I_rus[k]);
	  }
	  cout<<I_cun[0]<<endl;
	  cout<<I_cun[1]<<endl;
	  cout<<I_cun[2]<<endl;
 

}
void swap(int *a,int *b)
{
	int tem;
	tem =*a;
	*a  =*b;
	*b  =tem;
}
//产生随机数1-i
int  random(int min, int max)
{
    srand( (unsigned)time( NULL ) );
    return (min+rand() % (max-min+1))-1;
}

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