蓄水池抽样及实现
蓄水池抽样（Reservoir Sampling ）是一个很有趣的问题，它能够在o（n）时间内对n个数据进行等概率随机抽取，例如：从1000个数据中等概率随机抽取出100个。另外，如果数据集合的量特别大或者还在增长（相当于未知数据集合总量），该算法依然可以等概率抽样。
说蓄水池抽样之前，先说一下等概率随机抽取问题，等概率随机抽取是一个很有用的东西，因为在很多情况下，尤其是搞模式识别时，需要这个东西。比如，我们想从10000个样本中随机抽取5000个作为训练集，5000个作为测试集，那么等概率随机抽取便派上用场了。那么，究竟该如何做等概率随机抽取呢？一般的想法应该是，随机生成一个（0，n-1）之间的数x，然后抽取集合中第x个数据，如果本次生成的数x‘与之前某次生成的数x是相同的，那么继续随机生成，直到生成一个与之前所有生成过的数不同的数。并将这样的随机生成做m次。
这样做思路上是最简单的，但是问题却出现了，如果m比较小还好，比如n=10000，m=100，就是说10000个数里面随机挑100个，基本上没什么重合的概率，因为当做第100次随机生成时，之前才生成了99个，所以至多有99/10000≈1%的概率生成与之前重复的数。那么，我们可以很顺利的随机的、等概率的生成100个数，并且基本上只调用100次rand函数就行。但是如果m比较大呢？还是刚才那个例子，如果n=10000，m=5000呢？这涉及到一个求期望的问题，设sum为该算法调用rand函数的次数，我们来求一下E（sum）,由于每次随机生成的数都是独立的，因此E（sum）=E(num1)+E(num2)+...E(numm),所以我们求出生成第x个数需要调用多少次rand函数即可。
假设前面已经生成了x个数，我们要生成第x+1个数，那么，有如下概率：
调用1次rand函数就成功生成该数的概率为：（1-x/n）
调用2次rand函数就成功生成该数的概率为：（1-x/n）*（x/n）
。。。
调用k次rand函数就成功生成该数的概率为：（1-x/n）*（x/n）^(k-1)
则E（numx）=1*（1-x/n）+2*（1-x/n）*（x/n）+...+k*（1-x/n）*（x/n）^(k-1)...
这个计算并不难，在此我仅给出我计算出的期望：n/(n-x)（默认x从0开始，程序员的习惯）
最后，E（sum）=n/n+n/(n-1)+n/(n-2)+....+n/(n-m+1)
这个期望就比较难算了，复杂度大概是O（n*（lg(n）-lg(n-m)））级别的。在n比较大，m也比较大的时候，这个规模比O（n）可大多了。因此，蓄水池抽样在这个时候就有优势了，而且，对于另一种比较变态的情况，假设n非常大，以至于我们并不知道n的确切数量，而且n还在动态的增长，我们要不停的随机等概率的抽取n的一定比例（例如10%），这种情况下，上面所介绍的普通抽样方法就很难做到了。
蓄水池抽样：从N个元素中随机的等概率的抽取k个元素，其中N无法确定
先给出代码：
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
上述伪代码的意思是：先选中第1到k个元素，作为被选中的元素。然后依次对第k+1至第N个元素做如下操作：
每个元素都有k/x的概率被选中，然后等概率的（1/k）替换掉被选中的元素。其中x是元素的序号。
算法的成立是用数学归纳法证明的：
每次都是以 k/i 的概率来选择 例: k=1000的话， 从1001开始作选择，1001被选中的概率是1000/1001，1002被选中的概率是1000/1002，与我们直觉是相符的。 接下来证明： 假设当前是i+1, 按照我们的规定，i+1这个元素被选中的概率是k/i+1，也即第 i+1 这个元素在蓄水池中出现的概率是k/i+1 此时考虑前i个元素，如果前i个元素出现在蓄水池中的概率都是k/i+1的话，说明我们的算法是没有问题的。 对这个问题可以用归纳法来证明：k < i <=N 1.当i=k+1的时候，蓄水池的容量为k，第k+1个元素被选择的概率明显为k/(k+1), 此时前k个元素出现在蓄水池的概率为 k/(k+1), 很明显结论成立。 2.假设当 j=i 的时候结论成立，此时以 k/i 的概率来选择第i个元素，前i-1个元素出现在蓄水池的概率都为k/i。 证明当j=i+1的情况： 即需要证明当以 k/i+1 的概率来选择第i+1个元素的时候，此时任一前i个元素出现在蓄水池的概率都为k/(i+1). 前i个元素出现在蓄水池的概率有2部分组成, ①在第i+1次选择前得出现在蓄水池中，②得保证第i+1次选择的时候不被替换掉 ①.由2知道在第i+1次选择前，任一前i个元素出现在蓄水池的概率都为k/i ②.考虑被替换的概率： 首先要被替换得第 i+1 个元素被选中(不然不用替换了)概率为 k/i+1，其次是因为随机替换的池子中k个元素中任意一个，所以不幸被替换的概率是 1/k，故 前i个元素(池中元素)中任一被替换的概率 = k/(i+1) * 1/k = 1/i+1 则(池中元素中)没有被替换的概率为: 1 - 1/(i+1) = i/i+1 综合① ②,通过乘法规则 得到前i个元素出现在蓄水池的概率为 k/i * i/(i+1) = k/i+1 故证明成立
个人觉得这个算法在模式识别中，对数据集做随机划分（划分为训练集和测试集）的时候非常好用，因为一来我不需要预先知道数据的总量，二来对于数据量很大的情况速度还是比第一种方法快的。正好之前需要写一个随机划分的程序，小试了一下，感觉还不错，代码如下：
主要功能：
输入一个包含2类的libsvm样式的样本集合文件，并制定需要随机提取的样本数量
输出随机提取的样本集合
RandomSelect.h
#include<stdio.h> #include<windows.h> #include<vector> #include<time.h> #include<math.h> using namespace std; #define MAX_LENGTH 0x7fff void reservoirSampling(long select_num,long* pos_select,long* neg_select,char* input_path); void saveFile(long select_num,long* pos_select,long* neg_select,char* input_path,char* ouput_path); void usage();
RandomSelect.cpp
/* Author: SongQi Create Time:2012/9/19 Function: A program which can randomly select samples from the input sample file. It can be used to seperate training samples and testing samples.Now it only can select two-class samples. Input: 1.input sample file path 2.output sample file path 3.sample number need to select for each postive and negative samples. Output: 1.output sample file */ #include "stdafx.h" #include "RandomSelect.h" #include <fstream> using namespace std; int main(int argc, char** argv) { //input sample file path char* input_path = NULL; //output sample file path char* ouput_path = NULL; //selected sample number from the input file long select_num=0; if( argc != 4) { usage(); system("pause"); return 0; } else { input_path=argv[1]; ouput_path=argv[2]; select_num=atol(argv[3]); } long *pos_select=new long[select_num]; long *neg_select=new long[select_num]; //random select samples from the input file reservoirSampling(select_num,pos_select,neg_select,input_path); saveFile(select_num,pos_select,neg_select,input_path,ouput_path); system("pause"); return 0; } void reservoirSampling(long select_num,long* pos_select,long* neg_select,char* input_path) { FILE *input_file=fopen(input_path,"r"); if(input_file==NULL) { printf("the input sample file does not exist!\n"); usage(); system("pause"); return; } fseek(input_file,0,SEEK_SET); //set to the start of the file int label=0; long index=0; long pos_count =0; long neg_count =0; srand( time(NULL) ); char stuff; while(true) { if(fscanf(input_file,"%d ",&label)!=1) break; //printf("%d\n",label);RAND_MAX if(label==0) { if(neg_count<select_num) neg_select[neg_count]=index; else { long is_select=rand()%(neg_count+1)+1; if(is_select<=select_num) neg_select[rand()%select_num]=index; } neg_count++; } else { if(pos_count<select_num) pos_select[pos_count]=index; else { long is_select=rand()%(pos_count+1)+1; if(is_select<=select_num) pos_select[rand()%select_num]=index; } pos_count++; } //printf("pos:%dneg:%d\n",pos_count,neg_count); index++; char* detect_buffer = new char[2](); while(strcmp(detect_buffer,"\n")!=0) fread(detect_buffer,sizeof(char),1,input_file); } printf("%d\n",index); FILE *out_file1=fopen("pos_index.txt","w"); FILE *out_file2=fopen("neg_index.txt","w"); for(int i=0;i<select_num;i++) { fprintf(out_file1,"%d\n",pos_select[i]); fprintf(out_file2,"%d\n",neg_select[i]); } fclose(out_file1); fclose(out_file2); fclose(input_file); } void saveFile(long select_num,long* pos_select,long* neg_select,char* input_path,char* ouput_path) { ifstream fin(input_path); FILE *output_file=fopen(ouput_path,"w"); char line[MAX_LENGTH]; long index=0; while( fin.getline(line, MAX_LENGTH)) { //printf("%s",line); for(long count=0;count<select_num;count++) { if(index==pos_select[count]) { fprintf(output_file,"%s\n",line); } if(index==neg_select[count]) { fprintf(output_file,"%s\n",line); } } index++; } fclose(output_file); return; } void usage() { printf(" arg1:input sample file path\n"); printf(" arg2:output sample file path\n"); printf(" arg3:sample number need to select for each postive and negative samples\n"); }
