蒙特卡罗方法及其应用.doc

计算机处理之蒙特卡罗方法及其应用【标题】 蒙特卡罗方法及其应用 【摘要】 蒙特卡罗方法是一种随即抽样方法，建立一个与求解有关的概率模型或随即现象来求得所要研究的问题的解。这种利用计算机进行模拟的抽样方法以其精度高，受限少等优点广泛应用于数理计算，工程技术，医药卫生等领域。本文介绍蒙特卡罗方法的简要内容，起源，基本思路及应用优点，并简要介绍了一些蒙塔卡罗方法在相关医学方面的应用，并提出了一些今后发展与应用上的展望。【关键词】 蒙特卡罗方法 基本内容 应用 【正文】 一 蒙特卡罗方法简介 1 概述蒙特卡罗(Monte Carlo) 方法, 又称随机抽样法，统计试验法或随机模拟法。是一种用计算机模拟随机现象,通过仿真试验,得到实验数据,再进行分析推断,得到某些现象的规律或某些问题的求解的方法。蒙特卡罗方法的基本思想是,为了求解数学、物理、工程技术或生产管理等方面的问题,首先建立一个与求解有关的概率模型或随机过程,使它的参数等于所求问题的解,然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征,最后给出所求解的近似值。概率统计是蒙特卡罗方法的理论基础,其手段是随机抽样或随机变量抽样。对于那些难以进行的或条件不满足的试验而言,是一种极好的替代方法。蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，很少受几何条件限制，收敛速度与问题的维数无关。例如在许多工程、通讯、金融等技术问题中,所研究的控制过程往往不可避免地伴有随机因素,若要从理论上很好地揭示实际规律,必须把这些因素考虑进去。理想化的方法是在相同条件下进行大量重复试验,采集试验数据,再对数据进行统计分析,得出其规律。但是这样需要耗费大量的人力、物力、财力,尤其当一个试验周期很长,或是一个破坏性的试验时,通过试验采集数据几乎无法进行,此时蒙特卡罗方法就是最简单、经济、实用的方法。因此它广泛应用在粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。蒙特卡罗方法研究的问题大致可分为两种类型,一种是问题本身是随机的;另一种本身属于确定性问题,但可以建立它的解与特定随机变量或随机过程的数字特征或分布函数之间的联系,因而也可用随机模拟方法解决,如计算多重积分,求解积分方程、微分方程、非线性方程组,求矩阵的逆等。2 起源蒙特卡罗方法的起源可以追溯到18世纪著名的蒲丰问题, 1777年,法国科学家蒲丰(Buffon)提出用投针试验计算圆周率值的问题。2.1蒲丰问题蒲丰问题是一个古典概率问题,在平面上有彼此相距为2a的平行线,向平面任意投一长度为2l的针,假定la,容量算出针与平行线相交的概率为由(2)式可知,要求解就必须知道p,而采用常用的方法是无法得到p 的,然而,从统计学的角度却可以通过实验来得到p,这就是进行投针实验. 投针实验N次可能有n次针与平行线相交,当N充分大时,可以认为显然,实验次数N越多, p的近似程度越好。需要指出的是,上述由投针试验求得的近似值的方法,是进行真正的试验,并统计试验结果,要使获得的频率值与概率值偏差小,就要进行大量的试验,这在实际中,往往难以做到。所以,在现代计算机技术出现之前,用频率近似概率的方法抑或称为雏形时代的蒙特卡罗方法并没有得到实质上的应用。2.2蒙特卡罗方法与蒲丰问题随着计算机和计算机技术的迅速发展,可以非常方便地利用计算机模拟随机实验。用数值模拟方法代替上述真正的投针实验,是利用均匀分布于( 0, 1)之间的随机数序列,并构造出随机投针的数学模型,然后进行大量的随机统计并求得到的近似值。如图1建立坐标系,平面上一根针的位置可以用针中心Ml 的坐标x和针与平行线的夹角来决定,在y方向上的位置不影响相交性质,任意投针,意味着x与都是任意取的. 但的范围可限于 0, , x的范围可限于 0, a ,在这种情况下,针与平行线相交的数学条件是其次,怎样模拟投针呢? 亦即如何产生任意的 x, , x在 0, a 任意取值,意味着x在 0, a 上取哪一点的概率都一样,即x的概率密度函数为由此,产生任意( x, )的过程就变为由f1 ( x)抽样x,由f2 ()抽样的过程,容易得到式中, 1 , 2均为(0, 1)上均匀分布的随机数,只要随机数的均匀性和独立性良好,如此构造的数值模型就很好地模拟了实际试验中的一次投针,并用下式判断是否相交且记录统计结果:是相交几率p的估计值,这样就实现了用数值方法模拟真正投针实验。3蒙特卡罗方法的基本思路和特点用蒙特卡罗方法求解问题时,应建立一个概率模型,使待解问题与此概率模型相联系,然后通过随机试验求得某些统计特征值作为待解问题的近似解,与此相似,在一些物理问题,如核裂变、直流气体放电等过程中,粒子的输运过程及粒子输运总效应,也是可以与某些概率过程联系起来,例如,电子与原子、分子、离子的碰撞过程,实际上就是与碰撞截面有关的概率过程,这样,从数学物理特征来说,类似于用随机投针方法计算的近似值,确定条件下的核裂变、直流气体放电中粒子的输运过程及粒子输运的总效应可以用多次掷骰子的方法近似求出。随着现代计算机技术的出现和飞速发展,用计算机模拟概率过程,实现多次模拟试验并统计计算结果,进而可获得所求问题的近似结果,计算机的大存储量、高运算速度使得在短时间内,获得精度极高且内容丰富的模拟结果,在历史上,也正是原子弹工程研究初期阶段的工作,为模拟裂变物质的中子随机扩散,提出了运用大存储量、高运算速度计算机的要求,这也成为当时推动计算机技术发展的重要动力,也就是在第二次世界大战期间,冯诺依曼和乌拉姆两人把他们所从事的与研制原子弹有关的秘密工作对裂变物质的中子随机扩散进行直接模拟以摩纳哥国的世界闻名赌城蒙特卡罗(Monte Carlo)作为秘密代号来称呼。用赌城名比喻随机模拟,风趣又贴切,很快得到了广泛接受,此后,人们便把这种计算机随机模拟方法称为蒙特卡罗方法。二 蒙特卡罗方法在相关医学方面的应用（一） 蒙特卡罗方法在辐射防护中的应用 1 蒙特卡罗方法与MCNP程序蒙特卡罗方法利用已知的光子反应截面数据, 模拟各种微观物理过程, 通过概率抽样对源粒子的行为进行跟踪, 决定每次碰撞后次级粒子的运动方向和速度。根据需要对相应的物理量进行统计, 逐次跟踪下去, 就可以得到所需的结果。该方法相当于一种计算机模拟实验。由于射线与物体作用是一个随机过程, 所得到的宏观物理量又是一个统计值, 对复杂条件下辐射场的计算, 射线衰减与散射过程及空间物质的几何分布有关, 做准确的解析困难很大。在这种情况下, 蒙特卡罗方法是很有效的求解方法。在反应堆及实验装置中, 常用吸收性很强的物质作为中子和光子的屏蔽材料, 求中子或光子经过不同介质, 不同厚度的屏蔽层后的穿透概率和能量分布。当屏蔽物的形状复杂, 散射各向异性, 材料介质不均匀, 核反应截面与能量、位置有关时, 迁移方程难以用数值方法求解, 用蒙特卡罗方法能够得到满意的结果。因此, 可将蒙特卡罗方法应有于辐射屏蔽防护上。我们采用MCNP作为计算程序。MCNP是一个大型多功能的蒙特卡罗计算程序, 可处理复杂场所三维几何结构的中子- 光子耦合输运问题。MCN P 具有较强的通用性, 在源描述、空间物质的几何分布上具有很大的灵活性, 可处理任意三维几何结构问题, 适用面宽, 现已用于射线无损检测系统、辐射屏蔽、核仪器设计和保健物理等许多问题上。它可以很好的用于跟踪计算、决定辐射剂量、物理实验模拟、辐射屏蔽防护上。对光子的输运问题,MCNP详细处理了各种微观物理过程。MCNP程序通过一个输入文件IN P 和有关元素的截面数据文件对物理问题进行计算。输入文件包括描述问题所必需的全部信息, 由包含不同输入信息的数据卡片组成, 卡片具有指定的格式。在每一卡片中填写量化的数据信息。输入卡片按类主要分为栅元卡、曲面卡、数据卡三个部分。栅元卡和曲面卡描述物体分布的空间几何信息, 每一个几何体通过栅元由描述几何体各表面的曲面按一定关系构成, 空间几何越复杂, 需要的曲面卡和栅元卡就越多。数据卡包括问题(光子、中子) 类型、栅元物理参数、曲面物理参数、源描述、材料描述、结果计数描述、问题截断条件等。另外还有一些专门的数据卡片提供降低方差、减少计算所需时间的技巧方法。使用该程序主要是对它的输入卡IN P 的记录形式、计算结果输出、误差估计、减少相对误差技巧的掌握。减小相对误差的主要方法有: (1) 增加输运粒子数; (2)强迫碰撞; (3) 增加粒子的重要性; (4) 源偏倚; (5) 使用能量截断卡。运用这些方法可使在运行相同粒子数情况下的相对误差降低。 2 蒙特卡罗方法的优势 蒙特卡罗方法在辐射防护领域的应用是该方法的重要的应用领域之一。由于受物理条件的限制, 为了得到所求结果, 必须借助理论计算。蒙特卡罗方法具有逼真地描述真实的物理过程的特点。具体地说, 蒙特卡罗方法具有以下几个方面的优势:一、由于某些核实验或者工程项目耗资巨大且需要很长的时间, 而在理论上又无法进行推导, 利用蒙卡方法可有效节省资金和时间。二、由于蒙卡方法是利用计算机进行模拟实验, 可有效避免放射性物质对环境的污染。三、对于有些实验, 可通过对模拟结果和实验结果的比较,达到减小误差的目的。对于外照射防护, 一般采取下面三种方法中的一种, 或几种联合应用: 1、缩短受照时间; 2、增大与辐射源之间的距离; 3、在人与辐射源之间增加辐射屏蔽。辐射屏蔽在辐射防护领域占有相当重要的地位, 用MCN P 程序进行屏蔽计算, 最重要的就是要设置好粒子源和粒子探测器, 在输入文件中要设定好源项和计数卡, 在计数卡设定好粒子探测器的种类、位置。（二） 蒙特卡罗方法在光动力疗法中的应用 光动力疗法( photodynamic therapy, PDT)利用光敏剂受光激发,诱发化学反应产生光毒性物质,选择杀伤病变组织。蒙特卡罗方法计算大量光子迁移轨迹,统计复杂几何形状非均匀生物组织模型光子宏观能量分布,建立两种不同光学参数组织光动力剂量数学模型,可提供分析光动力疗法选择性损伤新途径。 蒙特卡罗方法认为光子与悬浮粒子碰撞,步长和方向改变遵循统计规律。逐个计算大量光子迁移轨迹,统计光子宏观能量分布。利用光动力剂量,光敏剂浓度和能量分布之间关系,建立光动力剂量数学模型。模拟过程中,假定氧是充足的,忽略它对组织光学性质影响;忽略光敏剂滞留对组织光学参数影响,有研究证明, 光敏剂滞留使吸收系数只增加0. 01cm- 1 cm;忽略光漂白效应对光敏剂的影响。光损伤是个阈值现象,每条光动力剂量曲线可被视为阈值曲线。增加光剂量(或光敏剂浓度)使损伤阈值向组织深处移动,表现为光动力剂量曲线向纵深发展。用组织坏死光动力剂量阈值评价模型中选择性损伤程度。光动力疗法对肿瘤组织损伤程度与光敏剂剂量和光子能量分布有关。 采用蒙特卡罗方法,利用光动力剂量、光敏剂浓度和能量分布之间关系,建立两种不同光学参数组织光动力剂量数学模型,研究光动力疗法中选择性光损伤。蒙特卡罗方法逐个计算大量光子迁移轨迹,统计嵌于正常组织的肿瘤中光子宏观能量分布。相比正常组织,肿瘤组织表现出较高吸收率。利用损伤阈值衡量光动力疗法选择性损伤程度。光动力剂量数学模型研究光剂量和光敏剂浓度相互制约关系,提供分析光动力疗法选择性损伤新途径。（三） 蒙特卡罗方法在肿瘤放射物理学中的应用粒子与物质相互作用时服从统计学规律,发生作用的位置、作用的形式(如对光子而言,有光电效应、康普顿效应、电子对效应),发生作用后粒子可能被吸收或散射,散射粒子的运动方向和能量、2次作用位置间的距离等参数均是随机变量。蒙特卡罗方法可以模拟粒子与物质相互作用的全过程,通过模拟10万甚至100万个粒子的输运过程,就可以比较精确地计算出粒子束与物质相互作用的宏观特征,如注量分布、吸收剂量分布。用蒙特卡罗方法解粒子输运问题一般包括3个过程:源分布抽样过程,产生粒子的初始状态;空间、能量和运动方向的随机游动过程,产生粒子的运动状态序列;记录贡献与分析结果,记录每个粒子对所求量的贡献并分析所求量的误差。蒙特卡罗方法的优点是可以处理粒子疏运的各种复杂情况。（四）蒙特卡罗方法确定生物组织体内的温度分布利用计算机技术进行生物组织体内温度场的实时模拟和重构是当前生物医学工程领域研究的热点,也是激光医学研究的热点。由于生物组织兼含固体、液体和气体,因此其热传递包括传导、对流和辐射等多种方式。有关生物组织热传递的研究,可以溯源到1948 年Pennes将人体手臂简化为圆柱模型,并建立了所谓“生物传热方程”。在激光医学及其技术的应用中,人们普遍感兴趣的是激光辐照靶组织诱发热传递的宏观表现,这时热传导起主要作用,组织体内热辐射和对流对温度场的影响可以忽略,即便这样,由于生物组织结构的复杂性,特别是边界条件的复杂性,基于Pennes 思想的热传递方程一般只能采用数值近似的有限差分法(FDM) 、有限元法(FEM)和边界元素法(BEM)等进行求解。近年来,这些方法已成功地被运用于求解激光辐照生物组织诱发温度分布等相关问题。然而这些数值分析方法都存在局限性,即必须同时求解整个研究区域,并还必须处理复杂的生物体几何特征与生物体的属性,编程相对复杂,计算耗时,从而影响了它们在临床上的实际应用。Haji 等率先提出基于蒙特卡罗思想的概率方法求热传导方程,随后Kowsary 等成功地把这种方法应用于各向异性热传导方程的求解,最近Deng 等人则把这种思想运用于生物热传递方程的求解。与传统的数值分析方法相比较,蒙特卡罗方法求解时对研究对象的维度和几何细节依赖小,如采用蒙特卡罗方法可以独立于空间内的其它点而单独获取其一期望点的解,这对于只要分析或求解一些孤立点温度的情形尤其有用,因此这种方法很适合运用于生物组织热传递的分析和求解,特别是联合红外热成像技术无损获取组织体表面的温度分布,可以实现由体表温度分布重构和反演生物组织体内的温度场。（五） 蒙特卡罗方法在X荧光无损分析分析中的应用同位素源激发X 射线荧光分析(Xrayfluorescence analysis)，XRF胙为微量元素分析工具，以其设备简单、操作方便、灵敏度高、多元素同时分析和对样品非破坏性分析等优点，在环境科学、生物医学、地质、冶金、考古等许多领域内都得到了广泛地应用。然而，受基体吸收和增强的影响，常规XRF无损分析(直接测试)因找不到准确可靠的、同类物质的标准参考物质进行比较，致使所测结果不准确；很多文献采用PIXE(Proton inducedXray emission)进行无损分析，适当提高了分析结果精度，但PIXE需启动运行费用很高的加速器，不经济也不能推广为常规实验方法。可见，研究提高XRF样品无损分析准确度问题具有重大的现实意义。XRF在样品无损分析中的关键技术就是如何配制与待测样品在化学组成和物理形状上相同或相近的标准样品。由于待测样品各种各样、成分复杂、形态各异，很难通过实验手段实现。蒙特卡罗方法具有能够逼真地描述具有随机性质的事物的特点及物理实验过程、受几何条件限制少、收敛速度与问题的维数无关、能够同时计算多个方案和多个未知量、误差容易确定、程序结构简单且易于实现等优点，可以很方便地解决这一难题。国内外许多研究者都曾经尝试过把蒙特卡罗方法应用于XRF中，且都对各模拟过程的散射增强影响给出了定量估计。但上述研究者都只是停留在方法的探讨上，还未有人将这一方法推广到珍贵样品的无损定量分析领域，用来解决XRF在样品无损定量分析时标准样品难制、直接测试精度不够的难题。（六） 蒙特卡罗方法模拟量表量表作为测量工具广泛应用于医疗卫生领域的研究。在量表设计过程中, 应答条目的级数以及条目数量是影响信度的重要因素。若以蒙特卡罗(Monte Garlo) 方法研究二者对内部一致性信度的影响, 首先需要模拟出符合一定要求的矩阵以表示各条目的测量得分。由MA THWORK 公司推出的MA TLAB 是建立在向量、数组和矩阵上的软件,利用它所提供的各种函数可以比较容易地产生出符合研究要求数据阵。下面介绍第一军医大学有关人员利用自编程序做的相关研究。1 对数据阵的要求1.1假设量表中的应答条目是一长100mm 的直线, 各个被测者在各条目上的得分是0 100 的正整数。1.2如果量表各条目均数能准确测量被测症状, 则每个条目的得分应与被测症状在人群中的分布相同或相似, 除非许多条目设计的不合理, 不能测量出实际被测特征。以被测特征呈正态分布为例, 假设所设计的量表所有条目在某被测人群中测量的平均得分服从均数为50、标准差为16.7的正态分布,则某个条目在该被测人群的平均得分亦应服从或近似服从正态分布, 且均数不会偏离50太远。1.3为使该量表的内部一致性信度系数足够大, 各列数据两两间应有一定的正相关性。2 产生数据阵2.1 用下式估计所需的样本容量该研究取用自编程序1 作图。若取内部一致性信度= 0. 8 0. 9 (以下称原始信度) ,则取样本容量n= 180。2.2 用e= normrand(50,16.7,180,1) 产生一列共180 个服从均数为50、标准差为16.7 的正态分布的数据。2.3设h 为正整数, 则用z = h3 *(- 0. 5+ rand (180, 1) ) 可产生一列- h/2 - h/2 的服从均匀分布的数据。其中, h 是满足特定的样本容量n、条目数量k 以及原始信度通过运行自编程序2 经过反复尝试而确定的。2.4将步骤2、3 产生出的两列数据c 和z 中的各对应元素两两相加得到一列数据, 通过取其中0 100 的数据, 并将该范围外的其他数据由0 100 内的服从均匀分布的随机数代替,最后对所有数据四舍五入所得的一列数据作为该量表中第一个条目对n 个被测者的测量得分。如此重复k 次, 所得数据阵x 即为k 个条目的n 个得分。相应的语句如下:关键语句:3模拟次数t 的确定从理论上讲, 模拟次数越多, 精度越高, 但所耗机时也越长; 模拟次数过少, 则误差过大。因此, 需权衡模拟次数t 和误差。设原始信度系数的总体均数为, 标准差为, 根据中心极限定理服从标准正态分布。用t 个系数的标准差s 作为的估计值, 则按95% 的置信水平可得此为确定t 的计算公式。通过运行自编程序3 得到图2。 可见, 当t 100 时, E对t 的变化非常敏感; 当100t400 时, 对t 的敏感程度逐渐降低并趋向于稳定; 当n400时, 增加t 对减小的影响已不明显。因此, 模拟次数取为400。4检验数据4.1分布检验若取条目数k= 10, 运行自编程序4 产生400个这样的数据阵, 各列数据(即各条目得分) 的均数分别为50. 0115, 50. 0362, 49. 9428, 49. 9320, 50. 1913, 50. 0072,50. 0326, 50. 1156, 50. 0567, 50. 0876, 平均为5010414。标准差分别为2210863, 2219496, 22. 6689, 23. 2700, 21. 3576,21. 5275, 23. 3341, 21. 9156, 22. 6159, 22. 3370, 平均为22.4062。随机抽取其中任一个数据阵, 对各列数据在统计软件SPSS 1010上用one- sample Kolmogo rov-Sm irnoy 检验, 结果P 值均远大于0.05, 一般在0.5 以上, 各列数据均符合正态分布。4.2相关性检验运行程序4 可得所有400 个数据阵的相关矩阵的平均值:可见各列数据两两间相关系数均在0.33 以上, 这就保证了内部一致性信度系数能达到要求, 事实上系数平均为0.8352。4.3随机性检验随机抽取上述任一数据阵, 对各列数据两两间的随机性检验在SPSS 1010 软件上用runs(游程检验)进行, 结果证明各列中的数据均是随机出现的。因此, 所产生的数据阵符合作者的要求。三 总结及展望本文主要讲了蒙特卡罗方法的概述起源和思路特点及其在医学领域的应用，其中包括：在辐射防护、光动力疗法、肿瘤放射物理学、生物阻止体内的温度分布、X荧光无损分析、模拟量表等医学方面的应用。蒙特卡罗方法是重要的数值计算方法,可以模拟许多大型的、难以实观的复杂实验或社会行为过程。由于有了蒙特卡罗方法,计算机已经不仅仅是数学和理论科学的重要工具,并正日益成为实验物理学、应用科学、社会科学及基础研究人员的第二实验场所。但是，究其起源及发展现状，我们可以有以下展望：1 集团抽样方法一种灵活而简便的蒙特卡罗技巧 集团抽样方法不仅具有灵活、简便和普适性强的特点，而且，通过选择台适的变量还可以明显地提高计算效率，真可谓是一种适于蒙特卡罗方法通用软件中使用的技巧。2 零方差技巧一实现零方差的重要途径：相关技巧 零方差技巧几乎是在蒙特卡罗方法出现的同时便被提了出来。众所周知，方差每减小一个数量级，蒙特卡罗抽样数便可以减小两个数量级。因零方差技巧如此鼓舞人心，有关的研究工作一直在进行中。不过，已见到的零方差技巧几乎全是建立在偏倚抽样基础上的，忽略了用其它蒙特卡罗技巧解决的重要性。 比较建立在相关原理基础上的上述零方差技巧和建立在偏倚抽样原理基础上的零方差技巧。前者不需要改变原随机游动过程，而且不需要增加任何附加条件。可是后者，不仅需要改变原随机游动过程，而且还需要增加不少附加条件，这些附加条件往往还无法实现。3 自改善方法一蒙特卡罗方法通用软件的智能化问题无论是哪一种零方