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

你给出的解是正确的

首先假设有编号为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块钱

蒙特卡洛模拟技术，是用随机抽样的方法抽取一组满足输入变量的概率分布特征的数值，输入这组变量计算项目评价指标，通过多次抽样计算可获得评价指标的概率分布及累计概率分布、期望值、方差、标准差，计算项目可行或不可行的概率，从而估计项目投资所承担的风险。

蒙特卡洛模拟的步骤如下：

第一步，通过敏感性分析，确定风险变量。

第二步，构造风险变量的概率分布模型。

第三步，为各输入风险变量抽取随机数。

第四步，将抽得的随机数转化为各输入变量的抽样值。

第五步，将抽样值组成一组项目评价基础数据。

第六步，根据基础数据计算出评价指标值。

第七步，整理模拟结果所得评价指标的期望值、方差、标准差和它的概率分布及累计概率，绘制累计概率图，计算项目可行或不可行的概率。

蒙特卡洛模拟程序如图7-26所示。

图7-26 蒙特卡洛模拟程序图

实训Ⅷ某项目建设投资为1亿元，流动资金1000 万元，项目两年建成，第三年投产，当年达产。不含增值税年销售收入为5000万元，经营成本2000万元，附加税及营业外支出每年为50万元，项目计算期12 a。项目要求达到的项目财务内部收益率为15%，求内部收益率低于15%的概率。

由于蒙特卡洛模拟的计算量非常大，必须借助计算机来进行。本案例通过手工计算，模拟20次，主要是演示模拟过程。

(1)确定风险变量。通过敏感性分析，得知建设投资、产品销售收入、经营成本为主要风险变量。流动资金需要量与经营成本线性相关，不作为独立的输入变量。

(2)构造概率分布模型。建设投资变化概率服从三角形分布，其悲观值为13亿元、最大可能值为1亿元、乐观值为9000万元，如图7-27所示。年销售收入服从期望值为5000万元、σ＝300万元的正态分布。年经营成本服从期望值为2000万元、σ＝100 万元的正态分布。

图7-27 投资三角形分布图

建设投资变化的三角形分布的累计概率，见表7-16及图7-27所示。

表7-16 投资额三角形分布累计概率表

(3)对投资、销售收入、经营成本分别抽取随机数，随机数可以由计算机产生，或从随机数表中任意确定起始数后，顺序抽取。本例从随机数表(表7-20)中抽取随机数。假定模拟次数定为k＝20，从随机数表中任意从不同地方抽取三个20 个一组的随机数，见表7-17。

表7-17 输入变量随机抽样取值

(4)将抽得的随机数转化为各随机变量的抽样值。

这里以第1组模拟随机变量产生做出说明。

1)服从三角形分布的随机变量产生方法。

根据随机数在累计概率表(表7-16)或累计概率图(图7-28)中查取。投资的第1个随机数为48867万元，查找累计概率048 867所对应的投资额，从表7-16中查得投资额在10300与10600之间，通过线性插值可得

第1个投资抽样值＝10300+300×(48867-39250)/(52000-39250)＝10526万元

2)服从正态分布的随机变量产生方法。

从标准正态分布表(表7-21)中查找累计概率与随机数相等的数值。例如销售收入第1个随机数06242，查标准正态分布表得销售收入的随机离差在-153与-154之间，经线性插值得-15348。

图7-28 投资的累计概率分布图

第1个销售收入抽样值＝5000-15348×300≈4540万元。

同样，经营成本第一个随机数66 903相应的随机变量离差为04328，第一个经营成本的抽样值＝2000+100×04328＝2043万元。

3)服从离散型分布的随机变量的抽样方法。

本例中没有离散型随机变量。另举例如下，据专家调查获得的某种产品售价的概率分布见表7-18。

表7-18 某种产品售价的概率分布

根据上表绘制累计概率如图7-29所示。

若抽取的随机数为43252，从累计概率图纵坐标上找到累计概率为043252，划一水平线与累计概率折线相交的交点的横坐标值125元，即是售价的抽样值。

(5)投资、销售收入、经营成本各20个抽样值组成20组项目评价基础数据。

(6)根据20组项目评价基础数据，计算出20 个计算项目评价指标值，即项目财务内部收益率。

(7)模拟结果达到预定次数后，整理模拟结果按内部收益率从小到大排列并计算累计概率，见表7-19所示。

从累计概率表可知内部收益率低于15%的概率为15%，内部收益率高于15%的概率为85%。

图7-29 售价累计概率曲线

表7-19 蒙特卡洛模拟法累积概率计算表

①每次模拟结果的概率＝1/模拟次数。

贴一个蒙特卡洛方法的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)

在业务中使用电子表格的三十多年中，我所见过的大多数电子表格模型和预测都使用了统计学家所说的确定性方法。

为了说明这一点，如果我要预测一个时期的利润，那么确定性模型将使用一个数字作为我的销售预测，使用另一个数字作为我的运营支出预测，依此类推。

相反，一种更有用的方法将采用概率方法，并得到蒙特卡洛方法的支持。

过去，我没有考虑太多，而是编写了宏来支持蒙特卡洛方法。

但是那时候我错了。有一种更简单的方法。

蒙特卡罗方法

蒙特卡洛方法支持“随机”或“概率”系统。

法线曲线

根据Wikipedia所说，“纯随机系统是状态被随机确定的系统，具有随机的概率分布或模式，可以进行统计分析，但可能无法精确预测。”

此图说明了正态概率分布，这可能是大多数业务使用的最佳方法。

要使用此方法，首先要设置要由该曲线定义的关键假设。然后，工作簿中的模型会从此正态分布中随机选择值，在其计算中使用这些值，记录关键结果，多次重复该过程，然后汇总结果记录。此过程称为蒙特卡洛方法。

在本文中，我将向您展示如何使用Excel Data Tables记录每次计算产生的结果。

但是，在详细介绍之前，我需要提到两件事。

首先，每当您打开使用数据表的蒙特卡洛分析时，请确保“蒙特卡洛”工作簿是唯一打开的工作簿。这是因为它将需要多次重新计算，并且如果您打开了其他工作簿，它们也将不必要地重新计算。这会使您的模拟非常缓慢。

其次，您可能会对我的示例感到不满。

如果您从事财务工作，您可能会抱怨以下示例大大简化了通常是复杂的财务建模过程的过程。在这种情况下，我会回答说我使模型非常简单，因此您可以了解我建议的对标准建模过程的更改。

而且，如果您不从事财务工作，您可能会抱怨我应该使用您自己的专业而不是财务方面的例子。在这种情况下，我会回答我以简单的损益表为例，因为即使您从事工程，运营，市场营销或任何地方的工作，您都可以理解简单的损益表。因此，您将能够了解模型的功能，并且能够将我的技术应用于您自己的模型和预测。

] - 排队论中的一个仿真程序，主要是用于仿真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程序问题。需要用到蒙特卡洛方法全部的内容，包括:matlab程序问题。需要用到蒙特卡洛方法、蒙特卡洛模拟法、怎样用matlab实现蒙特卡洛仿真等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！