蓄水池抽样

如何等概率的从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, tmp,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 <N; 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


3、Given a long log file ,pick 1000items evenly from them

已标记关键词 清除标记
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页