从随机数生成到随机采样的C++实现

随机数的妙用：上回书说到的关于用蒙特卡洛方法进行一些数值计算（比如pi的逼近、特殊几何图像面积（存在解析解）的计算），都还是随机数较为间接的应用（也就是我所说的，将随机数转换为一种均匀分布的概率表现），这节我们考察随机数更为直接的应用，那就是随机采样。

问题

生成随机整数的有序子集。程序的输入包括两个整数m和n，m < n，输出是[0, n)范围内的有序随机整数列表。从概率的角度要求，得到的每个选择都是等概率。

这里我们还可以对比matlab的随机采样的高层工具箱函数。randperm()和randsample。（对matlab不感兴趣的可以略过此节）

• randperm函数形式较为简单

p = randperm(n);        # 对n的全排列(permutation)
p = randperm(n, k);     # 返回全排列的k个值（不允许重复）

小应用（实现对原始数据的按行重新排列）

X = randn(n, m);
Y = X(randperm(size(X, 1)), :);

• randsample函数重载较多，功能也相对强大

y = randsample(n,k)
y = randsample(population,k)
y = randsample(n,k,replacement)
y = randsample(population,k,replacement)

• 重载1用法同randperm；
• 重载2或者写作y = randsample(ARRAY, k)，从数组中随机选出k个不同的数
• replacement是一个bool变量，为1的时候，允许重复，反之不允许，默认不允许。

实现类似的数据重新按行排列的语句：

X = randn(n, m);
Y = X(randsample(1:size(X, :), size(X, 1), 0), :)

当X是一个向量时：

X = randn(n, 1);
Y = randsample(X, length(X), 0);

解决方案之一

我们以考察0, 1, 2, ..., n-1整数序列为例，逐个考察每个整数的随机选择情况，不妨先取m = 2, n = 5，则第一个整数0被选中的概率为：

if (rand() % n < 2)

我们却不能按照此形式继续考察1, 2, ...

• 如果0被选中，则紧接着1被选中的概率为

if (rand()% (n-1) < (m-1)) 
//if (rand() % (5-1) < (2-1)）
// 也即是在余下的(n-1)个数中，仅有(m-1)次被选中的机会

• 如果0不被选中，则紧接着1被选中的概率为：

if (rand()%(n-1) < m)   // 这里的深层含义是：以条件判断表达概率
// if (rand()%(5-1) < 2)
// 也即是在余下的(n-1)个数中，仍有m次被选中的机会

也即是，如果要从$\n$个余下的整数中选中其中的$\m$个，我们以$\fracmn$选择下一个数

实现如下：

#include <iostream>
using namespace std;
void randsample(unsigned n, unsigned m)
{
    unsigned select = m;
    unsigned remaining = n;
    unsigned i;
    for (i = 0; i < n; ++i)
    {
        // 每一次循环，比如第i次循环，考察的是它被选中的概率
        // 考察到第i次时，可见之前已经历了共i次的判断（从0开始计数）
        // 也就是在余下的(n-i)个值中找到i出现的概率
        if (rand() % remaining < select)
        {
            cout << i << " ";    // 表示i被选中
            --select;
        }
        --remaining;
    }
    cout << endl;
}

我们还可以继续简化这个函数，并使用一个数组存放这些随机采样值：

#include <iostream>
using namespace std;
void randsample(unsigned n, unsigned m, unsigned arr[])
{
    unsigned idx = 0;
    for (unsigned i = 0; i < n; ++i)
    {
        if (rand() % (n-i) < m)
        {
            arr[idx++] = i;
            --m;
        }
    }
}

我们也可递归化（忽然想到一个话题，就是递归和数学归纳法的关系，这个问题的讨论留待以后了）这个实现。

void randsample(unsigned n, unsigned m)
{
    if (n <= m)
        return ;
    if (rand() % n < m)
    {
        cout << n - 1 << " ";
        randsample(n-1, m-1);
    }
    else
        randsample(n-1, m);
}

解决方案之二（集合实现）

我们利用集合自身的性质（即不允许重复），以及C++STL中set类模板的性质（自动被排序，因标准的STL set以RB-tree为底层机制）来实现随机采样值的生成。

#include <set>
using namespace std;
//typedef set<unsigned> SampleSet;
void randsample(unsigned n, unsigned m)
{
    set<unsigned> s;
    while (s.size() < m)
        s.insert(rand()%n);    // 集合的性质保证了最终得到的数据没有重复
    set<unsigned>::const_iterator = it;
    for (it == s.begin(); it != s.end(); ++it)
        cout << *it <<  " ";
    cout << endl;
}

解决方案之三（集合实现的优化）

我们来分析解决方案二可能存在的问题就是当m非常接近n时，进行

s.insert(rand()%n);         // 会存在大量无意义的值的添加操作

如何保证m的时间复杂度呢？

void randsample(unsigned n, unsigned m)
{
    set<unsigned> s;
    for (unsigned i = n-m; i < n; ++i)
    {
        unsigned t = rand()%(i+1);
        if (s.find(t) == s.end())
            // t <= i
            s.insert(t);
        else
            // i >= i,这样每次循环，循环m次都有一个随机值的产生，只不过是否等概率，就要打个问号了
            s.insert(i);
    }
}

解决方案之四

生成随机整数的有序子集的另一种方法是首先把0, 1, ..., n-1的数组顺序打乱，这称为shuffle，也是一种排列，permutation。

这里首先介绍一个randint(unsigned l, unsigned u)函数，

unsigned randint(unsigned l, unsigned u)
{
    return l + rand() % (u-l+1);
    // rand() % (u-l+1) -> [0, u-l]
    // [0, u-l] -> [l, u]
    // 也即是实现了[l, u]闭区间的随机整数采样
}

void genshuffle(unsigned n, unsigned m)
{
    unsigned* x = (unsigned*)malloc(sizeof(unsigned)*n);
    for (unsigned i = 0; i < n; ++i)
        x[i] = i;
    for (unsigned i = 0; i < m; ++i)
    {
        unsigned j = randint(i, n-1);
        unsigned t = x[i]; x[i] = x[j]; x[j] = t; // swap(x[i], x[j]);
    }
    sort(x, x+m); // 如果不进行排序，genshuffle的功能就类似于randperm(n, m)
    for (unsigned i = 0; i < m; ++i)
        cout << x[i] <<  " ";
    cout << endl;
}

x[0, ...,i-1]记录已选中的元素集合，x[i, ..., n-1]表示未选中的元素集合。以这种机制避免了对新元素是否已经选中的判断。