一个蒙特卡罗的小例子.txt

蒙特卡罗方法概述和它的主要引例：蒲丰投针问题

x4=34
x5=42
x6=18
x7=26
x8=2
如此可见x8=x0，可见这列数字得周期是8，若令x0=1，其余不变，则数列为
5，25，51，49，53，9，45，33，37，57，29，17，21，41，13，1
则周期为16
若按照上式产生得随机序列其值域为0～N-1，若要产生0～1之间的得随机数只要把每个随机数在除以N即可。
2、2随机性统计检验
一个好得随机数发生器或一个好得随机数生成程序要有两个条件
1、生成的随机数有足够长的周期。2、生成的随机数列应当有真正随机序列的统计性质。其周期的长短比较容易测试和判断，尔统计性质的优劣则不那么简单了。下面着重讨论统计性质的两种常用的检验方法，频数分布检验和行程频数检验。
先讨论频数分布检验。对于一个均匀分布的随机数发生器，设所产生的随机数序列的值域为[0,1]，则所产生的随机数字应与从0～1之间均匀的频数分布一致，为检验频数分布情况，可按画统计直方图的方法，将整个至于分成M个宽度相等的子区间，设xi是第i子区间内出现的随机数的个数，即第i个子区间的频数，则所有子区间中随机数的个数的平均频数
郁闷，没有公式编辑器，这后面的公式就不介绍了，大概了解下算了
就是每一个子区间中出现x个数字的几率要服从高斯分布规律
行程频数检验中的行程是指数字序列中单调上升或单调下降的连续数字串中数字的个数。由统计规律可知，对一个真正随机数序列其中行程长度k出现次数的期望值为E(K)
k=1,E(K)＝1/12*(5N+1)
k=2,E(K)=1/60*(11N-14)
……
k<N-1,E(K)=2/(k+3)!*[N(k2+3k+1)-(k3+3k2-k-4)]
行程长度k为任意值E(K)＝1/3(2N-1)
检验的方法就是输出给定随机数序列中出现的各种行程长度的次数，并与上述的期望值想比较，符合程度好的随机性好。
上述的两种检验方法都获得满意结果的随机数，一般 来说是可以信任的均匀分布的随机数序列。但仍然不是真正的随机数，因为并不能保证他们满足真正随机数的其他统计性质。
3
、用mc方法计算积分
已知函数f(x)在区间[a,b]上连续，A点x=a，B点x=b，其函数曲线f(x)－x，如图。
图片点击可在新窗口打开查看此主题相关图片如下：
下面介绍用mc方法来计算对f(x)的在[a,b]上的定积分
如图在图上做一个矩形，其宽度是b-a，高度是f(x)在[a,b]上的最大值f(m)，那么矩形的面积就是s1=(b-a)*f(m)，给出两个随机数xi和yi，且满足
a<=xi<=b,0<=yi<=f(m)，则所有的随机点都落在图中的矩形内。如果不等式yi<f(xi)成立，那么随即点（xi,yi）还落在矩形中的阴影区域内。
   总共设产生了N个随机点，落在阴影中的点数为M，那么当N足够大的时候对f(x)求区间[a,b]上的积分就换算成了求阴影的面积s2=M/N*(b-a)*f(m)。其实这和前面计算圆的面积没有本质的差别。对于比较复杂的问题，其基本处理方法还是相同的。如用mc方法计算一个具有空洞的球体质量和质心坐标等。
   随着电子计算机技术的发展，大容量、高速度新型计算机不断问世，蒙特卡罗方法越来越明显地发挥它的威力。上面已经说到，它首先是在研制核武器及核反应过程中发展起来的。如今，在各种新式武器的研制中，包括核武器、强激光武器、粒子束武器等都得到了广泛应用，而且在核电站设计、聚变堆研究中也发展了成套的蒙特卡罗程序软件。另外，在军事科学技术研究中也求助于蒙待卡罗方法：在实战前，对作战双方的军事实力、政治、经济、地理、气象等因素进行模拟，但这些因素可能随时发生变化，如果在计算机上进行“战斗”模拟，计算机就可以在很短时间内把一个很长的战斗过程模拟下来，告诉我们可能的结果。这样，军事指挥人员就可以进行成千上万次的战斗模拟，从中选择对自己一方量有利又最稳妥的作战方案，赢得战争的胜利。这相当于用计算机进行大规模的军事演习。现在世界上已有不少国家采用这种模拟方法，并在实际战斗中取得了成功。蒙特卡罗方法在国民经济、科学研究中亦已被广泛应用。在航天飞行器设计，稀薄气体分子动力学计算、高能物理、理论物理、凝聚态理论、材料科学、气象预报、可靠性工程及风险评价、地质探矿乃至医学科学中都大量使用了此种方法。我们可以毫不夸张地预言，随着计算机技术包括硬件、软件的新发展，蒙特卡罗方法应用的领域及深度还将有更大的发展。同时，人们也不会再感到它深奥莫测，而是得心应手的科学研究手段。
学习的话请进

3.另一个蒙特卡罗的小例子（关于MCMC）     -|killniu 发表于 2006-4-18 19:38:00  
 

    呵呵。就当是上课的笔记吧。抽样方法为Rejection Sampling。注意这个方法要用重尾分布来模拟轻尾分布。我们用Cauchy（密度函数为g）来模拟正态分布（密度函数为f）。设sup(f/g)=C。模拟的效率应该是1/C。也就是说n个g里面会有n/C个f。

    步骤：

0：计算出C。

1：产生Y~g.

2：Accept Y as X with prob=f(y)/{Cg(y)}。即产生i.i.d. 的 u~U(0,1)，Accept Y if u<=f(y)/{Cg(y)}。

3：repeat。

下面是一个例子，可以在R中运行。

x=NULL
y=rcauchy(100000,0,1)   #产生100000个Cauchy分布的rv
u=runif(100000,0,1)        #产生u
x=y[u<dnorm(y)/dcauchy(y)/1.52]   #以概率f(y)/{Cg(y)}让x=y；这儿C=1.52。

运行，length(x)=65784 几乎就是 100000/1.52。

qqnorm(x)后画的图几乎就是一条线。说明x确实是正态的。

中间可以range(x),range(y)看看。前面一个范围很小，轻尾；后面一个范围很大，重尾。

作业：用均匀分布U(0,1)产生标准正态的随机变量。

标准正态分布密度为f=1/sqrt(2*pi)*exp(-x^2/2))，f/g=f显然当x=0时达到最大为c=1/sqrt(2*pi)=0.3989423.
所以程序为：
c=0.3989423
y=runif(10000,0,1)
u=runif(10000,0,1)
x=y[u<dnorm(y)/dunif(y)/c] 
qqnorm(x)发现产生的效果不好。我觉得可能是因为c<1导致的。我们知道mc产生的效率是1/c。c<1的时候导致效率大于1了。有点问题。反过来，用正态产生均匀分布呢？

 

# 用U(0,1)均匀分布来产生指数分布(lambda=1).
# 显然 C=1.
u=runif(10000,0,1)
y=runif(10000,0,1)
x=y[u<dexp(y)/dunif(y)]
plot(density(x))   # 画出来的图线是效果不是很好，马马虎虎


# 用U(0,1)均匀分布来产生指数分布(lambda=1)
c=0.3989423
y=runif(10000,0,1)
u=runif(10000,0,1)
x=y[u<dnorm(y)/dunif(y)/c] 
plot(density(x))   # 画出来的图线是效果不是很好
qqnorm(x)  #另一个检验正态的图形方法，效果不好


4.
   用蒙特卡罗法求定积分和多重积分的问题
   举个例子吧
f(x,y,z)=x^2+y^2+z^2-3
其中 -sqrt(3)<=x<=sqrt(3)
-sqrt(3)<=y<=sqrt(3)
由于是半球，因此，0<=z<=sqrt(3)
编制M程序：
mtcl.m
function y=mtcl(fun,x,y,z,n)
if nargin<5 n=10000;end
n0=0;
xh=x(2)-x(1);
yh=y(2)-y(1);
zh=z(2)-z(1);
for i=1:n
a=xh*rand+x(1);
b=yh*rand+y(1);
c=zh*rand+z(1);
f=feval(fun,a,b,c);
if f<=0
n0=n0+1;
end
end
y=n0/n*xh*yh*zh;
球半球
f.m
function f=f(x,y,z);
f=x.^2+y.^2+z.^2-3;
积分计算
>>x=[-sqrt(3),sqrt(3)];
>>y=[-sqrt(3),sqrt(3)];
>>z=[-sqrt(3),sqrt(3)];
>>mtcl('f',x,y,z)
ans=
10.7124
当增加随机点的个数时
>>mtcl('f',x,y,z,100000)
ans=
10.9159