几种随机算法的实现原理

转自：http://www.androidstar.cn/%E5%87%A0%E7%A7%8D%E9%9A%8F%E6%9C%BA%E7%AE%97%E6%B3%95%E7%9A%84%E5%AE%9E%E7%8E%B0%E5%8E%9F%E7%90%86/

　　在日常工作中，经常需要使用随机算法。比如面对大量的数据， 需要从其中随机选取一些数据来做分析。 又如在得到某个分数后， 为了增加随机性， 需要在该分数的基础上， 添加一个扰动， 并使该扰动服从特定的概率分布（伪随机）。本文主要从这两个方面出发， 介绍一些算法， 供大家参考。

　　首先假设我们有一个使用的随机函数float frand()， 返回值在(0, 1)上均匀分布。大多数的程序语言库提供这样的函数。 在其他的语言如C/C++中， 可以通过间接方法得到。如 frand()= ((float)rand() ) / RAND_MAX;

1, 随机选取数据

　　假设我们有一个集合A(a_1,…,a_n), 对于数m，0≤m≤n, 如何从集合A中等概率地选取m个元素呢？
　　通过计算古典概率公式可以得到， 每个元素被选取的概率为m/n。 如果集合A里面的元素本来就具有随机性， 每个元素在各个位置上出现的概率相等， 并且只在A上选取一次数据，那么直接返回A的前面m个元素就可以了， 或者可以采取每隔k个元素取一个等类似的方法。这样的算法局限很大， 对集合A的要求很高， 因此下面介绍两种其他的算法。

1.1 假设集合A中的元素在各个位置上不具有随机性， 比如已经按某种方式排序了，那么我们可以遍历集合A中的每一个元素a_i, 0<=n 根据一定的概率选取ai。如何选择这个概率呢？

设m’为还需要从A中选取的元素个数， n’为元素a_i及其右边的元素个数， 也即n’=(n-i+1)。那么选取元素a_i的概率为 m’/n’。
由于该算法的证明比较繁琐， 这里就不再证明。
我们简单计算一下前面两个元素(2<=m<=n)各被选中的概率。
1) 设p(a_i=1)表示a_i被选中的概率。显而易见， p(a_1=1)=m/n, p(a_1=0)为(n-m)/n;
2)第二个元素被选中的概率为
p(a_2=1)= p(a_2=1,a_1=1)+p(a_2=1,a_1=0)
= p(a_1=1)*p(a_2=1│a_1=1)+ p(a_1=0)* p(a_2=1│a_1=0)
= m/n * (m-1)/(n-1) + (n-m)/n*m/(n-1)
= m/n

我们用c++语言， 实现了上述算法

template<class T>
bool getRand(const vector vecData, int m, vector& vecRand)
{
        int32_t nSize = vecData.size();
    if(nSize  < m || m < 0)
        return false;
    vecRand.clear();
        vecRand.reserve(m);
    for(int32_t i = 0, isize = nSize; i < isize ; i++){
            float fRand = frand();
            if(fRand <=(float)(m)/nSize){
                vecRand.push_back(vecData[i]);
                m--;
            }
        nSize --;
    }
    return true;
}

　　利用上述算法， 在m=4， n=10， 选取100w次的情况下， 统计了每个位置的数被选取的概率

位置 概率
1    0.399912
2    0.400493
3    0.401032
4    0.399447
5    0.399596
6    0.39975
7    0.4
8    0.399221
9    0.400353
10    0.400196

　　还有很多其他算法可以实现这个功能。比如对第i个数， 随机的从a_i, …, a_n中， 取一个数和a_i交换。这样就不单独介绍了。

1.2 在有些情况下，我们不能直接得到A的元素个数。比如我们需要从一个很大的数据文件中随机选取几条数据出来。在内存不充足的情况下，为了知道我们文件中数据的个数， 我们需要先遍历整个文件，然后再遍历一次文件利用上述的算法随机的选取m个元素。

　　又或者在类似hadoop的reduce方法中， 我们只能得到数据的迭代器。我们不能多次遍历集合， 只能将元素存放在内存中。 在这些情况下， 如果数据文件很大， 那么算法的速度会受到很大的影响， 而且对reduce机器的配置也有依赖。

这个时候，我们可以尝试一种只遍历一次集合的算法。

1. 取前m个元素放在集合A’中。
2. 对于第i个元素(i>m), 使i在 m/i的概率下， 等概率随机替换A’中的任意一个元素。直到遍历完集合。
3. 返回A’

下面证明在该算法中，每一个元素被选择的概率为m/n.

1. 当遍历到到m+1个元素时， 该元素被保存在A’中的概率为 m/(m+1), 前面m个元素被保存在A’中的概率为 1- (m/m+1 * 1/m) = m/m+1
2. 当遍历到第i个元素时，设前面i-1个元素被保存在A’中的概率为 m/(i-1)。根据算法， 第i个元素被保存在A’中的概率为m/i , 前面i-1各个元素留在A’中的概率为 m/(i-1) * (1-(m/i* 1/m) = m/i;
3. 通过归纳，即可得到每个元素留在A’中的概率为 m/n;

我们在类似 hadoop的reduce函数中， 用java实现该算法。

public void reduce(TextPair key, Iterator value, OutputCollector collector, int m)
{
    Text[] vecData = new Text[m];
    int nCurrentIndex = 0;
    while(value.hasNext()){
        Text tValue = value.next();
        if(nCurrentIndex < m){
           vecData[nCurrentIndex] = tValue;
        }
        else if(frand() < (float)m / (nCurrentIndex+1)) {
           int nReplaceIndex = (int)(frand() * m);
           vecData[nReplaceIndex] = tValue;
        }
        nCurrentIndex ++;
    }
 //collect data
…….
}

　　利用上述算法，在m=4， n=10， 经过100w次选取之后， 计算了每个位置被选择的选择的概率

位置 概率
1    0.400387
2    0.400161
3    0.399605
4    0.399716
5    0.400012
6    0.39985
7    0.399821
8    0.400871
9    0.400169
10    0.399408

2. 随机数的计算

在搜索排序中，有些时候我们需要给每个搜索文档的得分添加一个随机扰动， 并且让该扰动符合某种概率分布（这还算是随机扰动吗O(∩_∩)O）。假设我们有一个概率密度函数f(x), min<=x<=max, 并且有

那么可以利用f(x)和frand设计一个随机计算器r(frand()), 使得r(frand())返回的数据分布， 符合概率密度函数f(x)。
令

那么函数

符合密度函数为f(x)的分布。
下面对这个以上的公式进行简单的证明：

由于g(x)是单调函数， 并且x在[0,1]上均匀分布，那

由于上述公式太复杂， 计算运算量大， 在线上实时计算的时候通常采用线性差值的方法。
算法为：
1）在offline计算的时候， 设有数组double A[N+1]；对于所有的i, 0<=i<=N, 令

2）在线上实时计算的时候，令

f = frand()，
lindex = (int) (f* N);
rindex = lindex +1;

那么线性插值的结果为 A[lindex]*(A[rindex]-f) + A[rindex] * (f – A[lindex])
我们做了一组实验，令f(x)服从标准正太分布N(0,1)， N=10000， 并利用该算法取得了200*N个数。对这些数做了个简单的统计， 得到x轴上每个小区间的概率分布图。

3后记
　　在日常工作中， 还有其他一些有趣的算法。比如对于top 100w的query， 每个query出现的频率不一样， 需要从这100w个query， 按照频率越高， 概率越高的方式随机选择query。限于篇幅， 就不一一介绍了。

几种随机算法比较