fortran产生随机数方法介绍

《fortran产生随机数方法介绍》由会员分享，可在线阅读，更多相关《fortran产生随机数方法介绍（9页珍藏版）》请在金锄头文库上搜索。

1、fortran 产生随机数方法介绍(附代码)注意：现在计算机产生的随机数都是伪随机数。10-1之间均匀分布的随机数random_number(x)产生一个0到1之间的随机数(x可以是向量)，但是每次总是那几个数。 用了 random_seed ()后，系统根据日期和时间随机地提供种子，使得随机数更随机了。program randomimplicit nonereal : xcall random_seed ()!系统根据日期和时间随机地提供种子call random_number (x) !每次的随机数就都不一样了write(*,*) xstopend program random2任意区间均

2、匀分布的随机数function my_random (lbound,ubound)implicit nonereal : lbound,uboundreal : lenreal : my_randomreal : tlen二ubound-lbound !计算范围大小call random_number(t) !t 是 0-1 之间的随机数my_random=lbound+len*treturnend注意：在循环外call random_seed()3产生一个随机数数组，只需加一个循环即可function my_random (lbound,ubound)implicit nonereal :

3、lbound,uboundreal : leninteger sizereal : my_random(size) !size代表数组元素的个数real : tinteger ilen二ubound-lbound !计算范围大小do i=1,10call random_number(t) !t 是 0-1 之间的随机数my_random(i)=lbound+len*t !把 t 转换成 lbound-ubound 间的随机数 end doreturnend注意：同理在循环外call random_seed()4标准正态分布随机数/高斯分布随机数(1)徐士良的那本程序集里介绍了正态分布随机数产生

4、的原理，不过他的方法只能产生较为简 单的随机数，随机数的质量并不高，特别是随机数的数目较多时。(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

5、和蒙特卡罗模拟_随机数一、为什么选择 SAS 做蒙特卡罗模拟？为什么要用SAS做蒙卡？首先，对我来说，我只会用SAS，而且打算用SAS完成我所有的工作。当然， 其他一些通用的理由有(Fan, etc.,2002):1. 蒙卡是个计算密集的活，而SAS Base、SAS Macro、SAS/IML强大而灵活的编程能力能满足这一 要求；2. 做蒙卡时要用到大量的统计/数学技术，而SAS就内置了大量的统计/数学函数(在SAS/Stat和 SAS/ETS)；用Fortran或C+当然也是非常好的主意，只是他们缺少内置的统计函数，代码就要冗长复 杂很多。二、什么是蒙卡？一个启发性例子好，开始，什么是蒙卡

6、？ 了解它背景知识的最好办法当然是wiki-Monte_Carlo_method。蒙特卡罗是位 于摩洛哥的一家赌场，二战时，美国Los Alamos国家实验室把它作为核裂变计算机模拟的代码名称。 作为模拟方法，蒙卡以前就叫统计抽样(statistical sampling)，我们感兴趣的结果因为输入变量的不确定 而不可知，但如果能依概率产生输入变量的样本，我们就可以估计到结果变量的分布。跟蒙卡对应的， 还有一种模拟技术叫系统模拟，包括排队、库存等模型，这些模型都跟随时间推移而出现的事件序列有 关。下面举个蒙卡的例子，来自 Evans, etc.( 2 00 1)的超级简化版。 假设一家企业，利

7、润是其需求量的函数，需求是随机变量。为了简化讨论，假定利润就是需求的两倍。 这里输入变量就是不可控的需求，结果变量就是我们感兴趣的利润。假设需求以相同的概率取 10、 20、 30、 40、 50、 60 这六种情况。在这样的简化下，我们就可以投一枚均匀的骰子来产生需求的样本，如 果点数为 1，对应得需求就是 10，点数为 2，需求就是 20，以下类推。这样，我们的模拟过程就是:1.投骰子；2. 根据骰子的点数确定需求量；3. 根据需求量，求利润。掷10 次骰子，假设我们的模拟结果如下:重复次数骰子点数需求量利润15501002330603330604660120511020633060744

8、04085501009220201055050这样通过模拟需求的样本，我们对结果利润的分布也就有所知晓，比如平均利润可以算出就是 63。蒙卡 一个重要的步骤就是生成随机数，这里我们是用投骰子来完成。三、又一个例子：利用蒙特卡罗模拟方法求圆周率n(pi)再举个很有名的例子，就是估计圆周率口的值，来自Ross(2006)。这个试验的思路正好可以帮我们温习 一下几何概型的概念。我们知道概率的古典概型，就是把求概率的问题转化为计数：样本空间由 n 个样本点组成，事件A由k个样本点组成，则事件A的概率就是k/n。考虑到 概率和面积在测度上具有某种共性，几何概型的基本想法是把事件跟几何 区域相对应，用面积

9、来计算概率，其要点是：认为样本空间Q是平面上的某个区域，其面积记为u(Q)； 向区域Q随机投掷一点，该点落入Q内任何部分区域内的可能性只与这部分区域的面积成比例，而与这部分区域的位置和形状无关；设事件A是Q的某个区域，面积为u(A)，则向区域Q上随机投掷 一点，该点落在区域A的可能性称为事件A的概率，P(A)二u(A)/u(Q)。扯远了，回到用蒙卡估计圆周率n的实验思路。假设样本空间是一个边长 为2面积为4的正方形，我们感兴趣的事件是正方形内的一个半径为1面 积为n的圆，所以向正方形内随机投掷一点，落在圆里面的概率为n/4。实验的思路如下：1. 1)生成随机数生成n个均匀落在正方形内的点；2.

10、 2)对落在正方形内的n个点，数一数正好落在圆里面的点的个数，假设为k (另外n-k个点就落在 圆外面的正方形区域内)。3. 3)k/n就可以大致认为是圆的面积与正方形的面积之比，另其等于n/4，就可以求出圆周率n的估 计值。又，Forcode提供了利用电子表格Excel求解以上问题的详细过程，有兴趣可以看看。另外，这里也有 一个详细的说明，用 Mathematica 实现。四、随机数很重要(以及下期预告) 在第一个例子中，我强调了骰子要是均匀的。那里骰子就是我们的随机数生成器，骰子的均匀程度就是 我们随机数的“随机”成分。在上面的10次试验中，投出了 3次点数“3”和3次点数“5”，然后点数

11、“1”、“2”、 “4”、“6”各出现一次因为只是重复了 10 次，这样我们对骰子的均匀程度不好评估，但如果重复无数 次假如是十万次，如果还出现类似比例的结果，那么就有理由怀疑这粒骰子的均匀程度了。如果骰 子灌了铅，比如投出点数“6”的可能性从六分之一上升到6分之三，那么根据这粒骰子来做模拟，我们就 有高估需求以及利润的危险。下次就讲讲随机数的生成原理，当然结合SAS来讲，或者也会提一下Excel。五、其他软件包在商业领域，基于Excel的插件比较兴盛，比如Crystal Ball应用就较为普遍，有试用版可以下载。其他 竞争产品有Risk、Risk Solver等。现在据说Risk Simul

12、ator也开始兴盛起来，大有后来居上之势。他 的开发者 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

13、上均匀分布的随机变量。一般地，我们把0,1上均匀分布 随机变量的抽样值称为随机数，其他分布随机变量的抽样都是借助于随机数来实现的。以下谈的都是所 谓“伪随机数(Pseudo Random Numbe。产生随机数，可以通过物理方法取得(很久很久以前，兰德 公司就曾以随机脉冲源做信息源，利用电子旋转轮来产生随机数表)，但当今最为普遍的乃是在计算机 上利用数学方法产生随机数。这种随机数根据特定的迭代公式计算出来，初值确定后，序列就可以预测 出来，所以不能算是真正的随机数(就成为“伪随机数”) 。不过，在应用中，只要产生的伪随机数列能 通过一系列统计检验，就可以把它们当成“真”随机数来用。现在大多软件

14、包内置的随机数产生程序，都是使用同余法(Congruential Random Numbers Generato。 “同余”是数论中的概念。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同余。上面的定义可以读写成，对

15、整数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)=moda*R(i),m，其中初值