随机数的妙用:上回书说到的关于用蒙特卡洛方法进行一些数值计算(比如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]
表示未选中的元素集合。以这种机制避免了对新元素是否已经选中的判断。