【蓄水池抽樣應用之】如何等機率的從N個元素中選取出K個元素

如何等機率的從N個元素中選取出K個元素？這個問題就是一個蓄水池抽樣（Reservoir 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=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 
  故證明成立 

 

對於抽樣問題，最近看見了一些方法，做個總結：

問題：要求從1,2,3..n中，以等機率的方式，抽取m個元素

1、使用上面的蓄水池抽樣

void　sample_pool(const int N, const int m){       int i, rd;        int* x = new int[N];        for(i = 0 ; i < N ; i ++)            x[i] = i + 1;         for(i = m ; i < N; i ++ )        {             rd = rand()%i;             if(rd < m)             swap(x[i],x[rd]);         }         for(i = 0 ; i < m; i ++)             cout<<x[i]<<" ";         delete []x;         x = NULL;}//空間和時間均為O(N)

2 、從N個中選取m個， 可以先確定一個後，然後從身下的N-1個中選取m-1個出來。

void sample_rand(const int N,const int m){        int select = m,i,rd;        int remain = N;        for(i = 0; i < N ; i++)        {             rd = rand()%remain;             if(rd < select)            {                  cout<< i<<" ";                  select--;             }             remaining--;          }}

 上面這個方法非常經典，是Knuth在the art of computer programming中提出的。使用的額外空間為O(1)，時間為O(N)。其機率的證明也是非常簡單的。簡單推到可發現，是等機率選擇每個元素的。而且，最後肯定會選擇剛剛好m個元素，前面一直沒選擇的話，則當remaining == select時，就會都選擇。

3、將抽樣的看成是一個集合，則要從N中選擇出m個不同的元素，存入到集合中，可用set來完成

利用STL中的set來完成這個功能。

void sample_set(const int N,const int m){set<int>s;while(s.size()<m){s.insert(rand()%n);}for(set<int>::iterator it = s.begin();it!=s.end();it++)cout<<*it<<" ";}

4、擾亂一個遞增序列。

for i =[0,N)

swap(x[i],x[rand(i,n-1)];

有人證明，只要擾亂前m個就可以。

void sample_shuf(const int N,const int m){int i, j;int *x = new int[N];for(i = 0 ; i <Ｎ; i++)  x[i]=i+1;for(i = 0 ; i < m ; i ++){j = rand(i,n-1);swap(x[i],x[j]);}sort(x,x+m);Print(x,m);delete []x;x= NULL;}

關於採樣的幾個問題：

1、Given a random number generator which can generate the number in rang(1,5)uniformly, how can u use it to build a random number generator which can generate the number in range(1,7) uniformly?

解答：利用拒絕採樣定理

       首先，將（1，5）之間的隨機發生器使用兩次，按照五進位進行使用，拼成一個（1,25）的隨即發生器既：（[gen][gen]）5,每一[]為一個5進位上的位，換算為十進位為：x=gen*5+gen,在十進位上的範圍為：6-30，進行一個簡單的左移動，可換算成1-25範圍上的值; 然後將（1,25）平均分配到7中情況上面，考慮21是7的倍數，因此可以每三個做一個映射（當然，也可以不管，直接截斷7後面的數字，但是範圍太小，效率不高），1-3--》1,4-6--》2,19-21--》7，此時就是等機率的，如果產生了22-25之間的數字，可以有兩種方法決定結果：

     （1）拒絕採樣，重新再運算

     （2）如果得到了22-25之間的數字，則此次的隨即發生器結果，直接使用上一次得到的結果。這個方法有人證明過，是等機率的，演算法Metropolis Algorithm。

五進位表示：                     --> 對應的十進位：                減5平移

11 12 13 14 15                     6   7   8   9   10              1   2   3   4   5

21 22 23 24 25                     11  12  13  14  15              6   7   8   9   10

31 32 33 34 35                     16  17  18  19  20              11  12  13  14  15

41 42 43 44 45                     21  22  23  24  25              16  17  18  19  20

51 52 53 54 55                     26  27  28  29  30              21  22  23  24  25

2、Generate a random permutation for a deck of cards

解答：

        從後往前，第k步的時候，隨機產生一個1-k，之間的數字j，然後交換j和k處的數字，可以很容易的最後這個排列就是一個等機率得到的排列。

for k=N:1

    j = rand(1,k)

     swap(j,k)

end

      同樣的，也可以從前往後進行這個過程，不過產生的範圍就是變成k-N之間了。

for k = 1:N

     j = rand(k,N)

      swap(j,k)

end