随机数的生成
引言
能够产生伪随机数对于模拟事件、估计概率和其它数量、进行随机分配或选择以及对符号式结果进行数值型测试等是非常重要的. 这种应用程序可能需要均匀分布或非均匀分布的数,以及放回或无放回取样的元素.
函数 RandomReal、RandomInteger 和 RandomComplex 生成均匀分布的随机数. RandomVariate 生成内置分布的数. RandomPrime 生成一定范围内的素数. 函数 RandomChoice 和 RandomSample 从一列值中进行放回或无放回抽样. 元素的权重可以相同也可以不同. 此外,还有一个用于为随机数生成定义额外方法和分布的框架.
一系列非重复事件可以通过 RandomSample 进行模拟. 例如,通过按次序对从1到
的整数随机抽样的概率可以它来模拟.
这里估计的是

从2到8变化时,按次序得到

个元素的概率.
| Out[10]= |  |
| Out[11]= |  |
随机数生成在蒙特卡洛估计中居核心地位. 函数
的期望值的估计可以通过生成对应分布的函数值并将这些值代入
求其平均值得到.
| Out[12]= |  |
| Out[13]= |  |
随机过程的模拟可以通过生成具有所需属性的一系列数字实现. 随机漫步可以通过对伪随机数递归加和创建.
| Out[14]= |  |
随机数的替代可用于测试符号表达式的等价性. 例如,两个表达式之差的绝对值可以在随机生成的点上进行计算,以测试表达式的不相等.
这里没有提供证据表明

和

在实数域内不同.
| Out[15]= |  |
这里提供了证据表明

和

至少对于某些复数值是不等的.
| Out[16]= |  |
RandomPrime 等概率的选择素数. 它对于生成RSA加密大素数等操作是很有用的. 素数在该区域中的素数上均匀分布,但对于整个范围并不是如此,因为总体上,素数在正整数的范围内并不是均匀分布的.
| Out[17]= |  |
随机数生成函数
主要函数是 RandomReal、RandomInteger、RandomComplex、RandomVariate、RandomChoice 和 RandomSample. RandomReal、RandomInteger 和 RandomComplex 生成给定数值范围内的数值. RandomVariate 生成来自一种统计分布的数值. RandomChoice 和 RandomSample 生成的元素可以来自包括非数字值的有限集合.
随机数
RandomReal 在一个给定的实数域内生成伪随机实数. RandomInteger 在一个给定的整数域内生成伪随机整数. RandomComplex 在复数平面上一个给定的矩形区域上生成伪随机复数. RandomVariate 生成的伪随机数服从一定的统计分布. RandomPrime 生成在一定范围内概率相等的质数.
随机实数的生成.
随机整数的生成.
随机复数的生成.
服从某种分布的随机值的生成.
随机素数的生成.
如果定义域以
和
的形式指定,RandomReal 和 RandomInteger 在指定范围内生成均匀分布的数值. RandomVariate 采用的规则由指定的分布定义. 另外,定义新的方法和分布的机制也包括在内.
双参数接口提供了一种便捷的途径来一次获得多个随机数. 更重要的是,一次生成大量伪随机数有显著的效率优势.
生成

个0到1之间的数只需要几分之一秒.
| Out[18]= |  |
逐次生成

个数需要的时间大约是上面的5倍.
| Out[19]= |  |
对于维数为
到
的多维数组,先按照所需总数生成
个伪随机数,然后再进行分区. 这使生成多维数组的效率尽可能的达到最高,原因是生成全部随机数的效率达到了最高,而分区所用的时间是可以忽略不计的.
生成一个100×100×100×10数组与生成一个

个数的向量所需的时间大致相同.
| Out[20]= |  |
| Out[21]= |  |
一个维数相同的数组,如果一次生成10个数,所需的时间将是上面的好几倍.
| Out[22]= |  |
对于统计分布来说,一次生成多个数的速度优势可能更加明显. 除了继承所用的均匀数字生成程序的这种效率优势,许多统计分布还得益于初等和特殊函数向量化的计算. 例如,WeibullDistribution 受益于初等函数 Power、Times 和 Log 的向量计算.

个韦伯数在瞬间就能生成.
| Out[23]= |  |
如果逐次生成

个 Weibull 数则要花上好几秒.
| Out[24]= |  |
随机数生成在进行探索性研究时可能会非常有用. 例如,您可能会在一个较长的数字序列中寻找一个随机数字序列.
| Out[25]= |  |
将

的头一百万位数字转化成一个整数字符串.
给出5位小数组成的字符串出现在

的头一百万位数字中的位置.
| Out[27]= |  |
在进行分布估计时,如果封闭形式的结果未知或难以计算,随机数生成也是非常有用的. 随机矩阵属性提供了一个例子.
估计均匀分布的 5×5 实数矩阵有实数特征值的概率.
| Out[28]= |  |
| Out[29]= |  |
用于贝叶斯统计的 Gibbs 采样程序是一个模拟多元分布的例子[1]. Gibbs 采样器提供了一种手段来模拟多元分布的值,只要这种多元分布的每一个坐标相对于其它坐标的条件分布是已知的. 在某些限制条件下,通过对条件分布的反复抽样构造的随机向量的分布将收敛于真正的多元分布.
下面的例子将为 Casella 和 George 给出的例子构造一个 Gibbs 采样程序[2]. 所观察的分布是二元的.
已知时,
的条件分布是二项分布,而
已知时,
的条件分布是
分布. 正如 Casella 和 George 所提到的,有各种使用 Gibbs 采样程序来检测收敛性和采样的的对策可以推荐. 为简单起见,假设通过1000次迭代可以达到收敛. 该分布的大小为
的样本将从第1000次迭代后的
个值中采取. 应该注意的是这
个值是相关的.
定义一个二项分布和

条件分布的采样程序.
Gibbs 采样程序也可以被定义为随机数字生成的分布框架内的一个分布对象. 在 "定义分布生成程序"中提供了一个关于 Gibbs 采样程序作为分布对象的例子.

是一个长度为

的样本.
| Out[33]= |  |
| Out[34]= |  |
只要条件分布的数据足够,条件分布应该与所假设的二项分布和
分布非常一致. 最大的数据量发生在边缘分布的密度最高时,因此这些值可用于比较. 下面的图形比较了经验分布和假设的条件分布,其中利用了宽度为0.05的箱来估算连续值的概率.
在

时,比较

的经验和理论分布.
| Out[35]= |  |
在

时,比较

的经验和理论分布.
| Out[36]= |  |
任意精度的实数和复数
RandomReal 和 RandomComplex 默认生成机器精度数. 对于连续分布,RandomVariate 默认生成机器数. 通过设置 WorkingPrecision 选项可以获得任意精度的数.
RandomReal、RandomComplex 和 RandomVariate 的选项.
该选项对于均匀分布的实数、复数以及内置分布的实数是有效的. WorkingPrecision 也可以被纳入用户定义的分布.
这里是一个精确到25位有效数字的5到50之间的实数.
| Out[37]= |  |
这里给出一个精确到50位有效数字的

分布的数.
| Out[38]= |  |
在模拟过程中,如果有精度损失的可能并且需要高度精确的结果时,提高 WorkingPrecision 可能会很有用. 提高精度也可用于估计计算中的精度损失.
估计区间

上计算

时最大的精度损失.
| Out[39]= |  |
如果输入的精度小于指定的 WorkingPrecision,函数将会发出问题警告,然后输入的精度将被人为地增加来生成所需精度的伪随机数.
| Out[40]= |  |
WorkingPrecision 不是 RandomInteger 的选项. 整数具有无限的精度,因此精度完全由函数名指定.
| Out[41]= |  |
随机元素
RandomChoice 和 RandomSample 从一列可能元素中生成伪随机的选择. 该元素可以是数字或非数字.
| RandomChoice[{e1,e2,...}] | 给出一个 的伪随机选择 |
| RandomChoice[list,n] | 从 list 中给出一个由 n 个伪随机选择组成的列表 |
| RandomChoice[list,{n1,n2,...}] | 从 list 中给出一个 × ×... 的伪随机选择 |
| RandomChoice[{w1,w2,...}->{e1,e2,...}] |
| 给出一个权重为 的伪随机选择 |
| RandomChoice[wlist->elist,n] | 给出一个由 n 个加权选择组成的列表 |
| RandomChoice[wlist->elist,{n1,n2,...}] |
| 给出一个加权选择的 × ×... 数组 |
从列表中作出随机选择.
| RandomSample[{e1,e2,...},n] | 给出由 中 n 个元素组成的伪随机样本 |
| RandomSample[{w1,w2,...}->{e1,e2,...},n] |
| 给出使用权重 从 中选出的 n 个元素组成的伪随机样本 |
| RandomSample[{e1,e2,...}] | 给出 的一个伪随机排列 |
| RandomSample[wlist->elist] | 使用初始权重 wlist 给出 elist 的一个伪随机排列 |
从列表中进行随机采样.
RandomChoice 与 RandomSample 的主要区别在于 RandomChoice 在
中进行的是放回抽样,而RandomSample 进行的抽样是无放回的. 通过 RandomChoice 所选出的元素的数量不受 elist 中元素数量的限制,一个元素
可以被多次选择. 由 RandomSample 所返回的样本的大小受 elist 中元素数量的限制,样本中不同的元素出现的次数受该元素在 elist 中出现次数的限制.
如果 RandomChoice 或 RandomSample 的第一个参数是一个列表,则各元素被选择的概率相同. 权重的设定定义了
的集合上的一个分布. 权重必须是整数,但权重之和不必为1. 对于权重
,
在初始分布中的概率是
. 由于 RandomSample 进行的是无放回抽样,根据每次选择后剩余的总权重,权重进行内部更新.
RandomChoice 可以用于对可能结果为有限列表的独立同分布事件的模拟.
| Out[42]= |  |
| Out[43]= |  |
RandomChoice 可以用于生成有限支撑上任意离散分布的观察数据.
| Out[44]= |  |
| Out[46]= |  |
当结果为有限集合,结果列表中每个元素只出现一次时,RandomSample 可以用来模拟观察数据. 列表中不同元素出现的次数可以大于1.
这是从一个装有80个蓝色物体和45个红色物体的容器中,7次取物的模拟.
| Out[47]= |  |
对列表中所有元素进行随机采样的结果是一个随机排列.
| Out[48]= |  |
对元素的权重进行分配的结果是,值的权重越高,在随机排列中出现的往往越早.
| Out[49]= |  |
对于同一个加权或未加权的列表,RandomSample[#, 1]& 的分布与 RandomChoice 等价.
这将给出由样本量为1的

个随机样本确定的经验概率密度函数.
| Out[51]= |  |
| Out[52]= |  |
| Out[53]= |  |
这两个例子中的概率非常接近,并接近于理论值.
| Out[54]= |  |
RandomSample 也可用于组的随机分配,例如临床试验. 下面使用的是整数,但也可用其它识别值代替,例如姓名或身份证号码等.
| Out[55]= |  |
RandomChoice 和 RandomSample 受 SeedRandom 中 Method 选项的变化的影响. 内置方法的介绍见"方法". 另外,定义新方法的机制在"定义自己的生成程序"中介绍.
起点和局部化
伪随机数生成程序从算法上创造具有明显随机性水平的数字. 伪随机数的生成方法通常使用递推关系从当前状态产生一个数,并建立一个产生下一个数字的新状态. 状态的设定可以通过给生成程序指定一个整数作为起点,这个起点将用于算法中递推关系的初始化.
初始起点或称种子的给定就确定了伪随机数生成程序. 在许多情况下,往往希望局部地或全局地为一个随机数生成程序设置起点,以获取序列不变的随机值. 如果全局设置,起点将影响未来的伪随机数,除非新的起点被显式设置. 如果局部设置,起点只影响局部化代码内随机数和元素的生成.
局部化和起点函数.
SeedRandom 函数为随机生成程序提供了一种设置起点的方式. 单独使用时,SeedRandom 将全局地为随机生成程序设置起点. BlockRandom 函数为随机生成程序提供了一种局部设置或改变起点而不影响全局状态的方式.
| Out[56]= |  |
下面给出了两个不同的数,原因是第一个 RandomReal 是在 BlockRandom 内部生成,而第二个是在 BlockRandom 外部生成.
| Out[57]= |  |
SeedRandom 还提供了随机生成程序的切换机制.
用于 SeedRandom 的选项.
一个单独的生成程序起点的设置可以直接通过 Method 选项指定该生成程序来实现. 所有的生成程序都可以通过Method->All 来设置起点.
这里默认生成程序的起点为

,但

生成程序则不是.
| Out[58]= |  |
将

生成程序的起点设置为

,得到的是不同的随机数.
| Out[59]= |  |
并行计算中的 SeedRandom 和 BlockRandom
在并行计算中使用指令 SeedRandom 和 BlockRandom 有一些微妙之处. 在一个并行计算中,这些指令仅影响用于当前线程的生成程序. 而通常您希望在整个并行计算之前或之内使用.
对于并行计算而言,在每一线程都有一个生成程序,生成的随机数独立于其它线程上生成程序产生的随机数这一点是非常有利的. 在Mathematica 中,并行计算中所用的各线程将被给予一个从零开始的唯一编号值(通常会依次直到 $ProcessorCount),这将用于给出各线程上的不同种子及生成程序.
下表介绍了在序列及并行计算中使用的一些不同.
| | |
| SeedRandom[seed] | 用 seed 作为当前所有序列随机生成程序的种子,用 seed + i 作为当前所有并行生成程序的种子,其中 i 为并行线程的编号 | 仅将 seed 用作当前线程随机生成程序的种子 |
| SeedRandom[seed,Method->"ParallelGenerator"] | 用 seed + i 作为并行生成程序的种子,其中 i 为并行线程的编号 | 无影响 |
| SeedRandom[Method->method] | 将序列随机生成程序的方法变为 method | 仅将当前线程随机生成程序的方法变为 method |
| BlockRandom[expr] | 将所有伪随机生成程序局域化来计算 expr | 仅将当前线程的伪随机生成程序局域化来计算 expr |
序列与并行计算中的 SeedRandom 和 BlockRandom.
| Out[61]= |  |
| Out[62]= |  |
尽管具有相同的种子,结果是不同的. 主要的不同在于排序,原因在于并行调度程序在一个计算重复时,可能会在一个线程之前运行另一个线程.
| Out[63]= |  |
很多,但不是所有的相同结果都能在这两个计算中找到. 这是因为当一个计算重复时,不能保证一个给定的线程使用完全相同的次数.
由于前面的并行计算是在
BlockRandom 内部完成的,并行生成程序的状态已经恢复到使用之前的状态,因此再次运行实际上是一种重复.
| Out[65]= |  |
SeedRandom 和 BlockRandom 在并行计算内部的使用应该谨慎进行,因为通过同一线程完成的不同部分可能会产生相同的结果.
| Out[80]= |  |
您会注意到有些结果看上去是相同的. 这可以用 Union 来检查.
| Out[81]= |  |
这种情况下,在20个结果中仅有8个不同的和式. 如果运行,并集的长度通常等于您机器的处理器数目. 这是因为每个处理器的生成程序在每次使用之前都被重新给予种子,并且每种情形中 RandomReal 的使用是相同的,因此结果完全相同.
| Out[92]= |  |
| Out[93]= |  |
在并行计算内部通过 SeedRandom 可完成的一项任务是设置生成程序. 假设你希望将各线程的生成程序设置为不同种子的默认
生成程序.
下面定义一个经编译的函数,将随机生成程序变为

方法,并用
s 作为种子.
| Out[9]= |  |
通过与生成程序以相同方式设置的序列计算进行比较,可以验证使用的是这些生成程序.
进行序列计算,局域设置生成程序,设置方式与并行计算生成程序的设置方式相同.
| Out[10]= |  |
并行结果只是上面结果的一种置换.
| Out[31]= |  |
| Out[32]= |  |
以这种方式设置生成程序并不可取,因为使用相同的生成程序,仅改变种子不能保证生成的数字不以某种方式相关.
设置并行生成程序的一种较简单并更可靠的方法是通过 "ParallelGenerator" 方法提供,在方法中介绍.
方法
有5个伪随机生成程序方法存在于所有的系统中. 在这5种方法中,Mersenne Twister 法在序列或并列版本中均有提供. 一个依赖于平台的第六种方法存在于基于英特尔处理器的系统中. 方法名称用于处理并行计算时的生成程序. 在"定义自己的生成程序"中所介绍的定义新方法的框架也被纳入.
| "Congruential" | 线性同余生成程序(低质量的随机性) |
| "ExtendedCA" | 扩展的元胞自动生成程序(默认) |
| "Legacy" | Mathematica 6.0 之前的默认生成程序 |
| "MersenneTwister" | Mersenne Twister 移位寄存器生成程序 |
| "MKL" | 英特尔 MKL 生成程序(基于英特尔处理器的系统) |
| "ParallelGenerator" | 并行计算时用于生成程序的初始化和种子设置 |
| "ParallelMersenneTwister" | 周期为 的1024个 Mersenne Twister 生成程序集 |
| "Rule30CA" | Wolfram rule 30 生成程序 |
内置方法.
| Out[54]= |  |
| Out[55]= |  |
Congruential
使用一种线性同余生成程序. 这是最简单的伪随机数字生成程序之一,伪随机数在0 到 1 之间通过
得到,其中
根据模递推关系
给出,固定的整数
,
和
分别被称为乘数,增量和模量. 如果增量为0,该生成程序是一个乘法同余生成程序.
,
和
的值可以通过设置
方法的选项给出.
| | |
| "Bits" | Automatic | 为由位元构造的数指定位元的范围 |
| "Multiplier" | 1283839219676404755 | 乘数值 |
| "Increment" | 0 | 增量值 |
| "Modulus" | 2305843009213693951 | 模量值 |
| "ConvertToRealsDirectly" | True | 是否直接通过同余关系构建实数 |
用于 Method
的选项.
线性同余生成程序是周期性的,往往给出质量较低的随机性,尤其是当需要大量的随机数时. 如果直接通过同余关系生成实数,则周期小于或等于
.
默认选项值的选择原则是使周期很大并具有64位效率. 通过默认的选项设置,
生成程序通过了许多随机性的标准测试,尽管同余生成程序具有一些固有问题.
上述乘法同余生成程序的周期受到小于或等于模且相对于模为素数的正整数个数的限制. 上限是模的 Euler 函数.
| Out[57]= |  |
实际周期可以通过找到满足
模除
最小的整数
确定.
| Out[58]= |  |
| Out[59]= |  |
通过绘制生成数字序列的图形,也可以看到这些不同的数字.
| Out[60]= |  |
如果
设置为 False,实数可以通过序列元素中每次取8个位元来构造52位机器精度数而生成. 以这种方式产生的同余数字仍将循环,但循环过程将取决于位模式的重复,而不是初始的同余关系.
选项
可以是 Automatic,非零整数,或两个非零整数组成的数列,指定由位元构造数字所用的模
中位元的范围. 当
不是2的幂时,Automatic 使用
,否则使用
.
ExtendedCA
默认的
方法利用元胞自动机来生成高质量的伪随机数. 这种生成程序使用一种特别的5邻规则,因此每一个新单元取决于前一步中5个不相邻单元.
基于元胞自动机的随机数字生成程序根据一个确定规则来进行0和1组成的状态向量的演化. 对于一个给定的元胞自动化,新状态向量中给定位置处的一个元素(或元胞)由旧状态向量时该元胞的某些临近元胞确定. 状态向量中元胞的子集以随机位元的形式输出,并由此生成伪随机数.
用于
的元胞自动化产生的随机性水平非常高. 尽管一个元胞与前5个之间存在明显的相关性,输出中每一个单一元胞的使用都能够给出一个能通过许多随机性测试的位元流.
有用于修改状态向量大小、跳过的元胞及起始单元的选项. 默认值的选择原则着眼于质量和速度,而且通常没有必要修改这些选项.
| | |
| "Size" | 80 | 作为64的乘数的状态向量大小 |
| "Skip" | 4 | 要跳过的元胞数目 |
| "Start" | 0 | 从哪个元胞开始 |
用于 Method
的选项.
所用状态向量的长度默认时设置为
个元胞. 64 的倍数可由选项
控制. 一旦一个状态向量利用5邻规则通过元胞自动机的演化得到计算,随机数的位数将从 位
中选择.
将每一个状态向量的元胞四个一组进行分组,使用每一组中的第四个元胞,这种用法在实际应用中足以通过非常严格的随机测试. 这是
选项的默认用法. 为了更快地生成随机数,可以将
设置为2甚至是1,但随机数的质量会下降.
选项
与较大的
和
一起使用在设置一个用于并行计算的独立生成程序族时非常有用.

是默认的数字生成程序.
| Out[61]= |  |
| Out[62]= |  |
Legacy
方法
使用的生成程序在 Mathematica 6.0 之前的版本中被 Random 调用. 对于实数使用的是 Marsaglia-Zaman 带借位的减法生成程序(subtract-with-borrow generator). 整数生成程序基于Wolfram 规则30 元胞自动机生成程序. 规则30生成程序用于直接生成小整数以及生成大整数的一定数位.
| Out[63]= |  |
| Out[64]= |  |
为了保证与版本6.0之前产生的序列的一致性,为自动设置方法也适用于种子的"遗产"的方法. Automatic 方法的起点设置也应用于
方法.
方法中没有选项.
MersenneTwister
使用的是来自于 Matsumoto 和 Nishimura 的 Mersenne Twister 生成程序[3][4]. The Mersenne Twister 是一个周期为
的广义反馈移位寄存器生成程序.
通过 Mersenne Twister 生成程序生成5个随机数.
| Out[65]= |  |
方法中没有选项.
MKL
方法使用的是在英特尔的 MKL 库中所提供的随机数生成程序. MKL 程序库与平台有关. 在 微软Windows (32位,64位),Linux x86 (32位,64位) 以及 Linux Itanium 系统中可用均有方法
.
用于 Method
的选项.
| "MCG31" | 31位乘法同余生成程序 |
| "MCG59" | 59位乘法同余生成程序 |
| "MRG32K3A" | 带两个3阶组件的组合式多种递归生成程序 |
| "MersenneTwister" | Mersenne Twister 移位寄存器生成程序 |
| "R250" | 广义反馈移位寄存器生成程序 |
| "WichmannHill" | Wichmann-Hill 组合式乘法同余生成程序 |
| "Niederreiter" | Niederreiter 低差异序列 |
| "Sobol" | Sobol 低差异序列 |
方法.
前六个方法是均匀生成程序.
和
生成 Niederreiter 和 Sobol 序列. 这些序列是非均匀的,并且具有底层结构,有时在数值方法中是有益的. 例如,在多维蒙特卡洛积分中,这些序列通常能够带来较快的收敛.
下面显示的是一个2维 Niederreiter 序列的结构.
| Out[66]= |  |
| Out[67]= |  |
Rule30CA
方法
使用的是一种 Wolfram 规则30元胞自动机生成程序. 位元通过使用关系
将一个0和1组成的状态向量演化得到,其中
是
时刻元胞
的值.
Method
的选项.
所用状态向量的长度默认设置是
个元胞. 29的乘数可以通过选项
控制.
用

给出一个随机整数的2×3×4张量.
| Out[68]= |  |
方法
仅使用每一状态向量的第一个位元,这使得其速度比方法
慢,后者使用每一状态向量的多个位元.
ParallelMersenneTwister
使用 Matsumoto 和 Nishimura 所贡献的一组 Mersenne Twister 生成程序[3][4],参数的选择利用他们的"动态创建"程序 dcmt [19]. 该程序计算 Mersenne Twister 生成程序的互质参数,因此应该生成独立结果. 通过参数计算生成周期为
的广义反馈移位寄存器生成程序.
一个选项被包括在内,以选择所用的生成程序组.
| | |
| "Index" | 0 | 使用从 0 到 1023 的哪个生成程序 |
Method
的选项.
这里从不同的并行 Mersenne Twister 生成程序得到两组2500个随机数,并将每一数对作为一点制图.
| Out[10]= |  |
由这两个生成程序生成的数之间没有明显的相关性. 由于无相关性以及速度的缘故,这组生成程序用作并行计算的默认生成程序.
ParallelGenerator
是一个控制器方法,它允许您为用于并行计算的生成程序提供种子并进行改变.
一个选项被包括在内,以选择所用的生成程序组.
Method
选项.
赋给
方法的 Method 选项值可以是一个指定内置参数化方法的字符串,也可以是一个对非负整数给出随机生成程序指定的函数. 用于并行计算的各线程将被赋予一个从零开始的唯一编号(并且通常会依次直至$ProcessorCount),这将用于在各线程给出不同种子和生成程序.
| "ParallelMersenneTwister" | 周期为 的并行 Mersenne twister 生成程序 |
| "ExtendedCA" | 不同起始位置的扩展 CA 生成程序 |
| f | 用于第 i 个线程的生成程序 |
| "Default" | 还原默认方法 |
并行生成程序的方法.
字符串快捷方式作为一种方便的方式被提供,以得到两组高质量的独立生成程序.
的使用等价于使用函数f=Function[{i}, {"ParallelMersenneTwister", "Index"->i}]. 对于并行计算,这是默认的,因为生成程序速度很快并能产生高质量的随机数.
的使用通常等价于使用函数 f,该函数在下面根据用户机器上的处理器数目定义.
将方法重置为默认的
方法.
这里定义

选项的默认函数
Method->"ExtendedCA".
| Out[11]= |  |
参数的选项是这样的:如果您使用了机器上所有的 $ProcessorCount 处理器,您仍能得到与默认的序列
随机生成程序所产生的随机数质量一样好的随机数.
方法还以一种略微不同的方式对生成程序设置种子. 不同于仅在各处理器上使用相同的种子, SeedRandom[seed, Method->"ParallelGenerator"] 在各线程上使用
,其中 i 为该线程的独特编号. 这允许您在即使各线程的生成程序设置相同(即 Method->Function[{i}, "ExtendedCA"])的情况下,不同线程也可得到不同数字. 当然这种做法并不可取,因为即使种子不同,数字之间也可能会有不可预料的相关性.
一般地,针对不同线程给出生成程序方法的函数 f 可以返回随机生成程序方法所允许的任何结果.
| Out[3]= |  |
下面定义一个函数,给出编号在0到7之间的不同的生成程序方法.
下面将并行生成程序变为由函数给出的生成程序,并设置种子.
| Out[13]= |  |
下面进行序列式计算,局域地将生成程序设置为由函数给定的生成程序.
| Out[14]= |  |
| Out[15]= |  |
如要将并行生成程序恢复至默认方法,需要显式给出一个方法选项,否则,它仅改变种子.
定义自己的生成程序
只要方法遵循的是正确的模板,就可以插入随机的框架内. 生成程序对象的形式为
,其中 gsym 是识别生成程序及用于附加规则的符号. data 专用于与生成程序定义相关的顶层计算.
生成程序的初始化通过调用
来完成.
| Random`InitializeGenerator[gsym,opts] |
| 使用选项 opts 初始化生成程序 gsym |
生成程序初始化函数.
应该返回形如
的生成程序对象 gobj.
生成程序可以支持随机位流、随机整数和随机实数的生成. 如果生成程序支持位流,实数和整数的生成可以由位流转换. 在设定方法的时间,属性被查询以确定支持的对象和方式.
| GeneratesBitsQ | 如果方法生成位流,则设置为 True |
| GeneratesIntegersQ | 如果方法生成给定范围的整数,则设置为 True |
| GeneratesRealsQ | 如果方法生成给定范围和给定精确度的实数,则设置为 True |
生成程序的属性.
如果支持位流,则
应返回一个整数,该整数由 n 个随机位元,或由0或1组成的长度为 nbits 的列表来构成.
如果支持随机整数,则
应返回范围在
上的 n 个随机整数的列表. 如果结果超出范围则发出警告信息.
如果支持随机实数,则
应返回范围在
上、精确度为 prec 的 n 个随机实数的列表. 如果结果超出范围或精确度错误则发出警告信息.
任何生成函数都可以返回
,其中 res 是正确类型的结果,gobj 是一个新的生成程序对象(反映任何状态变化).
起点的设定通过设定
的一个整数值 seed 完成
预期返回的是一个新的生成程序对象.
例: 乘法同余生成程序
在下面的例子中,将定义一个乘法同余生成程序. 乘法同余生成程序遵循递推关系
下面所定义的生成程序将仅用于实数的生成.
这里将生成程序

的选项设置为默认.
生成程序的初始化将提取乘数和模的值. 如果这两个值中任何一个不是正整数,初始化将失败.
由内核对
的调用实际是被包装在 Catch 当中. Throw 可以用于初始化代码,一旦出现问题程序能够轻松退出.
这里证明

生成实数.
实数生成程序将返回所需数量的实数及一个新的
生成程序. 新生成程序的起点根据递推关系更新
使用

生成程序生成10个实数.
| Out[74]= |  |
| Out[75]= |  |
例:Blum-Blum-Shub 生成程序
Blum-Blum-Shub 生成程序是一个生成用于加密目的的伪随机位元的二次同余方法[5]. 同余式是指定素数
和
的模
×
.
将生成程序

设置默认选项.
在必要时,生成程序的初始化将在调用实际生成程序之前,提取选项值并发出错误信息.

是一个位元生成程序,并确定位宽.
这里使用

生成程序生成了5个整数和5个实数.
| Out[86]= |  |
统计分布
由非均匀统计分布生成随机变量的总体思路是,先生成一个介于0和1的均匀分布的随机变量,然后计算该随机值在所求分布中的逆累积分布函数. 但在实际中,如果需要大量的随机变量,直接应用这个方案计算量将会很大,特别是在逆累积分布函数复杂或不能以封闭形式表示时尤为如此.
在这种情况下,表查找,基于分布关系的直接构造,或接受拒收方法往往比累积分布函数的直接求逆更有效. 在某些层上,这些方法仍将依赖均匀分布的 RandomReal 值,均匀分布的 RandomInteger 值,加权的 RandomChoice 的观察值,或者是这些值的组合. 因此,通过 SeedRandom 设置方法,将对从统计分布得到的随机观察值产生作用.
所有内置统计分布的随机观测值都可以利用 RandomVariate 生成. 在 Mathematica 中,由 RandomVariate 所用的用于多种分布的方法服从 Gentle [6] 或其它文献中所建议或描述的方法.
统计分布随机值的生成.
统计分布观察值通过 RandomVariate 获得. 这包括所有内置分布和构件,包括一元或多元分布,连续或离散分布,参数或导出分布,以及由数据定义的分布.
| Out[2]= |  |
对于连续分布,WorkingPrecision 可用于得到较高精度的数值,正如它可以得到一定范围的较高精度均匀分布数值一样.
这是精确度为30位

分布的一个随机变数.
| Out[3]= |  |
多元分布的随机值可以通过相同方式生成.
| Out[90]= |  |
| Out[91]= |  |
这里生成的随机值服从的分布由它的概率密度函数定义.
| Out[9]= |  |
下面的章节讨论了生成随机变数的方法,并就这些方法如何在 Mathematica 中使用提供了具体范例.
连续分布
对于逆累积分布函数只包含初等函数的一元分布,通常直接计算逆累积分布函数得到随机均匀分布. 这可以被看作对一个均匀分布随机变量的直接计算. 属于这一类别的连续分布包括 CauchyDistribution、 ExponentialDistribution、ExtremeValueDistribution、GumbelDistribution、LaplaceDistribution、LogisticDistribution、ParetoDistribution、RayleighDistribution、TriangularDistribution 和 WeibullDistribution.
从多个均匀分布变数或从非均匀分布变数直接构造一个单一随机变数的方法也被使用. 正态变数是通过 Box-Müller 方法 由成对的随机均匀分布变数成对产生. 例如,随机变数 HalfNormalDistribution、LogNormalDistribution 和MultinormalDistribution 通过正态变数的直接变换获得.
InverseGaussianDistribution 使用一种涉及到正态和均匀变数的接受-补充方法. 该方法来自于 Michael,Schucany 和 Haas,并在 Gentle 中介绍 [6]. MaxwellDistribution 变数由 ChiDistribution 变数构造. chi 变数自身通过 GammaDistribution 变数的特殊情形 ChiSquareDistribution 变数获得.
在大多数情况下 FRatioDistribution 从一个单一随机
值构造每一个随机值. 而对于小自由度的情况, FRatioDistribution 变数则通过成对的伽玛变数得到,以避免在
构造中可能出现的0作除数的情况.
NoncentralChiSquareDistribution[
,
],
,变数的生成使用的是
分布的可加性,以避免对非整数
进行昂贵的逆累积分布函数的计算. 可加性在 Johnson, Kotz 和 Balakrishnan 中给出[7].
时,一个非中心的
变数可以以均值为
,方差为1的一个正态变数的平方的形式得到.
时,非中心
变数以一个中心和一个非中心的
随机变量之和的形式得到. 当
时,如果
,
,
的分布为
. 当
时,这种关系不成立. 这种时候构造的是
,其中
,
,且当
趋向于0时,
是非中心
的极限分布. 极限分布
是泊松和
变量的混合,在0处存在一个非零概率质量,对于正数值则存在一个连续密度. NoncentralFRatioDistribution 变数通过一个中心和一个非中心
变数获得.
对于多元统计软件包中的 WishartDistribution,矩阵是通过 Smith 和 Hocking 方法生成的[8]. 该方法通过 chi 分布的对角线的项和正态分布的偏离对角线的项组成的矩阵构造 Wishart 矩阵.
NoncentralStudentTDistribution、HotellingTSquareDistribution 和MultivariateTDistribution 均由一元随机变数直接构造.
GammaDistribution、BetaDistribution 和 StudentTDistribution 在一定程度上使用接受-拒绝方法.
对于 GammaDistribution[
,
],当
时,生成的是指数变数. 否则使用 Cheng 和 Feast [9] 及 Ahrens 和 Dieter [10] 所提出的方法.
变数根据
的参数
和
的值在多种方法之间切换而构造. 如果两个参数均为1,将生成均匀随机变数. 如果
参数中一个值为1,则使用一种封闭形式的逆累积分布函数. 否则,RandomVariate 在 Jöhnk [11]、Cheng [12] 及 Atkinson [13]的所提出的接受-拒绝方法之间切换. 从下面的一个例子中可以看到接受-拒绝方法比由两个伽马进行构造的方法的优越性. 直接使用接受-拒绝方法,其速度几乎是伽马对构造的速度的两倍.

变数的直接构造与接受-拒绝方法的比较.
| Out[92]= |  |
| Out[93]= |  |
对于 StudentTDistribution,RandomVariate 所使用的方法是由 Bailey 所提出的极拒绝方法[14]. 从下面的例子中可以看到这种方法比从正态和
变数直接构造的效率更高. 对于100万个学生
变数,直接构造所用的时间几乎是极法的1.5倍.
对于学生

,直接构造与 Bailey 的极拒绝方法的比较.
| Out[94]= |  |
| Out[95]= |  |
离散分布
GeometricDistribution、BetaBinomialDistribution 和 BetaNegativeBinomialDistribution 使用直接构造. GeometricDistribution 变数以
的形式生成,其中
由UniformDistribution[0, 1] 得到. BetaBinomialDistribution 和BetaNegativeBinomialDistribution 由 BinomialDistribution 和NegativeBinomialDistribution 变数构造,其中概率参数取随机变数 BetaDistribution 的值.
使用时,生成随机整数的表查询是通过 RandomChoice 使用分布的的概率质量函数作为权重来被执行的. 每当由于所需样本大小的原因使得构造一列权重值代价过大时,多数离散分布将切换到其它方法. 例如,当
趋向于1时,LogSeriesDistribution[p] 切换至直接构造
,其中
和
在区间
均匀分布[15]. 根据参数不同和样本的大小 NegativeBinomialDistribution[n, p] 可以切换至泊松-伽马混合分布的构造,其中泊松变量的均值服从伽马分布[6].
如果计算 PDF 值的计算量相对于所需随机值的数量较小,BinomialDistribution, HypergeometricDistribution 和 PoissonDistribution 将依赖于从密度函数中的直接取样,否则将切换至接受-拒绝方法. 当溢出或下溢会在直接计算PDF值时发生时,接受-拒绝方法也允许变数的生成,从而扩大了对这些随机数生成的参数值范围.
二项式分布和超几何分布切换到由 Kachitvichyanukul 和 Schmeiser 提出的稍经修正的接受-拒绝方法. 基于 BTPE (二项式,三角,平行四边形,指数)算法 [16] 的接受-拒绝部分的二项式方法,实际上使用一种有三个区域的分段 majorizing 函数和一个三角形 minorizing 函数,以进行快速接受测试. 该 majorizing 和 minorizing 函数围绕坐标调整后的二项式分布的中心创建一个双平行四边形外壳,majorizing 函数的尾部在坐标调整后的二项式分布的尾部形成一个指数外壳. 使用 BTPE 明显优于构造一个查询表的一种情形是当仅需要少量观察值而查询表将很大时.
基于 Kachitvichyanukul 和Schmeiser [17] 的 H2PE 算法的接受-拒绝部分的超几何方法使用一种 majorizing 函数,该函数围绕一个坐标调整后的超几何密度有三个区域. 密度的中间部分是由矩形区域包围,分布的尾部是由指数函数界定.
PoissonDistribution 所使用的接受-拒绝方法由 Ahrens 和 Dieter [18] 提出. 接受和拒绝是利用 PoissonDistribution[
] 随
增加趋向于 NormalDistribution[
,
] 的趋势,使用离散正态变数而实现的.
ZipfDistribution 的随机值通过 Devroye [15] 所介绍的接受-拒绝方法生成. 除了使用基本算术,该方法使用成对的均匀变数和仅涉及一个 Floor 和非整数幂的测试来有效地获得 Zipf 分布的值.
定义分布生成程序
Mathematica 中囊括了一定数目的分布构造程序,这使得定义新的分布对象,并将其像其它任何分布一样对待成为可能. 这包括随机数的生成. 然而,假定您只对服从某一分布的随机值的生成感兴趣,并有实现这一任务的算法. 在这种情况下,只定义随机数生成程序可能会很有益. 这种分布生成程序的定义通过
的规则得到支持.
分布的定义被
的规则所支持.
预期返回一个给定长度的向量,其精确度数值给定.
| Random`DistributionVector[dist,n,prec] |
| 定义由 dist 中生成精确度为 prec 的 n 个观察值的规则 |
定义由各种分布进行随机生成的函数.
预期返回一个给定长度、给定精度数字的向量. 由于表达式 expr 不是一个完全定义的分布对象,数字将通过 RandomReal 或 RandomInteger 生成,而不是通过 RandomVariate. 如果精度为 Infinity,值将通过 RandomInteger 生成. 否则,值将通过 RandomReal 生成.
由各种分布生成随机值的规则通常通过一个 TagSet 在分布的头部定义. 分布本身可以包含参数. 下面是一个简单的例子:为
定义规则,代表的是区间
上的一个均匀分布.
的随机数现在可以通过 RandomReal 生成.
下面从

中给出一个机器精度数与一个20位有效数字的数.
| Out[97]= |  |
矩阵和高维张量也可以直接通过 RandomReal 生成. RandomReal 使用用于
的定义来生成所需的随机值总数,并将该总数分配到指定维数上.
这是一个

数的 3×4 数组.
| Out[98]= |  |
离散分布生成程序可以以类似的方式定义. 主要区别是
的精确度参数现在将是 Infinity.
的离散版本提供了一个简单的例子.
的随机值现在可以通过 RandomInteger 获得.
这是10 个

数.
| Out[100]= |  |
虽然前面的例子说明了定义分布生成程序的基本框架,但分布本身并不是特别令人感兴趣. 事实上,这两个例子中,更简单的做法是只由 RandomVariate 生成值,并将最后的结果乘以
,而不是将定义附加到一个新的符号上. 下面的范例将展示更复杂一些的分布,在这种情况下,将定义附加到一个符号上将更加有用.
例:通过求逆获得正态分布
对于生成一个通用的单变量的统计分布的随机值,书本上的定义包括两个步骤:
- 在区间
上生成一个均匀随机数
例:在圆盘上的均匀分布
也可用于定义多维分布的生成程序. 例如,对于单位圆盘上均匀分布的随机点来说,满足
的实数点集
是理想的. 这样的一个随机点可以通过下述方式进行构造:
- 生成一个均匀分布在
上的一个随机角度
- 生成一个均匀分布在
上的一个随机向量
例:Gibbs 采样程序
Gibbs 采样程序也可定义为分布生成程序. 考察一个混合了
分布和二项分布的 Gibbs 采样程序作为例子. 该采样程序的一个特别情形在先前的例子中有介绍. 这里,该分布将用两个参数 m 和
来定义.
定义一个 Gibbs 采样程序

.
对于先前构造的具体 Gibbs 采样程序,m 为16,
为2.
这是

,

时,采样程序得到的5个向量.
| Out[115]= |  |
参考文献
[1] Geman, S. and D. Geman. "Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images." IEEE Transactions on Pattern Analysis and Machine Intelligence 6, no. 6 (1984): 721-741.
[2] Casella, G. and E. I. George. "Explaining the Gibbs Sampler." The American Statistician 46, no. 3 (1992): 167-174.
[3] Matsumoto, M. and T. Nishimura. "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator." ACM Transactions on Modeling and Computer Simulation 8, no. 1 (1998): 3-30.
[4] Nishimura, T. "Tables of 64-Bit Mersenne Twisters." ACM Transactions on Modeling and Computer Simulation 10, no. 4 (2000): 348-357.
[5] Junod, P. "Cryptographic Secure Pseudo-Random Bits Generation: The Blum-Blum-Shub Generator." August 1999. http://crypto.junod.info/bbs.pdf
[6] Gentle, J. E. Random Number Generation and Monte Carlo Methods, 2nd ed. Springer-Verlag, 2003.
[7] Johnson, N. L., S. Kotz, and N. Balakrishnan. Continuous Univariate Distributions, Volume 2, 2nd ed. John Wiley & Sons, 1995.
[8] Smith, W. B. and R. R. Hocking. "Algorithm AS 53: Wishart Variate Generator." Applied Statistics 21, no. 3 (1972): 341-345.
[9] Cheng, R. C. H. and G. M. Feast. "Some Simple Gamma Variate Generators." Applied Statistics 28, no. 3 (1979): 290-295.
[10] Johnson, M. E. Multivariate Statistical Simulation. John Wiley & Sons, 1987.
[11] Jöhnk, M. D. "Erzeugung von Betaverteilten und Gammaverteilten Zufallszahlen." Metrika 8 (1964): 5-15.
[12] Cheng, R. C. H. "Generating Beta Variables with Nonintegral Shape Parameters." Communications of the ACM 21, no. 4 (1978): 317-322.
[13] Atkinson, A. C. "A Family of Switching Algorithms for the Computer Generation of Beta Random Variables." Biometrika 66, no. 1 (1979): 141-145.
[14] Bailey, R. W. "Polar Generation of Random Variates with the t-Distribution." Mathematics of Computation 62, no. 206 (1994): 779-781.
[15] Devroye, L. Non-Uniform Random Variate Generation. Springer-Verlag, 1986.
[16] Kachitvichyanukul, V. and B. W. Schmeiser. "Binomial Random Variate Generation." Communications of the ACM 31, no. 2 (1988): 216-223.
[17] Kachitvichyanukul, V. and B. W. Schmeiser. "Computer Generation of Hypergeometric Random Variates." Journal of Statistical Computation and Simulation 22, no. 2 (1985): 127-145.
[18] Ahrens, J. H. and U. Dieter. "Computer Generation of Poisson Deviates from Modified Normal Distributions." ACM Transactions on Mathematical Software 8, no. 2 (1982): 163-179.
[19] Matsumoto, M. and T. Nishimura. "Dynamic Creation of Pseudorandom Number Generators." In Proceedings of the Third International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing: Monte Carlo and Quasi-Monte Carlo Methods 1998, 56-69, 2000.