fortran产生随机数方法介绍

fortrai产生随机数方法介绍(附代码)注意：现在计算机产生的随机数都是伪随机数。1.0-1之间均匀分布的随机数random_number(x)产生一个0到1之间的随机数(x可以是向量)，但是每次总是那几个数。用了random_seed()后，系统根据日期和时间随机地提供种子，使得随机数更随机了。programrandomimplicitnonereal::xcallrandom_seed() !系统根据日期和时间随机地提供种子callrandom_number(x) !每次的随机数就都不一样了write(*,*)xstopendprogramrandom产生1-100的随机整数subroutinemy_random(abound,ubound)integer::abound,ubound,len,randomreal::tlen=ubound-aboundcallrandom_number(t)random二abound+floor(t*(len+1))returnendsubroutine2.任意区间均匀分布的随机数functionmy_random(lbound,ubound)implicitnonereal::lbound,uboundreal::lenreal::my_randomreal::tlen二ubound-lbound!计算范围大小callrandom_number(t)!t是0-1之间的随机数my_random=lbound+len*treturnend注意：在循环外callrandom_seed()产生一个随机数数组，只需加一个循环即可functionmy_random(lbound,ubound)implicitnonereal::lbound,uboundreal::lenintegersizereal::my_random(size) !size代表数组元素的个数real::tintegerilen二ubound-lbound!计算范围大小doi=l,10callrandom_number(t)!t是0-1之间的随机数my_random(i)=lbound+len*t!把t转换成lbound-ubound间的随机数enddoreturnend注意：同理在循环外callrandom_seed()标准正态分布随机数/高斯分布随机数徐士良的那本程序集里介绍了正态分布随机数产生的原理,不过他的方法只能产生较为简单的随机数，随机数的质量并不高，特别是随机数的数目较多时。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);那么XI,X2服从N(O,1)分布，且相互独立。等于说我们用两个独立的U(O,1)随机数得到了两个独立的N(0,1)随机数。值得说明的是，该方法产生的随机数质量很高。嘻嘻，本人亲自验证过。SAS和蒙特卡罗模拟_随机数一、 为什么选择SAS做蒙特卡罗模拟？为什么要用SAS做蒙卡？首先，对我来说，我只会用SAS，而且打算用SAS完成我所有的工作。当然，其他一些通用的理由有(Fan,etc.,2002):蒙卡是个计算密集的活，而SASBase、SASMacro、SAS/IML强大而灵活的编程能力能满足这一要求；2•做蒙卡时要用到大量的统计/数学技术，而SAS就内置了大量的统计/数学函数(在SAS/Stat和SAS/ETS)；用Fortran或C++当然也是非常好的主意，只是他们缺少内置的统计函数，代码就要冗长复杂很多。二、 什么是蒙卡？一个启发性例子好，开始，什么是蒙卡？了解它背景知识的最好办法当然是wiki-Monte_Carlo_method。蒙特卡罗是位于摩洛哥的一家赌场，二战时，美国LosAlamos国家实验室把它作为核裂变计算机模拟的代码名称。作为模拟方法，蒙卡以前就叫统计抽样(statisticalsampling)，我们感兴趣的结果因为输入变量的不确定而不可知，但如果能依概率产生输入变量的样本，我们就可以估计到结果变量的分布。跟蒙卡对应的，还有一种模拟技术叫系统模拟，包括排队、库存等模型，这些模型都跟随时间推移而出现的事件序列有关。下面举个蒙卡的例子，来自Evans,etc.(2001)的超级简化版。假设一家企业，利润是其需求量的函数，需求是随机变量。为了简化讨论，假定利润就是需求的两倍。这里输入变量就是不可控的需求，结果变量就是我们感兴趣的利润。假设需求以相同的概率取10、20、30、40、50、60这六种情况。在这样的简化下，我们就可以投一枚均匀的骰子来产生需求的样本，如果点数为1，对应得需求就是10，点数为2，需求就是20，以下类推。这样，我们的模拟过程就是：投骰子；根据骰子的点数确定需求量；根据需求量，求利润。

掷10次骰子，假设我们的模拟结果如下：重复次数骰子点数需求量利润1550100233060333060466012051102063306074404085501009220201055050这样通过模拟需求的样本，我们对结果利润的分布也就有所知晓，比如平均利润可以算出就是63。蒙卡一个重要的步骤就是生成随机数，这里我们是用投骰子来完成。三、又一个例子：利用蒙特卡罗模拟方法求圆周率n(pi)再举个很有名的例子，就是估计圆周率n的值，来自Ross(2006)。这个试验的思路正好可以帮我们温习一下几何概型的概念。我们知道概率的古典概型，就是把求概率的问题转化为计数：样本空间由n个样本点组成，事件A由k个样本点组成，则事件A的概率就是k/n。考虑到概率和面积在测度上具有某种共性，几何概型的基本想法是把事件跟几何区域相对应，用面积来计算概率，其要点是：认为样本空间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)生成随机数——生成n个均匀落在正方形内的点；2)对落在正方形内的n个点，数一数正好落在圆里面的点的个数，假设为k(另外n-k个点就落在圆外面的正方形区域内)。3)k/n就可以大致认为是圆的面积与正方形的面积之比，另其等于n/4,就可以求出圆周率n的估计值。又，Forcode提供了利用电子表格Excel求解以上问题的详细过程，有兴趣可以看看。另外，这里也有一个详细的说明，用Mathematica实现。四、 随机数很重要（以及下期预告）在第一个例子中，我强调了骰子要是均匀的。那里骰子就是我们的随机数生成器，骰子的均匀程度就是我们随机数的“随机”成分。在上面的10次试验中，投出了3次点数“3”和3次点数“5”，然后点数“1”、“2”、“4”、“6”各出现一次——因为只是重复了10次，这样我们对骰子的均匀程度不好评估，但如果重复无数次——假如是十万次，如果还出现类似比例的结果，那么就有理由怀疑这粒骰子的均匀程度了。如果骰子灌了铅，比如投出点数“6”的可能性从六分之一上升到6分之三，那么根据这粒骰子来做模拟，我们就有高估需求以及利润的危险。下次就讲讲随机数的生成原理，当然结合SAS来讲，或者也会提一下Excel。五、 其他软件包在商业领域，基于Excel的插件比较兴盛，比如CrystalBall应用就较为普遍，有试用版可以下载。其他竞争产品有@Risk、RiskSolver等。现在据说RiskSimulator也开始兴盛起来，大有后来居上之势。他的开发者JohnathanMun还在WileyFinance系列有一本书，ModelingRisk:ApplyingMonteCarloSimulation,RealOptionsAnalysis,Forecasting,andOptimizationTechniques。S语言（R或者S-Plus）由于其面对对象的特性，加之丰富的内置函数和诸多用户提供的库，使得R或者S-Plus也是蒙卡研究的得力工具——或者比SAS更有前景。这个我不熟，略之。RandomNumbers最简单、最基本、最重要的随机变量是在［0,1］上均匀分布的随机变量。一般地，我们把［0,1］上均匀分布随机变量的抽样值称为随机数，其他分布随机变量的抽样都是借助于随机数来实现的。以下谈的都是所谓“伪随机数”（PseudoRandomNumbers）。产生随机数，可以通过物理方法取得（很久很久以前，兰德公司就曾以随机脉冲源做信息源，利用电子旋转轮来产生随机数表），但当今最为普遍的乃是在计算机上利用数学方法产生随机数。这种随机数根据特定的迭代公式计算出来，初值确定后，序列就可以预测出来，所以不能算是真正的随机数（就成为“伪随机数”）。不过，在应用中，只要产生的伪随机数列能通过一系列统计检验，就可以把它们当成“真”随机数来用。现在大多软件包内置的随机数产生程序，都是使用同余法(CongruentialRandomNumbersGenerators)。“同余"是数论中的概念。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)=l。现在我们可以开讲“同余”了。设m是正整数，用m去除整数a、b,得到的余数相同，则称“a与b关于模m同余”。上面的定义可以读写成，对整数a、b和正整数m,若mod(a,m)=mod(b,m)，则称“a与b关于模m有相同的余数(同余)”，记做a三b(modm)(这就是同余式)。举个例子，mod(13,4)=l,mod(l,4)=l,则读成13和1关于模4同余，记做13三l(mod4)。当然，同余具有对称性，上式还可以写成1=13(mod4)。a三b(modm)的一个充要条件是a=b+mt,t是整数，比如13=1+3*4。a=b+mt可以写成(a-b)/m二t,即m能整除(a-b)。这些就够了。更多基础性的介绍，可以参考《同余(数据基础)》。乘同余法同余法是一大类方法的统称，包括加同余法、乘同余法等。因为这些方法中的迭代公式都可以写成上面我们见过的同余式形式，故统称同余法。常用的就是下面的乘同余法(MultiplicativeCongruentialGenerator.)。符号不好敲，做些约定，如R(i)就是R加一个下标i。乘同余法随机数生成器的同余形式如下：R(i+1)=a*R(i)(modm)。这个迭代式可以写成更直观的形式，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：ABCDE1iR(i)a*R(i)mod[a*R(i),m]mod[a*R(i),m]/m201=B2*4=mod(C2,13)=D2/1331=D2=B3*4=mod(C3,13)=D3/1342=D3=B4*4=mod(C4,13)=D4/13上面显示的是公式(Excel中，公式与计算结果的转换，用快捷键ctrl+~实现)，结果看起来是这样的，E列就是我们想要的在［0,1］上均匀分布的随机数：ABCDE1iR(i)a*R(i)mod[a*R(i),m]mod[a*R(i),m]/m201440.3076923083141630.23076923142312120.923076923上表中，a*R(l)=16,R(2)=3,m=13,—个同余式就是R(2)=a*R(1)(modm),3=16(mod13)，或者说，16-3能够被13整除——同余法就是这意思了。说“乘同余法”，要点在于a*R(i)是乘法形式。在SAS系统中，a=397,204,094,m=2,1-1=2147483647(—个素数)，随机种子R(0)可以取1到m-1之间的任何整数。伪随机数的检验现在内置于各大软件包的随机数生成器都经过了彻底的检验，我们当然可以安心地使用这个黑箱。或则我们也可以多了解些。一个好的“伪”随机数列，应该看起来就跟从［0,1］上均匀分布的随机数列中随机抽取出来的一样。两个比较直观的检验有：均匀性检验——所有的数是否均匀地分布在［0,1］区间；独立性(或不相关)检验——序列中的数不存在序列相关，每个数都跟其前后出现的数独立无关。一个例子是如0.1、0.2、0.3、0.4、0.5，这里每一个相继的数都正好比它前面的一个数大0.1，这样这个数列就似乎有了某种“格局”。具体的检验方法就略之了，更多请参见徐钟济(1985)。生成其他概率分布的随机数前面提到过，我们把［0,1］上均匀分布随机变量的抽样值称为随机数，其他分布随机变量的抽样都是借助于随机数来实现的。对其他连续型分布，(累积)分布函数为F(x)，r是在［0,1］上均匀分布的随机变量，另r=F(x)，解出的x就是该连续型分布的随机抽样。由于x可以写成r的反函数，该变换就称作“逆变换"(InverseTransformationmethod)。对离散型分布，思路类似，不过繁琐些，具体可以见徐钟济(1985)、Ross(2006)和Evans等(2001)。大多数相关软件包都会提供各种分布的随机数生成函数，知道有这回事就可以。最重要的，是我们陈述一类问题时，知道需要用哪种概率分布来描述Evans等(2001)：常用的连续分布均匀分布(UniformDistribution)描述在某最小值和最大值之间所有结果等可能的随机变量的特性。正态分布(NormalDistribution)是对称的，具有中位数等于均值的性质，就是我们熟悉的钟形形状。各种类型的误差常常是正态分布的。最后，作为中心极限定理的推论，大量具有任意分布的随机变量的平均数也呈正态分布。三角形分布(TriangularDistribution)由三个参数来定义，最大值、最小值和最可能值。临近最可能值的结果比那些位于端点的结果有更大的出现机会。指数分布(ExponentialDistribution)常用来构建在时间上随机重现的事件的模型，如顾客到达服务系统的时间间隔、机器元件失效前的工作时间等。它的主要性质是“无记忆性"(Memoryless)，即当前时间对未来结果没有影响。例如，不论机器原先已经运转了多长时间，它继续运转直至出现故障的时间间隔总有同样的分布。对数正态分布(LognormalDistribution):若随机变量x的自然对数是正态的，则x就服从对数正态分布。对数正态分布的常见例子是股票价格。常用的离散分布贝努利分布(BernoulliDistribution)描述的是只有两个以常数概率出现的可能结果的随机变量的特性；二项分布(BinomialDistribution)给出每次试验成功概率为p的n次独立重复贝努利实验的模型；泊松分布(PoissonDistribution)用于建立某种度量单位内发生次数的模型的一种离散分布，如某时段内某事件发生的次数。其他有用的分布如伽玛分布(GammaDistribution)、威布尔分布(WeibullDistribution)、贝塔分布(BetaDistribution)略之。下期预告上次落了些东西。说到蒙卡与C++,在数量金融领域，不能不提的一本书就是MarkJoshi的C++DesignPatternsandDerivativesPricing(CambridgeUni.Press,2004),面试中一定要说自己读过的，这样人家以为你至少会用C++做一个蒙特卡罗模拟。这个系列只讲SAS,下次先讲如何用SAS生成随机数，然后就是具体的项目。SAS随机数函数及CALL子程序《SAS和蒙特卡罗模拟(1)：开篇》一一简介，通过例子建立起蒙卡的直观概念；参考软件包及书目SAS和蒙特卡罗模拟(2):随机数基础》——随机数入门；参考书目

*********************************************************************sl**X**X**X**X**X**X**X*一、SAS随机数函数和CALL子程序SAS系统产生随机数，两种方式，利用SAS函数(Functions)或者CALL子程序(CALLRoutines),它们的语法格式是(具体的区别容后讨论)：方式代码说明函数var二name(seed,<arg>)var为存储随机数列的变量，name为特定的分布函数形式,seed为随机数种子,<arg>为特定分布要求的参数(可选)CALL子程序callname(seed,<arg>,var)同上，记得seed=0,±1,±2,,土(2**31-2)SAS可用的随机数函数列表如下(可以参见SASHelpandDocumentation-SASProducts-BaseSAS-SASLanguageDictionary-FunctionsandCALLRoutines-FunctionsandCALLRoutinesbyCategory)：分布SAS函数说明二项分布(Binomial)ranBin(seed,n,p)n为独立实验的次数,p为成功概率指数分布(Exponential)ranExp(seed)正态分布(Normal)ranNor(seed),ornormal(seed)ranNor和normal是同质的，但normal没有相对应的CALL子程序泊松分布(Poisson)ranPoi(seed,m)m为均值均匀分布(Uniform)ranUni(seed),oruniform(seed)ranUni和uniform是同质的，但uniform没有相对应的CALL子程序柯西分布(Cauchy)ranCau(seed)伽玛分布(Gamma)ranGam(seed,a)a>0为形状参数由分布律表格决定的离散分布(tabledprobabilitydistribution)ranTbl(seed,pl,p2,...pn)Ep=1三角分布ranTri(seed,h)h为斜边(最可能值)

(Triangular)上面的随机函数，除了normal和uniform,都可以由直接相应的CALL子程序调用。二、SAS随机数函数和CALL子程序：比较用SAS随机数函数同时创建的多个随机数变量其实都属于同一个随机数列。这话费解，一个例子先，创建两个随机数变量，各包含3个记录，其中x1的种子为123，x2的种子为456：dataranuni(drop=i);retainseed1123seed2456;doi=1to3;x1=ranuni(seed1);x2=ranuni(seed2);output;end;run;procprintdata=ranuni;run;结果为：Obsseedlseed2xlx2ll234560.750400.3209l2l234560.l78390.906033l234560.357l20.22lll好像没什么异样。我们把上面的xl增加为6个记录:dataranuni2(drop=i);retainseedll23;doi=lto6;xl=ranuni(seedl);output;end;run;procprintdata=ranuni2;run;结果如下，把上下用红色标出的数字对照看一看：ObsseedlxlObsseedlxl1230.750401230.320911230.320913 123 0.178394 123 0.906031230.357121230.357121230.22111什么意思？在第一段代码中，x2的种子根本不起作用，把x2的记录安插到x1里，看起来就是用种子123产生的随机数列加长了而已。x2的第一个值并不是由种子456产生的，而是产生第一个x1后的下一个值，xl、x2其实属于同一个随机数列，尽管x2的种子被指定为456，而xl的被指定为123。现在就可以重复上面的一句话：用SAS随机数函数同时创建的多个随机数变量其实都属于同一个随机数列。用CALL子程序就能够同时产生独立的随机数列。dataranuni3(drop=i);retainseed3l23seed4456;doi=lto3;callranuni(seed3,x3);callranuni(seed4,x4);output;end;run;procprintdata=ranuni3;run;结果如下：Obs seed3 seed4x3x4Obs seed3 seed4x3x4l6ll463328689l533263830888547364405l6l6ll463328689l533263830888547364405l67740697946869447500.750400.3209l0.l78390.342930.360450.3l988以上x3就是xl。xl和x3的初始种子都是123,但x3那个结果显示的种子是当前种子值。要在SAS随机数函数语句中显示随机种子的当前值，可以由以下代码给出：dataranuni4(drop=i);retainseedll23;doi=lto3;xl=ranuni(seedl);seed=xl*(2**3l-l);output;end;run;procprintdata=ranuni4;varseed1seedx1;run;结果如下，可以跟上面由CALL子程序得出的结果对照:Obsseed1seedx1112316114633280.7504021236891533260.320911233830888540.17839 参考资料 XitaoFan,etc..MonteCarloStudies:AGuideforQuantitativeResearchers.SASInstituteInc.,2002朱世武《SAS编程技术与金融数据处理》，北京：清华大学出版社，2003Matlab()随机数生成方法：第一种方法是用random语句，其一般形式为y=random(分布的英文名',Al,A2,A3,m,n),表示生成m行n列的mxn个参数为(Al,A2,A3)的该分布的随机数。例如：R=random('Normal',0,1,2,4):生成期望为0,标准差为1的(2行4列)2x4个正态随机数R=random('Poisson',1:6,1,6):依次生成参数为1到6的(1行6列)6个Poisson随机数第二种方法是针对特殊的分布的语句：一．几何分布随机数(下面的P，m都可以是矩阵)R=geornd(P)(生成参数为P的几何随机数)R=geornd(P,m)(生成参数为P的xm个几何随机数)1R=geornd(P,m,n)(生成参数为P的m行n列的mxn个几何随机数)例如R=geornd(1./2.人(1:6))(生成参数依次为112,1122,到1/2人6的6个几何随机数)R=geornd(0.01,[15])(生成参数为0.01的(1行5列)5个几何随机数).Beta分布随机数R=betarnd(A,B)(生成参数为A,B的Beta随机数)R=betarnd(A,B,m)(生成xm个数为A,B的Beta随机数)1R=betarnd(A,B,m,n)(生成m行n列的mxn个数为A,B的Beta随机数).正态随机数R=normrnd(MU,SIGMA) (生成均值为MU,标准差为SIGMA的正态随机数)R=normrnd(MU，SIGMA,m) (生成xm个正态随机数)

R=normrnd(MU,SIGMA,m,n)(生成m行n列的mxn个正态随机数)例如(1)R=normrnd(0,1,[15])生成5个正态(0,1)随机数(2)R=normrnd([123;456],0.1,2,3) 生成期望依次为[1,2,3;4,5,6],方差为0.1的2x3个正态随机数．四．二项随机数：类似地有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=binor