蒙特卡洛法ppt课件.ppt

MonteCarlo 蒙特卡罗方法又称随机抽样技巧或统计试验方法 半个多世纪以来 由于科学技术的发展和电子计算机的发明 这种方法作为一种独立的方法被提出来 并首先在核武器的试验与研制中得到了应用 蒙特卡罗方法是一种计算方法 但与一般数值计算方法有很大区别 它是以概率统计理论为基础的一种方法 由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程 解决一些数值方法难以解决的问题 因而该方法的应用领域日趋广泛 利用随机模拟和统计实验的方法求某些数学 物理和工程问题近似解的数值方法 计算圆周率的问题 谱丰问题Buffon1777 任投一针的概率意义为 1 针的中点至距它最近的平行线的距离x是 0 a 上均匀分布的随机数 2 针与线的夹角 是 0 上均匀分布的随机数 3 随机数x 互相独立 针线相交的充要条件 若投针n次 相交m次 当n 时 x和 的概率密度为1 a和1 最后得到 值为 一些人进行了实验 其结果列于下表 为了得到具有一定精确度的近似解 所需试验的次数是很多的 通过人工方法作大量的试验相当困难 甚至是不可能的 因此 蒙特卡罗方法的基本思想虽然早已被人们提出 却很少被使用 本世纪四十年代以来 由于电子计算机的出现 使得人们可以通过电子计算机来模拟随机试验过程 把巨大数目的随机试验交由计算机完成 使得蒙特卡罗方法得以广泛地应用 在现代化的科学技术中发挥应有的作用 计算机模拟 通过一定的公式产生服从 0 1 均匀分布的随机数ri i 1 2 通过变换得到服从 0 和 0 a 均匀分布的随机数 如果满足 则针线相交 进行n模拟 其中m次相交 计算高维积分 解椭圆和抛物线型偏微分方程解 代数方程组 计算逆矩阵等 蒙特卡洛方法的收敛性 误差及特点 设某个随机变量的简单样本为 其算术平均为 由强大数定理 蒙特卡罗方法的近似值与真值的误差问题 概率论的中心极限定理给出了答案 该定理指出 如果随机变量序列X1 X2 XN独立同分布 且具有有限非零的方差 2 即f X 是X的分布密度函数 则 当N充分大时 有如下的近似式其中 称为置信度 1 称为置信水平 通常 蒙特卡罗方法的误差 定义为上式中与置信度 是一一对应的 根据问题的要求确定出置信水平后 查标准正态分布表 就可以确定出 下面给出几个常用的 与的数值 关于蒙特卡罗方法的误差需说明两点 第一 蒙特卡罗方法的误差为概率误差 这与其他数值计算方法是有区别的 第二 误差中的均方差 是未知的 工程结构的破坏概率可以表示为 用蒙特卡洛方法 当选取95 置信度时 抽样方差 用相对误差表示 由于一般pf是一个小量 可以近似表示为 如果e 0 1 减小方差的各种技巧 显然 当给定置信度 后 误差 由 和N决定 要减小 或者是增大N 或者是减小方差 2 在 固定的情况下 要把精度提高一个数量级 试验次数N需增加两个数量级 因此 单纯增大N不是一个有效的办法 另一方面 如能减小估计的均方差 比如降低一半 那误差就减小一半 这相当于N增大四倍的效果 因此降低方差的各种技巧 引起了人们的普遍注意 后面课程将会介绍一些降低方差的技巧 蒙特卡洛方法特点 程序结构简单 模拟过程灵活 限制条件少 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程 误差具有概率的特征 即是概率意义的误差误差只与样本标准差和样本大小有关 与样本元素所在的空间无关 收敛速度与问题的维数无关具有同时计算多个方案与多个未知量的能力 收敛速度慢 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程 蒙特卡罗方法可以部分代替物理实验 甚至可以得到物理实验难以得到的结果 用蒙特卡罗方法解决实际问题 可以直接从实际问题本身出发 而不从方程或数学表达式出发 它有直观 形象的特点 受几何条件限制小 在计算s维空间中的任一区域Ds上的积分时 无论区域Ds的形状多么特殊 只要能给出描述Ds的几何特征的条件 就可以从Ds中均匀产生N个点 得到积分的近似值 其中Ds为区域Ds的体积 这是其他数值方法难以作到的 另外 在具有随机性质的问题中 如考虑的系统形状很复杂 难以用一般数值方法求解 而使用蒙特卡罗方法 不会有原则上的困难 收敛速度与问题的维数无关 由误差定义可知 在给定置信水平情况下 蒙特卡罗方法的收敛速度为 与问题本身的维数无关 维数的变化 只引起抽样时间及估计量计算时间的变化 不影响误差 也就是说 使用蒙特卡罗方法时 抽取的子样总数N与维数s无关 维数的增加 除了增加相应的计算量外 不影响问题的误差 这一特点 决定了蒙特卡罗方法对多维问题的适应性 而一般数值方法 比如计算定积分时 计算时间随维数的幂次方而增加 而且 由于分点数与维数的幂次方成正比 需占用相当数量的计算机内存 这些都是一般数值方法计算高维积分时难以克服的问题 具有同时计算多个方案与多个未知量的能力 对于那些需要计算多个方案的问题 使用蒙特卡罗方法有时不需要像常规方法那样逐个计算 而可以同时计算所有的方案 其全部计算量几乎与计算一个方案的计算量相当 例如 对于屏蔽层为均匀介质的平板几何 要计算若干种厚度的穿透概率时 只需计算最厚的一种情况 其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到 另外 使用蒙特卡罗方法还可以同时得到若干个所求量 例如 在模拟粒子过程中 可以同时得到不同区域的通量 能谱 角分布等 而不像常规方法那样 需要逐一计算所求量 误差容易确定 对于一般计算方法 要给出计算结果与真值的误差并不是一件容易的事情 而蒙特卡罗方法则不然 根据蒙特卡罗方法的误差公式 可以在计算所求量的同时计算出误差 对干很复杂的蒙特卡罗方法计算问题 也是容易确定的 一般计算方法常存在着有效位数损失问题 而要解决这一问题有时相当困难 蒙特卡罗方法则不存在这一问题 程序结构简单 易于实现 在计算机上进行蒙特卡罗方法计算时 程序结构简单 分块性强 易于实现 收敛速度慢 如前所述 蒙特卡罗方法的收敛速度为 一般不容易得到精确度较高的近似结果 对于维数少 三维以下 的问题 不如其他方法好 误差具有概率性 由于蒙特卡罗方法的误差是在一定置信水平下估计的 所以它的误差具有概率性 而不是一般意义下的误差 计算结果与系统大小有关 对于大系统或小概率事件的计算问题 计算结果往往比真值偏低 中子穿透问题 已知中子垂直进入厚度为3d的铅壁 设每个中子在铅壁中每次走过d后才与铅原子碰撞 碰撞后随机弹射 走过d后再和第二个铅原子碰撞 如此反复 每个中子可能穿透铅壁 返回 若经10次碰撞后没有穿透或返回则被铅壁吸收 求穿透 返回和吸收的概率 是 0 2 均匀分布 计算机模拟5000个中子 结果穿透的占26 3 吸收的占22 返回的占51 0 1 均匀分布随机数的生成 乘同余法 混合同余法 若C取正奇数 M 2k 4q 1 x0取任意非负数 最大周期T 2k 由具有已知分布的总体中抽取简单子样 在蒙特卡罗方法中占有非常重要的地位 总体和子样的关系 属于一般和个别的关系 或者说属于共性和个性的关系 由具有已知分布的总体中产生简单子样 就是由简单子样中若干个性近似地反映总体的共性 随机数是实现由已知分布抽样的基本量 在由已知分布的抽样过程中 将随机数作为已知量 用适当的数学方法可以由它产生具有任意已知分布的简单子样 随机数的统计检验 参数检验 矩检验 检验随机数子样的均值与理论均值的差异是否显著 在 0 1 上均匀分布的随机数的期望和方差为 设随机变量R共有n个观察值r1 r2 rn 由中心极限定理得 其中 设随机变量R共有n个观察值r1 r2 rn 由中心极限定理得 其中 渐近服从标准正态分布N 0 1 所以可以进行u检验 当给定显著水平后 即可根据正态分布表确定临界值 检验时一般可取显著水平为 0 05 此时的临界值为1 96 当有显著差异 均匀性检验又称频率检验 检验随机数子样的经验频率与理论频率的差异是否显著 把 0 1 等分成k个区间 以 i k i k i 1 2 k 表示第k个小区间 如ri是 0 1 均匀分布随机数的一个样本 则落在任意一个小区间上的概率Pi为这些小区间的长度1 k 故n个样本落在任意小区间内的概率为mi nPi n k 设落在第i个小区间内有ni个 则统计量 均匀性检验 渐近服从 2 k 1 分布 如果给定显著水平 0 05 若 独立性检验 独立性检验 检验随机数子样中前后各数的相关性是否显著 两个随机变量的相关系数可以反映它们之间的线性相关程度 设给定随机数r1 r2 rn 计算前后距离为k的样本相关系数 其中 在统计假设 k 0成立时 当n充分大时 统计量 渐近服从标准正态分布N 0 1 故可以进行u检验 除了上述三种基本检验外 还有其他检验 如组合规律检验 连贯性检验等 任意分布随机数的生成 利用 0 1 均匀分布的随机数可以生成任意分布的随机数 1 反函数法 针对分布函数是显函数时的方法 反函数定理 如果随机数x的分布函数F x 连续 则R F x 是 0 1 上均匀分布的随机变量 直接抽样的基础 如果F x 是显函数 则有x F 1 R 如果已经得到 0 1 上均匀分布的随机数r1 r2 rn xi F 1 ri 则 x1 x2 xn 是 上服从F x 分布 产生 a b 上均匀分布的随机数 给定ri 则有 产生服从指数分布的随机数 给定ri 则有 由于ri为 0 1 上均匀分布的随机数 则有1 ri也为 0 1 上均匀分布的随机数 则有 2 舍选法 已知f x 并集中在有限区间 a b 内 即 步骤 选取 a b 为f x 的取值区间 选取常数c使其在 a b 有cf x r2为服从f x 分布的第一个随机数 如果cf z1 r2舍去 重复 得出一系列服从f x 分布的随机数 这种方法适合于形状比较平坦的f x 如果f x 在某一小区间内的取值变化比较大 该方法就不稳定 产生标准正态分布的随机数 已知标准正态分布的概率密度函数为 选择 并取区间 a a a根据精度取5或6 产生一对 0 1 上均匀分布的随机数r1 r2 并将r1变成 a a 上均匀分布的随机数 如果 如果不满足上述条件 则舍去 继续取一对 0 1 上均匀分布的随机数r1 r2进行舍选 最终得出一系列服从标准正态分布的随机数zi 作变换 的随机数 则取z1为一个服从f x 的随机数 3 离散逼近法 设已知f x 并集中在有限区间 a b 内 将 a b 分成n等份其分点为 a x0 x1 x2 xn b 计算每个子区间中的概率 显然有 又将 0 1 分成n等份 其分点为 0 r0 r1 r2 rn 1 同时有ri ri 1 pi 产生 0 1 均匀分布的随机数r 作变换 4 随机变量函数法 设随机变量Z是其它随机xi i 1 2 n的函数 如果随机变量xi i 1 2 n互相独立 各自均有不同的分布且已知 可以利用上述的方法产生服从各自分布的随机变量xi 代入函数关系得到zj j 1 2 Z落在区间 xi 1 xi 的概率为 只要n充分大 z就可以近似作为服从f x 分布的随机数 几种常见的抽样公式 已知R R1 R2为 0 1 上均匀分布的随机数 标准正态分布 指数分布 X1 x2为互相独立的随机变量 威布尔分布 正态分布 对数正态分布 多维随机变量舍选法 已知n维概率密度函数f x1 x2 xn 并集中在有限区间 a1 b1 a2 b2 an bn 内 选取常数c使其在定义域上有cf x1 x2 xn rn 1为服从f x1 x2 xn 分布的第一个随机数 如果cf z1 z2 zn rn 1舍去 加速收敛 蒙特卡洛方法是用随机变量的样本推断母体的性质 而样本总量是有限的 因此推断总有误差 已证明当n充分大时 估计值的标准误差与成反比 也可以说收敛速度与成反比 因此要将蒙特卡洛方法的计算精度提高一位 理论上必须成百倍地提高样本容量 实际上样本容量应提高一千倍 采用加速收敛的方法可以用同样的样本得到更精确的估计量 其主要的途径是改进抽样方法 对偶变量法 z z 为z的无偏估计 取za z z 2 E za E z za也是无偏估计量 如果z z 负相关 则 负相关随机变量的产生 由n个均匀分布的随机数ri i 1 2 n得到估计量z 由1 ri i 1 2 n得到z 则z z 负相关 重要性抽样法 若存在一个抽样密度函数h X 满足条件 分层抽样法 分层抽样法的概念与重要性抽样法相似 它们都是要使对pi贡献大的抽样出现的次数更多 分层抽样法并不改变原始的密度函数 而是将抽样区间分成一些区间 并使各个区间的抽样点数不同 在贡献大的区间抽取更多的样本 将积分区域G X 0分成M个互不相交的子区域Lj 在每个区间内取Nj个均匀分布在该区间的均匀分布的随机数 其它加速收敛的方法有 相关抽样法 控制变量法 复合重要性抽样法 V空间重要性抽样法等 例 计算积分 一般MC法 方差 重要性取样法0 13 对偶变量法0 05 一般MC法1 例 一简单刚架 受重力w 等效地震载荷kw的作用 设立柱的极限弯矩为M1 横梁的极限弯矩为M2 其相关系数为0 8 假设失效是由于产生塑性铰所致 各杆为等截面且抗弯刚度相同 除w外 k