09 Sep 2012

简单取样问题

最近在看概率分析的面试题。

首先，得知道两个最基本的离散分布：几何分布和二项分布，两者都基于独立Bernoulli试验。所谓Bernoulli试验即有两种结果的试验。几何分布描述的随机变量是取得一次成功要进行的试验次数，二项分布则描述n次试验中成功的次数。在考虑类似二项分布描述的问题时，可以用指示器随机变量求期望（这个看似没用的东西本质用了统计变量的数字特征将问题简化）。

有个比较经典的取样问题：n个数中随机取m个。

m=1时，只要random(1, n)。当m=2时，你可以两次random(1,n)。但当m比较大如m=n时，这种方法最后random出的几个数重复的可能性很大，效率较低，称算法1。算法1理想的情况是O(m)，实际效率与n/m有关，比值越大效率越高。

算法2与算法1等价，但确是算法1的泛化：先将结果初始化成1…m，然后把m次random出的结果当成优先级，根据这个优先级去排序。依然有重复问题，不过这时有新的办法了：办法1是如果重复继续对重复的优先级做random知道分出大小关系，办法2是增加random的范围，减少碰撞可能，如random(1，n^3)，这时所有元素唯一的概率至少是1-1/n。算法2时间耗在排序上，这里排序的key是优先级，value是1…m，总体是O(n+mlgm)，O(n)是初始化代价。

可以看到，算法2在做优先级的时候还是会有极端的大量重复情况。事实上，当m < n/2时，选第m个优先级时不重复的概率pm>1/2，而前面的所有m-1个做优先级时不重复的概率也都>1/2，因此有上界：m次成功概率>1/2的Bernoulli试验的期望。又几何分布的期望是1/p，所以这里的期望<2。也就是说，m<n/2时，用算法2最多要random2m次，显然没有之后排序复杂度高。

算法2还很自然地用在这样的问题上：rand5()模拟rand7()，依然有办法1和办法2去处理重复优先级的情况。

当m > n / 2时，算法2就不靠谱了，这样n/2后面元素不重复的期望上界奔着无穷大去了。这时可以用Knuth同学或者Bob Floyd同学的方法。下面来说说knuth的solution，叫它算法3，基本想法是比如m=2, n=5，选第一个是2/5，选第二个是1/4或者2/4，分别对应第一个数选到和没选到的random办法：

    for(i = 0; i < n; ++i)
    if(random() % (n - i) < m) 
        m--; // bang!

可以看到，分母(n-i)每个循环步都减1，最多减到0，分子只有在if里面才会减1。极端情况是全n-m个都没进if，那么剩下m个必定都进if。另外注意到算法3bang出来结果是大小有序的，生产环境中如果有乱序要求就需要再对结果的index做shuffle，但更多环境中要的是有序的，这也是算法3的优势。算法3的复杂度是O(n)。

Bob Floyd的算法4更进一步，但是它是从算法2的思路出发的，基于集合的，另外用到了搜索，复杂度是Ｏ(m)，我要解释又一大堆，网上有说的很好的。

当m=n时的问题叫shuffle，这时候肯定要求无序。经典的洗牌算法是基于交换，每次random一个(i,n)的index。这时也有需要初始化的额外开销，m！= n时也可以按这思路去做，复杂度是O(n+mlgm)。而对于初始化有优化的反案，即当无规律用到时初始化，可以用额外的from和to数字加一个top指针来做。在考虑shuffle的时候有两种错误：一是每次random(1,n)，这时因为前k个已经交换过了，如果还把它们的index放在里面概率就不是1/n了。二是随机一个offset，这时不能到1/n!，只能有1/(n!-n)。

有的环境n不知道，这时可以用所谓蓄水池方法去构造方案，然后证明概率相等，比如m=1时可以遍历所有行然后以1/行序号的概率替换已选结果。m>1时类似，证明都非常简单。

我们在学校或者书本学习了一类方法的理论，实际环境中往往有更加简单但不那么完美的方法。比如本文的问题当n为100w，m=100w-10时候可以生成一个包含10个元素的有序随机样本，然后bang出不再样本中的整数。再如n=2^31,m=1000w时候，可以random出一个1100w个整数然后去重得到结果。

写完文章两点体会是：一是像这种思路性的东西要当日记写，说这么细反而让人恶心。二是写前目标是不超过三段，每段几个字，因为牛文都那样，写完发现自己离写牛文差很远。

最后留个问题：random一个无结构稀疏矩阵，给定稀疏度，且无空行空列。

Hong Wu at 16:02