matlab程序问题。需要用到蒙特卡洛方法

matlab程序问题。需要用到蒙特卡洛方法,第1张

你给出的解是正确的

首先假设有编号为1~16的16个球,其中

编号1~8,8个球是红色,那么9~16,8个球是白色

n=1e6; %游戏100万次

A=0;B=0;C=0;D=0;E=0; %得奖统计清零

for i=1:n

examp=randperm(16); %随机打乱1~16,16个自然数

num=sum(examp(1:8)<=8); %examp(1:8)取出前8个,就是从16个球中随机抽8个

%sum(examp(1:8)<=8),就是统计编号小于8的球的数量

%也就是红球的个数

if num==0||num==8

A=A+1; %如果8个都是红球,或者没有红球,A类统计加以

elseif num==1||num==7

B=B+1; %以下的判断依次类推

elseif num==2||num==6

C=C+1;

elseif num==3||num==5

D=D+1;

else

E=E+1;

end

end %100万次抽球后,A,B,C,D,E的次数都一一统计记录

t=10A/n+1B/n+05C/n+02D/n-3E/n

%A/n就是得到A奖的概率,以此类推

%用每个将的奖金乘以奖的概率,再相加,得到了奖金的期望

%结果表明,按照概率统计,平局每局要亏掉差不多1块钱

蒙特卡罗法也称统计模拟法、统计试验法。是把概率现象作为研究对象的数值模拟方法。是按抽样调查法求取统计值来推定未知特性量的计算方法。蒙特卡罗是摩纳哥的著名赌城,该法为表明其随机抽样的本质而命名。故适用于对离散系统进行计算仿真试验。在计算仿真中,通过构造一个和系统性能相近似的概率模型,并在数字计算机上进行随机试验,可以模拟系统的随机特性

蒙特卡洛法(又称统计试验法)是描述装备运用过程中各种随机现象的基本方法,而且它特别适用于一些解析法难以求解甚至不可能求解的问题,因而在装备效能评估中具有重要地位。

用蒙特卡洛法来描述装备运用过程是1950年美国人约翰逊首先提出的。这种方法能充分体现随机因素对装备运用过程的影响和作用。更确切地反映运用活动的动态过程。在装备效能评估中,常用蒙特卡洛法来确定含有随机因素的效率指标,如发现概率、命中概率、平均毁伤目标数等;模拟随机服务系统中的随机现象并计算其数字特征;对一些复杂的装备运用行动,通过合理的分解,将其简化成一系列前后相连的事件,再对每一事件用随机抽样方法进行模拟,最后达到模拟装备运用活动或运用过程的目的。[2]

基本思路

蒙特卡洛法的基本思想是:为了求解问题,首先建立一个概率模型或随机过程,使它的参数或数字特征等于问题的解:然后通过对模型或过程的观察或抽样试验来计算这些参数或数字特征,最后给出所求解的近似值。解的精确度用估计值的标准误差来表示。蒙特卡洛法的主要理论基础是概率统计理论,主要手段是随机抽样、统计试验。用蒙特卡洛法求解实际问题的基本步骤为:

(1)根据实际问题的特点.构造简单而又便于实现的概率统计模型.使所求的解恰好是所求问题的概率分布或数学期望;

(2)给出模型中各种不同分布随机变量的抽样方法;

(3)统计处理模拟结果,给出问题解的统计估计值和精度估计值。[2]

优缺点

蒙特卡罗法的最大优点是:

1方法的误差与问题的维数无关。

2对于具有统计性质问题可以直接进行解决。

3对于连续性的问题不必进行离散化处理

蒙特卡罗法的缺点则是:

1对于确定性问题需要转化成随机性问题。

2误差是概率误差。

3通常需要较多的计算步数N

蒙特卡罗法作为一种计算方法,是由美国数学家乌拉姆(Ulam , S M)与美籍匈牙利数学家冯·诺伊曼(von Neumann,J)在20世纪40年代中叶,为研制核武器的需要而首先提出来的实际上,该方法的基本思想早就被统计学家所采用了例如,早在17世纪,人们就知道了依频数决定概率的方法。

步骤

蒙特卡洛法是一种用来模拟随机现象的数学方法,这种方法在作战模拟中能直接反映作战过程中的随机性。在作战模拟中能用解析法解决的问题虽然越来越多,但有些情况下却只能采用蒙特卡洛法。使用蒙特卡洛法的基本步骤如下:

(1)根据作战过程的特点构造模拟模型;

(2)确定所需要的各项基础数据;

(3)使用可提高模拟精度和收敛速度的方法;

(4)估计模拟次数;

(5)编制程序并在计算机上运行;

(6)统计处理数据,给出问题的模拟结果及其精度估计。

在蒙特卡洛法中,对同一个问题或现象可采用多种不同的模拟方法,它们有好有差,精度有高有低,计算量有大有小,收敛速度有快有慢,在方法的选择上有一定的技巧。[3]

应用举例

在我方某前沿防守地域,敌人以1个炮兵排(含两门火炮)为单位对我方进行干扰和破坏。为躲避我方打击,敌方对其指挥所进行了伪装并经常变换射击地点。经过长期观察发现,我方指挥所对敌方目标的指示有50%是准确的,而我方火力单位在指示正确时,有1/3的射击效果能毁伤敌人1门火炮,有1/6的射击效果能全部消灭敌人。

解:希望能用某种方法把我方将要对敌人实施的20次打击结果显示出来,确定有效射击的比率及毁伤敌方火炮的平均值。这是一个概率问题,可以通过理论计算得到相应的概率和期望值。但这样只能给出作战行动的最终静态结果,而显示不出作战行动的动态过程。

为了显示我方20次射击的过程,必须用某种方式模拟出以下两件事:一是观察所对目标的指示正确或不正确;二是当指示正确时,我方火力单位的射击结果。对第一件事进行模拟试验时有两种结果,每一种结果出现的概率都是1/2。因此,可用投掷1枚硬币的方式予以确定。当硬币出现正面时为指示正确,反之为不正确。对第二件事进行模拟试验时有3种结果,毁伤1门火炮的可能为1/3,毁伤2门火炮的可能为1/6,没能毁伤敌火炮的可能为1/2。这时,可用投掷骰子的办法来确定,如果出现的是1、2、3三个点则认为没能击中敌人,如果出现的是4、5点则认为毁伤敌1门火炮,如果出现6点则认为毁伤敌2门火炮。

通过上面的方式,就可把我方20次射击的过程动态地显现出来。

] - 排队论中的一个仿真程序,主要是用于仿真M/M/1、M/D/1模型。输入排队模型相关参量,返回计算结果。[SYschedulerar

] - 用仿真模型可以对企业生产、销售进行决策,从而采用最合理的生产或销售方案。本课题论文运用蒙特卡洛仿真方法对SY公

司生产订单排产系数进行仿真分析,探讨是否应进行改进。[

3rar] - 蒙特卡洛资料及算法···········[bank_tellerrar

] - c++写的实现数学排队论的算法程序,上研究生同志们用得上,不过我也是找了个model来稿的,算法很难,喜欢编程的大虾们都应该对算法感兴趣把[mcrar] - MATLAB编写的蒲丰氏问题的蒙特卡洛仿真[queurerar] - 针对M/PH/1(k)排队系统推导出该排队系统在任意时刻、到达时刻、退去时刻的队列长度状态概率分布、平均队列长度、以及平均等待时间,并编写程序进行数值计算,并对M/E2/1(k), M/M/1(k), M/H2/1(k)排队系统的性能进行数值比较。

贴一个蒙特卡洛方法的matlab程序,供大家使用。 {3 x& K/ i1 i( D8 C0 c$ O % Example Monte Carlo Simulation in Matlab 0 O5 \; P" t# t7 v8 c& @ % Function: y = x2^2/x1 5 Z0 W4 e9 q, d5 B+ c % % Generate n samples from a normal distribution 4 s! c6 y, I6 H" d) K+ v Y; X: Q % r = ( randn(n,1) sd ) + mu 4 U F Q) t, T# q w/ K' Q % mu : mean / E( P8 U" c o! G8 s/ x % sd : standard deviation % % Generate n samples from a uniform distribution 2 u# ^& K [0 z% F) @1 y % r = a + rand(n,1) (b-a) - D+ }& U$ w- M9 @& Q9 W, Z % a : minimum % b : maximum n = 100000; % The number of function evaluations 7 x5 a" @- F& O- Z; w5 j % --- Generate vectors of random inputs ! K& x0 ^# X+ q( V6 { % x1 ~ Normal distribution N(mean=100,sd=5) % x2 ~ Uniform distribution U(a=5,b=15) x1 = ( randn(n,1) 5 ) + 100; 2 B' l3 n) V) D$ ~ x2 = 5 + rand(n,1) ( 15 - 5 ); \: O: Y( w3 [9 d: V4 r( k4 { % --- Run the simulation % Note the use of element-wise multiplication - ~% x$ `7 A6 v9 R F y = x2^2 / x1; ' g$ O7 U; R F% ` % --- Create a histogram of the results (50 bins) hist(y,50); / M9 m+ s( [ w" J2 I% s/ X % --- Calculate summary statistics y_mean = mean(y) y_std = std(y) ; R7 A2 y M/ T" p, h m y_median = median(y)

下面是在Excel中模拟一只股票价格的例子。假设股票价格

的对数收益率服从正态分布,均值为0,每日变动标准差为01,

模拟股票价格1年的路径,过程如下:

用到两个内置函数,即用rand()来产生0到1之间的随机数,然后用norminv()来获得服从既定分布的随机数,即收益率样本=norminv(rand(), 0, 01)。假定股票价格的初始值是100元,那么模拟的价格就是 S=100 exp(cumsum(收益率样本))。

其中的cumsum()不是Excel的内置函数,其意思就是收益率样本的累积,每个时刻的值都是当前样本及此前所有样本的和,如,收益率样本从单元格C3开始,当前计算C15对应的模拟价格,则模拟价格计算公式是:100 exp(sum($C$3:C15))。

由此可以得到股票价格的一条模拟路径。

其他非正态分布也可以通过类似方式得到分布的抽样,即分布函数的逆函数,这些函数Excel都内置了。所以,做蒙特卡洛模拟的时候,关键是先确定所需模拟的分布,然后进行抽样,然后应用层面的各种公式就可以在抽样的基础上进行计算了。

--------以下是补充的--------

根据上面提到的思路,其实可以很便捷地为期权做定价。下面就用蒙特卡洛方法为一个普通的欧式看涨期权定价(蒙特卡洛在为普通期权plain vanilla option定价时不占优势,因为相对于解析法而言计算量很大。但是,如果要给结构比较复杂的奇异期权定价时,可能蒙特卡洛法就比较实用,有时可能成为唯一的方法)。

1)假设这个期权是欧式看涨期权,行权价格为50元,标的股票当前的价格也是50元,期权剩余时间是1天。

2)假设标的股票的价格服从对数正态分布,即股票的每日收益率服从正态分布,均值为0,每日标准差为1%。

根据分布假设,首先用rand()函数产生在0到1之间的均匀分布样本。为了提高精确度,这里抽样的数量为1000个(其实1000个是很少的了,通常需要10万个甚至50万个,但是在Excel表格中 *** 作这么多数字,不方便,这是Excel的不足之处)。

下一步,用norminv(probability, mean, std)函数来获得股票收益率分布的1000个抽样,其中的probability参数由rand()产生的抽样逐个代入,mean=00, std = 001。注意这里抽样得到的日度收益率。也就是说,这个样本对应的下一个交易日股票价格的收益率分布。

下一步,股票价格=50×exp(收益率样本),得到股票价格分布的抽样,有1000个样本。

根据我做的实验,这1000个样本的分布图形(histogram)跟对数正态分布是比较接近的,如下图所示:

图的横轴是股票价格,纵轴是样本中出现的频率。

得到了股票价格未来一天分布的样本之后,就可以以此样本来计算期权的价格了。

欧式看涨期权的定义为:

C=max(S-K,0)

所以,根据这个计算公式可以计算出在到期那天在特定的价格下期权的价值。在Excel中,相当于 期权价值=max(股票价格样本 - 50,0)。由此就可以得到了该期权未来1天价值的样本。

然后,将未来价值贴现回来(用无风险利率贴现,假设无风险利率为005,则贴现公式是=exp(-005/360)×期权价值,得到期权价格的1000个样本。

最后,对期权价格的1000个样本求平均,Excel函数average(期权价格样本),就可以得到期权的价格了。

我这里算出来的是:02015元。

而根据Black-Scholes期权定价公式算出来的理论价格则是02103元。二者比较接近,但是还是有差距。

而且,每次刷新Excel表格,就重新做一次模拟,得到的模拟价格变动比较大,有时是02043元,有时是01989元。由于这个抽样的数量比较小(1000个样本),所以估算的结果受到样本的影响会比较大。如果把抽样数量提高100倍甚至500倍,那么样本变动的影响可能会小一个或者两个数量级。但是计算量就大了,如果计算机性能不够高,那么利用Excel来做的话,比较困难。

这就是我的工作台:

------ 再来一个 --------

看到有人提到利用蒙特卡洛方法来估计圆周率Pi,挺有意思,也简单,所以就在Excel中做了一个实验。

基本原理在于在直角坐标系中的第一个象限中的一个单位圆,如下图所示:

在这个面积为1的正方形中,有四分之一的圆,圆的半径与正方向的边长都是1。那么根据圆的面积公式,这个图形中阴影部分的面积应该是 Pi/4。

下面开始进入蒙特卡洛的解法。

即,如果我们对这个正方形平面中的点进行均匀地抽样,随着抽样点的增多,那么落入阴影内的点的数量与总抽样数量的比,应该基本上等于阴影的面积Pi/4与整个正方形面积1的比,即Pi/4。用数学表示,就是

阴影内的样本点数量 ÷ 总数量 = Pi/4

所以,Pi = 4 × 阴影内的样本点数量 ÷ 总数量。

下面就在Excel中进行实验。

用rand()函数生成2000个随机数,作为随机样本点的X轴坐标,

再用rand()函数生成2000个随机数,作为随机样本点的Y轴坐标。

如此就得到了2000个随机样本点,这些点的X轴坐标和Y轴坐标都大于零且小于1,所以是在前面所说的正方形之中的点。

下一步,判断样本点是否处于阴影之内,由于这个阴影就是单位圆在直角坐标系第一想象的四分之一,所以圆阴影内的点都符合如下不等式:

翻译到Excel中,就是用IF函数来判断,例如:

IF(A2^2 + B2^2 <=1, 1, 0)

即,如果样本点在阴影中,得到1,否则得到0。这样就把样本点区分开来了。

最后,把所有得到的1和0加总,就知道所有样本点中处于阴影中样本点的数量了。

最后根据

Pi = 4 × 阴影内的样本点数量 ÷ 总数量

就可以算出Pi来了。

我这个试验中算出来的 Pi=3142。

以下是样本点的散点图:

由于样本数量有限,所以计算出来的Pi的精度并不高。

以下是工作界面,挺简单的。

来源:知乎

贴一个蒙特卡洛方法的matlab程序,供大家使用。

{3 x& K/ i1 i( D8 C0 c$ O

% Example Monte Carlo Simulation in Matlab 0 O5 \; P" t# t7 v8 c& @

% Function: y = x2^2/x1 5 Z0 W4 e9 q, d5 B+ c

%

% Generate n samples from a normal distribution 4 s! c6 y, I6 H" d) K+ v Y; X: Q

% r = ( randn(n,1) sd ) + mu 4 U F Q) t, T# q w/ K' Q

% mu : mean / E( P8 U" c o! G8 s/ x

% sd : standard deviation

%

% Generate n samples from a uniform distribution 2 u# ^& K [0 z% F) @1 y

% r = a + rand(n,1) (b-a) - D+ }& U$ w- M9 @& Q9 W, Z

% a : minimum

% b : maximum

n = 100000; % The number of function evaluations 7 x5 a" @- F& O- Z; w5 j

% --- Generate vectors of random inputs ! K& x0 ^# X+ q( V6 {

% x1 ~ Normal distribution N(mean=100,sd=5)

% x2 ~ Uniform distribution U(a=5,b=15)

x1 = ( randn(n,1) 5 ) + 100; 2 B' l3 n) V) D$ ~

x2 = 5 + rand(n,1) ( 15 - 5 ); \: O: Y( w3 [9 d: V4 r( k4 {

% --- Run the simulation

% Note the use of element-wise multiplication - ~% x$ `7 A6 v9 R F

y = x2^2 / x1; ' g$ O7 U; R F% `

% --- Create a histogram of the results (50 bins)

hist(y,50); / M9 m+ s( [ w" J2 I% s/ X

% --- Calculate summary statistics

y_mean = mean(y)

y_std = std(y) ; R7 A2 y M/ T" p, h m

y_median = median(y)

以上就是关于matlab程序问题。需要用到蒙特卡洛方法全部的内容,包括:matlab程序问题。需要用到蒙特卡洛方法、蒙特卡洛方法原理、详细说明:蒙特卡洛方法在排队论中的应用,和在Excel中如何应用蒙特卡洛-Meng te ka luo excel等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!

欢迎分享,转载请注明来源:内存溢出

原文地址: https://outofmemory.cn/zz/9284483.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-04-26
下一篇 2023-04-26

发表评论

登录后才能评论

评论列表(0条)

保存