MATLAB中的随机数发生函数

MATLAB中有两个基本的(伪)随机数发生函数, rand与randn. 前者产生0和1之间均匀分布的随机数, 后者产生均值为0, 方差为1的正态分布的随机数. 统计工具箱(Statistics Toolbox)中那些较为复杂的随机数发生函数通过调用它们来实现, 比如均匀分布函数unifrnd要调用rand, 正态分布函数normrnd要调用randn, 等等.

每次启动MATLAB时rand与randn的初始状态都会被自动重置, 从而会产生完全相同的伪随机序列. 这一点可以通过一个小实验验证. 启动MATLAB(我使用的版本是Version 7.6.0.324 (R2008a)), 输入命令rand, 得到0.8147. 紧接着输入randn, 得到-0.4326. 重启MATLAB后再次运行rand和randn, 得到的结果与上述结果一样.

为了得到不同的伪随机序列, 可以自行重置随机数发生函数的状态, 相应的语句是rand(method, s)和randn(method, s). 这里的两个参数都需要稍加说明.method是一个字符串, 用来指明随机数发生器使用那种算法. rand为method的取值提供三种选择:'twister', 'state', 'seed', 其中最新的'twister'方法产生的伪随机序列重复周期最长, 也是MATLAB 7.4及更高版本的默认方法.randn中只包括后两种method, 默认方法是'state'.

另一个参数s设置所选定方法的状态. 可以用rand(method)和randn(method)查看method方法当前的状态(这一命令并不改变随机数发生函数当前使用的算法):'seed'方法的状态是一个数, 而'state'与'twister'方法的状态是一个数组. 用来设置状态的参数s应该是一个自然数, 它将随机数发生器所选定的方法的状态设置为``第s个''可能的状态. 需要特别指出的是, 无论使用何种方法,s的取值最好限制在0到2^31-1之间. 超出这个范围时, 不同的s可能对应同一个状态. 可以验证, 所有不小于2^32-1的s所对应的状态都与s = 2^31-1所对应的状态相同.s取2^31到2^32-2之间的不同值时对应的状态可能不一样, 但这些状态都已在s取0到2^31-1时出现过.

有时我们希望每次运行得到的伪随机序列都不一样, 这样看起来更像随机的. rand的帮助文档中给了一个示例方案:

rand('twister',sum(100*clock))

这一思路是将状态值设置为时变的. clock函数返回反映当前时间的1x6数组[年, 月, 日, 时, 分, 秒](注意sum(100*clock)一般不是整数. 这个例子似乎不太合乎规范, 但运行起来没问题). 但是这一方法效率并不高, 主要原因是总状态数为2^31, 约等于2.15e9; 而sum(100*clock)的上界仅约为2.1e5, 而且并非小于这个上界的每个值都能取到.这个帖子对此进行了详细分析, 并且给出了一个效率较高的方法:

rand('twister', fix(mod(1e11*(sum(clock)-2009), 2^31)))

除了rand和randn这两个基本函数之外, MATLAB的统计工具箱还提供一个函数randg,用来产生规一化的服从Gamma分布的随机数. Gamma分布随机数产生函数gamrnd和泊松分布随机数产生函数poissrnd就需要通过调用randg实现功能. 看一下randg的帮助文档不难发现, 这个随机数产生器的状态是由rand和randn二者的状态共同决定的.

另外, 根据mathworks网站上最新的说明文档, 上述手动设置随机数产生器状态的方法已经过时了, 保留这一方式只是为了保证后向兼容. 新的方法我还没试过(原文: For version 7.7, use the default stream as described in the @RandStream reference documentation.), 毕竟我的MATLAB版本才到7.6.