蒙卡方法初步（讲稿）

蒙卡方法初步蒙卡方法介绍：粒子衰变盒中粒子蒙卡积分：中心极限定理，ImportanceSampling，VegasMC模拟：随机行走Markov链Metropolis高能物理事例产生子：权重，Unweight，直方图，Tree，bin，error拟合，Chi-2，Likelihood

蒙卡方法初步 蒙卡方法介绍： 粒子衰变 盒中粒子 蒙卡积分： 中心极限定理, Importance Sampling , Vegas MC模拟： 随机行走 Markov链 Metropolis 高能物理事例产生子: 权重, Unweight, 直方图, Tree, bin, error 拟合，Chi-2, Likelihood

在很多物理问题中，可以直接模拟物理过程，而不引言：粒子衰变用写下比如描述系统的微分方程。唯一的要求是需要知道几率分布函数。假设在初始时刻t=0，我们有N(O)个x核子，该种核子会辐射衰变，单位时间衰变几率为の，则有dN(t) = -oN(t)dt,1解为可定义平均寿命为=N(t) = N(O)e-ot0假设X核子辐射衰变到Y，Y也会衰变，则有dNx(t)= -0xNx(t),dtdNy(t)-0yNy (t) + OxNx(t)dt可以取大量样本，每个样本包含多个粒子，模拟这些粒子的演进：产生随机数与衰变几率比较

引言: 粒子衰变 假设在初始时刻t=0，我们有N(0)个X核子， 该种核子会辐射衰变，单位时间 衰变几率为ω，则有 解为 可定义平均寿命为 假设 X核子辐射衰变到Y，Y也会衰变，则有 可以取大量样本，每个样本包含多个粒子，模拟这些粒子的演进： 产生随机数 与 衰变几率 比较 在很多物理问题中，可以直接模拟物理过程，而不 用写下比如描述系统的微分方程。唯一的要求是需 要知道几率分布函数

引言：粒子衰变voidmc sampling(intinitial n_particles,intmax time,doubledecay_probability,int *ncumulative)(inttime,np,n_unstable,particle_limit90000n unstable=initial n particles;80000// accumulatethenumberof particlespertime70000step per trial60000ncumulative[o]=initial_n_particles;50000// loopovereachtimestep40000for (time=1;time<= max_time; time++)particlelimit=nunstable30000for(np=1; np<=particle_limit; np++)(20000if(double(rand()/RAND_MAX<=10000decay_probability)(n unstable=n unstable-1;101001人Time1//endofloopoverparticlesncumulative[time]=n_unstable;1//endof loopovertimesteps//endmcsamplingfunction

引言: 粒子衰变 void mc_sampling(int initial_n_particles, int max_time, double decay_probability, int *ncumulative) { int time, np, n_unstable, particle_limit; n_unstable = initial_n_particles; // accumulate the number of particles per time step per trial ncumulative[0] = initial_n_particles; // loop over each time step for (time=1; time <= max_time; time++){ particle_limit = n_unstable; for ( np = 1; np <= particle_limit; np++) { if( double(rand())/RAND_MAX <= decay_probability) { n_unstable=n_unstable-1; } } // end of loop over particles ncumulative[time] = n_unstable; } // end of loop over time steps } // end mc_sampling function

在很多物理问题中，可以直接模拟物理过程，而不引言：盒中粒子用写下比如描述系统的微分方程。唯一的要求是需要知道几率分布函数。t=o单位时间内，左边粒子移动到右边的几率为/N.Choosethenumberof particlesN.Make a loop over time,where the maximum time (or maximum number of steps)should belarger than thenumber of particlesN· For everytime step At there is a probabilityni/N for a move to the right.Comparethisprobabilitywith arandomnumberx..If x<n/N,decrease thenumber of particles inthe lefthalfbyone, i.e.,n=n-1.Else,moveaparticlefromtheright half tothe left,i.e.,n=n+l.Increasethetimebyoneunit(theexternal loop)

引言: 盒中粒子 在很多物理问题中，可以直接模拟物理过程，而不 用写下比如描述系统的微分方程。唯一的要求是需 要知道几率分布函数。 单位时间内，左边粒子移动到 右边的几率为 / l n N

引言：盒中粒子int main(int argc, char* argvll)(srand(unsigned(time(0)));// seedtherandomizerint initial_n_particles,max_time;inttime,random_n,nleft;Particle In Boxofstreamofile("PinBox.dat")10000cout>initial n_particles;9000nleft=initial_n_particles8500max time=2*initial n particles8000for(time=0;time<=max_time;time++)N7500randomn=(int)(initial nparticles7000*double(rand()/RAND_MAX);6500if (random_n <= nleft)6000nleft -= 1;550015000101001000110000100000else (Timenleft += 1;1ofile<<setiosflags(ios::showpoint/ios::uppercase)ofile << setw(15)<< time;ofile<<setw(15)<<nleft;return O;1//endmainfunction

引言: 盒中粒子 int main(int argc, char* argv[]) { srand(unsigned(time(0))); // seed the randomizer int initial_n_particles, max_time; int time, random_n, nleft; ofstream ofile("PinBox.dat"); cout > initial_n_particles; nleft = initial_n_particles; max_time = 2*initial_n_particles; for( time=0; time <= max_time; time++){ random_n = (int)(initial_n_particles *double(rand())/RAND_MAX); if ( random_n <= nleft){ nleft -= 1; } else { nleft += 1; } ofile << setiosflags(ios::showpoint | ios::uppercase); ofile << setw(15) << time; ofile << setw(15) << nleft; } return 0; } // end main function

MC积分

MC积分

MC积分。蒙特卡洛方法最重要的应用在于对高维定积分的计算。采用蒙卡模拟的方法来求解这种确定性的数学物理问题时，必须先选择一个合适的概率模型，利用此概率模型试验所得的随机事件的统计结果等价于待求问题的解。。设想一个物理系统由相互作用的多个粒子组成，如凝聚态物质中的原子、原子中的电子，而每个粒子都可以有几个自由度，要描述这个系统就要涉及高维的积分。●如对m个原子组成的气体，经典配分函数是：=dr...drmeap-β2>U(rii)i<i这里的β=（kBT）-1代表温度，U（rii）是两体相互作用势。除非m相当小时，这个式子要用其他任意一种数值积分计算方法都是不可能算出的

MC积分

。假设我们要计算地图上某个区域的面积，可以均匀的在包围此区域的一投点法个方框里撒上小石子，待计算面积与正方形面积之比即为该区域中的石子数目与总数之比。显然，石子越小，撒的数自越多越均匀，则求得的面积越准确。这就是用MonteCarlo方法计算定积分的原理。因为f（x）的定积分值即为曲线下的面积值Sf(r)dr = So(n/N)这里的S.为方形区域的面积，N是总点数，n是郑入f（x）下面积区域中的点数，该方法从原理上与舍选抽样法是一致的。例如，用MonteCarlo方法求元值。产生一对在一1，1区间中均匀分布的随机数作为点的坐标x，y)值，判断条件2+y2<1是否成立，成立则计数n值，当总点数N足够大时，有：Tlim22=N8

投点法

平均值法平均值法是建立在大数定理的基础之上：。如果函数f(x在[a.bl区间，以均匀的概率分布密度随机地取N个数xi，对每个x；计算出函数值f(x）。大数法则告诉我们这些函数值之和除以N所得的值将收敛于函数的期望值，即：N1)●故有：Var(f)=v2Var(f)N(b-aN2Nf(ci)f()drNi=1。在简单抽样中，我们由均匀分布中选取随机数，并不考虑被积函数的具体情况。因此，被积函数的极大处和极小处有相同的抽样权重。。但是对积分贡献较大的更多在函数值较大处，故直观上就可以看出，非权重的简单抽样效率较低。即为了获得较高计算结果的精度，需要大量的抽样

平均值法 𝑉𝑎𝑟() = 𝑉 2 𝑁2 𝑉𝑎𝑟(𝑓) = 𝑁 1 𝑉 2 𝑉𝑎𝑟(𝑓) 𝑁

中心极限定理与大数法则一样，概率论中的中心极限定理同样是MonteCarlo方法应用于统计计算的基础。它是对积分的蒙特卡洛估计值收敛于该积分的正确结果的收敛程度的研究中心极限定理告诉我们：在有足够大，但又有限的抽样数N的情况下蒙特卡洛估计值的误差分布。即：→ d(β)af/vA。其中μ=一aJf(α)dc，f是函数f(x)的标准偏差，Φ(β)是Gauss正态分布，那么对于积分值得标准偏差。s，就有：ag=f)-μlα，从上式可以看到蒙特卡洛积分的误差和与N有关。在方差固定时，增加模拟次数可以有效地减小误差

中心极限定理