Question

蒙特卡罗随机模拟

我想通过matlab用蒙特卡罗方法模拟200年的灾害损失数据，灾害损失程度服从参数为u=9.1388，σ=1.3577的对数正态分布，灾害损失次数服从λ=8的泊松分布，然后从每年的数据中选出最大值。急用！请写请楚每一步的注释，谢谢！如能让我满意可以再给你追加一百分。
非常感谢你的回答，这次我把数据通过插图传上来了。因为是从每年当中选出一个极大值，最后选出的极大值应该有200个。还好，用的时候说没有“normrnd”这个函数？麻烦你再帮我弄一次，再次感谢……

Answer 1

楼主这个题目有点问题，灾害损失数据不知道，所以不知道你给出的数据是否符合你所说的那种特定分布（比如正态分布、泊松分布），但你要从某种分布（我用你给出的分布：正态、泊松分布做出模拟）的数据中找出最大值还是有办法的，希望能对你具体计算有所启发。如下：

    （我用随机数模拟的方法，假设灾害损失数据个数为200个，灾害损失次数为200个）Matlab程序：

     %%%%%%%%%%

     %蒙特卡罗方法模拟灾害问题

     clc;clear;

     r=normrnd(9.1388,1.3577,[1 200])  %建立灾害损失程度随机数组（u为 9.1388, sigama为1.3577的正态分布，数据为200个）

     p=poissrnd(8,[1 200])  %建立灾害次数随机数组（lamda等于8的poisson分布，数据为200个）

     rmax=max(max(r))  %输出正态矩阵中最大值

     pmax=max(max(p))  %输出泊松矩阵中最大值

     %%%%%%%%%%

结果如下:

rmax =12.6008

pmax = 21

这个很简单的。掌握了这两个函数就行了。有什么问题欢迎探讨QQ：182335530

6月1日·回答补充

    首先，normrnd这个函数是存在的。其次，看了你发的数据图，我对你问的问题有点迷惑，请问你要的“从每年当中选出一个极大值”这个“极大值”是否是指的：“每年最有可能遭到的最大经济损失”？如果是的话，答案很明显，是1亿元乘以每年发生的灾害次数，并不需要matlab计算，直接观察就可以得到答案；如果指的是：“每年最有可能遭受的经济损失”，那么这是一个二维随机变量的计算，容易求出。 

6月2日·回答补充

    首先说声对不起，问题反复研究，我发现你给出的条件不够，问题在于给出的灾害经济损失数据是确定的，而又说其服从对数正态分布，其实此处服的数据服从哪种分布已经不重要了，因为数据是确定的。如果这个数据需要用对数正态分布的蒙特卡罗方法模拟的话或许会有其他办法，因为你问题的题目是用的蒙特卡罗方法。除此之外，我能力有限，暂时想不到其他方法来有效解决你所提出的问题，再次再次说声抱歉。

    现在按照我的理解解这道题，我求的最终结果是每年最有可能遭受的经济损失，而不是求其中某次最大的。蒙特卡罗方法的体现在于对于其中年份数据（年灾害次数）的模拟，使用了lamda=8的泊松分布，函数为poissrnd；然后求出每年最有可能遭受的经济损失（期望）；再给出年灾害次数与各类灾害造成经济损失（亿元）的二维随机变量分布；最后计算出这个分布列中期望与频数的数量积，此值就为：每年最有可能遭受的经济损失。

%%%%%%%%%%%%%Matlab程序如下：%%%%%%%%%%%%%%%

format long;clc;clear;

%相应灾害等级的损失频数：

f1=[142 50 26 13 6 4 1 1 0 1];

%1990年到2006年的灾害次数数据

fa=[8 8 5 8 7 5 9 9 9 9 7 10 4 18 10 11 7];

%模拟1990年之前对应183年间每年发生的灾害次数，服从参数为8的泊松分布如下：

fb=poissrnd(8,[1 183]) ;

%连接前后两组数据

f2=cat(2,fa,fb);         %实际数据

subplot(2,1,1);

plot(1807:2006,f2,'*r'); %自1807年至2006年每年发生灾害的次数图像

title('每年灾害次数服从amda为8的泊松分布');

subplot(2,1,2);

plot(1:10,f1,'*b');

title('灾害经济损失服从对数正态分布')

%每次灾害时发生相应经济损失（亿元）的概率：

for i=1:10

    p1(i)=f1(i)/sum(f1);

end

%200年间每年发生灾害的概率

for j=1:200

    p2(j)=f2(j)/sum(f2);

end

%关于经济损失和年份的二维随机变量

for j=1:200

    for i=1:10

        c(j,i)=p2(j)*p1(i);

    end

end

%计算每年造成的经济损失期望emax

for i=1:10

    e(i)=p1(i)*i;

end

emax=sum(e);

%损失期望指数；参考：周俊华《中国台风灾害风险评估》p52-p55

%计算相应年份发生灾害的损失

f3=f2' ;%转置每年灾害次数

money=f3*emax;

disp('相应年份最可能造成的经济损失');

money

fp=fopen('d:\money.txt','wt');

fprintf(fp,'%4.2f ',money);

fclose(fp);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

计算结果如下：

14.79 14.79 9.24 14.79 12.94 9.24 16.64 16.64 16.64 16.64 12.94 18.48 7.39 33.27 18.48 20.33 12.94 16.64 7.39 11.09 11.09 12.94 12.94 14.79 14.79 12.94 7.39 11.09 12.94 12.94 12.94 11.09 18.48 16.64 12.94 12.94 14.79 9.24 25.88 11.09 29.57 14.79 20.33 14.79 16.64 14.79 9.24 18.48 16.64 14.79 11.09 20.33 9.24 18.48 7.39 11.09 24.03 20.33 14.79 5.55 14.79 16.64 5.55 11.09 14.79 12.94 9.24 7.39 14.79 24.03 3.70 33.27 11.09 5.55 9.24 7.39 11.09 16.64 18.48 18.48 12.94 14.79 20.33 11.09 7.39 12.94 9.24 14.79 20.33 20.33 9.24 18.48 7.39 18.48 14.79 14.79 18.48 11.09 9.24 9.24 12.94 35.12 29.57 22.18 9.24 16.64 5.55 12.94 9.24 20.33 12.94 20.33 14.79 16.64 18.48 9.24 12.94 25.88 7.39 12.94 11.09 9.24 16.64 9.24 11.09 14.79 14.79 16.64 16.64 7.39 11.09 14.79 11.09 16.64 16.64 18.48 5.55 18.48 14.79 16.64 14.79 11.09 9.24 11.09 16.64 11.09 12.94 22.18 12.94 16.64 20.33 12.94 11.09 12.94 11.09 11.09 16.64 11.09 7.39 18.48 16.64 12.94 12.94 18.48 20.33 7.39 14.79 3.70 9.24 24.03 25.88 14.79 20.33 12.94 9.24 20.33 29.57 9.24 12.94 20.33 16.64 18.48 14.79 18.48 7.39 29.57 9.24 24.03 22.18 9.24 7.39 12.94 12.94 11.09 24.03 18.48 14.79 11.09 14.79 11.09

%%%%%%%%文件存储于D:/money.txt，是short型数据%%%%%%%%%%%% 

模拟点如下：