统计计算课程设计题目基于蒙特卡洛方法求数值积分中文摘要蒙特卡洛方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。

传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。

1利用随机投点法，平均值法，重要性采样法，分层抽样法，控制变量法，对偶变量法，运用R软件求1d xe xθ-=⎰，42d x e xθ-=⎰和12d1xexxθ-=+⎰数值积分。

计算以上各种估计的方差，给出精度与样本量的关系，比较各种方法的效率，关键字蒙特卡洛随机投点法平均值法 R软件21 绪论蒙特卡洛的基本思想是，当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。

蒙特卡洛方法解题过程的三个主要步骤：（1）构造或描述概率过程对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。

即要将不具有随机性质的问题转化为随机性质的问题。

（2）实现从已知概率分布抽样构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡洛方法模拟实验的基本手段，这也是蒙特卡洛方法被称为随机抽样的原因。

最简单、最基本、最重要的一个概率分布是（0，1）上的均匀分布（或称矩形分布）。

随机数就是具有这种均匀分布的随机变量。

随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。

产生随机数的问题，就是从这个分布的抽样问题。

在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。

另一种方法是用数学递推公式产生。

这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。

不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。

由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。

由此可见，随机数是我们实现蒙特卡洛模拟的基本工具。

（3）建立各种估计量一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。

建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。

122方法介绍随机投点法随机投点法是进行n 次试验，当n 充分大的时候，以随机变量k/n 作为期望值E(X)的近似估计值，即n k p X E /)(=≈其中k 是n 次实验中成功的次数。

若一次投点试验的成功概率为p ，并以⎩⎨⎧=，表明试验失败表明试验成功0,1i X则一次试验成功的均值与方差为p p p X E i =-⋅+⋅=)1(01)()1()1(01)(222p p p p p X Var i -=--⋅+⋅=若进行n 次试验，其中k 次试验成功，则k 为具有参数为（n,p ）的二项分布，此时，随机变量k 的估计为n k p /=显然，随机变量的均值和方差满足()()p k E nn k E p E ==⎪⎭⎫ ⎝⎛=1()()np p p Var -=1dd3设计算的定积分为()dx x f I ba⎰=，其中a,b 为有限数，被积函数f(x)是连续随机变量ξ的概率密度函数，因此f(x)满足如下条件:()⎰+∞∞=-1)(dx x f x f 非负，且显然I 是一个概率积分，其积分值等于概率)(b a P <≤ξ。

下面按给定分布f(x)随机投点的办法，给出如下Monte Carlo 近似求积算法：(1)产生服从给定分布的随机变量值,i=1,2,…,N;(2)检查i x 是否落入积分区间。

如果条件b x a i <≤满足，则记录i x 落入积分区间一次。

假设在N 次实验以后，落入积分区间的总次数为n ，那么用 Nn I = 作为概率积分的近似值，即 Nn I ≈ 平均值法任取一组相互独立、同分布的随机变量{}i ξ，i ξ在[a,b]内服从分布率p(x),令,则也是一组相互独立、同分布的随机变量，而且(){}()()()I dx x f dx x p x p p E bab ai ===⎰⎰**ξ由强大数定理()1*1lim 1=⎥⎦⎤⎢⎣⎡=∑=∞→N i i N I p N p ξ4若记 ()I p N i Ni =∑=ξ1*1,则依概率1收敛到I ，平均值法就是用I 作为I 的近似值。

假如所需计算积分为()dx x f I ba ⎰=，其中被积函数在[a,b]内可积，任意选择一个有简单办法可以进行抽样的概率密度函数p(x)，使其满足条件：● ()()()b x a x f x p <≤≠≠时当0,0●1)(=⎰badx x p记 ()()()()()⎪⎩⎪⎨⎧=≠=0,00,*x p x p x p x f x p 则所求积分为()()dx x p x p I ba⎰=*因而Monte Carlo 近似求积算法为：(1) 产生服从分布率p(x)的随机数 N i x i ,...,2,1,=(2) 计算均值()∑==Ni i p N I 1*1ξ,即有I I ≈重要性采样法从数学角度上看定积分可以看成()()()dx x g x f x g I ⎰=10 其中g(x)是某个随机变量X 的密度函数，因此积分值I 可看成随机变量5Z=f(x)/g(x)的数学期望值()()()∑∑===≈=Ni i i Ni ix g x f N z NZ E I 1111为了减少模拟实验的方差应适当选取g(x),使Var （I ）尽可能小，如果被积函数f(x)>0,可取g(x)=cf(x),当c=1/I 时就有Var(I)=0.一般应选取和f(x)相似的密度函数g(x)，使f(x)/g(x)接近于常数，故而Var(I)接近于0，以达到降低模拟实验的方差，这种减少方差的模拟试验法为重要抽样法。

分层抽样法分层抽样法是利用贡献率大小来降低估计方差的方法。

它首先把样本空间D 分成一些不交的小区间mi U D 1==，然后在各个小区间内的抽样数由其贡献大小决定。

即，定义()dx x f p iD i ⎰=，则iD 内的抽样数in 应与ip 成正比。

考虑积分()dx x f ⎰=1θ将[0,1]分成m 个小区间：1...010=<<<=m a a a 则()()∑∑⎰⎰=====-mi imi a a Idx x f dx x f ii 1111θ记1--=i i i a a l 为第i 个小区间的长度，i=1,2,...,m ，在每个小区间上的积分值可用均值法估计出来，然后将其相加即可给出θ的一个估计。

具体步骤为：1）独立产生U(0，1)随机数{,...1,=j u ij 2）计算,...1,1=+=-j u l a x ij i i ij63）计算()∑==il j ijii x f n l I 1ˆ于是θ的估计为∑==mi iI 1ˆˆθ，其方差为()∑==mi i i i n l Var 122ˆσθ 其中，()2221⎥⎦⎤⎢⎣⎡-=⎰-i i a a i i l I dx l x f ii σ对偶变量法控制变量法利用数学上积分运算的线性特性:()()()[]()dx x g dx x g x f dx x f ⎰⎰⎰+-=选择函数g(x)时要考虑到：g(x)在整个积分区间都是容易精确算出, 并且在上式右边第一项的运算中对(f-g)积分的方差应当要比第二项对f 积分的方差小。

在应用这种方法时，在重要抽样法中所遇到的，当g(x)趋于零时，被积函数(f-g)趋于无穷大的困难就不再存在，因而计算出的结果稳定性比较好。

该方法也不需要从分布密度函数g(x)， 解析求出分布函数G(x)。

由此我们可以看出选择g(x)所受到的限制比重要抽样法要小些。

模拟过程：1)独立产生U(0，1)随机数{,...1,=j u ij2)计算)(1ˆ15∑==n i i x f n θ ()∑=-=ni i x f n 1511~θ73)计算()()()∑=-+=+=ni iix f x f n1552112~ˆˆθθθ控制变量法通常在蒙特卡洛计算中采用互相独立的随机点来进行计算。

对偶变量法中却使用相关联的点来进行计算。

它利用相关点间的关系可以是正关联的，也可以是负关联的这个特点。

两个函数值1f 和2f 之和的方差为{}{}{}{}(){}(){}221121212f E f f E f E f V f V f f V --++=+如果我们选择一些点，它们使1f 和2f 是负关联的。

这样就可以使上 式所示的方差减小。

当然这需要对具体的函数1f 和2f 有充分的了解 1）独立产生U(0，1)随机数{,...1,=j u ij2）计算()∑==ni i x f n 121ˆθ，找g(x),f(x)是相关的，且E[g(x)]=μ 3）计算()∑=-+=ni i x g n 126)(1ˆˆμλθθ 3 程序及实现结果10d x e x θ-=⎰的求解随机投点法先利用R 软件产生服从[0,1]上均匀分布的随机数 X,Y, ，计算)(x f y <的个数，即事件发生的频数，求出频率，即为积分的近似值。

R程序s1<-function(n){f<-function(x) exp(-x)a<-0b<-1x<-runif(n)y<-runif(n)m<-sum(y<f(x))j=m/nvar<-1/n*var(y<f(x))lis<-list(j,var)return(lis)}s1(10^4)s1(10^5)s1(10^6)89s1(10^7)对模拟次数n 调试了4次，分别为7654,,,n n n n ，得到精确值和模拟值。

表 随机投点法的模拟次数和模拟值精确值为 平均值法先用R 软件产生n 个服从[0,1]上均匀分布的随机数i x ，计算)(i x f ,再计算)(i x f 的平均值，即为定积分的近似值 R 程序p1<-function(n) {f<-function(x) exp(-x) a<-0 b<-1 x<-runif(n) y<-mean(f(x))n 410 510 610 710模拟值 方差10var<-1/n*var(f(x)) lis<-list(y,var) return(lis) }p1(10^4) p1(10^5) p1(10^6) p1(10^7)对模拟次数n 调试了4次，分别为7654,,,n n n n ，得到精确值和模拟值。