采用Fortran 90程序实现生成任意随机数列的方法
——谷辰先生
在科学计算/数值计算中的很多场合（比如蒙特卡洛模拟）都需要用到随机数。当然首先需注意，通过计算机产生的随机数都是伪随机数，并不是真正的随机数，因为真正意义上的随机数在某次产生过程中是按照物理实验过程（如掷骰子）中表现的分布概率随机产生的，其结果是不可预测的。而计算机中的随机函数是按照一定算法模拟产生的，其结果是确定的，是可见的。所以用计算机随机函数所产生的“随机数”并不随机，是伪随机数。但是，这个一般来说不影响我们在数值模拟中的使用。
Fortran自带一个random_number(x)函数可以产生一个0~1之间的随机数（x可以是向量），但是如果你要让每一次运行产生的随机数都不一样，则还应该调用random_seed()函数，然后系统就会根据日期和时间随机产生种子，从而得到“真正的”随机数。
下面将从最简单生成均匀分布随机数开始，逐步到复杂的瑞利分布，通过Fortran代码范例实战的方式讲解如何生成满足任一分布的随机数。
(1) 以下代码可以生成10组0~1之间的随机数（每组），而且每次运行结果都不一样：
!*****************************************************************
program tutorial01_random_number
implicit none
integer(4) :: k
real(8) :: x(3)
call random_seed()
do k = 1,10
call random_number(x)
write(*,*) x
write(*,*)
end do
stop
end program tutorial01_random_number
!*****************************************************************
运行结果如下：
请务必注意，call random_seed()应该在循环之外，如果放到循环内的话，即：!*****************************************************************
program tutorial01_random_number
implicit none
integer(4) :: k
real(8) :: x(3)
do k = 1,10
call random_seed()
call random_number(x)
write(*,*) x
write(*,*)
end do
stop
end program tutorial01_random_number
!*****************************************************************
在我windows下intel fortran编译器上运行的结果会是这样子的：
所以请务必注意。
matlab生成随机数
(2) 生成一组[xmin,xmax]范围内均匀分布的数组(对应子程序random_vec)
以下直接给出子程序
!*****************************************************************
subroutine random_vec(xmin,xmax,x,n)
implicit none
integer(4) :: n
real(8) :: xmin,xmax,x(n)
!----------------------------------------------------------------
!  **To generate a series of random numbers bounded in [xmin,xmax]**
!  Note: You have to use "call random_seed()" outside this subroutine,
!    otherwise you may get the SAME result every time you run it.
!
----------------------------------------------------------------
call random_number(x)
x = xmin + (xmax-xmin)*x
return
end subroutine
!*****************************************************************
正如其中所强调的，为了每次运行都出现不同的结果，应该在调用random_vec之前先call random_seed()一下。
(3) 生成一组标准正态分布随机数组(对应子程序GaussDistribution)
这里我们采用Box和Muller在1958年给出的由均匀分布的随机变量生成正态分布随机变量的算法，即：
X1 = sqrt(-2*log(U1)) * cos(2*PI*U2);
X2 = sqrt(-2*log(U1)) * sin(2*PI*U2);
其中U1, U2 是区间(0, 1) 上均匀分布的随机变量，且相互独立。
程序这里不再给出，请见文后下载地址。
以下是我们生成10000组正态分布随机数后，将结果在Matlab中用hist函数做成柱状图之后的结果：
(4) 生成一组任意正态分布随机数组(对应子程序GaussDistributionx)，均值为miu，标准差为sigma
在(3)的结果的基础上我们可以得到一组任意正态分布随机数组，均值为miu，标准差为sigma。下面给出的是一组均值miu为10，标准差sigma为20情况下的结果用Matlab画出来的分布图。
(5)对于如何生成满足其他任意分布的的随机数，可以参考文章《各种分布的随机数的生成》，这里不再赘述。其主要思想就是通过计算累积分布函数的反函数的方法来得到。例如，对于瑞利分布，我们也可以得到该分布函数的程序。瑞利分布的密度函数为：
利用Fortran通过模拟10000组满足瑞利分布的随机数组，将结果通过Matlab表示可以得到如下分布图（取miu=2）：

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。