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

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

1。
抽样方法为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)  #另一个检验正态的图形方法，效果不好
蒙特卡罗直接抽样程序：
function mcs_1
syms i x s mm k
mm=1000;
nn=mm^2;
s=0;
i=1;
u=[450 162 135];%矩阵U第一列代表抗力R,第二列代表恒载G,第三列代表活载L
p=[90.4 12.96 21.6];%矩阵p的三列分别为R,G,L的标准差
m=pi/(p(1,3)*sqrt(6));%m为极值I型当量正态化过程中的一个常数
n=u(1,1)-0.5772/m;%n为极值I型当量正态化过程中的另一个常数
uu=log(u(1,1)/sqrt(1+(p(1,1)/u(1,1))^2));
pp=sqrt(log(1+(p(1,1)/u(1,1))^2));
a=lognrnd(u(1,1),p(1,1),mm,mm);%对数正态分布随机抽样
b=normrnd(u(1,2),p(1,2),mm,mm);%正态分布随机抽样
c=norminv(rand(mm))*pp+uu*ones(mm);%极值I型分布随机抽样
for i=1:nn
    if a(i)-b(i)-c(i)<0
        s=s+1;
    end
        i=i+1;
end
bb=norminv(1-s/nn)%可靠度指标
disp(s)
2.
蒙特卡罗方法概述

序言：为了丰富知识和学习，以后将陆续介绍一些新兴学科和交叉学科，让大家对近代科学有大致的了解，起到科普知识学习的目的!由于本人能力有限，介绍的偏重数学，物理，计算机再加少部分化学！希望大家喜欢！这些大部分摘自网上，恕不能一一写出摘自哪里，如有侵权请告诉我！

Monte Carlo method
   计算方法的名称大多是与数学术语相关的，唯有蒙特卡罗方法名称特殊。大家知道，蒙特卡罗是地中海北岸的一个风景优美的城市，它是一座世界上有名的赌城。饶有兴趣的是科学家们竞用它命名了一种数学计算方法。那么，赌博与随机模拟(亦称统计试验)有什么关系呢?其实，赌博本身就可以看作是一种最简单的统计试验。例如掷骰子，一个均匀正方形的六面体，它在每一面刻有1到6的六个数码，如果要想知道掷骰子时某一面向上机会的多少，既可以用数学理论来证明，也可以用多次投掷方法来验证。只要投掷的数目足够多，就一定能证明每一面朝上的机会是1／6。类似扔硬币正反面朝上机会显然各都是1／2。有人做了30000次投掷试验，统计的结果正面向上的次数是14994次。也就是说，十分接近1／2了。这二个例子告诉我们，少量的试验结果可随是无规律的．数学上称随机性，而大量的试验就能检验或发现某些规律，也就是说有必然性，数学上叫统计性。对赌博来说，既有随机性也有统计性。
    当然，实际问题要比掷骰子和硬币复杂得多。但是，统计数学的现象和规律是类似的，由于这种方法受到赌博的启示，所以人们就形象地用赌博城蒙特卡罗来命名这种计算方法。独具风格的计算方法。
    MC（Monte Carlo）是上世纪四十年代，Los Alamos实验室（LANL）曼哈顿计划中由Stanislaw Ulam和Von Neumann两位数学家倡导和创建的一种独具特色的数值模拟方法。为体现方法的内涵，他们用世界著名赌城的名字将其命名为Monte Carlo方法。
    MC方法也被称为统计实验方法，Monte Carlo借助计算机，发展功能强大和计算高效的程序是应用的前提。MCNP即是蒙特卡罗大型科学计算程序，用于处理连续能量、时间相关、三维几何、中子－光子－电子辐射输运问题。
蒙特卡罗方法不仅名称独特，它的计算模拟过程也独具风格，因而它常给人们带上一层神秘的色彩。但是，分析这种方法的基本特点就会发现，它的规律并不难掌握。用蒙特卡罗求解问题首先要建立一个概率模型或随机过程，使它的参数等于问题的解，然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征，最后再输出所求解的近似值。解的精度可以用估计值的标准误差来表示。根据概率论的定理，只要重复抽样随次数N趋向无穷大时，则其随机变量的算术平均将精确等于它的数学期望，即所求问题的解。当然，N实际上不可能取无穷大，但只要取到与问题相应的足够大便可十分逼近于所求的真实解。反之，N小了就得不到正确的结果。这就决定了只有在高速电子计算机上才能有效地应用这种方法。蒙特卡罗方法既可求解随机性问题，又可以求解确定性问题。所谓随机问题，是指问题的过程或参量受到随机性的影响，当然最后结果遵从统计规律， 是确定的。例如，在核反应或核武器爆炸中，中子或伽玛射线在介质中的输运问题是随机性问题。因为粒子和介质元素发生什么相互作用，碰撞后如何运动都有随机性。我们可以用蒙待卡罗方法直接模拟这种输运过程，只要足够多的抽样，便可获得精确解。1946年，美国科学家首先在电子计算机上对中子连锁反应进行了模拟，并且把第一个实验程序命名为蒙特卡罗程序。
    因此，蒙特卡罗方法首先在核科学研究中得到应用。其他如运筹学中的库存问题，随机服务系统中的排队问题等都属随机性问题；另一类是确定性问题，如计算定积分、解偏微分方程、线性代数方程等，它们可以用差分法、数值积分法等方法求解，但同样也可以用蒙特卡罗方法求解。在某些情况下，蒙特卡罗方法还更有优越性，如统计物理中常遇到的计算某物理系统的平均量，象内能、自由能、熵等均是高维积分，用通常数值积分法十分困难，而用Metroplis提出的特殊蒙特卡罗方法则能精确高效地求解出结果，从而使其在统计物理、凝聚态物理研究中成为一种标准算法。
MC方法主要应用于粒子输运问题上，它的适用范围分为实验核物理、反应堆核物理以及粒子辐射效应、抗辐射加固等方面。
    MC方法求解粒子输运问题时，应该遵循四个主要的步骤：弄清楚粒子输运的全部物理过程；确定适用的MC技巧；确定描述粒子运动的状态参数和状态序列；确定粒子输运过程中有关分布的抽样方法。
    MC模拟粒子输运过程包括以下几个物理过程，它对于我们理解MCNP程序的运算是很有用的。模拟粒子输运过程大致分为：源分布抽样，确定输运粒子的初始状态；空间、能量和运动的随机游动过程，确定粒子输运过程中的状态；记录贡献与分析结果。
    MC统计模拟方法可以采用一些基本技巧，包括直接模拟方法、简单加权方法、统计估计方法、指数变换方法等，通过程序设计者对模拟程序的优化，可以极大地提高MCNP的能力。
蒙特卡罗方法有其独特的优点，第一，与所求解问题的几何维数及问题条件关系不大，几何越复杂，它相对优点越明显。例如，在粒子输运问题中，用差分法解二维问题比一线问题几乎要多4倍以上的时间，而采用蒙特卡罗方法则几乎不受影响。第二，适应性强。例如，积分域形状特殊时，用一般积分法求解困难大，而蒙特卡罗法则不受影响，对问题也不一定要进行离散化，可连续处理，如粒子能量可以连续跟踪模拟等。第三，程序结构简单，所需计算机存贮单元比其他数值方法少，容易建立通用性很强的应用软件。应当指出，蒙特卡罗计算结果的收敛性与抽样数N的平方根成反比，这就决定了要得到高精度解就必须用高速电子计算机实现大N抽样。但是它的误差只依赖于标准方差σ及 。因此，要节省时间就既要研究各种模拟方法以减少方差、提高效率，所以降低方差的各种抽样技巧便成为研究该方法的重要内容。
蒙特卡罗方法有很强的适应性，问题的几何形状的复杂性对它的影响不大。该方法的收敛性是指概率意义下的收敛，因此问题维数的增加不会影响它的收敛速度，而且存贮单元也很省，这些是用该方法处理大型复杂问题时的优势。因此，随着电子计算机的发展和科学技术问题的日趋复杂，蒙特卡罗方法的应用也越来越广泛。它不仅较好地解决了多重积分计算、微分方程求解、积分方程求解、特征值计算和非线性方程组求解等高难度和复杂的数学计算问题，而且在统计物理、核物理、真空技术、系统科学 、信息科学 、公用事业、地质、医学，可靠性及计算机科学等广泛的领域都得到成功的应用。
为对mc方法有个初步的认识先看个用mc方法求圆周率的例子。
   设园的半径是r，圆心位于xoy平面的（r,r）处，且内切于边长为2r的正方形中，那么显然正方形的面积是S1＝4r×r，用mc方法计算圆的面积的基本思想是随机的正方形范围内部画点，若共画了N个点，而落在圆内部的有M个点，当N足够大的时候，圆的面积可为s=M/N*S1,基本方法是首先产生两个随机数x和y其值域为[0，2r]，然后进行判断（x,y）是否落在圆内
（x－r）^2＋（y-r）^2<r^2
记录下总的点数N和落在圆中的点数M，则圆的面积为S＝4r^2M/N，从而得pi＝4×M/N
这里要注意两点：第一，两个随机数要是在矩形范围内均匀分布得。大部分计算机语言提供了随机数发生器，可产生0－1之间得均匀随机数，经过适当得变换不难变成任意值域得均匀随机数。第二，随机数得个头N必须要足够大。已确保一定得精度。负责由于统计涨落过大而使得误差很大，这是统计方法得自身得规律。
mc方法不仅能用于物理方程得数值计算，也可以用于物理过程得数字模拟。氢原子的电子云模拟就是一个简单的例证，由量子力学可以得知氢原子的s态波函数ψ＝ψ（r)只是半径的函数，于θ，φ无关，而氢原子中的电子沿半径r的分布密度即电子在半径r处单位厚度球壳内出现的几率
D=4×pi×r^2×ψ2
习惯上把这种分布称作为电子云。
氢的基态1s态（n=1,l=0,m=0），有
D=1/(a1*a1*a1)4×r^2×e-2r/a1
Dmax=1.1
r0=0.25nm
其中a1＝5。29×10－2 nm，是D的最大值Dmax处的r值，其值是与波尔半径相同，r0是Dmax处的r值，其值与波尔半径相同。r0是D收敛处的r值，即D的收敛点。
氢的2s态（n=2,l=0,m=0）有
D=1/(8×a1*a1*a1)×（2－r/a1)2×e-r/a1
Dmax=0.14
r0=1.0nm
氢原子的3s态（n=3,l=0,m=0）有
D=1/(81×81×3×a1*a1*a1)[27-18*r/a1+2*(r/a)2]×e-2r/3a1
Dmax=0.2
r0=2.0nm
氢的电子云的模拟是依据上述的分布函数用绘图点的密度来描绘电子的概率分布函数，在氢原子例子中设随机数发生器可产生0－1的随机数，记做rand(k)（k=1,2,3……），首先利用一个随机数rand(1)来产生一个随机的电子半径轨道
r=r0*rand(1)
显然0<r<r0，由r计算出D（r)。再产生一个随机的概率判据D0＝Dmax* rand(2)进行判断，如果D(r）<D0，则从新开始，否则再产生一个 随机的角度值
θ＝2*pi*rand(3)。最后计算要描绘的点的坐标x=srcosθ y=srsinθ。其中s是控制图形的尺寸因子，与每纳米点数相对应。绘出点（x,y)后再从新上述的过程，经过足够多的次数电子云的图貌就由点的疏密反应出来了。

随机数的真与伪
   蒙特卡罗方法模拟问题都离不开随机抽样，于是就需要产生各种概率分布的随机变量，其中最简单也最重要的随机变量是在［0，1］上均匀分布的随机变量。我们把［0，1］上均匀分布随机变量的抽样值称之为随机数。因此在电子计算机上如何产生这种“真正”的随机变量，就变得特别重要。
   随机抽样的一种方法是物理方法，它应用某种物理现象的随机性，在计算机上设有专门的附件实现，称为“随机数发生器”。它所产大的随机数倒是“真”的随机数，但其缺点是过程一去不复返，不能进行重复检查，而且设备费用昂贵，不实用。另一种是目前广泛使用的数学方法，它用迭代过程实现，每个相邻数 由其前一个数 ，或前一组的数通过某些算术或逻辑运算求得。显然，这样的一系列数就不是真正的随机数，但它们只要能通过一系列局部随机性检验，就可当作随机数来使用。科学家们定义这种数为“伪随机数”，其优点是在计算机中只要贮存几个初值，其他数就可通过迭代算得，速度快，费用低，可以重复。这正适合蒙特卡罗模拟计算需要，“伪”的就比“真”的有用。
   产生伪随机数的方法又是一门特殊学科，要涉及数论等理论，这里就不愿详述了。选择方法必须注意使产生的伪随机数序列有随机性好、在计算机上容易实现、省时、周期长等特点。现在，各种大、中、小型电子计算机的软件系统中一般都有产生伪随机数的专用于程序，连微机及计算器中也有简单的伪随机数产生子程序，可见用途之广。日益广泛的应用
2、假随机数的产生
   真的随机数如同掷骰子一样，产生1－6范围内的随机整数，抽奖用的摇号码机产生的0－9范围内的随机整数，这些真的随机数除了统计规律外无其他规律可循。假随机数就是用某种算法给出的似乎随机的数，既然数的给出是按照某种算法的那么就有某种规律，那么这些数必然有一定的周期，设周期是n那么n+1就一定是第一个数了，此后均要一次重复出现。当然如果周期足够大，可使在整个使用过程中不至于呈现周期性，那么假随机数也是使用的。例如，让计算机的假随机数发生器其周期大于机器的记忆单元数。此外假随机数的统计性质也是表征随机数品质的又一个重要的指标。均匀分布的随机数，即要求数的出现是随机均匀的，也要求数的分布是均匀的。至于评价数的均匀，这里不作介绍。
2、1 均匀随机数的产生
产生均匀随机数常用幂剩余法，例如可按如下的公式产生随机数
xn=cx n-1(modN)
或者写成xn=cx n-1－Nint（cx n-1/N)
其中c、N为给定的常数，给出x0后，就可以用上式依次给出x1，x2……等等，如何确定常数c、N、X0是十分关键的问题，人们仍然在不断的研究中，现面给出c、N和X0的一般原则。N的取值一般取N＝2 m-1，其中m为计算机的二进制数的字长，N-1一般取计算机能表现出来的最大整数，如字长16位的N=2^15=32768;字长32位的N=2^31=2147483648,c的取值一般取c=8*M+-3，其中M位任一正整数，如有取c=16897，c=65539，c=397204099等，建议取c~N1/2，这样统计性比较好，关于x0的取值，一般取x0为奇数，如x0=13。可以检验其他条件不变，x0是奇数时，周期是T,偶数时周期时T/2。如N=64,c＝5，x0＝2则得
x1=10
x2=50
x3=58