6.2 随机数与蒙特卡罗模拟

6.2 随机数与蒙特卡罗模拟

6.2.1 随机数的生成

MATLAB统计工具箱中函数名以md三个字符结尾的函数用来生成常见分布的随机数,如表6-4所列。


以上函数直接或间接调用了rand函数或randn函数。下面以案例形式介绍normrnd和random函数的用法。

【例6-3】调用normrnd函数生成1000×3的正态分布随机数矩阵,其中均值μ=75,标准差σ=8,并作出各列的频数直方图。

以上命令生成的随机数矩阵比较长,此处略去,作出的频数直方图如图6-5所示。

【例6-4】调用normrnd函数生成1000×3的正态分布随机数矩阵,其中各列均值μ分别为0、15、40,标准差σ分别为1、2、3,并作出各列的颇数直方图。

以上命令生成的随机数矩阵略去,作出的频数直方图如图6-6所示。

【例6-5】调用random函数生成1000×1的二项分布随机数向量,然后作出频率直方图。其中二项分布的参数为n=10,p=0.3。

以上命令生成的随机数略去,作出的频率直方图如图6-7所示。

【例6-6】调用random函数生成10000×1的卡方分布随机数向量,然后作出频率直方图,并与自由度为1 0的卡方分布的密度函数曲线作比较,其中卡方分布的参数(自由度)为n=10。

以上命令生成的随机数略去,作出的频率直方图及自由度为10的卡方分布的密度函数曲线如图6-8所示。从图中可以看出,由random函数生成的卡方分布随机数的频率直方图与真正的卡方分布密度曲线附和得很好。

6.2.2 蒙特卡罗模拟

兹特卡罗方法是20世纪40年代美国的S.M.乌拉姆和J.冯●诺伊曼坠首先提出的。数学家冯●诺伊曼用驰名世界的赌城一一摩纳哥的Monte Carlo一一来命名这种方法,为它必上了一层神秘色彩。

蒙特卡罗方法的仰题过科可以分为三个上要步骤:构造或描述概率过程,实现从已知概率分布抽样和建立各种估计量。

1. 构造或描述概率过程

对于本身就具有随机性质的问题,如粒子输运问题,主要是正确描述和模拟这个概率过程。对于本来不是随机性质的确定性问题,如计算定积分,就必须事先构造一个人为的概率过程,它的某些参量正好是所要求问题的解,即要将不具有随机性质的问题转化为随机性质的问题。

2. 实现从已知概率分布抽样

构造了概率模型以后,由于各种概率模型都可以看作是由各种各样的概率分布构成的,因此产生己知概率分布的随机变量(或随机向量)就成为实现蒙特卡罗方法模拟实验的基本手段,这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布(或称矩形分布),随机数就是具有这种均匀分布的随机变量,随机数序列就是具有这种分布的一个简单子样,也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题,就是从这个分布抽样的问题。在计算机上可以用物理方法产生随机数,但价格昂贵,不能重复,使用不便。另一种方法是用数学递推公式产生,这样产生的序列与真正的随机数序列不同,所以称为伪随机数或伪随机数序列。不过经过多种统计检验表明,它与真正的随机数或随机数序列具有相近的性质,因此可把它作为真正的随机数来使用。己知分布随机抽样有很多方法,与从(0,1 )上均匀分布抽样不同,这些方法都是借助随机序列实现的,也就是说都是以产生随机数为前提的。由此可见,随机数是实现蒙特卡罗模拟的基本工具。

3. 建立各种估计量

一般来说,构造了概率模型并从中抽样后,需要确定一个随机变量作为所要求的问题的解,我们称它为无偏估计。建立各种估计量,相当于对模拟实验的结果进行考察和登记,从中得到问题的解。

【例6-7】使用蒙特卡罗方法计算圆周率,程序为simulation_pi.m。

当模拟投点次数为10万次时,计算的圆周率存在一定的误差,我们将测试投点次数从10万次逐渐提升至1000万次时,圆周率的波动越来越小并向真实的圆周率收敛。程序为simulation_piN.m,结果如图6-9所示。

使用上述程序的结果可以画出投点效果图,如图6-10所示。