蒙特卡洛方法初探.ppt

随机数的产生及蒙特卡洛随机模拟方法 实验目的 实验内容 学习随机数的产生及蒙特卡洛随机模拟方法的基本过程与方法 1 数学模拟的方法 4 实验作业 3 蒙特卡洛随机模拟实例 2 产生随机数的计算机命令 数学模拟的方法 在实际问题中 面对一些带随机因素的复杂系统 用分析方法建模常常需要作许多简化假设 与面临的实际问题可能相差甚远 以致解答根本无法应用 这时 计算机模拟几乎成为唯一的选择 在一定的假设条件下 运用数学运算模拟系统的运行 称为数学模拟 现代的数学模拟都是在计算机上进行的 称为计算机模拟 计算机模拟可以反复进行 改变系统的结构和系数都比较容易 蒙特卡洛 MonteCarlo 方法是一种应用随机数来进行计算机模拟的方法 此方法对研究的系统进行随机观察抽样 通过对样本值的统计分析 求得所研究系统的某些参数 用蒙特卡洛方法进行计算机模拟的步骤 1 设计一个逻辑框图 即模拟模型 这个框图要正确反映系统各部分运行时的逻辑关系 2 模拟随机现象 可通过具有各种概率分布的模拟随机数来模拟随机现象 一 产生模拟随机数的计算机命令 在Matlab软件中 可以直接产生满足各种分布的随机数 命令如下 1 产生m n阶 a b 均匀分布U a b 的随机数矩阵 unifrnd a b m n 产生一个 a b 均匀分布的随机数 unifrnd a b 当只知道一个随机变量取值在 a b 内 但不知道 也没理由假设 它在何处取值的概率大 在何处取值的概率小 就只好用U a b 来模拟它 2 产生mm nn阶离散均匀分布的随机数矩阵 R unidrnd N R unidrnd N mm nn 当研究对象视为大量相互独立的随机变量之和 且其中每一种变量对总和的影响都很小时 可以认为该对象服从正态分布 若连续型随机变量X的概率密度函数为其中 0为常数 则称X服从参数为的指数分布 指数分布的期望值为 排队服务系统中顾客到达间隔 质量与可靠性中电子元件的寿命通常服从指数分布 例顾客到达某商店的间隔时间服从参数为10 分钟 的指数分布 指数分布的均值为10 指两个顾客到达商店的平均间隔时间是10分钟 即平均10分钟到达1个顾客 顾客到达的间隔时间可用exprnd 10 模拟 设离散型随机变量X的所有可能取值为0 1 2 且取各个值的概率为其中 0为常数 则称X服从参数为的泊松分布 泊松分布在排队系统 产品检验 天文 物理等领域有广泛应用 泊松分布的期望值为 6产生1个参数为n p的二项分布的随机数binornd n p 产生m n个参数为n p的二项分布的随机数binornd n p m n 掷一枚均匀硬币 正面朝上的次数X服从参数为 p的二项分布 X B 1 p 总结 常见分布的随机数产生语句 补充 随机数的产生命令MATLAB可以直接产生满足各种分布的随机数具体命令如下 产生m n阶 0 1 上均匀分布的随机数矩阵rand m n 产生一个 0 1 上均匀分布的随机数rand 产生m n阶 a b 上均匀分布的随机数矩阵unifrnd a b m n 产生一个 a b 上均匀分布的随机数unifrnd a b 产生一个1 n的随机排列 元素均出现且不重复 p randperm n 注意 randperm 6 与unifrnd 1 6 1 6 的区别 产生m n阶均值为mu方差为sigma的正态分布的随机数矩阵normrnd mu sigma m n 产生一个均值为mu方差为sigma的正态分布的随机数normrnd mu sigma 产生m n阶期望值为mu mu 1 的指数分布的随机数矩阵exprnd mu m n 产生一个期望值为mu的指数分布的随机数exprnd mu 注意 产生一个参数为 的指数分布的随机数应输入exprnd 1 产生m n阶参数为A1 A2 A3的指定分布 name 的随机数矩阵random name A1 A2 A3 m n 产生一个参数为为A1 A2 A3的指定分布 name 的随机数random name A1 A2 A3 举例 产生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 二 频率的稳定性模拟 1 事件的频率在一组不变的条件下 重复作n次试验 记m是n次试验中事件A发生的次数 频率f m n 2 频率的稳定性 掷一枚均匀硬币 记录掷硬币试验中频率P 的波动情况 functionliti1 p mm pro zeros 1 mm randnum binornd 1 p 1 mm a 0 fori 1 mma a randnum 1 i pro i a i endpro pronum 1 mm plot num pro 在Matlab中编辑 m文件输入以下命令 在Matlab命令行中输入以下命令 liti1 0 5 1000 在Matlab命令行中输入以下命令 liti1 0 5 10000 练习掷一枚不均匀硬币 正面出现概率为0 3 记录前1000次掷硬币试验中正面频率的波动情况 并画图 在Matlab命令行中输入以下命令 liti1 0 3 1000 例2掷两枚不均匀硬币 每枚正面出现概率为0 4 记录前1000次掷硬币试验中两枚都为正面频率的波动情况 并画图 在Matlab中编辑 m文件输入以下命令 functionliti2 p mm pro zeros 1 mm randnum binornd 1 p 2 mm a 0 fori 1 mma a randnum 1 i randnum 2 i pro i a i endpro pro num 1 mm plot num pro liti2 0 4 100 liti2 0 4 10000 二 几何概率模拟 1 定义 向任一可度量区域G内投一点 如果所投的点落在G中任意可度量区域g内的可能性与g的度量成正比 而与g的位置和形状无关 则称这个随机试验为几何型随机试验 或简称为几何概型 2 概率计算 P A A的度量 S的度量 例5两人约定于12点到1点到某地会面 先到者等20分钟后离去 试求两人能会面的概率 解 设x y分别为甲 乙到达时刻 分钟 令A 两人能会面 x y x y 20 x 60 y 60 P A A的面积 S的面积 602 402 602 5 9 0 5556 functionproguji liti5 mm mm是随机实验次数frq 0 randnum1 unifrnd 0 60 mm 1 randnum2 unifrnd 0 60 mm 1 randnum randnum1 randnum2 proguji 0 forii 1 mmifabs randnum ii 1 20frq frq 1 endendproguji frq mm liti5 10000 proguji 0 5557 三 蒲丰投针实验 法国科学家蒲丰 Buffon 在1777年提出的蒲丰投针实验是早期几何概率一个非常著名的例子 蒲丰投针实验的重要性并非是为了求得比其它方法更精确的 值 而是它开创了使用随机数处理确定性数学问题的先河 是用偶然性方法去解决确定性计算的前导 由此可以领略到从 概率土壤 上开出的一朵瑰丽的鲜花 蒙特卡罗方法 MC 蒲丰投针实验可归结为下面的数学问题 平面上画有距离为a的一些平行线 向平面上任意投一根长为l l a 的针 假设针落在任意位置的可能性相同 试求针与平行线相交的概率P 从而求 蒲丰投针实验 如右图所示 以M表示针落下后的中点 以x表示M到最近一条平行线的距离 以 表示针与此线的交角 针落地的所有可能结果满足 其样本空间视作矩形区域 面积是 针与平行线相交的条件 它是样本空间 子集A 面积是 积分计算symslphi int l 2 sin phi phi 0 pi ans l因此 针与平行线相交的概率为 从而有 特别当时p为统计频率 蒲丰投针实验的计算机模拟 formatlong 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 frame i getframe 描点并取帧title CurrentPoint num2str i Total num2str n counter counter 1 统计针与线相交的次数endendfren counter n pihat 2 l a fren 用频率近似计算 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 例在我方某前沿防守地域 敌人以一个炮排 含两门火炮 为单位对我方进行干扰和破坏 为躲避我方打击 敌方对其阵地进行了伪装并经常变换射击地点 经过长期观察发现 我方指挥所对敌方目标的指示有50 是准确的 而我方火力单位 在指示正确时 有1 3的概率能毁伤敌人一门火炮 有1 6的概率能全部消灭敌人 现在希望能用某种方式把我方将要对敌人实施的1次打击结果显现出来 利用频率稳定性 确定有效射击 毁伤一门炮或全部消灭 的概率 五 复杂概率模拟 选学 分析 这是一个复杂概率问题 可以通过理论计算得到相应的概率 为了直观地显示我方射击的过程 现采用模拟的方式 需要模拟出以下两件事 1 问题分析 1 观察所对目标的指示正确与否 模拟试验有两种结果 每一种结果出现的概率都是1 2 因此 可用投掷一枚硬币的方式予以确定 当硬币出现正面时为指示正确 反之为不正确 2 当指示正确时 我方火力单位的射击结果情况 模拟试验有三种结果 毁伤一门火炮的可能性为1 3 即2 6 毁伤两门的可能性为1 6 没能毁伤敌火炮的可能性为1 2 即3 6 这时可用投掷骰子的方法来确定 如果出现的是 三个点 则认为没能击中敌人 如果出现的是 点 则认为毁伤敌人一门火炮 若出现的是 点 则认为毁伤敌人两门火炮 2 符号假设 i 要模拟的打击次数 k1 没击中敌人火炮的射击总数 k2 击中敌人一门火炮的射击总数 k3 击中敌人两门火炮的射击总数 E 有效射击 毁伤一门炮或两门炮 的概率 3 模拟框图 functionliti6 p mm efreq zeros 1 mm randnum1 binornd 1 p 1 mm randnum2 unidrnd 6 1 mm k1 0 k2 0 k3 0 fori 1 mmifrandnum1 i 0k1 k1 1 elseifrandnum2 i 3k1 k1 1 elseifrandnum2 i 6k3 k3 1 elsek2 k2 1 endendefreq i k2 k3 i endnum 1 mm plot num efreq 在Matlab中编辑 m文件输入以下命令 在Matlab命令行中输入以下命令 liti6 0 5 2000 在Matlab命令行中输入以下命令 liti6 0 5 20000 5 理论计算 6 结果比较 模拟结果与理论计算近似一致 能更加真实地表达实际战斗动态过程 作业 1 掷两枚不均匀硬币 每枚正面出现概率为0 4 记录前1000次掷硬币试验中两枚都为正面频率的波动情况 并画图 2 两船欲停靠同一个码头 设两船到达码头的时间各不相干 而且到达码头的时间在一昼夜内是等可能的 如果两船到达码头后需在码头停留的时间分别是1小时与2小时 试求在一昼夜内 任一船到达时 需要等待空出码头的概率 频率估计概率 functionliti2 p mm pro zeros 1 mm randnum binornd 1 p 2 mm a 0 fori 1 mma a randnum 1 i randnum 2 i pro i a i endpro pro num 1 mm plot num pro 在Matlab中编辑 m文件输入以下命令 习题1解答 liti2 0 4 100 liti2 0 4 10000 3 某厂生产的灯泡能用1000小时的概率为0 8 能用1500小时的概率为0 4 求已用1000小时的灯泡能用到1500小时的概率 频率估计概率 2 在一袋中有10个相同的球 分别标有号码1 2 10 今有放回任取两个球 求取得的第一个球号码为奇数 第二个球的号码为偶数