fortran产生随机数方法介绍

fortran产生随机数方法介绍（附代码）
注意：现在计算机产生的随机数都是伪随机数。

1.0-1之间均匀分布的随机数
random_number(x) 产生一个0到1之间的随机数（x可以是向量），但是每次总是那几个数。

用了random_seed ()后，系统根据日期和时间随机地提供种子，使得随机数更随机了。

program random
implicit none
real :: x
call random_seed () ! 系统根据日期和时间随机地提供种子
call random_number (x) ! 每次的随机数就都不一样了
write(*,*) x
stop
end program random
2.任意区间均匀分布的随机数
function my_random (lbound,ubound)
implicit none
real :: lbound,ubound
real :: len
real :: my_random
real :: t
len=ubound-lbound !计算范围大小
call random_number(t) !t是0-1之间的随机数
my_random=lbound+len*t
return
end
注意：在循环外call random_seed()
3.产生一个随机数数组，只需加一个循环即可
function my_random (lbound,ubound)
implicit none
real :: lbound,ubound
real :: len
integer size
real :: my_random(size) !size代表数组元素的个数
real :: t
integer i
len=ubound-lbound !计算范围大小
do i=1,10
call random_number(t) !t是0-1之间的随机数
my_random(i)=lbound+len*t !把t转换成lbound-ubound间的随机数
end do
return
end
注意：同理在循环外call random_seed()
4.标准正态分布随机数/高斯分布随机数
(1)徐士良的那本程序集里介绍了正态分布随机数产生的原理，不过他的方法只能产生较为简
单的随机数，随机数的质量并不高，特别是随机数的数目较多时。

(2)Box 和 Muller 在 1958 年给出了由均匀分布的随机变量生成正态分布的随机变量的算法。

设 U1, U2 是区间 (0, 1) 上均匀分布的随机变量，且相互独立。

令
X1 = sqrt(-2*log(U1)) * cos(2*PI*U2);
X2 = sqrt(-2*log(U1)) * sin(2*PI*U2);
那么 X1, X2 服从 N(0,1) 分布，且相互独立。

等于说我们用两个独立的 U(0,1) 随机数得到了
两个独立的 N(0,1)随机数。

值得说明的是，该方法产生的随机数质量很高。

嘻嘻，本人亲自验证过。

SAS和蒙特卡罗模拟_随机数
一、为什么选择SAS做蒙特卡罗模拟？
为什么要用SAS做蒙卡？首先，对我来说，我只会用SAS，而且打算用SAS完成我所有的工作。

当然，其他一些通用的理由有(Fan, etc.,2002)：
1.蒙卡是个计算密集的活，而SAS Base、SAS Macro、SAS/IML强大而灵活
的编程能力能满足这一要求；
2.做蒙卡时要用到大量的统计/数学技术，而SAS就内置了大量的统计/数
学函数(在 SAS/Stat和SAS/ETS)；用Fortran或C++当然也是非常好的
主意，只是他们缺少内置的统计函数，代码就要冗长复杂很多。

二、什么是蒙卡？一个启发性例子
好，开始，什么是蒙卡？了解它背景知识的最好办法当然是
wiki-Monte_Carlo_method。

蒙特卡罗是位于摩洛哥的一家赌场，二战时，美国Los Alamos国家实验室把它作为核裂变计算机模拟的代码名称。

作为模拟方法，蒙卡以前就叫统计抽样(statistical sampling)，我们感兴趣的结果因为输入变量的不确定而不可知，但如果能依概率产生输入变量的样本，我们就可以估计到结果变量的分布。

跟蒙卡对应的，还有一种模拟技术叫系统模拟，包括排队、库存等模型，这些模型都跟随时间推移而出现的事件序列有关。

下面举个蒙卡的例子，来自Evans, etc.( 2001)的超级简化版。

假设一家企业，利润是其需求量的函数，需求是随机变量。

为了简化讨论，假定利润就是需求的两倍。

这里输入变量就是不可控的需求，结果变量就是我们感兴趣的利润。

假设需求以相同的概率取10、20、30、40、50、60这六种情况。

在这样的简化下，我们就可以投一枚均匀的骰子来产生需求的样本，如果点数为1，对应得需求就是10，点数为2，需求就是20，以下类推。

这样，我们的模拟过程就是：
1.投骰子；
2.根据骰子的点数确定需求量；
3.根据需求量，求利润。

掷10次骰子，假设我们的模拟结果如下：
重复次数骰子点数需求量利润
1 5 50 100
2 3 30 60
3 3 30 60
4 6 60 120
5 1 10 20
6 3 30 60
7 4 40 40
8 5 50 100
9 2 20 20
10 5 50 50
这样通过模拟需求的样本，我们对结果利润的分布也就有所知晓，比如平均利润可以算出就是63。

蒙卡一个重要的步骤就是生成随机数，这里我们是用投骰子来完成。

三、又一个例子：利用蒙特卡罗模拟方法求圆周率∏(pi)
再举个很有名的例子，就是估计圆周率∏的值，来自Ross(2006)。

这个试验的思路正好可以帮我们温习一下几何概型的概念。

我们知
道概率的古典概型，就是把求概率的问题转化为计数：样本空间由
n个样本点组成，事件A由k个样本点组成，则事件A的概率就是
k/n。

考虑到概率和面积在测度上具有某种共性，几何概型的基本想
法是把事件跟几何区域相对应，用面积来计算概率，其要点是：
1.认为样本空间Ω是平面上的某个区域，其面积记为υ(Ω)；
2.向区域Ω随机投掷一点，该点落入Ω内任何部分区域内的
可能性只与这部分区域的面积成比例，而与这部分区域的位置和形状无
关；
3.设事件A是Ω的某个区域，面积为υ(A)，则向区域Ω上随机投掷一
点，该点落在区域A的可能性称为事件A的概率，P(A)=υ(A)/υ(Ω)。

扯远了，回到用蒙卡估计圆周率∏的实验思路。

假设样本空间是一个边长为2
面积为4的正方形，我们感兴趣的事件是正方形内的一个半径为1面积为∏的圆，所以向正方形内随机投掷一点，落在圆里面的概率为∏/4。

实验的思路如下：
1.1）生成随机数——生成n个均匀落在正方形内的点；
2.2）对落在正方形内的n个点，数一数正好落在圆里面的点的个数，假设
为k（另外n-k个点就落在圆外面的正方形区域内）。

3.3）k/n就可以大致认为是圆的面积与正方形的面积之比，另其等于∏/4，
就可以求出圆周率∏的估计值。

又，Forcode提供了利用电子表格Excel求解以上问题的详细过程，有兴趣可以看看。

另外，这里也有一个详细的说明，用Mathematica实现。

四、随机数很重要（以及下期预告）
在第一个例子中，我强调了骰子要是均匀的。

那里骰子就是我们的随机数生成器，骰子的均匀程度就是我们随机数的“随机”成分。

在上面的10次试验中，投出了3次点数“3”和3次点数“5”，然后点数“1”、“2”、“4”、“6”各出现一次——因为只是重复了10次，这样我们对骰子的均匀程度不好评估，但如果重复无数次——假如是十万次，如果还出现类似比例的结果，那么就有理由怀疑这粒骰子的均匀程度了。

如果骰子灌了铅，比如投出点数“6”的可能性从六分之一上升到6分之三，那么根据这粒骰子来做模拟，我们就有高估需求以及利润的危险。

下次就讲讲随机数的生成原理，当然结合SAS来讲，或者也会提一下Excel。

五、其他软件包
在商业领域，基于Excel的插件比较兴盛，比如Crystal Ball应用就较为普遍，有试用版可以下载。

其他竞争产品有@Risk、Risk Solver等。

现在据说Risk Simulator也开始兴盛起来，大有后来居上之势。

他的开发者Johnathan Mun还在Wiley Finance系列有一本书，Modeling Risk: Applying Monte Carlo Simulation, Real Options Analysis, Forecasting, and Optimization Techniques。

S语言（R或者S-Plus）由于其面对对象的特性，加之丰富的内置函数和诸多用户提供的库，使得R或者S-Plus也是蒙卡研究的得力工具——或者比SAS更有前景。

这个我不熟，略之。

Random Numbers
最简单、最基本、最重要的随机变量是在[0,1]上均匀分布的随机变量。

一般地，我们把[0,1]上均匀分布随机变量的抽样值称为随机数，其他分布随机变量的抽样都是借助于随机数来实现的。

以下谈的都是所谓“伪随机数”(Pseudo Random Numbers)。

产生随机数，可以通过物理方法取得（很久很久以前，兰德公司就曾以随机脉冲源做信息源，利用电子旋转轮来产生随机数表），但当今最为普遍的乃是在计算机上利用数学方法产生随机数。

这种随机数根据特定的迭代公式计算出来，初值确定后，序列就可以预测出来，所以不能算是真正的随机数（就成为“伪随机数”）。

不过，在应用中，只要产生的伪随机数列能通过一系列统计检验，就可以把它们当成“真”随机数来用。

现在大多软件包内置的随机数产生程序，都是使用同余法(Congruential Random Numbers Generators)。

“同余”是数论中的概念。

0.预备知识：同余
捡回小学一年级的东西："4/2=2"读作“4除以2等于2”，或者，“2除4等于2”。

还有求模的符号mod(number,divisor)，其中，number是被除数（在上式中，为4），divisor是除数（上式中的2）。

这样的约定对SAS和Excel都通用，如mod(4,13)=4，mod(13,4)=1。

现在我们可以开讲“同余”了。

设m是正整数，用m去除整数a、b，得到的余数相同，则称“a与b关于模m同余”。

上面的定义可以读写成，对整数a、b 和正整数m，若mod(a,m)=mod(b,m)，则称“a与b关于模m有相同的余数（同余）”，记做a≡b(mod m)（这就是同余式）。

举个例子，mod(13,4)=1，mod(1,4)=1，则读成13和1关于模4同余，记做13≡1(mod 4)。

当然，同余具有对称性，上式还可以写成1=13(mod 4)。

a≡b(mod m)的一个充要条件是a=b+mt，t是整数，比如13=1+3*4。

a=b+mt可以写成(a-b)/m=t，即m能整除(a-b)。

这些就够了。

更多基础性的介绍，可以参考《同余（数据基础）》。

1.乘同余法
同余法是一大类方法的统称，包括加同余法、乘同余法等。

因为这些方法中的迭代公式都可以写成上面我们见过的同余式形式，故统称同余法。

常用的就是下面的乘同余法(Multiplicative Congruential Generator.)。

符号不好敲，做些约定，如R(i)就是R加一个下标i。

乘同余法随机数生成器的同余形式如下：R(i+1)=a*R(i) (mod m)。

这个迭代式可以写成更直观的形式，R(i+1)=mod[a*R(i),m]，其中初值R(0)称为随机数种子。

因为mod(x,m)总是等于0到m-1的一个整数，所以最后把R(i+1)这个随机序列都除以m，就可以得到在[0,1]上均匀分布的随机数。

下面用电子表格演示一遍，假设随机数种子R(0)=1，a=4，m=13：
上面显示的是公式（Excel中，公式与计算结果的转换，用快捷键ctrl+~实现），结果看起来是这样的，E列就是我们想要的在[0,1]上均匀分布的随机数：
上表中，a*R(1)=16，R(2)=3，m=13，一个同余式就是R(2)=a*R(1) (mod m)，3=16(mod 13)，或者说，16-3能够被13整除——同余法就是这意思了。

说“乘同余法”，要点在于a*R(i)是乘法形式。

在SAS系统中，a=397,204,094，m = 2^31-1=2147483647（一个素数），随机种子R(0)可以取1到m-1之间的任何整数。

2.伪随机数的检验
现在内置于各大软件包的随机数生成器都经过了彻底的检验，我们当然可以安心地使用这个黑箱。

或则我们也可以多了解些。

一个好的“伪”随机数列，应该看起来就跟从[0,1]上均匀分布的随机数列中随机抽取出来的一样。

两个比较直观的检验有：
1.均匀性检验——所有的数是否均匀地分布在[0,1]区间；
2.独立性（或不相关）检验——序列中的数不存在序列相关，每个数都跟
其前后出现的数独立无关。

一个例子是如0.1、0.2、0.3、0.4、0.5，
这里每一个相继的数都正好比它前面的一个数大0.1，这样这个数列就
似乎有了某种“格局”。

具体的检验方法就略之了，更多请参见徐钟济（1985）。

3.生成其他概率分布的随机数
前面提到过，我们把[0,1]上均匀分布随机变量的抽样值称为随机数，其他分布随机变量的抽样都是借助于随机数来实现的。

对其他连续型分布，（累积）分布函数为F(x)，r是在[0,1]上均匀分布的随机变量，另r=F(x)，解出的x就是该连续型分布的随机抽样。

由于x可以写成r的反函数，该变换就称作“逆变换”(Inverse Transformation method)。

对离散型分布，思路类似，不过繁琐些，具体可以见徐钟济（1985）、Ross（2006）和Evans等（2001）。

大多数相关软件包都会提供各种分布的随机数生成函数，知道有这回事就可以。

最重要的，是我们陈述一类问题时，知道需要用哪种概率分布来描述Evans等（2001）：
常用的连续分布
1.均匀分布(Uniform Distribution)描述在某最小值和最大值之间所有结
果等可能的随机变量的特性。

2.正态分布(Normal Distribution)是对称的，具有中位数等于均值的性
质，就是我们熟悉的钟形形状。

各种类型的误差常常是正态分布的。

最
后，作为中心极限定理的推论，大量具有任意分布的随机变量的平均数
也呈正态分布。

3.三角形分布(Triangular Distribution)由三个参数来定义，最大值、最
小值和最可能值。

临近最可能值的结果比那些位于端点的结果有更大的
出现机会。

4.指数分布(Exponential Distribution)常用来构建在时间上随机重现的
事件的模型，如顾客到达服务系统的时间间隔、机器元件失效前的工作
时间等。

它的主要性质是“无记忆性”(Memoryless) ，即当前时间对
未来结果没有影响。

例如，不论机器原先已经运转了多长时间，它继续
运转直至出现故障的时间间隔总有同样的分布。

5.对数正态分布(Lognormal Distribution)：若随机变量x的自然对数是
正态的，则x就服从对数正态分布。

对数正态分布的常见例子是股票价
格。

常用的离散分布
1.贝努利分布(Bernoulli Distribution)描述的是只有两个以常数概率出
现的可能结果的随机变量的特性；
2.二项分布(Binomial Distribution)给出每次试验成功概率为p的n次独
立重复贝努利实验的模型；
3.泊松分布(Poisson Distribution)用于建立某种度量单位内发生次数的
模型的一种离散分布，如某时段内某事件发生的次数。

其他有用的分布如伽玛分布(Gamma Distribution)、威布尔分布(Weibull Distribution)、贝塔分布(Beta Distribution)略之。

4.下期预告
上次落了些东西。

说到蒙卡与C++，在数量金融领域，不能不提的一本书就是Mark Joshi的C++ Design Patterns and Derivatives Pricing (Cambridge Uni. Press, 2004)，面试中一定要说自己读过的，这样人家以为你至少会用C++做一个蒙特
卡罗模拟。

这个系列只讲SAS，下次先讲如何用SAS生成随机数，然后就是具体的项目。

SAS随机数函数及CALL子程序
********************************************************************
******************************************************
《SAS和蒙特卡罗模拟（1）：开篇》——简介，通过例子建立起蒙卡的直观概念；参考软件包及书目
《SAS和蒙特卡罗模拟（2）：随机数基础》——随机数入门；参考书目
********************************************************************* ********************************************************************* ********
一、SAS随机数函数和CALL子程序
SAS系统产生随机数，两种方式，利用SAS函数(Functions)或者CALL子程序(CALL Routines)，它们的语法格式是（具体的区别容后讨论）：
SAS可用的随机数函数列表如下（可以参见SAS Help and Documentation-SAS Products-Base SAS-SAS Language Dictionary-Functions and CALL Routines-Functions and CALL Routines by Category）：
上面的随机函数，除了normal和uniform，都可以由直接相应的CALL子程序调用。

二、SAS随机数函数和CALL子程序：比较
用SAS随机数函数同时创建的多个随机数变量其实都属于同一个随机数列。

这话费解，一个例子先，创建两个随机数变量，各包含3个记录，其中x1的种子为123，x2的种子为456：
data ranuni(drop=i);
retain seed1 123 seed2 456;
do i=1 to 3;
x1=ranuni(seed1);
x2=ranuni(seed2);
output;
end;
run;
proc print data=ranuni;run;
结果为：
Obs seed1 seed2 x1 x2
1 123 456 0.75040 0.32091
2 12
3 456 0.17839 0.90603
3 123 456 0.35712 0.22111
好像没什么异样。

我们把上面的x1增加为6个记录：
data ranuni2(drop=i);
retain seed1 123;
do i=1 to 6;
x1=ranuni(seed1);
output;
end;
run;
proc print data=ranuni2;run;
结果如下，把上下用红色标出的数字对照看一看：
Obs seed1 x1
1 123 0.75040
2 12
3 0.32091
3 123 0.17839
4 123 0.90603
5 123 0.35712
6 123 0.22111
什么意思？在第一段代码中，x2的种子根本不起作用，把x2的记录安插到x1里，看起来就是用种子123产生的随机数列加长了而已。

x2的第一个值并不是由种子456产生的，而是产生第一个x1后的下一个值，x1、x2其实属于同一个随机数列，尽管x2的种子被指定为456，而x1的被指定为123。

现在就可以重复上面的一句话：用SAS随机数函数同时创建的多个随机数变量其实都属于同一个随机数列。

用CALL子程序就能够同时产生独立的随机数列。

data ranuni3(drop=i);
retain seed3 123 seed4 456;
do i=1 to 3;
call ranuni(seed3,x3);
call ranuni(seed4,x4);
output;
end;
run;
proc print data=ranuni3;run;
结果如下：
Obs seed3 seed4 x3 x4
1 1611463328 736440516 0.75040 0.34293
2 689153326 774069794 0.32091 0.36045
3 38308885
4 686944750 0.17839 0.31988
以上x3就是x1。

x1和x3的初始种子都是123，但x3那个结果显示的种子是当前种子值。

要在SAS随机数函数语句中显示随机种子的当前值，可以由以下代码给出：
data ranuni4(drop=i);
retain seed1 123;
do i=1 to 3;
x1=ranuni(seed1);
seed=x1*(2**31-1);
output;
end;
run;
proc print data=ranuni4;
var seed1 seed x1;
run;
结果如下，可以跟上面由CALL子程序得出的结果对照：
Obs seed1 seed x1
1 123 1611463328 0.75040
2 12
3 689153326 0.32091
3 123 38308885
4 0.17839
---------参考资料---------
1.Xitao Fan, etc..Monte Carlo Studies: A Guide for Quantitative
Researchers. SAS Institute Inc.,2002
2.朱世武《SAS编程技术与金融数据处理》，北京：清华大学出版社，2003 Matlab() 随机数生成方法：
第一种方法是用random 语句，其一般形式为
y = random('分布的英文名',A1,A2,A3,m,n)，
表示生成m 行n 列的m × n 个参数为( A1 , A2 , A3 ) 的该分布的随机数。

例如:
(1) R = random('Normal',0,1,2,4): 生成期望为0,标准差为1 的(2 行4 列)2× 4 个正态随机
数
(2) R = random('Poisson',1:6,1,6):依次生成参数为1 到6 的(1 行6 列)6 个Poisson 随机数
第二种方法是针对特殊的分布的语句：
一．几何分布随机数（下面的P，m 都可以是矩阵）
R = geornd(P) （生成参数为P 的几何随机数）
R = geornd(P,m)（生成参数为P 的× m 个几何随机数）
１
R = geornd(P,m,n)（生成参数为P 的m 行n 列的m × n 个几何随机数）例如
(1)R = geornd(1./ 2.^(1:6)) ( 生成参数依次为1/2,1/2^2,到1/2^6 的6 个几何随机数)
(2)R = geornd(0.01,[1 5]) (生成参数为0.01 的（１行５列）5 个几何随机数).
二．Beta 分布随机数
R = betarnd(A,B)（生成参数为A,B 的Beta 随机数）
R = betarnd(A,B,m)（生成× m 个数为A,B 的Beta 随机数）
１
R = betarnd(A,B,m,n)（生成m 行n 列的m × n 个数为A,B 的Beta 随机数）.
三．正态随机数
R = normrnd(MU，SIGMA)（生成均值为MU，标准差为SIGMA 的正态随机数）
R = normrnd(MU，SIGMA,m)（生成× m 个正态随机数）
１
R = normrnd(MU，SIGMA,m,n) （生成m 行n 列的m × n 个正态随机数）
例如
(1) R = normrnd(0,1,[1 5]) 生成5 个正态(0,1) 随机数
(2) R = normrnd([1 2 3;4 5 6],0.1,2,3)生成期望依次为[1,2,3;4,5,6], 方差为0.1 的2× 3 个
正态
随机数．
四．二项随机数：类似地有
R = binornd(N,P)R = binornd(N,P,m) R = binornd(N,p,m,n)
例如
n = 10:10:60; r1 = binornd(n,1./n)或r2 = binornd(n,1./n,[1 6]) （都生成参数分别为
1 1
), L, ( 60, ) 的６个二项随机数．
(10,
10 60
五．自由度为V 的χ 2 随机数：
R = chi2rnd(V)R = chi2rnd(V R = chi2rnd(V
,m) ,m,n)
六．期望为MU 的指数随机数（即Exp 随机数）：
1
MU
R = exprnd(MU)R = exprnd(MU,m)R = exprnd(MU,m,n)
七．自由度为V1，V2 的 F 分布随机数：
R = frnd(V1，V2)R = frnd(V1，V2,m)R = frnd(V1，V2,m,n)
八．Γ ( A, λ ) 随机数：
R = gamrnd（A,lambda）R = gamrnd（A,lambda,m)R = gamrnd（A,lambda,m,n) 九．超几何分布随机数：
R = hygernd(N,K,M)R = hygernd(N,K,M,m)R = hygernd(N,K,M,m,n)
十．对数正态分布随机数
R = lognrnd(MU，SIGMA)R = lognrnd(MU，SIGMA,m)R = lognrnd(MU，SIGMA,m,n) 十一．负二项随机数：
R = nbinrnd(r,p)R = nbinrnd(r,p,m)R = nbinrnd(r,p,m,n)
十二．Poisson 随机数：
R = poissrnd(lambda) R = poissrnd(lambda,m)R = poissrnd(lambda,m,n) 例如，以下3 种表达有相同的含义：lambda = 2;R = poissrnd(lambda,1,10)
（或R = poissrnd(lambda,[1 10])或R = poissrnd(lambda(ones(1,10)))
十三．Rayleigh 随机数：
R = raylrnd(B) R = raylrnd(B,m)R = raylrnd(B,m,n)
十四．V 个自由度的t 分布的随机数：
R = trnd(V) R = trnd(V,m)R = trnd(V,m,n)
42
十五．离散的均匀随机数：
R = unidrnd(N) R = unidrnd(N,m)R = unidrnd(N,m,n)
十六．[A,B] 上均匀随机数
R = unifrnd(A,B) R = unifrnd(A,B,m)R = unifrnd(A,B,m,n)
例如unifrnd(0,1:6)与unifrnd(0,1:6,[1 6]) 都依次生成[0,1] 到[0,6]的６个均匀随机数．: 十七．Weibull 随机数
R = weibrnd(A,B) R = weibrnd(A,B,m)R = weibrnd(A,B,m,n)。