蓄水池抽样及实现

蓄水池抽样(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");
}

 

 

posted @ 2012-11-27 13:21  handspeaker  阅读(...)  评论(...编辑  收藏