蓄水池抽樣及實現

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