数学建模 计算机仿真-蒙特卡洛方法
蒙特卡罗算法
xi 1 6xi mod11, ui+1 xi 1 /11 （a 6, m 11 ）
如果令 a 3, x0 1 ,得到序列: 1,3,9,5, 4,1,3,9........ 如果令 a 3, x0 2, 得到序列： 2, 6, 7,10,8, 2, 6.......
14
3. 指数分布 指数分布随机变量X的概率密度为
e x , x 0 f ( x) 0, x 0
指数分布常用来描述寿命问题.
15
二.随机数的生成
1.蒙特卡罗模拟的关键是生成优良的随机数。 2.在计算机实现中,我们是通过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模仿了随机数的性质(如可以通过 统计检验)。我们通常称之为伪随机数(pseudorandom numbers)。 3.在模拟中，我们需要产生各种概率分布的随机数， 而大多数概率分布的随机数产生均基于均匀分布 U(0,1)的随机数。
2
基本思想很早以前就被人们所发现和 利用。17世纪，人们就知道用事件发 生的“频率”来决定事件的“概率”。 19世纪人们用投针试验的方法来决定π。 高速计算机的出现，使得用数学方法 在计算机上大量模拟这样的试验成为 可能。
3
从Buffon(蒲丰)投针问题谈起
Buffon 投针问题:平面上画很多平行线，间距为 a．向此平面投掷长为 l ( l < a) 的 针,此针与任一平行线相交的概率 p。
Monte Carlo Simulation Methods (蒙特卡罗模拟方法)
主要内容: 一. M-C方法概述. 二. 随机数的生成. 三. 模拟训练. 四. 实验题目.
成信院数学与信息科学系 李胜坤
系统建模与仿真第12讲 Monte Carlo蒙特卡洛方法
Nicholas Metropolis (1915-1999)
Monte-Carlo, Monaco
引言(Introduction)
Monte Carlo模拟的应用： 自然现象的模拟： 宇宙射线在地球大气中的传输过程； 高能物理实验中的核相互作用过程； 实验探测器的模拟 数值分析： 利用Monte Carlo方法求积分
2
3.141528 3.141528 3.141509 3.141553 3.141506
3
3.141527 3.141521 3.141537 3.141527 3.141538
n
(i )2
si
i1
n 1
0.000012
0.0000032
s si / n
ua s t(0.683, n 1) 0.0000033
引言(Introduction)
Monte Carlo模拟在实际研究中的作用
引言(Introduction)
Monte Carlo模拟的步骤： 1. 根据欲研究的系统的性质，建立能够描述该系统特性的理 论模型，导出该模型的某些特征量的概率密度函数； 2. 从概率密度函数出发进行随机抽样，得到特征量的一些模 拟结果； 3. 对模拟结果进行分析总结，预言系统的某些特性。
k n 1
3.1415279
14
例1 在我方某前沿防守地域，敌人以一个炮排（含两 门火炮）为单位对我方进行干扰和破坏．为躲避我方 打击，敌方对其阵地进行了伪装并经常变换射击地 点．
经过长期观察发现，我方指挥所对敌方目标的指 示有50％是准确的，而我方火力单位，在指示正确 时，有1/3的射击效果能毁伤敌人一门火炮，有1/6 的射击效果能全部消灭敌人．
数学建模——蒙特卡洛简介
数学建模——蒙特卡洛简介-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN数学建模——蒙特卡洛方法（案例）一、概述蒙特卡罗方法是一种计算方法。
原理是通过大量随机样本，去了解一个系统，进而得到所要计算的值。
它非常强大和灵活，又相当简单易懂，很容易实现。
对于许多问题来说，它往往是最简单的计算方法，有时甚至是唯一可行的方法。
它诞生于上个世纪40年代美国的"曼哈顿计划"，名字来源于赌城蒙特卡罗，象征概率。
二、π的计算第一个例子是，如何用蒙特卡罗方法计算圆周率π。
正方形内部有一个相切的圆，它们的面积之比是π/4。
现在，在这个正方形内部，随机产生10000个点（即10000个坐标对 (x, y)），计算它们与中心点的距离，从而判断是否落在圆的内部。
如果这些点均匀分布，那么圆内的点应该占到所有点的π/4，因此将这个比值乘以4，就是π的值。
通过R语言脚本随机模拟30000个点，π的估算值与真实值相差0.07%。
三、积分的计算上面的方法加以推广，就可以计算任意一个积分的值。
比如，计算函数 y = x2 在 [0, 1] 区间的积分，就是求出下图红色部分的面积。
这个函数在 (1,1) 点的取值为1，所以整个红色区域在一个面积为1的正方形里面。
在该正方形内部，产生大量随机点，可以计算出有多少点落在红色区域（判断条件 y < x2）。
这个比重就是所要求的积分值。
用Matlab模拟100万个随机点，结果为0.3328。
四、交通堵塞蒙特卡罗方法不仅可以用于计算，还可以用于模拟系统内部的随机运动。
下面的例子模拟单车道的交通堵塞。
根据 Nagel-Schreckenberg 模型，车辆的运动满足以下规则。
▪当前速度是 v 。
▪如果前面没车，它在下一秒的速度会提高到 v + 1 ，直到达到规定的最高限速。
▪如果前面有车，距离为d，且 d < v，那么它在下一秒的速度会降低到 d -1 。
蒙特卡洛模拟
2013年9月2日
29
蒲丰投针问题
2013年9月2日
30
蒙特卡洛与21点
▪ 大多数赌场使用6副牌或8副牌玩这种游戏，以 防止“数牌点”，在你的模拟中使用两副牌(共 104张)。只有2位参与者，你和庄家。游戏开始 时每人得到两张牌，对于牌面为2~10的牌，点 数和面数相同；对于为人脸(J、Q、K)的牌， 点数为10；牌面为A的牌，点数为1或者11.游 戏的目的是得到总数尽量接近21点的牌，不得 超过(超过称“爆了”)，并使你得到的总点数 多于庄家。
2013年9月2日
3
Monte Carlo方法的发展历史
▪ 1777年，古稀之年的蒲丰在家中请来好些 客人玩投针游戏（针长是线距之半），他事 先没有给客人讲与π有关的事。客人们虽然 不知道主人的用意，但是都参加了游戏。他 们共投针2212次，其中704次相交。蒲丰说， 2212/704=3.142，这就是π值。这着实让人 们惊喜不已。
39
误差具有概率性
▪ 由于蒙特卡罗方法的误差是在一定置信水平下 估计的，所以它的误差具有概率性，而不是一 般意义下的误差。
2013年9月2日
40
蒙特卡罗方法的主要应用范围
▪ 蒙特卡罗方法所特有的优点，使得它的应用范 围越来越广。它的主要应用范围包括：粒子输 运问题，统计物理，典型数学问题，真空技术， 激光技术以及医学，生物，探矿等方面，特别 适用于在计算机上对大型项目、新产品项目和 其他含有大量不确定因素的复杂决策系统进行 风险模拟分析。随着科学技术的发展，其应用 范围将更加广泛。
1不满足相互独立的要求2不可避免的出现重复问题因此我们将计算机产生的随机数称为伪随机数大连大学数学建模工作室产生伪随机数的方法大连大学数学建模工作室matlab中生成随机数的函数大连大学数学建模工作室matlab中生成随机数的函数exprnd指数分布的随机数生成器geornd几何分布的随机数生成器poissrnd泊松分布的随机数生成器unidrnd离散均匀分布的随机数生成器unifrnd连续均匀分布的随机数生成器betarnd贝塔分布的随机数生成器binornd二项分布的随机数生成器matlab中生成随机数的函数大连大学数学建模工作室其它函数大连大学数学建模工作室蒙特卡洛方法的实例讲解计算圆周率在平面上画一个半径r的圆和边长为2r的正方形让他们的中心重合
数学建模十大经典算法( 数学建模必备资料)
建模十大经典算法1、蒙特卡罗算法。
该算法又称随机性模拟算法，是通过计算机仿真来解决问题的算法，同时通过模拟可以来检验自己模型的正确性。
2、数据拟合、参数估计、插值等数据处理算法。
比赛中通常会遇到大量的数据需要处理，而处理数据的关键就在于这些算法，通常使用Matlab作为工具。
3、线性规划、整数规划、多元规划、二次规划等规划类问题。
建模竞赛大多数问题属于最优化问题，很多时候这些问题可以用数学规划算法来描述，通常使用Lindo、Lingo、MATLAB软件实现。
4、图论算法。
这类算法可以分为很多种，包括最短路、网络流、二分图等算法，涉及到图论的问题可以用这些方法解决，需要认真准备。
5、动态规划、回溯搜索、分治算法、分支定界等计算机算法。
这些算法是算法设计中比较常用的方法，很多场合可以用到竞赛中。
6、最优化理论的三大非经典算法：模拟退火法、神经网络、遗传算法。
这些问题是用来解决一些较困难的最优化问题的算法，对于有些问题非常有帮助，但是算法的实现比较困难，需慎重使用。
7、网格算法和穷举法。
网格算法和穷举法都是暴力搜索最优点的算法，在很多竞赛题中有应用，当重点讨论模型本身而轻视算法的时候，可以使用这种暴力方案，最好使用一些高级语言作为编程工具。
8、一些连续离散化方法。
很多问题都是实际来的，数据可以是连续的，而计算机只认的是离散的数据，因此将其离散化后进行差分代替微分、求和代替积分等思想是非常重要的。
9、数值分析算法。
如果在比赛中采用高级语言进行编程的话，那一些数值分析中常用的算法比如方程组求解、矩阵运算、函数积分等算法就需要额外编写库函数进行调用。
10、图象处理算法。
赛题中有一类问题与图形有关，即使与图形无关，论文中也应该要不乏图片的，这些图形如何展示以及如何处理就是需要解决的问题，通常使用Matlab进行处理。
历年全国数学建模试题及解法赛题解法93A非线性交调的频率设计拟合、规划93B足球队排名图论、层次分析、整数规划94A逢山开路图论、插值、动态规划94B锁具装箱问题图论、组合数学95A飞行管理问题非线性规划、线性规划95B天车与冶炼炉的作业调度动态规划、排队论、图论96A最优捕鱼策略微分方程、优化96B节水洗衣机非线性规划97A零件的参数设计非线性规划97B截断切割的最优排列随机模拟、图论98A一类投资组合问题多目标优化、非线性规划98B灾情巡视的最佳路线图论、组合优化99A自动化车床管理随机优化、计算机模拟99B钻井布局0-1规划、图论00A DNA序列分类模式识别、Fisher判别、人工神经网络00B钢管订购和运输组合优化、运输问题01A血管三维重建曲线拟合、曲面重建01B 公交车调度问题多目标规划02A车灯线光源的优化非线性规划02B彩票问题单目标决策03A SARS的传播微分方程、差分方程03B 露天矿生产的车辆安排整数规划、运输问题04A奥运会临时超市网点设计统计分析、数据处理、优化04B电力市场的输电阻塞管理数据拟合、优化05A长江水质的评价和预测预测评价、数据处理05B DVD在线租赁随机规划、整数规划06A 出版资源配置06B 艾滋病疗法的评价及疗效的预测 07A 中国人口增长预测 07B 乘公交，看奥运 多目标规划 数据处理 图论 08A 数码相机定位 08B 高等教育学费标准探讨09A 制动器试验台的控制方法分析 09B 眼科病床的合理安排 动态规划 10A 10B赛题发展的特点：1.对选手的计算机能力提出了更高的要求：赛题的解决依赖计算机，题目的数据较多，手工计算不能完成，如03B ，某些问题需要使用计算机软件，01A 。
数学建模蒙特卡罗方法
• 求积分
b
f(x)dx 〔2.1〕 a
蒙特卡罗方法步骤如下：
• 1、在区间【a,b】上利用计算机均匀产生n个随 机数x1，x2·····xn，这个可以在MATLAB软件中 用unifrnd命令实现。
• 2、计算每一个随机数相应的被积函数值f〔x1〕，
z=16*data(:,3);
II=find(x>=y.^2&x<=y+2&z<=x.*(y.^2));
M=length(II);
V=192*M/10000
例5.15 用蒙特卡罗方法计算 (x2y2z2)dxdydz
其中,积分区域是由 z x2 y和2 z = 1 所围成。
被积函数在求积区域上的最大 值为2。所以有四维超立方体
蒙特卡洛计算方法具有随机性、不确定性.
N=10000; 递推公式和初始值确定后，整个随机数序列便被唯一确定。
V=8*M/10000
实验参考程序 蒙特卡罗方法计算体积
for k=1:L y=-1+3*data(:,2);
实验参考程序 蒙特卡罗方法计算体积
P=rand(N,3); 第二条语句实现了在积分区间上均匀产生N个随机数.
什么叫蒙特卡罗方法？
• 蒙特卡罗方法又称统计模拟法、随机抽样 技术，是一种随机模拟方法，以概率和统 计理论方法为根底的一种计算方法，是使 用随机数〔或伪随机数〕来解决很多计算 问题的方法。将所求解的问题同一定的概 率模型相联系，用电子计算机实现统计模 拟或抽样，以获得问题的近似解。为象征 性地说明这一方法的概率统计特征，故借 用赌城蒙特卡罗命名。
大小有关
• 所以在使用蒙特卡罗方法时，要 “扬长避短〞，只对问题中难以用解 析〔或数值〕方法处理的局部，使用
Monte-Carlo(蒙特卡洛方法)解析
常用的线性同余生成器
Modulus m 2^31-1
=2147483647
2147483399 2147483563
Multiplier a 16807
在 n 次中出现的频率。假如我们取 fn ( A) 作为 p P(A) 的估计，即 pˆ fn ( A) 。
然后取 ˆ
2l afn ( A)
作为
的估计。根据大数定律，当 n 时，
pˆ
fn ( A) a.s.
p.
从而有ˆ 2l P 。这样可以用随机试验的方法求得 的估计。历史上 afn ( A)
(2) 计算 X F -1(U ) ,则 X 为来自 F(x) 分布的随机数.
例 1 ：设 X ~ U (a,b) ，则其分布函数为
0
F
(
x)
x b
a a
1,
xa a xb
xb
F -1( y) a (b a) y ， 0 y 1
生成 U (0,1) 随机数 U，则 a (b - a)U 是来自
算法实现
许多程序语言中都自带生成随机数的方法, 如 c 中的 random() 函数, Matlab中的rand()函数等。 但这些生成器生成的随机数效果很不一样, 比如 c 中的函数生成的随机数性质就比较差, 如果用 c , 最好自己再编一个程序。Matlab 中的 rand() 函数, 经过了很多优化。可以产生性质很好的随 机数, 可以直接利用。
U (a,b) 的随机数。
例 2：
设 X ~ exp( ) 服从指数分布，则 X 的分布函数为：
数学建模蒙特卡罗模拟方法课件2
四
现值系数（i=10%)
五
净现值（税后）
累计净现值（税后）
合计
45306 45306 41353 26583 2514 2292 9964 3953
915
建设经营期
2005 0 0
16277 16277
0 0 0 -16277 -16277 1 -16277 -16277
2006 18064 18064 12329 8502 1003
表5-2 乙方案的主要经济技术指标
序号
项目
一
现金流入
1
销售收入
二
现金流出
1
开发建设投资
2
营业税金及附加
3
土地增值税
4
所得税
三
净现金流量（税后）
累计净现金流量（税后）
四
现值系数（i=10%)
五
净现值（税后）
累计净现值（税后）
合计
54660 54660 49215 30626 3034 4190 11365 5445
（二）、确定每个风险变量 的概率分布
同样采用“德尔菲法”估计 出以上 7 个风险变量概率分布 和其分布函数中的具体参数,如 下表所示：
住宅类销售收入 商业类销售收入 土地费用 前期费用 开发建设费用 营销费用 其他费用
住宅类销售收入 商业类销售收入 土地费用 前期费用 开发建设费用 营销费用 其他费用
三角分布
a：4007 b：4555 m：4218
三角分布
a：258 b：413 m：368
正态分布
均值：265 方差：30
第二年
分布
参数
三角分布
a：3996 b：5328 m：4440
数学建模算法之蒙特卡罗方法——原理、编程及应用
数学建模算法之蒙特卡罗方法——原理、编程及应用一、前言1946年，美国拉斯阿莫斯国家实验室的三位科学家John von Neumann,Stan Ulam和Nick Metropolis共同发明了蒙特卡罗方法。
此算法被评为20世纪最伟大的十大算法之一。
蒙特卡罗方法（Monte Carlo method），又称随机抽样或统计模拟方法，是一种以概率统计理论为指导的一类非常重要的数值计算方法。
此方法使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。
由于传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。
二、蒙特卡罗方法的基本原理以及思想1、蒲丰投针实验其基本思想源于法国数学家蒲丰提出著名的蒲丰投针实验，并以该方法求圆周率。
为了求得圆周率π值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（l＜a）的平行线相交的频率代替概率P，再利用准确的关系式：求出π值。
其中Ｎ为投针次数，n为针与平行线相交次数。
这就是古典概率论中著名的蒲丰氏问题。
2、射击问题设r表示射击运动员的弹着点到靶心的距离，ｇ(r)表示击中r处相应的得分数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射击水平。
则该运动员的射击成绩为用概率语言来说，<g>是随机变量ｇ(r)的数学期望，即当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。
有一个例子可以使你比较直观地了解蒙特卡洛方法：假设我们要计算一个不规则图形的面积，那么图形的不规则程度和分析性计算（比如，积分）的复杂程度是成正比的。
蒙特卡洛方法是怎么计算的呢？假想你有一袋豆子，把豆子均匀地朝这个图形上撒，然后数这个图形之中有多少颗豆子，这个豆子的数目就是图形的面积。
蒙特卡罗法在计算机仿真中的应用研究
本科毕业论文（设计、创作）题目：蒙特卡罗法在计算机仿真中的应用研究学生姓名：学号：0321002020所在院系：信息与通信技术系专业：电子信息工程入学时间：2010 年9 月导师姓名：傅有亮//朱亮职称/学位：副教授/硕士//讲师/硕士导师所在单位：完成时间：2014 年 5 月安徽三联学院教务处制蒙特卡罗法在计算机仿真中的应用研究摘要:在运用蒙特卡罗法计算求解问题的过程中会遇到一系列的问题：比如如何构造或描述概率过程、并且如何从已知概率分布抽样和建立估计量。
其中，构造或描述概率过程实际上就是建立随机试验模型，构造概率过程是对确定性的问题而言的，描述概率过程是对随机性的问题而言的，不同的问题所需要建立的随机试验模型各不相同。
此问题将是本论文的重点之所在。
所谓的从已知概率分布抽样指的是随机试验过程，随机模拟中必要包含某些已知概率分布的随机变量或随机过程作为输入，进行随机试验过程就是对随机变量的样本或随机过程的样本函数作为输入相应的输出过程，因此通常被称之为对已知概率分布的抽样。
如何产生已知分布的随机变量或随机过程是蒙特卡罗法中的一个关键问题，亦是本论文的关键。
总之，本论文所要阐述的主要问题包括如何产生随机数，如何描述概率过程以及如何使用计算机C语言程序来对蒙特卡罗法进行仿真。
关键词：蒙特卡罗法；仿真；概率；随机数；定积分The Research of the Monte carlo method in theapplication of computer simulationAbstract:Process calculation problem in using the Monte Carlo method will encounter a series of problems: such as how to structure or a probabilistic description of process, and from the known probability distribution of sampling and estimation. Among them, structure or describing the probability is actually a process of a random test model, construct probabilistic process is to the deterministic problem,describing the probability of random process is the problem in terms of the model of random test, different problems need to establish each are not identical. This problem will be the key point of the paper. The so-called from the known probability distribution of sample is a random process, the necessary simulation contains some known probability distribution of random variables or random process as input, the sample function of the random testing is the process of the samples of random variables or random process as the input and output process, it is often referred to as the known probability the sampling distribution. How to produce a known distribution of the random variables or random process is a key problem of Monte Carlo method, also is the key of this paper. Always, the main issues in this paper to set including how to generate random numbers, how to describe the probability process and how to use the computer Clanguage program to simulate the Monte Carlo method.Keywords:Monte Carlo method；simulation；probability ；random number；definite integral目录第一章绪论 (1)1.1 研究背景 (1)1.2 研究现状分析 (1)1.3 研究思路和方法 (2)第二章计算机仿真 (3)2.1 计算机仿真技术的概述 (3)2.2 计算机仿真技术的发展 (3)2.3 计算机仿真技术的发展现状及前景 (3)第三章定积分及其应用 (6)3.1 定积分的概念 (6)3.2 定积分的基本计算方法 (6)第四章蒙特卡罗法 (10)4.1 蒙特卡罗法的来源和概述 (10)4.2 概率模型和蒙特卡罗法 (11)4.3 基于蒙特卡罗法的定积分计算 (13)结语 (15)致谢 (16)参考文献 (17)第一章绪论1.1 研究背景蒙特卡罗方法在科学上又称统计模拟法、随机抽样技术，是随机模拟方法的一种，它的理论基础是以概率论和统计理论方法为前提的，或者说是通过使用随机数来对一些问题进行求解。
数学建模十大经典算法之蒙特卡罗方法讲义课件
f(X)是X的分布密度函数。则
lim P N
N
X N E( X ) x
1 x et2 / 2dt
2 x
当N充分大时，有如下的近似式
P
XN
E(X )
N
2 et2 / 2dt 1 2 0
其中α称为置信度，1－α称为置信水平。
这表明，不等式
XN
E(X )
N
近似地以概率
1－α成立，且误差收敛速度的阶为 O(N 1/ 2 )。
5) 误差容易确定。
6) 程序结构简单，易于实现。
1) 能够比较逼真地描述具有随机性质 的事物的特点及物理实验过程
从这个意义上讲，蒙特卡罗方法可以部分代替物 理实验，甚至可以得到物理实验难以得到的结果。用 蒙特卡罗方法解决实际问题，可以直接从实际问题本 身出发，而不从方程或数学表达式出发。它有直观、 形象的特点。
1901 3408
3.1415929
例2. 射击问题（打靶游戏）
设r表示射击运动员的弹着点到靶心的距离，ｇ(r)
表示击中r处相应的得分数（环数），f(r)为该运动员的 弹着点的分布密度函数，它反映运动员的射击水平。 该运动员的射击成绩为
g 0 g(r) f (r)dr
用概率语言来说，<g>是随机变量ｇ(r)的数学期
为数学期望<g>的估计值（积分近似值）。
➢ 基本思想
由以上两个例子可以看出，当所求问题的解是某 个事件的概率，或者是某个随机变量的数学期望，或 者是与概率、数学期望有关的量时，通过某种试验的 方法，得出该事件发生的频率，或者该随机变量若干 个具体观察值的算术平均值，通过它得到问题的解。 这就是蒙特卡罗方法的基本思想。
➢ 作业
蒙特卡罗法(Monte Carlo method)
实验十二计算机仿真实验实验目的：1. 掌握全概率公式与贝叶斯公式；2. 了解计算机仿真方法；3. 了解蒙特卡罗法(Monte Carlo method), 具有初级编程能力.实验原理：全概率公式: 设A 1, A 2, …, A n 为两两互斥事件，B 是A 1 + A 2 + … + A n 的子事件，则P(B )=P(A 1)P(B |A 1) + P(A 2)P(B |A 2) + … + P(A n )P(B |A n ).贝叶斯公式: P(A k |B )= P(A k )P(B | A k )/P(B ).计算机仿真: 就是在计算机上模拟各种实际系统的运行过程. 计算机仿真通常用来产生规定分布的随机变量.对于任意随机变量ξ ，其分布函数为F (x )，设 η = F (ξ )的分布函数为G (y )，则G (y ) = P{η ≤y }= P{F (ξ ) ≤y }= P{ξ ≤F -1( y )}= y ，这说明η 服从[0，1]的均匀分布.一般的编程语言都提供了均匀分布随机数发生器.应用随机数模拟试验的方法通常称为蒙特卡罗法(Monte Carlo method). 蒙特卡罗法不仅适用于处理随机性问题, 如存贮、排队、质量检验、市场营销、社会救急、生态竞争和传染病等问题；也可处理定性问题, 如计算多重积分、解积分方程及微分方程、解整数规划(特别是非线性整数规划)等.应用蒙特卡罗法解规划问题的基本思想是：先估计各个变量的大致取值范围，每次试验从中随机取出一个样本点，然后判断它是否为可行点. 若是则将其目标函数值与上一次的目标函数值相比较，记录下较优目标函数值与其样本点；否则重新抽样。
直到试验次数达到指定值或可行点数达到指定值为止.实验内容：1. 设有两个口袋，甲袋中盛有两个白球，一个黑球，乙袋中盛有一个白球，两个黑球.由甲袋任取一个球放入乙袋，再从乙袋中取出一个球.若从乙袋中取出的球是白球，那么从甲袋中取出放入乙袋的球是白球还是黑球? 用计算机模拟上述过程1000次，问理论判断是正确的有多少次?2. 用计算机模拟随机变量ξ ～ϕ (x ) =⎩⎨⎧≤>-0,0,0,e 55x x 的取值200次. 3. 用计算机模拟随机变量ξ ～N (120,102 )的取值800次，并画出统计直方图.4. 应用蒙特卡罗法解非线性规划问题:max z = - 2x 2 - y 2 + xy + 8x + 3ys.t. 3 x + y = 10x ≥0, y ≥0.C 语言简介(仅介绍本实验所用到的)1.标识符标识符是由程序员定义的单词，如函数名、变量名等. 标识符是由大小写字母、数字和下划线组成的，并以字母和下划线开始.2.关键字void (无值型) char (字符型) int (整数型) long (长整数型) float (浮点型) double (双精度浮点型)if (如果) else (或者)for (循环) while (循环) break (满足一定的条件终止循环)return (返回函数值)3.函数形式类型函数名(参数){}4.库函数简介double sin(double x) double cos(double x) double exp(double x) double log(double x) double sqrt(double x) double pow(double x, double y)double fabs(double x) int abs(int x)int printf(const char *format, ...); 屏幕格式化输出函数FILE *fopen(const char *filename, const char *mode); 文件打开int fprintf(FILE *fp, const char *format, ...); 文件格式化输出函数int fgetc(FILE *fp); 从文件中读出一字符int fclose(FILE *fp); 文件关闭void far setcolor(int color); 设置输出颜色void far line(int x1, int y1, int x2, int y2); 画直线void far rectangle(int left, int top, int right, int bottom); 画矩形unsigned far getpixel(int x, int y); 读出点(x, y)的颜色void far putpixel(int x, int y, int pixelcolor); 画点int random(int Num), 均匀产生0到Num-1中的一个随机数5.示例计算9！#include<stdio.h>void main(void){int i;long n;n=1;for(i=1;i<=9;i++)n*=i;printf("\n9!=%ld\n",n);getch();}6. Turbo C(2.0) 编辑命令F3 录入文件F9 编译Ctrl+F9 运行Ctrl+KB 定义块首Ctrl+KK 定义块尾Ctrl+KC 块粘贴Ctrl+KV 块移动Ctrl+Y 删除当前行7.部分源程序程序LAB1_1.C 求出方程sin x- x = 1在( -2，2)内的近似根#include<stdio.h>#include<math.h>float f(float x){ return sin(x)-x-1; }void main(void){float r=0.618,x0=-2,x1=2,x;int n=0;while(1){n++;x=(1-r)*x0+r*x1;if(f(x)*f(x0)<0.0)x1=x;else if(f(x)*f(x1)<0.0)x0=x;else break;if(fabs(x1-x0)<0.001)break;}printf("n=%d, x=%f\n",n,x);getch();}程序LAB2_3.C 给出正态分布函数表#include<stdio.h>#include<math.h>float f(float x){ return exp(-x*x/2); }void main(void){float x,x0,x1=0.0,F=0,h=0.0001;long n=0;FILE *fp;int i=0,p;fp=fopen("x.c","w");for(x1=0.0;x1<0.04;x1+=0.001){n=0; F=0; x0=-10;for(x=x0+h;x<x1;x+=h){if(n%2)F+=2*f(x);else F+=4*f(x);n++;}F+=f(x0)+f(x1);F/=3; F*=h; F*=0.39894;printf("x=%5.3f, F=%6.4f\n",x1,F);i++;p=F*10000;fprintf(fp,"%d,",p);if(i==10){i=0;fprintf(fp,"\n ");}if(p==9999)break;}fclose(fp);getch();}程序LAB3_2_1.C 解下列微分方程y - y tan x = sec x ，y (0) = 0, 并画出其图形: #include <stdio.h>#include <math.h>#include <graphics.h>float f(float x,float y){ return y*sin(x)/cos(x)+1/cos(x); }void main(void){int i=DETECT,j;float x=0.0,y=0.0,h=0.005;char *str="0.00";initgraph(&i,&j," ");setviewport(0,0,639,479,1);cleardevice();setbkcolor(BLUE);setcolor(WHITE);line(20,200,620,200);for(i=0;i<10;i++){line(20+i*60,195,20+i*60,200);str[0]=48+3*i/10;str[2]=48+3*i%10;outtextxy(20+i*60,205,str);}for(i=0;i<600;i++){y=y+f(x,y)*h;x+=h;putpixel(i+20,200-y*10,GREEN);}getch();closegraph();}程序LAB3_2_2.C 解下列微分方程组⎪⎪⎩⎪⎪⎨⎧+-=-=.01.0d d ,25.02d d xy y t y xy x t x x (0) = 100, y (0) = 8, 并画出其图形:#include <stdio.h>#include <math.h>#include <graphics.h>float f1(float x,float y){ return 2*x-0.25*x*y; }float f2(float x,float y){ return -y+0.01*x*y; }void main(void){int i=DETECT,j;float x=99.0,y=7.9,h=0.015;initgraph(&i,&j," ");setviewport(0,0,639,479,1);cleardevice();for(i=0;i<600;i++){x=x+f1(x,y)*h;y=y+f2(x,y)*h;putpixel(i+20,(102.5-x)*60,GREEN);putpixel(i+20,(8.6-y)*600,WHITE);}getch();closegraph();}程序LAB4.C 将矩阵化为行阶梯型,化为行最简型#include<stdio.h>#include<math.h>#define MAXR 20#define MAXC 40/*解线性方程组,以下是其增广矩阵*/float M_B[MAXR][MAXC]={{1,-2,2,1,-3},{2,1,1,-2,-1},{3,4,0,-5,1}}; void f1(int m,int n){int i,j,r=0,c=0;float x0;printf("\n以下将矩阵化为行阶梯型\n");for(i=0;i<m;i++){for(j=0;j<n;j++)printf("%8.2f",M_B[i][j]);printf("\n");}printf("按任一健继续...\n");getch();while(c<n){for(i=r;i<m;i++)if(fabs(M_B[i][c])>=0.0001)break;if(i<m){if(i!=r)for(j=0;j<n;j++){x0=M_B[i][j];M_B[i][j]=M_B[r][j];M_B[r][j]=x0;}x0=M_B[r][c];for(j=0;j<n;j++)M_B[r][j]/=x0;for(i=r+1;i<m;i++){x0=M_B[i][c];for(j=0;j<n;j++)M_B[i][j]-=x0*M_B[r][j];}r++;for(i=0;i<m;i++){for(j=0;j<n;j++)printf("%8.2f",M_B[i][j]);printf("\n");}printf("按任一健继续...\n");getch();}c++;if(r==m)break;}printf("\n以下将行阶梯型化为行最简型\n");while(r){r--;for(j=0;j<n-1;j++)if(fabs(M_B[r][j])>=0.0001)break; c=j;for(i=0;i<r;i++){x0=M_B[i][c];for(j=0;j<n;j++)M_B[i][j]-=x0*M_B[r][j];}for(i=0;i<m;i++){for(j=0;j<n;j++)printf("%8.2f",M_B[i][j]);printf("\n");}printf("按任一健继续...\n");getch();}printf("完毕,按任一健退出...\n");getch();}void f2(int n){int i,j;for(i=0;i<n;i++)for(j=0;j<n;j++){if(i!=j)M_B[i][n+j]=0;else M_B[i][n+j]=1;}f1(n,2*n);}void main(void){ f1(3,5); }程序LAB6_3.C 模拟产生服从N(120,400)分布的随机变量800次,并画出统计直方图#include <stdio.h>#include <stdlib.h>#include <math.h>#include <graphics.h>int Np[4000]={/*下列数据为标准正态分布函数值×10000, 由程序LAB2_3.C 产生*/ 5000,5003,5007,5011,5015,5019,5023,5027,5031,5035,…};int f(int p){int i;for(i=0;i<4000;i++)if(Np[i]>=p)break;return i;}void main(void){int i,j,x;float F,X[800];randomize();printf("\n\n");for(i=0;i<800;i++){x=random(10000);if(x<5000)F=-f(10000-x);else F=f(x);F*=0.001;X[i]=120+20*F;printf("%8.2f",X[i]);}getch();i=DETECT;initgraph(&i,&j," ");setviewport(0,0,639,479,1);cleardevice();setbkcolor(BLUE);setcolor(WHITE);line(20,400,620,400);j=20;for(F=90;F<150;F+=3){x=0;for(i=0;i<800;i++)if(X[i]>=F&&X[i]<F+5)x+=4;rectangle(j,400-x,j+30,400);j+=30;}getch();closegraph();}程序2000A01.C 统计文件2000A1.txt中A TCG的个数#include<stdio.h>void main(void){int n=-1,ATCG[40][4]={0,0};char c;FILE *fp;fp=fopen("2000A1.txt","r");while((c=fgetc(fp))!=EOF){if(c=='.')n;else if(c=='a')ATCG[n][0];else if(c=='t')ATCG[n][1];else if(c=='c')ATCG[n][2];else if(c=='g')ATCG[n][3];}fclose(fp);for(n=0;n<40;n)printf("{%d,%d,%d,%d},\n",ATCG[n][0], ATCG[n][1],ATCG[n][2],ATCG[n][3]);getch();}。
数学建模-蒙特卡罗方法
蒙特卡罗方法的特点 • 优点：
• 1、能够比较逼真地描述具有随机性质的事物的 特点及物理实验过程 • 2、受几何条件限制小 • 3、收敛速度与问题的维数无关 • 4、具有同时计算多个方案与多个未知量的能力 • 5、误差容易确定 • 6、程序结构简单，易于实现
• 缺点：
• 1收敛速度慢 • 2误差具有概率性 • 3在粒子输运问题中，计算结果与系统大小有关
什么是随机数？
• 在连续型随机变量的分布中，最简单而且最基本的 分布是单位均匀分布。由该分布抽取的简单子样称 为随机数序列，其中每一个体称为随机数 • 符号： 1 2
，
• 两个特点：独立性，均匀性
产生随机数 • 随机数表方法 • 物理方法
随机数表
• 随机数表是由0,1,2,3,4,5,6,7,8,9十个数字 组成，每个数字以0.1的等概率出现，数字 之间相互独立，这些数字序列叫作随机数 字序列。 • （如果要得到n位有效数字的随机数，只需 将表中每n个相邻的随机数字合并在一起， 且在最高位的前边加上小数点即可。例如， 某随机数表的第一行数字为7 6 3 4 2 5 8 9 1...，要想得到三位有效数字的随机数一次 为0.763，0.425，0.891...）
基本思想
• 当所求问题的解是某个事件的概率，或者 是某个随机变量的数学期望，或者是与概 率，数学期望有关的量时，通过某种试验 的方法，得出该事件发生的概率，或者该 随机变量若干个具体观察值的算术平均值， 通过它得到问题的解。 • 当随机变量的取值仅为 1 或 0 时，它的数学 期望就是某个事件的概率。或者说，某种 事件的概率也是随机变量（仅取值为 1 或 0 ） 的数学期望。
蒙特卡罗方法实验
面积、体积计算问题
数学建模_第3讲蒙特卡洛方法初探(2007)
几何概率模拟
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
2. 符号假设
i：要模拟的打击次数； k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数； k3：击中敌人两门火炮的射击总数． E：有效射击(毁伤一门炮或两门炮)的概率
3. 模拟框图
初始化:i=iti2(0.4,10000)
古典概率模拟
例3: 在一袋中有10 个相同的球，分别标有号码 1,2,…,10。每次任取一个球，记录其号码后放 回袋中，再任取下一个。这种取法叫做“有放 回抽取”。今有放回抽取3个球，求这3个球的 号码均为偶数的概率。（用频率估计概率）
function proguji=liti5(mm) %mm 是随机实验次数 frq=0; randnum1=unifrnd(0,60,mm,1); randnum2=unifrnd(0,60,mm,1); randnum=randnum1-randnum2; proguji=0; for ii=1:mm if abs(randnum(ii,1))<=20 frq=frq+1; end end proguji=frq/mm
蒙特卡罗方法
其中 >0为常数，则称X服从参数为 的指数 分布。 •指数分布的期望值为
•排队服务系统中顾客到达间隔、质量与可靠性 中电子元件的寿命通常服从指数分布。
例 顾客到达某商店的间隔时间服从参数为 10(分钟)的指数分布(指数分布的均值为10)
指两个顾客到达商店的平均间隔时间 是10分钟.即平均10分钟到达1个顾客. 顾客到 达的间隔时间可用exprnd(10)模拟。
2．产生mm*nn阶离散均匀分布的随机数 矩阵： R = unidrnd(N) R = unidrnd(N,mm,nn)
3.产生 m n 阶均值为 ， 标准差为 的正态分布的随机数矩阵： normrnd ( , ,m, n)
产生一个均值为 ，标准差的正态分布的随机数： normrnd ( , )
•当研究对象视为大量相互独立的随机变量之和， 且其中每一种变量对总和的影响都很小时，可以 认为该对象服从正态分布。
4．产生 m n 阶期望值为的指数分布的随机数矩阵： exprnd (,m, n )
•若连续型随机变量X的概率密度函数为
1 x / e f ( x ) 0 x0 x0
proguji = 0.5557
复杂概率模拟 例6 在我方某前沿防守地域，敌人以一个炮排 （含两门火炮）为单位对我方进行干扰和破 坏．为躲避我方打击，敌方对其阵地进行了伪装 并经常变换射击地点． 经过长期观察发现，我方指挥所对敌方目标的 指示有50％是准确的，而我方火力单位，在指示 正确时，有1/3的概率能毁伤敌人一门火炮，有 1/6的概率能全部消灭敌人． 现在希望能用某种方式把我方将要对敌人实施 的1次打击结果显现出来，利用频率稳定性，确定 有效射击(毁伤一门炮或全部消灭)的概率.