概率统计中的MonteCarlo方法及其建模应用PPT课件.ppt

概率统计中的Monte Carlo方法及其建模应用 1 主要内容 蒙特卡洛方法应用实例 2 作业内容 3 蒙特卡洛方法介绍 1 2020 3 1910 30 南京信息工程大学 蒙特卡洛方法介绍 蒙特卡洛起源与发展 1 随机数的产生原理 2 2020 3 1910 30 南京信息工程大学 蒙特卡洛起源与发展 1 2020 3 1910 30 南京信息工程大学 MonteCarlo的起源 MonteCarlo方法 又称随机模拟方法 对研究的系统进行随机观察抽样 通过对样本值的统计分析 求得所研究系统的某些参数它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来 曼哈顿计划 主持人之一 数学家 冯 诺伊曼用驰名世界的赌城 摩纳哥最大的城市MonteCarlo 来命名这种方法 Monte Carlo Monaco 2020 3 1910 30 南京信息工程大学 MonteCarlo方法的应用 物理 核物理 热力学与统计物理 粒子输运问题等数学 多重积分 解微分方程 非线性方程组求解等工程领域 真空技术 水力学 激光技术等经济学领域 期权定价 项目管理 投资风险决策等其他领域 化学 医学 生物 生产管理 系统科学 公用事业等方面 随着科学技术的发展 其应用范围将更加广泛 2020 3 1910 30 南京信息工程大学 MC的起源和发展 Buffon试验 假设平面上有无数条距离为1的等距平行线 现向该平面随机投掷一根长度为l的针 l 1 则我们可计算该针与任一平行线相交的概率 这里 随机投针指的是 针的中心点与最近的平行线间的距离X均匀地分布在区间 0 1 2 上 针与平行线的夹角 不管相交与否 均匀的分布在区间 0 上 此时 针与线相交的充要条件是 2020 3 1910 30 南京信息工程大学 从而针线相交的概率为根据上式 若我们做大量的投针试验并记录针与线相交的次数 则由大数定理可以估计出针线相交的概率p 从而得到 的估计值 Buffon试验 2020 3 1910 30 南京信息工程大学 functionpi estimation buffon llength n llength是针的长度 n是随机实验次数frq 0 xrandnum unifrnd 0 0 5 1 n phi unifrnd 0 pi 1 n forii 1 nif xrandnum 1 ii llength sin phi 1 ii 2 frq frq 1 endendpi estimation 2 llength frq n 2020 3 1910 30 南京信息工程大学 buffon 6 1000 piguji 3 1662 buffon 6 10000 piguji 3 1072 buffon 6 100000 piguji 3 1522 buffon 6 1000000 piguji 3 1386 buffon 6 1000000 piguji 3 1451 buffon 6 1000000 piguji 3 1418 buffon 6 1000000 piguji 3 1448 buffon 6 1000000 piguji 3 1405 buffon 6 1000000 piguji 3 1394 2020 3 1910 30 南京信息工程大学 建立统计模型 主要特征参量方面要与实际问题或系统相一致 问题的解对应于模型中随机变量的概率分布或其某些数字特征根据模型中各个随机变量的分布 在计算机上产生随机数 实现一次模拟过程所需的足够数量的随机数 进而进行随机模拟实验根据概率模型的特点和随机变量的分布特性 设计和选取合适的抽样方法 并对每个随机变量进行抽样 包括直接抽样 分层抽样 相关抽样 重要抽样等 按照所建立模型进行仿真试验 计算 求出问题的随机解统计分析模拟试验结果 给出问题的估计以及其精度估计 必要时 还应改进模型以降低估计方差和减少试验费用 提高模拟计算的效率 用蒙特卡洛方法进行计算机模拟的步骤 2020 3 1910 30 南京信息工程大学 蒙特卡洛模拟的理论基础 大数定律 贝努里 Bernoulli 大数定律中心极限定理 2020 3 1910 30 南京信息工程大学 蒙特卡洛模拟的误差分析 由中心极限定理可知 这表明 不等式近似地以概率1 成立 上式也表明 收敛到 的阶为O n 1 2 通常 蒙特卡罗方法的误差 定义为 2020 3 1910 30 南京信息工程大学 蒙特卡洛方法的特点 MonteCarlo方法及其程序结构简单产生随机数 通过大量简单重复抽样和简单计算计算相应的值收敛速度与问题维数无关MonteCarlo方法的收敛速度为O n 1 2 与一般数值方法相比很慢 因此 用MonteCarlo方法不能解决精确度要求很高的问题MonteCarlo方法误差 只与标准差 和样本容量n有关 而与样本所在空间无关 即MonteCarlo方法的收敛速度与问题维数无关 而其他数值方法则不然 MonteCarlo方法的适用性强MonteCarlo方法对多维问题的适用性在解题时受问题条件限制的影响较小例如 要计算s维空间中的任一区域Ds上的积分 2020 3 1910 30 南京信息工程大学 随机数的产生原理 2 2020 3 1910 30 南京信息工程大学 常用分布的随机数生成 1 均匀分布U a b 产生m n阶 a b 均匀分布的随机数矩阵R unifrnd a b m n 产生m n阶 0 N 离散均匀分布的随机数矩阵R unidrnd N R unidrnd N mm nn 2 正态分布N 2 产生m n阶均值为 标准差为 的正态分布的随机数矩阵 R normrnd m n 3 指数分布E 产生m n阶均值为 的指数分布的随机数矩阵 R exprnd m n 2020 3 1910 30 南京信息工程大学 常用分布的随机数生成 4 泊松分布 产生m n阶均值为 的泊松分布的随机数矩阵R poissrnd m n 5 二项分布B n p 产生m n个参数为n p的二项分布的随机数R binornd n p m n 2020 3 1910 30 南京信息工程大学 一般分布随机数产生方法 基本方法有如下三种 逆变换法复合抽样方法筛选法 2020 3 1910 30 南京信息工程大学 逆变换法 直接抽样方法 设随机变量X的分布函数为F x 定义F 1 y inf x F x y 0 y 1定理设随机变量U服从 0 1 上的均匀分布 则X F 1 U 的分布函数为F x 因此 要产生来自F x 的随机数 只要先产生来自U 0 1 的随机数 然后计算F 1 u 即可 其步骤为 2020 3 1910 30 南京信息工程大学 1 离散型分布 即 inf x F x u 其中令 1时 为了实现由任意离散型分布的随机抽样 直接抽样方法是非常理想的 2020 3 1910 30 南京信息工程大学 1 离散型分布 例1 掷骰子点数的抽样按照离散分布的直接抽样 1 由U 0 1 抽取u即 等价于 也可使用如下更简单的方法 2020 3 1910 30 南京信息工程大学 functiondiscreterandom liti11 mm Random unifrnd 0 1 1 mm fori 1 mmif floor 6 Random 1 i 6 Random 1 i Random 1 i 6 Random 1 i elseRandom 1 i floor 6 Random 1 i 1 endendcdfplot Random 1 离散型分布 2020 3 1910 30 南京信息工程大学 2 连续分布 对于连续型分布 如果分布函数F x 的反函数F 1 x 能够解析表示 则直接抽样方法是 2020 3 1910 30 南京信息工程大学 指数分布为连续型分布 其一般形式如下 其分布函数为 则 1 由U 0 1 抽取u因为1 u也是 0 1 上均匀随机数 可将上式简化为 例2 指数分布 2020 3 1910 30 南京信息工程大学 Randnum 2 log unifrnd 0 1 1 1000 cdfplot Randnum 例2 产生指数分布的随机数 2020 3 1910 30 南京信息工程大学 逆变换法 直接抽样方法 连续性分布函数的直接抽样方法对于分布函数的反函数容易实现的情况 使用起来是很方便的 但是对于以下几种情况 直接抽样法是不合适的 分布函数无法用解析形式表达 因而无法给出反函数的解析形式分布函数有解析形式 但是反函数的解析形式给不出来反函数有解析形式 但运算量很大下面叙述的抽样方法是能够克服这些困难的比较好的方法 2020 3 1910 30 南京信息工程大学 复合抽样方法 复合抽样方法的基本思想是由kahn提出的 考虑如下复合分布 其中f2 x y 为给定Y y时X的条件密度 F1 y 为Y的分布函数如果X密度函数f x 难于抽样 而X关于Y的条件密度函数f2 x y 以及Y的分布F1 y 均易于抽样 则X的随机数抽样 首先从分布F1 y 中抽样YF1 然后再从密度函数f2 x YF1 中抽样确定Xf2 x YF 2020 3 1910 30 南京信息工程大学 指数函数分布的一般形式为 例3指数函数分布的抽样 则 使用复合抽样方法 抽取服从该分布的样本 生成10000个随机数 画经验分布函数 n 5 引入如下两个密度函数 2020 3 1910 30 南京信息工程大学 对应分布函数为 使用复合抽样方法 首先从f1 y 中抽取y 从U 0 1 中抽取u 令 2020 3 1910 30 南京信息工程大学 指数分布 均值为1 y 再由f2 x yf1 中抽取xf 2020 3 1910 30 南京信息工程大学 functionliti110 n mm R1 unifrnd 0 1 mm 1 R2 unifrnd 0 1 mm 1 x zeros mm 1 y 1 R1 1 n x log R2 ycdfplot x functionliti110 n mm R1 unifrnd 0 1 mm 1 x zeros mm 1 y 1 R1 1 n x exprnd 1 y cdfplot x 使用复合抽样方法 首先从f1 y 中抽取y 再由f2 x yf1 中抽取X Liti110 5 10000 2020 3 1910 30 南京信息工程大学 筛选抽样 定理 设X的密度函数f x 且可将其表示成f x ch x g x 其中0 g x 1 c 1是常数 h 是一个密度函数 令U和Y分别服从U 0 1 和h y 则在U g Y 的条件下 Y的条件密度为 依据上述定理 若h 易于抽样 则X的抽样可如下进行 2020 3 1910 30 南京信息工程大学 例4令圆半径为R0 该圆上的点到圆心的距离为r r的密度函数和分布函数分别为 生成10000个随机数 画经验分布函数 1 直接抽样方法 缺点 开方运算在计算机上很费时间 functionliti111 1 R0 mm R unifrnd 0 1 mm 1 x R0 sqrt R cdfplot x Liti111 1 3 10000 2020 3 1910 30 南京信息工程大学 2 筛选抽样方法 取 则抽样框图为 显然 没有必要舍弃u1 u2的情况 此时 只需取 亦即 2020 3 1910 30 南京信息工程大学 functionliti111 2 R0 mm R1 unifrnd 0 1 mm 1 R2 unifrnd 0 1 mm 1 x zeros mm 1 forii 1 mmifR1 ii 1 R2 ii 1 x ii 1 R0 R2 ii 1 elsex ii 1 R0 R1 ii 1 endendcdfplot x Liti111 2 3 10000 2020 3 1910 30 南京信息工程大学 随机向量的抽样 直接抽样 1 分量X与Y相互独立 随机向量 X Y 的抽样若X Y相互独立 分布函数分别为FX x FY y 则从FX x 中抽样x 从FY y 中抽样y 得到二维随机变量 X Y 的抽样 x y 2 按照条件分布 抽取随机向量 X Y 的样本二维随机向量 X Y 的密度函数f x y 将其表示如下 其中fl x f2 y x 分别为X的边缘密度函数和给定X x条件下Y的条件密度函数 即 2020 3 1910 30 南京信息工程大学 根据上述边缘密度函数和条件密度函数 二维分布f x y 的抽样方法为 首先由fl x 中抽取xf1 再由f2 y xf1 中抽样yf2 xf1 yf2 就是该二维分布的一个抽样 例5对下面二维分布进行抽样 将f x y 写为 其中 2020 3 1910 30 南京信息工程大学 直接抽样方法 1 首先由fl x 中抽取xf1 2 再由f2 y xf1 中抽样yf2 2020 3 1910 30 南京信息工程大学 functionliti112 mm R1 unifrnd 0 1 mm 1 R2 unifrnd 0 1 mm 1 x 1 R1 y R1 log R2 plot x y axis 0 010010 holdonholdoff liti112 1000 2020 3 1910 30 南京信息工程大学 随机向量的抽样 筛选抽样 定理 设 X Y 的密度函数为f x y 且f x y ch x y g x y 其中0 g x y 1 c 1 h x y 是一个密度函数 令U和Z Xh Yh 分别服从U 0 1 和h x y 则在U g Z 的条件下 Z的条件密度为 若h 易于抽样 则从f x y 中的抽样 X Y 可如下进行 2020 3 1910 30 南京信息工程大学 例6生成单位圆内均匀分布的10000个随机数 并画散点图 相当于 1 1 1 1 上均匀分布 2020 3 1910 30 南京信息工程大学 mm 10000 xRandnum zeros 1 mm yRandnum zeros 1 mm ii 1 whileii mmRandnum1 unifrnd 1 1 Randnum2 unifrnd 1 1 s Randnum1 2 Randnum2 2 ifs 1 xRandnum 1 ii Randnum1 yRandnum 1 ii Randnum2 ii ii 1 endendplot xRandnum yRandnum 2020 3 1910 30 南京信息工程大学 练习 1 生成单位球内均匀分布的1行10000列随机数 并画散点图 2 设密度函数为生成 n 2 期望 5 的随机数10000个 并绘制经验分布图 为常数 2020 3 1910 30 南京信息工程大学 练习 3 设球壳内半径为R0 外半径为R1 点到球心的距离为r 则r的分布密度函数为生成10000个随机数 并画散点图 2020 3 1910 30 南京信息工程大学 蒙特卡洛方法应用实例 概率 计算模拟分析 1 定积分的MC计算 2 系统可靠性模拟计算 3 2020 3 1910 30 南京信息工程大学 概率 计算模拟分析 1 2020 3 1910 30 南京信息工程大学 复杂概率模拟 例6 在我方某前沿防守地域 敌人以一个炮排 含两门火炮 为单位对我方进行干扰和破坏 为躲避我方打击 敌方对其阵地进行了伪装并经常变换射击地点 经过长期观察发现 我方指挥所对敌方目标的指示有50 是准确的 而我方火力单位 在指示正确时 有1 3的概率能毁伤敌人一门火炮 有1 6的概率能全部消灭敌人 现在希望能用某种方式把我方将要对敌人实施的1次打击结果显现出来 利用频率稳定性 确定有效射击 毁伤一门炮或全部消灭 的概率 2020 3 1910 30 南京信息工程大学 2020 3 19 48 复杂概率模拟 分析 这是一个复杂概率问题 可以通过理论计算得到相应的概率 为了直观地显示我方射击的过程 现采用模拟的方式 1 问题分析需要模拟出以下两件事 1 观察所对目标的指示正确与否模拟试验有两种结果 每一种结果出现的概率都是1 2 因此 可用投掷一枚硬币的方式予以确定 当硬币出现正面时为指示正确 反之为不正确 2020 3 1910 30 南京信息工程大学 复杂概率模拟 2 当指示正确时 我方火力单位的射击结果情况模拟试验有三种结果 毁伤一门火炮的可能性为1 3 即2 6 毁伤两门的可能性为1 6没能毁伤敌火炮的可能性为1 2 即3 6 这时可用投掷骰子的方法来确定 如果出现的是 三个点 则认为没能击中敌人 如果出现的是 点 则认为毁伤敌人一门火炮 若出现的是 点 则认为毁伤敌人两门火炮 2020 3 1910 30 南京信息工程大学 复杂概率模拟 2 符号假设i 要模拟的打击次数 k1 没击中敌人火炮的射击总数 k2 击中敌人一门火炮的射击总数 k3 击中敌人两门火炮的射击总数 E 有效射击 毁伤一门炮或两门炮 的概率 2020 3 1910 30 南京信息工程大学 复杂概率模拟 3 模拟框图 2020 3 1910 30 南京信息工程大学 复杂概率模拟 functionliti26 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 2020 3 1910 30 南京信息工程大学 复杂概率模拟 liti26 0 5 2000 liti26 0 5 20000 2020 3 1910 30 南京信息工程大学 复杂概率模拟 5 理论计算模拟结果与理论计算近似一致 能更加真实地表达实际战斗动态过程 2020 3 1910 30 南京信息工程大学 定积分的MC计算 2 2020 3 1910 30 南京信息工程大学 定积分的MC计算 事实上 不少的统计问题最后都归结为定积分的近似计算问题 相对于其它方法 用MC方法比一般的数值方法有优点 主要体现在它的误差与维数m无关 下面考虑一个简单的定积分为了说明问题 我们首先介绍两种求 的简单的MC方法 然后给出几种较为复杂而更有效的MC方法 2020 3 1910 30 南京信息工程大学 方法简述 设a b有限 0 f x M x y a x b 0 y M 并设 X Y 是在 上均匀分布的二维随机变量 其联合密度函数为 则易见是 中y f x 曲线下方的面积假设我们向 中进行随机投点 则点落在y f x 下方的概率p 随机投点法 2020 3 1910 30 南京信息工程大学 若我们进行了n次投点 其中n0次点落入y f x 曲线下方 则用频率n0 n来估计概率p 即 那么我们可以得到 的一个估计 2020 3 1910 30 南京信息工程大学 注1随机投点法的思想简单明了 且每n次投点结果服从二项分布 故 其中注2可证是 的无偏估计 若用估计的标准差来衡量其精度 则估计的精度的阶为 2020 3 1910 30 南京信息工程大学 求解定积分的算例 例7计算定积分事实上 其精确解为用随机投点法求解 liti27 0 4 4 1000000 result 7 2336注 增加样本数目 可提高计算精度 但计算时间也会提高 functionresult liti27 a b m mm a是积分的下限 b是积分的上限 m是函数的上界 mm是随机实验次数frq 0 xrandnum unifrnd a b 1 mm yrandnum unifrnd 0 m 1 mm forii 1 mmif cos xrandnum 1 ii 2 yrandnum 1 ii frq frq 1 endendresult frq m b a mm 2020 3 1910 30 南京信息工程大学 设g x 是 a b 上的一个密度函数 改写 基本原理 对积分 其中 X是服从g x 的随机变量 可见 积分可以表示为X的函数的期望 由矩法 若有n个来自g x 的观测值x1 xn 则可给出 的一个矩估计 样本平均值法 2020 3 1910 30 南京信息工程大学 特别地 若a b有限 可取g x 为 a b 上均匀分布 此时 设x1 xn是来自U a b 的随机数 则 的一个估计为 具体步骤为 注可证是 的无偏估计 一般而言 样本均值法要比随机投点法更有效 样本平均值法 2020 3 1910 30 南京信息工程大学 例9计算定积分事实上 其精确解为样本平均值法求解 liti29 0 4 1000 result 7 1854liti29 0 4 10000 result 7 2153litti29 0 4 100000 result 7 2419注 增加样本数目 可提高计算精度 但计算时间也会提高 求解定积分的算例 functionresult liti29 a b mm a是积分的下限 b是积分的上限 积分函数cos x 2 mm是随机实验次数sum 0 xrandnum unifrnd a b 1 mm forii 1 mmsum sum cos xrandnum 1 ii 2 endresult sum b a mm 2020 3 1910 30 南京信息工程大学 几种降低估计方差的MC方法 样本均值法 样本均值法是假设g x 为均匀分布 采用均匀抽样 各xi是均匀分布的随机数 各xi对的贡献是不同 f xi 大则贡献大 但在抽样时 这种差别未能体现出来 重要抽样法 希望贡献率大的随机数出现的概率大 贡献小的随机数出现概率小 从而提高抽样的效率关键因素在于g x 的选取 使得估计的方差较小通过选取与f x 形状接近的密度函数g x 来降低估计的方差 2020 3 1910 30 南京信息工程大学 几种降低估计方差的MC方法 关联抽样法 将需要估计的积分分解成两个积分之差 对 的估计转化为对I1 I2的估计的差 即由于所以 估计方差的大小与I1 I2的估计的相关度有关 若两者的正相关程度越高 则 的估计方差越小 这便是关联抽样法的基本出发点 2020 3 1910 30 南京信息工程大学 例11利用MonteCarlo方法计算一个简单的积分 1 样本平均值法 产生n个U 0 1 随机数x1 xn 则 g x 1 0 x 1 为U 0 1 对应的概率密度 由此 1 7183 2020 3 1910 30 南京信息工程大学 functionresult liti211 a b mm a是积分的下限 b是积分的上限 积分函数 mm是随机实验次数sum 0 xrandnum unifrnd a b 1 mm forii 1 mmsum sum exp xrandnum 1 ii endresult sum mm result liti211 0 1 1000 1 7267result liti211 0 1 10000 1 7199result liti211 0 1 100000 1 7171 2020 3 1910 30 南京信息工程大学 例12利用MonteCarlo方法计算一个简单的积分 2 重要抽样法 利用线性近似 取 0 1 上密度函数 由重要抽样法思想 要选择一个与ex相似的密度函数 我们知道 ex的Taylor展开为 1 7183 设x1 xn是来自g x 的随机数 则 的估计为 2020 3 1910 30 南京信息工程大学 估计步骤 g x 的随机数对应分布函数为 1 产生n个U 0 1 随机数u1 un 则 2 xi 3 2020 3 1910 30 南京信息工程大学 result liti212 0 1 1000 1 7222result liti212 0 1 10000 1 7174result liti212 0 1 100000 1 7185functionresult liti212 a b mm a是积分的下限 b是积分的上限 积分函数 mm是随机实验次数sum 0 urandnum unifrnd a b 1 mm xrandnum 1 sqrt 1 3 unifrnd forii 1 mmsum sum exp xrandnum 1 ii 1 xrandnum 1 ii endresult 1 5 sum mm 2020 3 1910 30 南京信息工程大学 系统可靠性模拟计算 3 2020 3 1910 30 南京信息工程大学 系统的可靠性计算问题 一个元件 或系统 能正常工作的概率称为元件 或系统 的可靠性系统由元件组成 常见的元件连接方式 串联 并联 2020 3 1910 30 南京信息工程大学 例13设两系统都是由4个元件组成 每个元件正常工作的概率为p 每个元件是否正常工作相互独立 两系统的连接方式如下图所示 比较两系统的可靠性 S1 2020 3 1910 30 南京信息工程大学 例14S2 2020 3 1910 30 南京信息工程大学 例13 14 设两系统都是由4个元件组成 每个元件的寿命服从平均寿命为 a1 a2 b1 b2的指数分布 每个元件是否正常工作相互独立 两系统的连接方式如下图所示 求两系统寿命大于T 100的概率 functionRguji liti213 t thetaa1 thetaa2 thetab1 thetab2 mm t是要求系统生存的寿命 thetaa1是元件A1的数学期望 thetaa2是元件A2的数学期望 thetab1是元件B1的数学期望 thetab2是元件B2的数学期望 mm是随机实验次数frq 0 randnuma1 exprnd thetaa1 1 mm randnuma2 exprnd thetaa2 1 mm randnumb1 exprnd thetab1 1 mm randnumb2 exprnd thetab2 1 mm forii 1 mmif randnuma1 1 ii t endend Rguji frq mm 系统1 2020 3 1910 30 南京信息工程大学 functionRguji liti214 t thetaa1 thetaa2 thetab1 thetab2 mm t是要求系统生存的寿命 thetaa1是元件A1的数学期望 thetaa2是元件A2的数学期望 thetab1是元件B1的数学期望 thetab2是元件B2的数学期望 mm是随机实验次数frq 0 randnuma1 exprnd thetaa1 1 mm randnuma2 exprnd thetaa2 1 mm randnumb1 exprnd thetab1 1 mm randnumb2 exprnd thetab2 1 mm forii 1 mmif randnuma1 1 ii t randnumb1 1 ii t pass1 1 elsepass1 0 endif randnuma2 1 ii t randnumb2 1 ii t pass2 1 elsepass2 0 endif pass1 pass2 1frq frq 1 endendRguji frq mm 系统2 2020 3 1910 30 南京信息工程大学 优化问题求解 4 2020 3 1910 30 南京信息工程大学 2020 3 1910 30 南京信息工程大学 基本假设 试验点的第j个分量xj服从 aj bj 内的均匀分布 符号假设 求解过程 先产生一个随机数作为初始试验点 以后则将上一个试验点的第j个分量随机产生 其它分量不变而产生一新的试验点 这样 每产生一个新试验点只需一个新的随机数分量 当K MAXK或P MAXP时停止迭代 2020 3 1910 30 南京信息工程大学 2020 3 1910 30 南京信息工程大学 在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 转化为求最小值问题 2020 3 1910 30 南京信息工程大学 randlp mfunction sol r1 r2 randlp a b n 随机模拟解非线性规划debug 1 a 0 试验点下界b 10 试验点上界n 1000 试验点个数r1 unifrnd a b n 1 n 1阶的 a b 均匀分布随机数矩阵r2 unifrnd a b n 1 sol r1 1 r2 1 z0 inf fori 1 nx1 r1 i x2 r2 i lpc lpconst x1x2 iflpc 1z mylp x1x2 ifz z0z0 z sol x1x2 endendend ToMatlab randlp 返回 2020 3 1910 30 南京信息工程大学 案例分析 4 2020 3 1910 30 南京信息工程大学 坎雷渔业公司问题 克林特坎雷经营着Massachusetts一家拥有50条鳕鱼捕捉船的渔业公司 每个工作日 渔船早上离港 中午作业完毕 每次每条船能捕鱼3500单位 有许多港口都可以停靠并出售鳕鱼 每个港口每条的价格是不确定的 并且变化很大 而且港口之间价格也不一样 另外 每个港口的需求量是有限的 如果一条船比别的船晚到一个港口 那么它的鱼就卖不出去 要倒进海洋中 2020 3 1910 30 南京信息工程大学 1 坎雷渔业公司问题简化 简化问题假设渔业公司只有一条船 每次出海的成本为10 000美元 每次出海捕鱼3500单位 两个港口 格洛斯特 岩石港 可以停靠 格洛斯特是鳕鱼的集散地 价格一直稳定在每单位3 25美元 需求几乎是无限的岩石港比较小 价格较高但波动比较大 岩石港的价格服从均值为3 65标准差为0 20的正态分布 需求量服从表1的离散分布 并且我们假设两个港口之间的价格 需求量之间是相互独立的 每天渔船只能在一个港口停靠并出售它的鳕鱼 而岩石港每天的价格事先并不知道 坎雷想挣得尽可能大的利润 哪一个港口停靠更好 表1岩石港鳕鱼日需求分布表 2020 3 1910 30 南京信息工程大学 1 坎雷渔业公司问题简化 渔船在格洛斯特港停靠的利润G为 但是 停靠在岩石港的利润计算出P没这简单 因为价格和需求量都是不确定的 每天的利润是一个随机变量 为了决定选择哪个港口 下面的问题将是很有帮助的 a 使用岩石港日利润的概率分布大概是什么形状 b 使用岩石港利润高于使用格洛斯特港利润的概率是多少 c 使用岩石港亏本的概率是多少 d 使用岩石港日利润的期望值是多少 e 使用岩石港日利润的标准差是多少 2020 3 1910 30 南京信息工程大学 2 坎雷渔业公司初步分析 定义两个随机变量 PR 岩石港的鳕鱼价格 PR N 3 65 0 202 D 停靠岩石港坎雷面临的需求量 D的分布如表1记F为停靠岩石港的