23蒙特—卡罗

2.3蒙特—卡罗（M-C）模拟模型2.3.1蒙特-卡罗模拟的概要介绍A.用M-C方法求解确定性问题用统计试验方法去解决具有统计性质的数学—物理问题是理所当然的，然而能否用统计试验法或者叫做随机模拟方法去解决确定性的数学—物理问题呢？答案是肯定的。蒙特-卡罗（Monte-Carlo）方法就是一种通过随机变量的统计试验去求解数学—物理问题或者工程问题的一种数值计算方法，下面我们举几个简单的例子，予以说明。（1）问题一，求圆周率（π）的问题。众所周知，半径为R的圆，其面积S与半径R之间存在如下关系式。S=πR2则SR2外接圆的方形面积为S＇=4R2442RRss，这样求取值的问题，便转化为求取两个面积比例（ss）的问题，而这个面积比例值是可以通过随机试验获得。圆周的方程为２２２＝＋Ryx，经标准化后xy221。我们可以在（-1，1）区间内随机地取值（xi,yi），如果xyii221，则取随机变量ni=1，否则取“0”。这样，41lim1ssnNNiiN（2）问题二，求取任意函数f(x)的积分值问题。设为f（x）为定义在（0，1）区间的任意函数，如右图所示，求取其积分值。)1)(0()(102xfdxxf右图曲线下的面积即为f(x)的积分值θ2如果在正方形内随机投点（xi,yi），则（xi,yi）位于曲线f(x)下面的概率P（yf(x)）就是积分值θ210)(0))((xfdydxxfyP设随机变量η叫投点失败叫投点成功)(0)(1iiiiixfyifxfyif经过大量的相互统计独立的随机投点试验，则根据统计理论中的中心极限定律：NiiNN121lim样本的标准差，就是积分值2的统计无偏估计。SNniiN12112()由此可见，通过统计试验不仅可以得到积分值2，而且可以随时给出它的精度统计估计值。M-C方法的特点是对被积函数f(x)没有任何要求，如果我们用解析法去求其积分值，那么我们只能对一些特定的函数形态，才能给出其相应的积分值，否则就得应用各种近似计算方法去求取其值，所以M-C模拟方法是一种适应性极大的数值计算方法。这是M-C方法的第一个显著优点，但同时我应该注意到，M-C方法要求的独立试验次数是相当大的，否则计算精度就不高，换言之，M-C方法只有在当今计算机高度发达的今天，才能充分施展其才能。从另一个角度分析，如何提高计算效率，降低方差，加速收敛速度成为M-C方法成功与否的重要内容。事实上，我们可以重新设计，求取f(x)积分值的M-C算法。设xi为在（0，1）区间内均匀分布的随机变量。则NiiNxfN11)(1lim，两种不同的算法1与2都可以给出函数f(x)的积分值，这说明M-C模拟具有相当好的灵活性，其关键问题是如何将一个具体问题设计成一个概率试验模型，这可算作M-C方法的第二个显著特点吧。如果把上述例子推广到更一般的情况，设随机变量η是一个复杂的多变量函数（）xxxm12,,通过随机模拟，可以得到抽样值12，，N，统计地处理这些抽样数据，给出它的概率分布和各阶矩的估计值，通常条件下，求取的数学期望值，）（E与方差)(估计值已经足够，则问题便得到了解决，根据统计理论NdtePENNNtNii2%5.95211212则为成立，如果取置信水平以概率）（则由此可见：①M-C的自然收敛速度相当慢，为1N的数量级，换言之，要提高一个计算精度量级（即ε减小一个量级），则试验次数要成百倍地增长。②在模拟过程中，可随时估计其计算精度，便可避免盲目计算。③模拟精度与概率型的维数（m）无关，因此M-C方法特别适用于多维问题的数值求,解。④降低的方法，将会极大地改善M-C方法的收敛速度。B.加速收敛的方法我们在上节中给出了两种计算f(x)积分值的M-C模拟计算方法。实际上这两种算法的效率是不同的，为了能客观地比较它的效率，如用T1，T2与21,分别代表两种算法完成一次模拟计算所需的平均时间和它们的方差，则：211222TTC为它们的效率比，实验结果表明T2T1，现在需要讨论21和22的来源及其比较1022111)()(1dxxfxfNNii则此处θ代表f(x)数学期望值。1算法的方差是由f(x)对它的数学期望的偏差所引起的，故称它为“期望值估算法”。对2的模拟计算方法，实际上我们必须选定两个随机变量r1与r2NiiirrNrfrrfrrr1,2,122212121),(1)(0)(1),(统计量当当从概率类型角度看，这是一个典型的二项式概率分布类型，若设η取“1”的事件叫做p，η取“0”的事件叫做q，则p+q为一个必然事件在离散实验中二项式概率型的方差,为Npq2，此处N为总实验次数，当离散实验推广至连续型实验，N则代表积分区间，当今的积分区间为（0，1），所以N=1，此时pq2，所以1010102210210212222)(2)()(1.)()(11dxxfdxxfdxxfdxxfdxxf）（）（10)(1)()()()()()(1)()()()()(1)()()(2)(21122210102101021010101010102101021221010221022TTCdxxfxfdxxfdxxfdxxfdxxfdxxfdxxfdxxfdxxfdxxfdxxfxfdxxfdxxfdxxf由此证明了算法1比算法2更有效（收敛更快）其实102121)(,xfdxxx）（，所以在1的计算中已经用了）（21,rr对r1的数学期望值f(r2)，相对于）（21,rr的统计计算，它减少了一个变量的方差值，所以1222就很自然了，我们称2的算法“随机投点”法。结论为，如果在模拟计算过程中，在一处能用理论分析得到的数学期望值代替该处的统计模拟，即可减少结果的误差。下面讨论几种降低方差的具体抽样方法。①重要抽样法在上面所举的例子中，无论是x还是y均在（0，1）区间内均匀随机取样，但是由于f(x)不是均匀分布于x轴上，所以相同的△x对dxxf10)(的贡献率是不相等的，那么如果能在权重大的x部分加密采样，在权重小的部分减少采样密度，则可达到减少估计方差的目的。.)()()()(1010dxxgxgxfdxxf如果把g(x)看作被积函数fxgx()()的抽样概率密度函数，则其数学期望值仍为.如令hxfxgx()()().统计量NiihxhN1)(10)()(lim222210222dxxgxfnNh这只是从理论上证明了，对X采用不均匀采样，可以使方差等于“0”，但此式没有实用价值，因为还没有具体解决不均匀采样的方式问题。如果gx()与fx()具有相似性，换言之101)(Cdxxh,C为常数，则101)(dxxCh如设RN为一个在（0，1）区间内均匀采样的随机数，并令ixdxxChRN0)(，并由此反求出xi的值，显然xi在（0，1）区间内为非均匀分布，对积分值贡献大的区域xi密集，对积分值贡献小的区域xi稀疏。101010)()()()()(dxCxgxChdxxgxhdxxf由于假定了)(xf与)(xg具有相似性，所以h(x)与)(xf及)(xg在统计意义上均为不相关的函数。NiNcxigNdxcxgdxxchdxxBdxxAdxxBxABABdxAxBxABdxABAdxBBAAdxxBxA110)(1)(1)()()()()(0)()()(时，则当令则）不相关，则（）与（只要）（②相关抽样法利用随机变量之间的相关性，可以达到减少模拟方差的目的。若设fx()在（0，1）区间上是单调函数，则可采用对称化的处理方法，取gxfxfx()()()121，则随机变量)1()(21)(iiixfxfxg仍以θ为它的数学期望值，由于f(x)的单调性质，所以f(xi)与)()1(iixgxf具有负相关性质，所以的方差为)1(2)1()(41212iigxfxfD。)1()(,)(221iiixfxfxf与为其中的相关系数，由于是负相关性质，所以0，所以g212。2.3.2.辐射传输积分方程一般三维辐射传输方程可表达为:),(),(),()(),()(),('4rQdrIrprrIrrIse:为散射系数成。的吸收与散射过程所构：为消光系数，由介质。：为光子传输方向矢量：为坐标矢量。其中ser在此处我们暂时用I代表辐射亮度，因为我们在此处用符号L表示一种微分算子。在下节中我们用符号L代表LAI重取代坐标z，）（,rp为散射相位函数，并满足归一化条件，4),(1),,(rQdrp为源函数，它可以代表介质自身的热辐射，亦可以代表由外界进入的直射辐射，公式*可改变为：4),(),,()()(),(1eeseQdrIrPrrrIrPI）（如果把此等式的左边视为一种微分算子，用L表示，把等式的右边视为一种积分算子，用S表示，则公式可表示为eQSILI。IfQLISLIee令）（）（11，其物理意义为光子与介质的碰撞密度，则)()(11eeeQLISLf。此处L1为L的逆算子，即11LL，如果L为微分算子，那么1L必定是积分算子，其形式可以通过假定公式右边为零时，求解公式*左边的一阶微分方程的积分解而获得。因此：**),(),()()(),,(}),(exp{)('00*4rddrfrrrPdttrrrfesete）（其中t*为从r出发沿到达空间边界的距离。**为积分形式的辐射传输方程，以上积分辐射传输方程的获得主要靠数学符号，算子获得。比较抽象，不容易被理解，下面我们借助对物理过程的分析，同样可以获得辐射传输积分方程。如果把辐射在介质中的传输过程视为光子与介质的碰撞事件所组成的无后效的马尔柯夫链，其转移密度为Kxx(,)，换言之在相空间中（包括位置空间和辐射传输矢量空间）x处发生了碰撞，则在（xdxx,）处发生第二次碰撞的概率为),(xxKdx，设想在相空间中xr（）,处发生碰撞，试问光子被散射后穿过面元ds的概率是什么？（见图）rrrrdsrrrrrPes)('cos1)()(),(当光子穿越圆柱体时，圆柱体内的所有介质将对光子产生衰减作用，其衰减的表达式为：dhrrre),(exp)(，其中(,)rr为由rr到的光学距离，则发生在圆柱体内的绝对概率为dsdhrrrrYrPrrees2'),(exp)(),()()(，除以体积元dv=ds.dh，即为概率密度函数。）（rrrrrrrrrPrrxxkees2),(exp)()()()(),(令fn为第n次碰撞概率。)(),()(:1.1xdxfxxKxffkfxnnnn即。因为00021nnonffkfffff则有……。根据泛函数论中的不动点原理''')(),()(xdxfxxKxfx将上面的),('xxK表达式代入上式，即可获得公式**。把积分——微分形式的辐射传输