Monte-Carlo模拟教程.ppt

MonteCarlo模拟 内容提纲 1 引言2 MonteCarlo模拟基本思想3 随机数生成函数4 应用实例举例5 排队论模拟6 MonteCarlo模拟求解规划问题 引言 Introduction MonteCarlo方法 蒙特卡罗方法 又称随机模拟方法 属于计算数学的一个分支 它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的 亦称统计模拟方法 statisticalsimulationmethod 利用随机数进行数值模拟的方法 MonteCarlo名字的由来 MonteCarlo是摩纳哥 monaco 的首都 该城以赌博闻名 NicholasMetropolis 1915 1999 Monte Carlo Monaco MonteCarlo方法的基本思想 蒙特卡罗方法 或称计算机随机模拟方法 是一种基于 随机数 的计算方法 源于美国在第二次世界大战研制原子弹的 曼哈顿计划 该计划的主持人之一数学家冯 诺伊曼用驰名世界的赌城 摩纳哥的MonteCarlo 来命名这种方法 为它蒙上了一层神秘色彩 蒙特卡罗方法的基本思想很早以前就被人们所发现和利用 早在17世纪 人们就知道用事件发生的 频率 来决定事件的 概率 19世纪人们用蒲丰投针的方法来计算圆周率 上世纪40年代电子计算机的出现 特别是近年来高速电子计算机的出现 使得用数学方法在计算机上大量 快速地模拟这样的试验成为可能 蒲丰投针实验 法国科学家蒲丰 Buffon 在1777年提出的蒲丰投针实验是早期几何概率一个非常著名的例子 蒲丰投针实验的重要性并非是为了求得比其它方法更精确的 值 而是它开创了使用随机数处理确定性数学问题的先河 是用偶然性方法去解决确定性计算的前导 由此可以领略到从 概率土壤 上开出的一朵瑰丽的鲜花 蒙特卡罗方法 MC 蒲丰投针实验可归结为下面的数学问题 平面上画有距离为a的一些平行线 向平面上任意投一根长为l l a 的针 假设针落在任意位置的可能性相同 试求针与平行线相交的概率P 从而求 蒲丰投针实验 如右图所示 以M表示针落下后的中点 以x表示M到最近一条平行线的距离 以 表示针与此线的交角 针落地的所有可能结果满足 其样本空间视作矩形区域 面积是 针与平行线相交的条件 它是样本空间 子集A 面积是 symslphi int l 2 sin phi phi 0 pi ans l因此 针与平行线相交的概率为 从而有 蒲丰投针实验的计算机模拟 formatlong 设置15位显示精度a 1 l 0 6 两平行线间的宽度和针长figure axis 0 pi 0 a 2 初始化绘图板set gca nextplot add 初始化绘图方式为叠加counter 0 n 2010 初始化计数器和设定投针次数x unifrnd 0 a 2 1 n phi unifrnd 0 pi 1 n 样本空间 fori 1 nifx i l sin phi i 2 满足此条件表示针与线的相交plot phi i x i r counter counter 1 统计针与线相交的次数frame counter getframe 描点并取帧endendfren counter n pihat 2 l a fren 用频率近似计算 figure 2 movie frame 1 播放帧动画1次 一些人进行了实验 其结果列于下表 蒙特卡罗投点法是蒲丰投针实验的推广 在一个边长为a的正方形内随机投点 该点落在此正方形的内切圆中的概率应为该内切圆与正方形的面积比值 即n 10000 a 2 m 0 fori 1 nx rand 1 a y rand 1 a if x a 2 2 y a 2 2 a 2 2 m m 1 endenddisp 投点法近似计算的 为 num2str 4 m n x y o a 2 a 2 基本思想 由上面的例子可以看出 当所求问题的解是某个事件的概率 或者是某个随机变量的数学期望 或者是与之有关的量时 通过某种试验的方法 得出该事件发生的频率 再通过它得到问题的解 这就是蒙特卡罗方法的基本思想 蒙特卡罗方法的关键步骤在于随机数的产生 计算机产生的随机数都不是真正的随机数 由算法确定的缘故 如果伪随机数能够通过一系列统计检验 我们也可以将其当作真正的随机数使用 rand seed 0 1 rand 1 每次运行程序产生的值是相同的 rand state sum 100 clock rand rand 1 每次重新启动matlab时 输出的随机数不一样 注意 产生一个参数为 的指数分布的随机数应输入exprnd 1 产生m n阶参数为A1 A2 A3的指定分布 name 的随机数矩阵random name A1 A2 A3 m n 举例 产生2 4阶的均值为0方差为1的正态分布的随机数矩阵random Normal 0 1 2 4 name 的取值可以是 详情参见helprandom norm or Normal unif or Uniform poiss or Poisson beta or Beta exp or Exponential gam or Gamma geo or Geometric unid or DiscreteUniform 非常见分布的随机数的产生 逆变换方法 Acceptance Rejection方法最早由VonNeumann提出 现在已经广泛应用于各种随机数的生成 基本思路 通过一个容易生成的概率分布g和一个取舍准则生成另一个与g相近的概率分布f 为要模拟服从给定分布的随机变量 用生成一个易于实现的不可约遍历链作为随机样本 使其平稳分布为给定分布的方法 称为马氏链蒙特卡罗方法 马氏链蒙特卡罗方法 1 连续型随机变量 以指数分布为例 symstxlambda Fx int lambda exp lambda t t 0 x 分布函数symsr Fxinv finverse Fx x 求反函数Fxinv subs Fxinv x r 替换反函数变量x为rFxinv inline Fxinv x Fxinv 3 rand 产生参数lambda 3指数分布的随机数 指数分布随机数产生函数已经提供exprnd 1 3 1 1 2 离散型随机变量 以离散分布为例 x 2 4 6 8 px 0 1 0 4 0 3 0 2 以下为程序片段Fx 0 forn 1 length px Fx Fx sum px 1 n endr rand index find r Fx x index 1 1 通用离散分布随机数产生程序 例1在我方某前沿防守地域 敌人以一个炮排 含两门火炮 为单位对我方进行干扰和破坏 为躲避我方打击 敌方对其阵地进行了伪装并经常变换射击地点 经过长期观察发现 我方指挥所对敌方目标的指示有50 是准确的 而我方火力单位 在指示正确时 有1 3的射击效果能毁伤敌人一门火炮 有1 6的射击效果能全部毁伤敌人火炮 现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来 确定有效射击的比率及毁伤敌方火炮的平均值 分析 这是一个概率问题 可以通过理论计算得到相应的概率和期望值 但这样只能给出作战行动的最终静态结果 而显示不出作战行动的动态过程 为了能显示我方20次射击的过程 现采用模拟的方式 举例 需要模拟出以下两件事 2 当指示正确时 我方火力单位的射击结果情况 1 观察所对目标的指示正确与否 模拟试验有两种结果 每一种结果出现的概率都是1 2 因此 可用投掷一枚硬币的方式予以确定 当硬币出现正面时为指示正确 反之为不正确 模拟试验有三种结果 毁伤一门火炮的可能性为1 3 即2 6 毁伤两门的可能性为1 6 没能毁伤敌火炮的可能性为1 2 即3 6 这时可用投掷骰子的方法来确定 如果出现的是 三个点 则认为没能击中敌人 如果出现的是 点 则认为毁伤敌人一门火炮 若出现的是 点 则认为毁伤敌人两门火炮 问题分析 i 要模拟的打击次数 k1 没击中敌人火炮的射击总数 k2 击中敌人一门火炮的射击总数 k3 击中敌人两门火炮的射击总数 E 有效射击比率 E1 20次射击平均每次毁伤敌人的火炮数 符号说明 模拟框图 模拟结果 理论计算 结果比较 虽然模拟结果与理论计算不完全一致 但它却能更加真实地表达实际战斗动态过程 用蒙特卡洛方法进行计算机模拟的步骤 1 设计一个逻辑框图 即模拟模型 2 根据流程图编写程序 模拟随机现象 可通过具有各种概率分布的模拟随机数来模拟随机现象 3 分析模拟结果 计算所需要结果 随机投掷均匀硬币 验证国徽朝上与朝下的概率是否都是1 2 n 10000 给定试验次数m 0 fori 1 nx randperm 2 1 randperm 2 returnsa随机排列oftheintegers1 2 y x 1 ify 0 0表示国徽朝上 1表示国徽朝下m m 1 endendfprintf 国徽朝上的频率为 f n m n 小实例一 投掷硬币 随机投掷骰子 验证各点出现的概率是否为1 6 n 10000 m1 0 m2 0 m3 0 m4 0 m5 0 m6 0 fori 1 nx randperm 6 y x 1 switchycase1 m1 m1 1 case2 m2 m2 1 case3 m3 m3 1 case4 m4 m4 1 case5 m5 m5 1 otherwise m6 m6 1 endend 输出结果 小实例二 投掷骰子 设某班有m个学生 则该班至少有两人同一天生日的概率是多少 小实例五 生日问题 n 1000 p 0 m 50 设该班的人数为50fort 1 na q 0 fork 1 mb randperm 365 a a b 1 endc unique a 去掉重复数字从小到大重新排序iflength a length c p p 1 endendfprintf 至少两人同一天生日的概率为 f n p n 试验五源程序 clear m 50 p1 1 365 p2 1 365 m 365 ones 1 m p p1 p2 p 1 prod p fprintf 至少两人同一天生日的概率为 f n p 小实例五的理论值计算 排队问题随机模拟 排队论主要研究随机服务系统的工作过程 在排队系统中 服务对象的到达时间和服务时间都是随机的 排队论通过对随机服务现象的统计研究 找出反映这些随机现象平均特性的规律指标 如排队的长度 等待的时间及服务利用率等 1 系统的假设 1 顾客源是无穷的 2 排队的长度没有限制 3 到达系统的顾客按先后顺序依次进入服务 先到先服务 在某商店有一个售货员 顾客陆续来到 售货员逐个地接待顾客 当到来的顾客较多时 一部分顾客便须排队等待 被接待后的顾客便离开商店 设 1 顾客到来间隔时间服从参数为0 1的指数分布 对顾客的服务时间服从 上的均匀分布 排队按先到先服务规则 队长无限制 假定一个工作日为8小时 时间以分钟为单位 1 模拟一个工作日内完成服务的个数及顾客平均等待时间t 2 模拟100个工作日 求出平均每日完成服务的个数及每日顾客的平均等待时间 单服务员的排队模型模拟 w 总等待时间 ci 第i个顾客的到达时刻 bi 第i个顾客开始服务时刻 ei 第i个顾客服务结束时刻 xi 第i 1个顾客与第i个顾客之间到达的间隔时间 yi 对第i个顾客的服务时间 符号说明 c1 b1 c3 c4 c5 c2 e1 b2 e2 b3 e3 b4 e4 b5 ci ci 1 xiei bi yibi max ci ei 1 t 思路分析 流程框图 用蒙特卡洛法解非线性规划问题 先产生一个随机数作为初始试验点 以后则将上一个试验点的第j个分量随机产生 其它分量不变而产生一新的试验点 这样 每产生一个新试验点只需一个新的随机数分量 当K MAXK或P MAXP时停止迭代 在Matlab软件包中编程 共需三个 文件 randlp m mylp m lpconst m 主程序为randlp m mylp mfunctionz mylp x 目标函数z 2 x 1 2 x 2 2 x 1 x 2 8 x 1 3 x 2 转化为求最小值问题 randlp mfunction sol r1 r2 randlp a b n 随机模拟解非线性规划debug 1 a 0 试验点下界b 10 试验点上界n 1000 试验点个数r1 unifrnd a b n 1 a b 均匀分布随机数矩阵r2 unifrnd a b n 1 sol r1 1 r