写在前面:本篇博客将介绍经典的伪随机数生成算法,我们将 重点讲解 LCG(线性同余发生器) 算法与马特赛特旋转算法,在此基础上顺带介绍 Python 的 random 模块。 本篇博客还带有练习,无聊到喷水的练习,咳咳…… 学完前面的内容你就会了解到 Python 的 Random 模块的随机数生成的实现,是基于马特赛特旋转算法的,比如 random_uniform 函数。而本篇博客提供的练习会让你实现一个基于 LCG 算法的random_uniform,个人认为还是比较有意思的。练习题的环境为 Google Colaboratory(K80 GPU)Jupyter Notebook:https://colab.research.google.com/
Ⅰ. 蒙特卡洛模拟(Monte Carlo simulation)
0x00 引入:概念说明
引入:在 非确定性模型(deterministic model)的 概率模型(stochastic model)中,用分析法难以下手,这时我们就可以利用 “蒙特卡洛模拟” (或称蒙特卡洛仿真) 去解决问题。
概念:“蒙特卡洛模拟” 是一种反复使用随机采样技术进行模拟,并从总体概率分布计算出所需数值结果的一种计算算法。
命名的由来
蒙特卡洛模拟是在二战期间,在原子弹研制的项目中用来模拟裂变物质的中子随机扩散现象,由美国数学家冯·诺伊曼和乌拉姆研发的模拟算法。蒙特卡洛是在欧洲摩纳哥的一个城市,这个城市在当时是非常著名的一个赌城。因为赌博的本质是算概率,而蒙特卡洛模拟正是以概率为基础的一种方法,所以用赌城的名字为这种方法命名。
蒙特卡洛模拟可作用于数值积分、概率分布计算、基于概率的优化等领域。例如:
Simulation:为了模拟抛硬币,在
Monte Carlo method:倒一盒硬币,数
Monte Carlo simulation:在
百度百科:百科举了个容易理解的例子
举个例子,假如筐里有100个苹果,让我每次闭眼拿1个,挑出最大的。
于是我随机拿1个,再随机拿1个跟它比,留下大的,再随机拿1个……我每拿一次,留下的苹果都至少不比上次的小。拿的次数越多,挑出的苹果就越大,但我除非拿100次(即把筐中的苹果都拿了一遍),否则无法肯定是否挑出了最大的那一个。
这个挑苹果的算法,就属于蒙特卡罗算法。
告诉我们样本容量足够大,则最接近所要求解的概率。
0x01 蒙特卡洛仿真常规模式
统共分为如下四步:
- STEP1:定义要进行采样的区域(Range)
- STEP2:对定义的区域进行随机采样(Random sampling)
- STEP3:对收集的样本进行确定型计算(Deterministic computation)
- STEP4:统计结果并导出近似值(Approximate)
举个例子:两个简单的蒙特卡罗模拟的实例
① 掷两次骰子点数之和为
- STEP4:与数学计算的概率进行比较
② 求圆周率:
- 在提取出的点中,统计圆内部点的个数。
Ⅱ. 伪随机数发生器(PRNG)
0x00 引入:随机数生成问题
为了给蒙特卡洛模拟生成随机数,我们需要置随机数 “种子”。
这里我们选择采用伪随机数发生器 (Pseudo Random Number Genarator) ,简称
伪随机数并非完全随机
伪随机数发生器用于在系统需要随机数的时候,通过一系列种子值计算出来的伪随机数。因为生成一个真正意义上的“随机数”对于计算机来说是不可能的,伪随机数也只是尽可能地接近其应具有的随机性,但是因为有“种子值”,所以伪随机数在一定程度上是可控可预测的。
本章我们将介绍以下几种经典的伪随机数生成算法:(我们将重点讲解前两种算法)
- 线性同余发生器(Linear congruential generator)
- 马特赛特旋转演算法(Mersenne Twister)
- 冯·诺伊曼平方取中法(John von Neumann mid-square method)
- 乘法取中法(mid-product method)
- 常数乘子法(constant multiplier method)
- 斐波那契法(Fibonacci Method)
- 时滞斐波那契法(Lagged Fibonacci method)
- 移位法(Shift Method)
0x01 线性同余发生器(LCG)
线性同余发生器(Linear Congruential Generator),简称
它能产生具有不连续计算的伪随机序列的分段线性方程,生成器由循环关系定义如下:
为保持
其中参数 a, C, m 比较敏感,它们直接影响了伪随机数产生的质量,所以参数的取值非常重要!
0x02 马特赛特旋转算法(MersenneTwister)
“马特赛特旋转算法是当今生成随机数质量最好的算法”
马特赛特旋转算法(MersenneTwister),简称
补充:梅森素数
梅森素数的意思是形如的素数。利用梅森素数的性质可以设计出周期长度为梅森素数长度的随机数周期。比如目前Python、C++11等语言当中用的随机数计算包都是用的这种算法。目前常用的版本周期是,这是一个巨大的天文数字。
即使你从来没有听过这个算法,但你一定用过 Python 中的 random 模块:
Python 中的 random 模块就是采用了马特赛特旋转算法实现随机数的 ,不仅 Python 用,PHP、Perl 等流行编程语言内置的
- 第一阶段:获得基础的马特赛特旋转链
- 第二阶段:对于旋转链进行旋转算法
- 第三阶段:对于旋转算法所得的结果进行处理
- 马特赛特旋转能更快的生成高质量的伪随机数,随机数的反复周期取决于马特赛特的因数。
- 只需位运算就能实现的算法,速度异常迅速。
- 从密码学的角度来看,马特赛特旋转算法仍然是不安全的,当你知道随机数的特性时,仅利用有限的随机数(624个)就可以知道当前生成器所处状态并预测出即将出现的随机数。
0x03 冯·诺伊曼平方取中法(Mid-square Method)
平方取中法(Mid-square Method)是最早出现的伪随机数发生器,最早由冯诺依曼提出。
平方取中法的意思如同其名:
- 平方:是将一个 2s 位十进制随机数开平方后得到的一个 4s 位数(不足 4s 时高位补 0)。
- 取中:去头截尾,取中间 2s 位数作为一个新的随机数。
并将此数规范化(化为小于1的2s位数值),即第一个
产生伪随机数数列:
平方取中法的优点:
- 易于实现,内存占用少,计算简单。
平方取中法的缺点:
- 存对小数目偏倚现象,均匀性差,对初始数据依赖大,很难说明取什么样的种子值可以保证有足够长的周期。
- 很容易出现重复元素的段循环,而且容易退化为一常数,甚至退化为零,如果某个元素退化为0,那么后面所有元素都将会是0。
平方取中法的变形和推广:
- 乘法取中法(mid-product method):
若要产生具有10进制 2s 位的伪随机数序列,任意选取两个初始随机数
产生伪随机数数列:
- 常数乘子法(constant multiplier method):
递推公式:
产生伪随机数数列:
0x04 斐波那契法(Fibonacci Method)
遗憾的是,该方法并不是一个好的伪随机数发生器,该方法基于 Fibonacci 数列,其递推公式为:
该方法有
斐波那契法的变形和推广(了解):
- 时滞斐波那契法(Lagging Fibonacci Method):
时滞斐波那契生成器,简称
其中,新项由两个老项计算生成。
0x05 移位法(Shift Method)
计算机善于进行移位等逻辑运算,移位法就是利用计算机的这个特点去实现的。移位法是关于平方取中法的另一种改进形式,它是移位运算于指令相加运算相结合的迭代过程。方法如下:
取一个初始值
产生伪随机数数列:
移位法运算速度快,但是对初始值的依赖性很大,一般初始值不能取的太小,初始值选的不好会使伪随机数序列长度较短。同时独立性交叉,随后随机数序列周期
Ⅲ. Python 的随机数生成函数(Random 模块)
0x00 引入:Random 模块介绍
我们刚才讲解马特塞特旋转算法时提到了,Python 中的 Random 模块就是采用该算法实现的。
其生成
提供了可提取分布的函数:Uniform,Normal,lognormal,negative exponential,gamma, beta
使用前需引入 random 模块:
0x01 random – 生成 0.0 ~ 1.0 间的随机数
作用:随机生成一个
代码演示:生成实数
运行结果如下:
0x02 uniform – 生成指定范围的随机数
作用:随机生成一个指定范围的实数。
代码演示:生成一个 100~300 之间的数
运行结果如下:
0x03 randint – 生成指定范围的整数随机数
作用:随机生成一个指定范围的整数。
代码演示:生成一个100~300之间的整数
运行结果如下:
0x04 choice – 在样本集中随机选择
作用:在样本集中随机选择一个样本
代码演示:从样本集
运行结果如下:
0x05 sample – 在样本集中随机选择多次
作用:执行 n 次 “在样本集中随机选择一个样本” 。
代码演示:从样本集
运行结果如下:
Ⅳ. 练习(Assignment)
练习1:绘制对数直方图(利用LCG实现)
- 在 Colaboratory 中使用 Python 实现。
- 利用线性同余发生器(LCG)实现。
-
- 使用 Python 的 Matplotlib 模块,以生成直方图的形式展现分布,如下所示:
声明:本站部分文章及图片源自用户投稿,如本站任何资料有侵权请您尽早请联系jinwei@zod.com.cn进行处理,非常感谢!