蒲丰投针及蒙特卡罗模拟电子教案模拟
蒲丰投针――MonteCarlo算法
蒲丰投针 ―― Monte Carlo 算法背景：蒙特卡罗方法（Monte Carlo ），也称统计模拟方法，是在二次世界大战期间随着科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为基础的一类非常重要的数值计算方法。
蒙特卡罗方法在应用物理、原子能、固体物理、化学、生态学、社会学以及经济行为等领域中得到广泛利用。
蒙特卡罗方法的名字来源于世界著名的赌城 —— 摩纳哥的蒙特卡罗。
其历史起源可追溯到1777年法国科学家蒲丰提出的一种计算圆周的方法 —— 随机投针法，即著名的蒲丰投针问题。
问题：设在平面上有一组平行线，间距为d ，把一根长L 的针随机投上去，则这根针和平行线相交的概率是多少？(其中 L < d )分析：由于 L < d ，所以这根针至多只能与一条平行线相交。
设针的中点与最近的平行线之间的距离为 y ，针与平行线的夹角为 θ (0 ≤ θ ≤ π)。
相交情形 不相交情形易知针与平行线相交的充要条件是：sin 2Ly x θ≤=由于1[0,], [0, ]2y d θπ∈∈，且它们的取值均满足平均分布。
建立直角坐标系，则针与平行线的相交条件在坐标系下就是曲线所围成的曲边梯形区域（见右图）。
所以有几何概率可知针与平行线相交的概率是sin d 2212LL p d d πθθππ==⎰Monte Carlo 方法：随机产生满足平均分布的 y 和 θ，其中1[0,], [0, ]2y d θπ∈∈，判断 y 是否在曲边梯形内。
重复上述试验，并统计 y 在曲边梯形内的次数 m ，其与试验次数 n 的比值即为针与平行线相交的概率的近似值。
clear;n = 100000; L = 1; d = 2; m = 0;for k = 1 : ntheta = rand(1)*pi; y = rand(1)*d/2;if y < sin(theta)*L/2m = m + 1; end endfprintf('针与平行线相交的概率大约为 %f\n', m/n)计算π的近似值利用该方法可以计算 π 的近似值：sin d 22 22 1n LL m p d m d L d n πθθπππ⇒≈==≈⎰下面是一些通过蒲丰投针实验计算出来的 π 的近似值：蒲丰投针问题的重要性并非是为了求得比其它方法更精确的π值，而是在于它是第一个用几何形式表达概率问题的例子。
清华大学电力系统不确定性分析-----05蒙特卡罗模拟法_108801540
例题——利用变换方法生成标准 正态分布的随机变量
设对随机数 U1，U2 做下列变换：
N1 2 ln U 1 cos 2U 2 N 2 2 ln U 1 sin 2U 2
证明 N1， N2 为相互独立的服从 N(0， 1)分布的随机变量。
u1 1 2 n1 exp n12 n2 n1 2 u2 1 n 2 2 2 n1 2 n1 n2
舍选抽样方法
根据满足一定的检验条件进行舍选，以得到所需随机 变量的抽样值
设随机变量 η的密度函数 f x 可以表示 为 f x L
h x
g x, y dy
产生随机数（x，y）
其中 g x, y 是随机变量 X , Y 的联合密 度函数， h x 在 Y 的定义域上取常值，L 为一常量。则随机变量η的抽样过程为：
9
蒙特卡洛法的一般性原理
确定性问题
y
1
在正方形 0 x 1 ， 0 y 1 中随机地均匀地 投掷一点(ξ,η)，试求该点落在曲线 y f ( x) 下的
y=f(x)
概率 p。用 S 表示曲线 y f ( x) 下的区域；设ξ和η
S
相互独立，在[0,1]上满足均匀分布。于是有
dx 1
成立，其中α为置信度，1-α为置信水平。于是，可以根据问题的要求确定出 置信水平，按照正态分布表确定 cα，从而得到估计 GN 与真值 G 之间的误差
GN G
15
c N
蒙特卡洛模拟法计算指标的收敛过程示例
失稳概率 0.06 0.05 0.04 0.03 0.02 0.01 0 0
求出π值
Monte-Carlo模拟教程
举例
例1 在我方某前沿防守地域，敌人以一个炮排（含两门火炮） 为单位对我方进行干扰和破坏．为躲避我方打击，敌方对其阵地 进行了伪装并经常变换射击地点．
经过长期观察发现，我方指挥所对敌方目标的指示有50％是准 确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁 伤敌人一门火炮，有1/6的射击效果能全部毁伤敌人火炮．
蒙特卡罗方法的关键步骤在于随机数的产生， 计算机产生的随机数都不是真正的随机数(由算 法确定的缘故)，如果伪随机数能够通过一系列 统计检验，我们也可以将其当作真正的随机数 使用。
rand('seed',0.1);
rand(1) ％每次运ra行nd程('s序tat产e',s生um的(1值00*是clo相ck同)*r的and);
E = P(A0) = P(j=0)P(A0∣j=0) + P(j=1)P(A0∣j=1)
= 1 0 1 1 0.25 2 22
P(A1) = P(j=0)P(A1∣j=0) + P(j=1)P(A1∣j=1)
= 10 11 1 2 23 6
P(A2) = P(j=0)P(A2∣j=0) + P(j=1)P(A2∣j=1)
非常见分布的随机数的产生
• 逆变换方法
由定理 1 ，要产生来自 F(x) 的随机数，只要先 产生来自U (0,1) 随机数 u ，然后计算 F 1(u) 即 可。具体步骤如下：
(1) 生成 (0,1)上 均匀分布的随机数U。
(2) 计算 X F -1(U ) ,则 X 为来自 F(x) 分布的随机数.
蒙特卡罗方法的基本思想很早以前就被人们所发现和 利用。早在17世纪，人们就知道用事件发生的“频率” 来决定事件的“概率”。19世纪人们用蒲丰投针的方法 来计算圆周率π，上世纪40年代电子计算机的出现，特别 是近年来高速电子计算机的出现，使得用数学方法在计算 机上大量、快速地模拟这样的试验成为可能。
蒲丰投针实验模拟
一、蒲丰投针问题在平面上画有等距离的一些平行线，平行线间的距离为a(a>0) ，向平面上随机投一长为l(l<a)的针，针与平行线相交的概率p，结果发现π =2*l/(a*p).二、试验方法能够采纳MATLAB软件进行模拟实验，即用MATLAB编写程序来进行“蒲丰投针实验”。
1、基来源理因为针投到纸上的时候，有各样不一样方向和地点，但是，每一次投针时，其地点和方向都能够由两个量独一确立，那就是针的中点和偏离水平的角度。
以 x 表示针的中点到近来的一条平行线的距离，β表示针与平行线的交角。
明显有0<=x<=a/2 ，0<=β <=Pi 。
用边长为 a/2 及 Pi 的长方形表示样本空间。
为使针与平行线相交，一定x<=l*sinβ * ，知足这个关系的地区面积是从0 到Pi的l*sinβ对β的积分，可计算出这个概率值是(2l)/(Pi*a)。
只需随机生成n 对这样的x 和β，就能够模拟 n 次的投针实验，而后统计知足 x<=l*sin β * 的 x 的个数，就能够以为这是订交的次数。
而后利用公式求得π值。
2、MATLAB编程clear ('n')clear('a')clear('x')clear('f')clear ('y')clear ('m')disp(' 本程序用来进行投针实验的演示， a 代表两线间的宽度，针的长度 l=a/2 ，n 代表实验次数 '); a=input(' 请输入 a：');n=input(' 请输入 n：');x=unifrnd(0,a/2,[n,1]);f=unifrnd(0,pi,[n,1]);y=x<*a*sin(f);m=sum(y);PI=vpa(a*n/(a*m))三、实验数据 ( 部分程序截屏见后 )a n PI第一次310000第二次310000第三次3100000第四次3100000第五次31000000第六次31000000第七次3第八次3第九次3第十次3四、实验结论从上述数据剖析可知，跟着模拟次数的愈来愈多， PI 的值渐渐稳固在π值邻近，即愈来愈趋近于π，故蒲丰投针实验的确能够模拟出π的值。
数学教案：模拟方法——概率的应用
§3模拟方法-—概率的应用错误!教学分析这部分是新增加的内容．介绍几何概型主要是为了更广泛地满足随机模拟的需要,但是对几何概型的要求仅限于初步体会几何概型的意义，所以教科书中选的例题都是比较简单的．随机模拟部分是本节的重点内容．几何概型是另一类等可能概型，它与古典概型的区别在于试验的结果不是有限个，利用几何概型可以很容易举出概率为0的事件不是不可能事件的例子,概率为1的事件不是必然事件的例子．本节的教学需要一些实物模型为教具，教学中应当注意让学生实际动手操作,以使学生相信模拟结果的真实性，然后再通过计算机或计算器产生均匀随机数进行模拟试验，得到模拟的结果．在这个过程中，要让学生体会结果的随机性与规律性，体会随着试验次数的增加，结果的精度会越来越高．几何概型也是一种概率模型，它与古典概型的区别是试验的可能结果不是有限个．它的特点是在一个区域内均匀分布，所以随机事件的概率大小与随机事件所在区域的形状、位置无关，只与该区域的大小有关．如果随机事件所在区域是一个单点，由于单点的长度、面积、体积均为0，则它出现的概率为0，但它不是不可能事件；如果一个随机事件所在区域是全部区域扣除一个单点,则它出现的概率为1，但它不是必然事件．三维目标1．通过师生共同探究，体会数学知识的形成，正确理解几何概型的概念；掌握几何概型的概率公式：P（A）＝错误!，学会应用数学知识来解决问题，体会数学知识与现实世界的联系，培养逻辑推理能力．2．本节课的主要特点是随机试验多，学习时养成勤学严谨的学习习惯，会根据古典概型与几何概型的区别与联系来判别某种概型是古典概型还是几何概型，会进行简单的几何概率计算,培养学生从有限向无限探究的意识．重点难点教学重点:理解几何概型的定义、特点，会用公式计算几何概率．教学难点:等可能性的判断与几何概型和古典概型的区别．课时安排1课时错误!导入新课思路1。
复习古典概型的两个基本特点：(1)所有的基本事件只有有限个;(2)每个基本事件发生都是等可能的．那么对于有无限多个试验结果的情况相应的概率应如何求呢?为此我们学习几何概型．思路2。
Monte-Carlo(蒙特卡洛方法)解析
常用的线性同余生成器
Modulus m 2^31-1
=2147483647
2147483399 2147483563
Multiplier a 16807
在 n 次中出现的频率。假如我们取 fn ( A) 作为 p P(A) 的估计，即 pˆ fn ( A) 。
然后取 ˆ
2l afn ( A)
作为
的估计。根据大数定律，当 n 时，
pˆ
fn ( A) a.s.
p.
从而有ˆ 2l P 。这样可以用随机试验的方法求得 的估计。历史上 afn ( A)
(2) 计算 X F -1(U ) ,则 X 为来自 F(x) 分布的随机数.
例 1 ：设 X ~ U (a,b) ，则其分布函数为
0
F
(
x)
x b
a a
1,
xa a xb
xb
F -1( y) a (b a) y ， 0 y 1
生成 U (0,1) 随机数 U，则 a (b - a)U 是来自
算法实现
许多程序语言中都自带生成随机数的方法, 如 c 中的 random() 函数, Matlab中的rand()函数等。 但这些生成器生成的随机数效果很不一样, 比如 c 中的函数生成的随机数性质就比较差, 如果用 c , 最好自己再编一个程序。Matlab 中的 rand() 函数, 经过了很多优化。可以产生性质很好的随 机数, 可以直接利用。
U (a,b) 的随机数。
例 2：
设 X ~ exp( ) 服从指数分布，则 X 的分布函数为：
蒲丰氏投针问题的模拟过程
蒲丰氏投针问题的模拟过程，随机数发生器也是自编的，以供大家参考和提出建议。
谢谢。
（seed1和seed2最好选择3和5，为了使投针次数达到1000000，CVF进行如下设置Project->settings->link->output，将stack allocations reserve:设为1000000000）program getpiimplicit nonereal,parameter::a=5,L=4,pi=3.14159integer::n1,i,counter=0real,allocatable::R1(:),R2(:)real::theta,x,pi1write(*,*) 'input the size of the array:'read(*,*) n1allocate(R1(n1))allocate(R2(n1))call random(n1,R1,R2)do i=1,n1x=a*(2*R1(i)-1)theta=pi*R2(i)if(abs(x)<L*sin(theta)) counter=counter+1end dopi1=(1.0*n1)/(counter*1.0)*(2.0*L/a)write(*,*) 'this is PI:'write(*,*) piwrite(*,"('this is ratio of counter to total number',F10.6)")1.0&*counter/n1stopend programsubroutine random(n,R1,R2)implicit noneinteger n,seed1,seed2,i,little,mreal::R1(n),R2(n)integer::temp1(n),temp2(n)write(*,*) 'input seed1'read(*,*) seed1write(*,*) 'input seed2'read(*,*) seed2m=2**30m=m*2-1temp1(1)=397204094temp2(1)=temp1(1)R1(1)=(1.0*temp1(1))/(1.0*m)R2(1)=(1.0*temp2(1))/(1.0*m)little=0if(R1(1)<0.5) little=little+1do i=1,n-1temp1(i+1)=mod(seed1*temp1(i),m)R1(i+1)=(1.0*temp1(i+1))/(1.0*m)temp2(i+1)=mod(seed2*temp2(i),m)R2(i+1)=(1.0*temp2(i+1))/(1.0*m)if(R1(i+1)<0.5) little=little+1 .end do ;write(*,*) 'ratio of number which is little than 0.5'write(*,*) 1.0*little/nreturnend subroutine拟蒙特卡洛方法(Quasi-Monte Carlo)积分实例%使用Matlab提供的函数求积分,exp(-1/2*x^2)在(0,1)间积分 format long; syms xa = sym(1/2);f = exp(-a*x^2);ezplot(f)disp(int(f,-1,1));fprintf('integral result:%1.18f.\n',double(int(f,0,1)));%disp(double(int(f,0,1)));复制代码%使用拟蒙特卡洛方法积分%得到拟蒙特卡洛序列,即低偏差序列,halton法%如果有相关的工具箱的话,可以用Matlab里面的haltonset,faureset,sobolset 函数实现,x=halton(10000,2,5577);n=length(x);mju=0;for i=1:nmju=mju + exp(-0.5*x(i)^2);endmju=mju/n;fprintf('Quasi-Monte Carlo result:%1.18f.\n',mju);%disp(mju);%使用蒙特卡洛方法积分%得到Uniform序列,x=random('unif',0,1,10000,1);n=length(x);mju=0;for i=1:nmju=mju + exp(-0.5*x(i)^2);endmju=mju/n;fprintf('Monte Carlo result:%1.18f.\n',mju);%=============生成HALTON序列======================== function result = halton( m,base,seeder )%生成HALTON序列% Check inputsif nargin < 3seeder = 0;if nargin < 2error('MATLAB:Halton:NotEnoughInputs',...'Not enough input arguments. See Halton.'); endendres=0;n=length(base);for i=1:mfor j=1:nelement=0;temp=seeder+i;k=1;while temp>0element(k)=rem(temp,base(j));temp=fix(temp/base(j));k=k+1;endres(i,j)= 0;for k=1:length(element)res(i,j)=res(i,j)+element(k)/(base(j)^k); endendendresult=res;。
综合实验三 蒲丰投针问题实验
综合实验三 蒲丰投针问题实验一、实验目的1. 掌握几何概型、熟悉Monte Carlo 方法的基本思想；3．会用MATLAB 实现简单的计算机模拟二、实验内容在用传统方法难以解决的问题中，有很大一部分可以用概率模型进行描述．由于这类模型含有不确定的随机因素，分析起来通常比确定性的模型困难．有的模型难以作定量分析，得不到解析的结果，或者是虽有解析结果，但计算代价太大以至不能使用．在这种情况下，可以考虑采用Monte Carlo 方法。
下面通过例子简单介绍Monte Carlo 方法的基本思想．Monte Carlo 方法是计算机模拟的基础，它的名字来源于世界著名的赌城——摩纳哥的蒙特卡洛，其历史起源于1777年法国科学家蒲丰提出的一种计算圆周π的方法——随机投针法，即著名的蒲丰投针问题。
这一方法的步骤是：1) 取一张白纸，在上面画上许多条间距为d 的平行线，见图8.1(1)2) 取一根长度为()l l d <的针，随机地向画有平行直线的纸上掷n 次，观察针与直线相交的次数，记为m3）计算针与直线相交的概率．由分析知针与平行线相交的充要条件是 ϕs i n 21≤x 其中πϕ≤≤≤≤0,20d x 建立直角坐标系),(x ϕ,上述条件在坐标系下将是曲线所围成的曲边梯形区域，见图 8．l （2）．由几何概率知(*)22s i n 210d l d d G g p ππϕϕπ===⎰的面积的面积 4）经统计实验估计出概率,n m P ≈由(*)式即?2=⇒=ππd l n m Monte Carlo 方法的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型的参数或其他有关的特征量．然后通过模拟一统计试验，即多次随机抽样试验（确定m 和n ），统计出某事件发生的百分比．只要试验次数很大，该百分比便近似于事件发生的概率．这实际上就是概率的统计定义．利用建立的概率模型，求出要估计的参数．蒙特卡洛方法属于试验数学的一个分支．问题：(1) 经过n次试验后圆周率估计与的圆周 之间的差的绝对值的规律是？其中n分别取100，1000，2000，5000，10000，20000，50000(2) 参数l，d的不同选择，会对圆周率的估计有什么影响？可以选择d为l.5倍，2倍，3倍，4倍，5倍，8倍，10倍，20倍，50倍三、实验要求写出实验步骤、结果显示及分析四、实验分析以x 表示针的中点与最近一条平行线的距离，以j表示针与此线间的交角．显然0≤x≤a/20≤j≤p针与平行线相交的充要条件是x≤lsin(j)/2因(x,j)在图(4)中下面的矩形中等可能地取点，可见针与平行线相交的概率p 为图(4)正弦曲线线段与横轴围成的面积同图(4)中矩形面积的比．经计算得p= 另一方面得到如大量得投针实验，利用大数定理知：随着实验次数的增加，针与平行线相交的频率依概率收敛到概率p．那么在上式中以频率代替相应的概率p，则可以获得圆周率p的近似值．下面的程序是用matlab语言编写的计算机模拟投针以计算p 的近似值的程序．五、实验步骤1.编写MATLAB程序cleard=2l=0.5counter=0n=100x=unifrnd(0,d/2,1,n)fi=unifrnd(0,pi,1,n)for i=1:nif x(i)<1*sin(fi(i))/2counter=counter+1endendfren=counter/npihat=2*1/(d*fren)sqrt((pihat-pi)^2)结果显示:fren = 0.3300pihat =3.0303ans =0.1113以此类推：将n=1000，2000，5000，10000，20000，50000分别代入，可得：当n=1000时，fren =0.3240pihat =3.0864ans =0.0552当n=2000时，fren =0.3230pihat =3.0960ans =0.0456当n=5000时，fren =0.3204pihat =3.1211ans =0.0205当n=10000时，fren =0.3190pihat =3.1348ans =0.0068当n=20000时，fren =0.3172pihat =3.1521ans =0.0105当n=50000时，fren =0.3177pihat =3.1478ans =0.00622.改变d的取值，分别为1.5，2 ，3 ，4，5，8，10，20，50倍仍用1中的程序：cleard=3l=0.5counter=0n=100x=unifrnd(0,d/2,1,n)fi=unifrnd(0,pi,1,n)for i=1:nif x(i)<1*sin(fi(i))/2counter=counter+1endendfren=counter/npihat=2*1/(d*fren)sqrt((pihat-pi)^2)结果显示：d为1.5倍时fren =0.2300pihat =2.8986ans =0.2430d为2倍时fren =0.1700pihat =2.9412ans =0.2004d为3倍时fren =0.1100pihat =3.0303ans =0.1113d为4倍时fren =0.0800pihat =3.1250ans =0.0166d为5倍时fren =0.0600pihat =3.3333ans =0.1872d为8倍时fren =0.0400pihat =3.1250ans =0.0211d为10倍时fren =0.0300pihat =3.3333ans =0.1872d为20倍时fren =0.0100pihat =5ans =1.8539d为50倍时fren =0pihat =Infans =Inf六、结果分析1．经过n次试验后圆周率估计与的圆周π之间的差的绝对值的规律是：n的次数取值越多，圆周率估计与的圆周π之间的差的绝对值越小：圆周率越接近真值。
清华大学电力系统不确定性分析-----05蒙特卡罗模拟法_108801540
蒙特卡洛法的一般性原理
确定性问题
y
1
在正方形 0 x 1 ， 0 y 1 中随机地均匀地 投掷一点(ξ,η)，试求该点落在曲线 y f ( x) 下的
y=f(x)
概率 p。用 S 表示曲线 y f ( x) 下的区域；设ξ和η
S
相互独立，在[0,1]上满足均匀分布。于是有
u1 1 2 n2 exp n12 n2 n2 2 u2 1 n 2 1 2 n2 2 n1 n2
g n1 , n2 J
n12 exp 2 2 1
21
2 1 n2 2 exp 2
0.06 0.05
500000
0.04 0.03 0.02 0.01
1000000
1500000
2000000
> >
0 10000
30000
16
50000
70000
90000
蒙特卡洛法的基本特点
优点： ① 收敛速度与问题的维数无关 ② 受问题的几何条件影响不大 ③ 具有处理连续问题的能力 ④ 具有直接处理随机性问题的能力 缺点： a) 对于维数少的问题，一般是一维和二维问题，其 收敛速度慢 ，效率不及其它方法 b) 大的几何系统问题和小概率计算问题的求解偏差 c) 误差是概率误差而不是一般意义下的误差
0.00021
13787 33840 34366 20503 8164 2429 544
101
458166.0 699144.9 506523.4 236566.1 76105.1 18995.7 3744.1
蒲丰投针与蒙特卡洛(MonteCarlo)方法
蒲丰投针与蒙特卡洛（Monte —Carlo）方法1777年法国科学家蒲丰（Buffon ）提出并解决了如下的投针问题：桌面上画有一些平行线，它们之间的距离都是，一根长为a )(a l l ≤的针随机地投在桌面上。
问：此针与任一直线相交的概率是多少？设表示针的中点到最近的一条平行线的距离，Y 表示针与平行线的夹角（如图），如果X 2sin l Y X <， 或Y lX sin 2<时，针与一条直线相交。
由于向桌面投针是随机的，所以用来确定针在桌面上位置的是二维随机向量。
并且在),(Y X X ⎟⎠⎞⎜⎝⎛2,0a 上服从均匀分布，在Y ⎟⎠⎞⎜⎝⎛2,0π上服从均匀分布,与Y 相互独立。
由此可以写出的联合概率密度函数：X ),(Y X⎪⎩⎪⎨⎧<<<<=其它20,204),(ππy ax ay x f 于是，所求概率为：∫∫∫∫===⎭⎬⎫⎩⎨⎧<<20sin 20sin 224),(sin 2πππal dxdy adxdy y x f Y l X P y ly lx ①由于最后的结果与π有关，因此有些人想利用它来计算π的值。
其方法是向桌面投针次，若针与直线相交次，则针与直线相交的频率为n k n k ，以频率代替概率，则有al n k π2=，所以aknl2=π。
下表列举了这些试验的有关资料。
投针试验的历史资料（折算为1）a 试验者 年份 针长投针次数n 相交次数k π的试验值Wolf 1850 0.85000 2532 3.1596 Smith1855 0.63204 1219 3.1554 De.Morgan 1860 1600 383 3.137 Fox 1884 0.751030 489 3.1595 Lazzerini 1901 0.833408 1801 3.1415929 Reina1925 0.5425208593.1795这个思路已被人们发展成为统计学的一个分支—随机试验法或称为蒙特卡洛（Monte—Carlo ）方法，其中随机试验可借助计算机大量重复，以致结果更接近真值。
Monte Carlo(蒙特卡洛方法)解析
于是有： l p P( X sin ) 2 0
l sin 2
0
2 2l dxd a a
2l ap
若我们独立重复地作 n 次投针试验，记 n ( A) 为 A 发生的次数。 fn ( A) 为 A
U(0,1)随机数的生成
乘同余法：
xi 1 axi
mod m
ui 1 xi 1 / m 其中 xi , a, m 均为整数, x0 可以任意选取。
x0称为种子,a 是乘因子,m是模数
一个简单的例子
当 x0 1 时,得到序列: 1,6,3,7,9,10,5,8,4,2,1,6,3......
1 确定行为的模拟
例：曲线下的面积
本节以曲线下的面 积为例说明蒙特卡罗 模拟在确定行为建模 中的应用．
下面的算法给出了用蒙特卡罗方法求曲线下面积 的计算机模拟的计算格式．
在给定区间上曲线y=cosx下面积的真值是2．注意到即使对 于产生的相当多的点数，误差也是可观的．对单变量函数，一般 说来，蒙特卡罗方法无法与在数值分析中学到的积分方法相比， 没有误差界以及难以求出函数的上界M也是它的缺点．然而，蒙 特卡罗方法可以推广到多变量函数，在那里它变得更加实用．
ˆ f n ( A) 。 在 n 次中出现的频率。假如我们取 fn ( A) 作为 p P( A) 的估计，即 p
ˆ 然后取 2l a.s. ˆ fn ( A) 作为 的估计。根据大数定律，当 n 时， p p. af n ( A) 2l P 。这样可以用随机试验的方法求得 的估计。历史上 af n ( A)
蒙特卡罗随机模拟投点法
蒙特卡罗随机模拟投点法在数字积分中的应用数学与应用数学0901班：张瑞宸指导老师：任明慧摘要：本文首先介绍了蒙特卡罗方法的产生和发展，然后分析了蒙特卡罗方法计算数值积分的理论原理，最后给出了蒙特卡罗方法计算数值积分的MATLAB编程实现，全文主要是讨论了蒙特卡罗方法在定积分计算的应用。
而蒙特卡罗的优点：可以计算被积函数非常复杂的定积分、重积分，并且维数没有限制，这是别的数值积分方法还未达到的。
蒙特卡罗的缺点：收敛速度慢，误差一般较大，且是概率的误差，不是真正的误差。
关键词：蒙特卡罗方法，均值估计法，数值积分，Matlab编程Abstract：This paper first introduces the emergence and development of the Monte Carlo method, and then analyze the theoretical principles of Monte Carlo numerical integration method, Full-text mainly discussed the application of the Monte Carlo method in the definite integral. The advantages of Monte Carlo: can be calculated the integrable functions very complex definite integral, Multiple integrals, and dimension no limit, other numerical integration methods have not yet reached. Monte Carlo Disadvantages: slow convergence speed, error generally higher, and the probability of error, not a real error.Keywords: Monte Carlo method，Mean estimation method，numerical integral，Matlab programming0 引言历史上有记载的蒙特卡罗试验始于十八世纪末期（约1777年），当时布丰（Buffon）为了计算圆周率，设计了一个“投针试验”，后文会给出。
蒲丰氏问题
x ≤ l ⋅ sin θ
针在平行线间的位置
蒲丰氏问题的算法
如何产生任意的（x,θ）？x 在［0，a］上任意取值，表示 ， ］上任意取值，表示x 在［0，a］上是均匀分布的 ， ］上是均匀分布的， 其分布密度函数为： 类似地，θ的分布 的分布密度函数 的分布 为： 因此，产生任意的（x,θ） 产生任意的（ ） 产生任意的 的过程就变成了由f1(x)抽样 及 的过程就变成了由 抽样x及 抽样 抽样θ的过程了 由f2(θ)抽样 的过程了 抽样 的过程了。由此得 到：
当 −1 ≤ x ≤ 1 其它
令φ=2θ，则θ在［0,π］上均匀分布 在 ］上均匀分布，作变换 x = ρ cos θ y = ρ sin θ 其中0≤ρ≤1，0≤ρ≤π，则 x cos θ = x2 + y2 y sin θ = x2 + y2 (x,y) 表示上半个单位圆内的点。如果 (x,y) 在上半个单 位圆内均匀分布，则θ在［0,π］上均匀分布，由于 x2 − y2 2 2 cos ϕ = cos 2θ = cos θ − sin θ = 2 x + y2 2 xy sin ϕ = sin 2θ = 2 sin θ cos θ = 2 x + y2
2l 2l ≈ π= aP as N
请输入平行线相间的半距离（默认值 请输入平行线相间的半距离（默认值a=4.0） a=4 ） 请输入针的半长度（默认值 请输入针的半长度（默认值l=3.0） l=3 ） 请输入模拟试验次数（默认值 请输入模拟试验次数（默认值N=100000） N=100000 ）
s=0; for n=1:1:N; x2=mod(aa*x1,MM); x1=x2; randomx=x2*1.0/MM; x=a*randomx; x2=mod(aa*x1,MM); x1=x2; randomx=x2*1.0/MM; phi=pi*randomx; if (x<=l*sin(phi)) s=s+1; end end
第六章 第2节 投针试验
第六章第2节投针试验教学目标1经历试验、统计等活动过程，在活动中进一步发展学生合作交流的意识和能力。
2.能用试验的方法估计一些复杂的随机事件发生的概率。
教学重点:能用试验方法估计一些复杂的随机事件发生的概率.教学难点:借助大量的重复试验去感悟试验频率稳定于理论概率.课前准备：牙签（用来替代针）、直尺、计算器、画有平行线的纸。
多媒体辅助。
教学方法：本节课基本上是一节活动课，因而要体现学生的自主性，试验活动以及试验数据的汇总等都可以由学生自行组织完成。
让学生真正体验到了“当实验次数较大时，实验频率稳定于理论概率，并可据此估计某一事件发生的概率”。
力图让学生通过亲身的试验、统计过程获得用试验方法估计复杂事件发生的概率的体验。
教学过程：第一环节创设情境引入新课师：（出示多媒体）历史上法国有位数学家布丰做了一个有趣的试验————投针试验。
他在地上以相同的距离a画了一些平行线，站在一定距离远处，从一定高度随意投掷一定长度l ( l＜a )的针，该针可能与其中一条平行线相交，也可能与它们都不相交。
你能不能通过列表或画树状图求出该针与平行线相交的概率？生：不能，因为相交与不相交的发生可能性不同。
师：那么，怎样才能得到它们与平行线相交的概率呢？生：作试验得到概率. 因为我们前面学过，当试验次数较大时，试验频率稳定于理论概率”。
师：好的，那么我们现在就通过试验，解决这个问题。
板书课题（6.2.投针试验）设计意图：通过问题激发学生，学习的兴趣。
为试验教学打下基础。
第二环节合作交流探求新知试验活动1：投针试验试验道具：牙签（用来替代针）、直尺、计算器、画有平行线的纸.试验准备：在大的白纸上画出间距等于牙签长度的2倍的一组平行线.试验要求：（1）一定距离外（2）一定高度（3）随意抛掷（4）投针100次为一次试验试验步骤：（1）同桌两人组成一个合作小组，一人投牙签，一人记录相交次数.（2）计算本组试验数据中牙签与平行线相交的频率.每小组数据统计表（3）汇总各组试验数据，计算牙签与平行线相交的频率，从而估计相交的概率。
蒙特卡罗(Monte Carlo)模拟
如 对32位字长的计算机
real procedure random((xi)) integer array (li)n real array (xi)n l0 any integer that 1< l0 <231-1 for i=1 to n do
li =(231-1)除以16807 li-1的余数
1 b 1 n f ( x )dx f ( xi ) a ba n i 1
d b 1 1 n f ( x, y)dxdy f ( xi , yi ) c a (d c)(b a) n i 1
一般地，Af ( x )d ( measure of A)
作业
u
1 sin 2
中子屏蔽问题 Neutron Shielding problem 铅墙（长为5）
入 口
假设铅墙长为5，中子在铅中的平 均自由程为1，中子与铅原子碰撞 后各向同性散射。令碰撞8次后中 子能量耗尽，试求穿透铅墙的中 子的比例。 暂不考虑垂直纸面的运动，则中 子的水平位移是 1 cos cos ... cos 。
1
误差约
1 n
，它并不能和一些高级的数值积分算法比拟，
但对多维情况，MC方法却很有吸引力。
1 1 1
0 0 0
1 n f ( x, y, z )dxdydz f ( xi , yi , zi ) n i 1
我们可产生一系列随机数 1, 2 , 3 , 4 , 5 , 6 ,....... 可简单取3个随机数构成一个随机点，即 (1, 2 , 3 ), (4 , 5 , 6 ),....... 相应地，
算法 （常微分）
蛙 跳 如Mori et al. 1998, ApJ, 494, 430 如Liu, W. J. et al. 2007, Adv. Space Res., in press
学生的头脑是一个待被点燃的火种——《投针试验求π值的计算机模拟》案例的启示
学生的头脑是一个待被点燃的火种——《投针试验求π值的
计算机模拟》案例的启示
佚名
【期刊名称】《计算机教育》
【年(卷),期】2004(000)006
【摘要】《投针试验求π值的计算机模拟》一文是我校信息学院三年级本科生靳欣的论文，这是她在完成我开设的《电子商务》课程作业的基础上改写而成的。
【总页数】2页(P58,61)
【正文语种】中文
【相关文献】
1.点燃学生头脑中的火种——当代文学教学探索 [J], 吴梅菊
2.点燃学生心灵深处探究的火种——综合实践活动课案例与反思 [J], 李新
3.计算机模拟蒲丰投针试验近似计算圆周率研究 [J], 程登彪
4.用爱点燃火种——关爱留守儿童案例 [J], 周莹
5.用爱点燃火种——关爱留守儿童案例 [J], 周莹;
因版权原因，仅展示原文概要，查看原文内容请购买。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
概率模型的随机模拟与蒲丰投针实验第1章模拟1.1 模拟的概念每一个现实系统外部环境之间都存在着一定的数学的或者逻辑的关系，这些关系在系统内部的各个组成部分之间也存在。
对数学、逻辑关系并不复杂的模型，人们一般都可用解析论证和数值计算求解。
但是，许多现实系统的这种数学、逻辑模型十分复杂，例如大多数具有随机因素的复杂系统。
这些系统中的随机性因素很多，一些因素很难甚至不可以用准确的数学公式表述，从而无法对整个系统采用数学解析法求解。
这类实际问题往往可以用模拟的方法解决。
模拟主要针对随机系统进行。
当然，也可以用于确定性系统。
本文讨论的重点是其中的随机模拟。
采用模拟技术求解随机模型，往往需要处理大批量的数据。
因此，为了加速模拟过程，减少模拟误差，通常借助于计算机进行模拟，因此又称为计算机模拟。
计算机模拟就是在已经建立起的数学、逻辑模型的基础之上，通过计算机试验，对一个系统按照一定的决策原则或作业规则，由一个状态变换为另一个状态的行为进行描述和分析。
1.2 模拟的步骤整个模拟过程可以划分为一定的阶段，分步骤进行。
(1)明确问题，建立模型。
在进行模拟之前，首先必须正确地描述待研究的问题，明确规定模拟的目的和任务。
确定衡量系统性能或模拟输出结果的目标函数，然后根据系统的结构及作业规则，分析系统各状态变量之间的关系，以此为基础建立所研究的系统模型。
为了能够正确反映实际问题的本质，可先以影响系统状态发生变化的主要因素建立较为简单的模型，以后再逐步补充和完善。
(2)收集和整理数据资料。
模拟技术的正确运用，往往要大量的输入数据。
在随机模拟中，随机数据仅靠一些观察值是不够的。
应当对具体收集到的随机性数据资料进行认真分析。
确定系统中随机性因素的概率分布特性，以此为依据产生模拟过程所必需的抽样数据。
(3)编制程序，模拟运行。
选择适当的计算机语言。
按照系统的数学、逻辑模型编写计算机程序。
然后可以进行调试性模拟，分析模拟结果是否能够正确地反映现实系统的本质。
并相应修改模型和程序，反复调试之后确定模型。
根据系统的特点及问题的性质设定初始状态，在确定了模拟运行的时间、随机样本量的多少及独立模拟运行次数等之后，就可以开始进行大量的计算机模拟运行，从而获得模拟输出资料。
（4）分析模拟输出的结果。
对模拟的输出结果进行分析，一般包括一下几个方面的内容：模拟结果额统计特性，包括样本均值、方差及置信区间等；灵敏性分析；如果某个输入参数作微小的变化却导致输出结果巨大的改变，则应花费更多的时间和精力对它们求得更精确的估计。
根据前期确定的目标函数，在多种实现方案中选择最优的方案。
1.3 模拟的作用模拟技术广泛的应用于工业企业管理、物资管理、军事后勤系统、排队系统、交通运输系统等众多领域中。
模拟主要有以下作用：(1)可以处理难以用解析法解决的问题。
对于那些很难用解析方法加以处理的问题，模拟是一种有效的技术。
计算机模拟是一种数值试验手段，它可以模仿真实系统随机运行的特性。
将由此特性产生的抽样数据输入到事先编制好的程序中去，让计算机计算出相应的结果，从而得出实际问题的数值解。
(2)对假设和结论进行检验。
在建立数学模型的过程中，常常要对客观问题作一些假设。
这些假设并不一定符合实际，需要通过鉴定。
另外，理论研究的一些结论也往往需要加以检验，如果按照系统实际运行的进行试验，一般需要花费大量的人力、物力、财力和时间，有时甚至是不可能的。
使用计算机模拟技术代替实际系统的反复试验，可以大幅度地减少资源的耗费，快速地积累数据。
同时，经过认真分析计算机模拟出来的输出结果，可对系统的性能作进一步的分析与评价，从而能够确定假设的正确性，改进理论研究的一些成果，改善系统运行的方案等。
(3)对各种方案的评估，选择最优方案。
在模型建立后，人们通过收集、整理和分析有关数据资料，往往可能对该系统制订出多种不同的方案。
对这些不同的方案，可用计算机进行多次模拟，按照既定的目标函数对不同方案进行比较，从中选择最优的实现方案。
1.4 模拟的分类模拟分为动态模拟与静态模拟。
动态模拟又可以分为连续系统模拟和离散事件系统模拟，前者研究系统的状态随时间连续变化的情形，其模型一般是微分方程；后者讨论的是系统状态只在一些离散时间点上，由于事件驱动而变化，其模型一般只能用流程图或网络来表示。
而以蒙特卡罗（Monte Carlo）思想建立起来的蒙特卡罗模拟是典型的静态模拟。
第2章蒙特卡罗模拟2.1 蒙特卡罗模拟的起源蒙特卡罗是摩纳哥的一个城市，以赌博闻名于世界。
蒙特卡罗法借用这一城市的名称是为了象征性地表明该方法的概率统计的特点。
蒙特卡罗模拟是一种随机模拟方法。
以概率和统计理论方法为基础的一种计算方法。
将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的近似解。
由S.M.乌拉姆和J.冯·诺伊曼在20世纪40年代为研制核武器而首先提出。
在这之前，蒙特卡罗方法就已经存在。
1777年，法国科学家蒲丰提出用投针实验的方法求圆周率π，这就是蒙特卡罗模拟最早的雏形。
2.2 蒲丰抛针实验抛针问题主要内容是：在平面上等距离的画出一些平行线，向这个平面抛出某一特定长度的针，试求针与任一平行线相交的概率。
显然只有两种可能性：要么针与其中至少一条平行线相交；要么针与所有的平行线都不相交。
但是，古典概率无法解决这一个表面上很简单的问题。
不过，可以利用几何概率模型算出要求的概率值，正是随处可见的π。
图1 蒲丰投针试验方法假设平行线距离为a ，针长度为l ，针落下后中点为M ，x 则表示中点M 到最近的一条平行线的距离，ϕ表示针与平行线的交角。
那么，不难得到：020a x ϕπ⎧≤≤⎪⎨⎪≤≤⎩ 从几何的角度看，上述方程组代表了一个矩形区域Ω。
图2 蒲丰投针试验原理 这个矩形的面积为：()2a S πΩ= 由于针与平行线相交，则M 必定与最近的一条平行线相交，那么相交的充分必要条件是：0sin 20l x A ϕϕπ⎧≤≤⎪=⎨⎪≤≤⎩A 在图2中显示为Ω中的一个小区域，它的面积为：()01sin 2S A l d l πϕϕ==⎰ 由概率的几何意义，针与任一平行线相交的概率是：()()212S A l l P S a a ππ===Ω 从上面概率的表达式可以看出，只要l 跟a 的比值不变，概率P 的值也不变。
由于P 表达式中含有π，可以利用抛针试验得到的结果去估算π的值。
原理是尽量加大试验次数N ，记录针与平行线相交的次数n ，则P 的近似值是n N。
代入的P 表达式： 2lN anπ≈ 表1 抽针实验数据表蒲丰抛针试验揭示了蒙特卡罗方法的基本思想，即在大数定理的保证下，利用事件发生的“频率”作为事件的“概率”的近似值。
只要设计一个随机试验，使一个事件的概率与某未知数有关，然后通过重复试验，以频率近似表示概率，即可求得该未知数的近似解。
2.3 蒙特卡罗模拟的基本思想及其发展当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以 通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。
这就是蒙特卡罗方法的基本思想。
蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。
它是以一个概率模型为基础，按照这个模型所描绘的过程，通过计算机直接进行抽样试验得出结果，作为问题的近似解。
从理论上来说，蒙特卡罗方法需要大量的实验。
实验次数越多，所得到的结果才越精确。
计算机技术的发展，使得蒙特卡罗方法在最近10年得到快速的普及。
现代的蒙特卡罗方法，已经不必亲自动手做实验，而是借助计算机的高速运转能力，使得原本费时费力的实验过程，变成了快速和轻而易举的事情。
蒙特卡罗方法有很强的适应性，问题的几何形状的复杂性对它的影响不大。
蒙特卡罗方法的收敛性是指概率意义下的收敛，因此问题维数的增加不会影响它的收敛速度，而且存贮单元也很省，这些是用该方法处理大型复杂问题时的优势。
它不仅较好地解决了多重积分计算、微分方程求解、积分方程求解、特征值计算和非线性方程组求解等高难度和复杂的数学计算问题，而且在统计物理、核物理、真空技术、系统科学、信息科学、公用事业、地质、医学，可靠性及计算机科学等广泛的领域都得到成功的应用。
借助计算机技术，蒙特卡罗方法实现了两大优点：一是简单。
省却了繁复的数学报导和演算过程，使得一般人也能够理解和掌握；二是快速。
简单和快速，是蒙特卡罗方法在现代项目管理中获得应用的技术基础。
2.4蒙特卡罗模拟一般步骤(1)构造或描述概率过程：对于本身就具有随机性质的问题，主要是正确描述和模拟这个概率过程。
对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解，即要将不具有随机性质的问题转化为随机性质的问题。
(2)实现从已知概率分布抽样：构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量，就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。
最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）。
随机数就是具有这种均匀分布的随机变量。
产生随机数的问题，就是从这个分布的抽样问题。
在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。
另一种方法是用数学递推公式产生。
这样产生的序列，与真正的随机数序列不同，所以称为伪随机数。
不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。
由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。
由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。
(3)建立各种估计量：一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，并称它为无偏估计。
建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。