蒙特卡洛方法模拟小例子

       经过长期观察发现，我方指挥所对敌方目标的指示有50％是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部毁伤敌人火炮．

      现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。

     使用蒙特卡洛方法模拟50次打击结果：

function [out1 out2 out3 out4]=Msc(N)

% N开炮次数

% out1射中概率

% out2平均每次击中次数

% out3击中敌人一门火炮的射击总数

% out4击中敌人2门火炮的射击总数

k1=0;

k2=0;

k3=0;

for i=1:N

  x0=randperm(2)-1;

  y0=x0(1);

  if y0==1

 第%d次：指示正确||',i);

 第%d次：击中0炮

 第%d次：击中1炮

 第%d次：击中2炮

  else

 第%d次：指示错误，击中0炮

  end

end

out1=(k2+k3)/N;

out2=(0*k1+k2+2*k3)/20;

out3=k2/N;

out4=k3/N;

运行： 

1.[out1 out2 out3 out4]=Msc(50)

结果： 

1.第1次：指示正确||第1次：击中2炮||

2.第2次：指示错误，击中0炮||

3.第3次：指示错误，击中0炮||

4.第4次：指示正确||第4次：击中0炮||

5.第5次：指示错误，击中0炮||

6.第6次：指示正确||第6次：击中1炮||

7.第7次：指示正确||第7次：击中0炮||

8.第8次：指示错误，击中0炮||

9.第9次：指示正确||第9次：击中2炮||

10.第10次：指示正确||第10次：击中1炮||

11.第11次：指示正确||第11次：击中1炮||

12.第12次：指示正确||第12次：击中2炮||

13.第13次：指示错误，击中0炮||

14.第14次：指示正确||第14次：击中1炮||

15.第15次：指示错误，击中0炮||

16.第16次：指示错误，击中0炮||

17.第17次：指示正确||第17次：击中0炮||

18.第18次：指示错误，击中0炮||

19.第19次：指示正确||第19次：击中1炮||

20.第20次：指示错误，击中0炮||

21.第21次：指示正确||第21次：击中0炮||

22.第22次：指示正确||第22次：击中1炮||

23.第23次：指示正确||第23次：击中0炮||

24.第24次：指示错误，击中0炮||

25.第25次：指示正确||第25次：击中1炮||

26.第26次：指示错误，击中0炮||

27.第27次：指示正确||第27次：击中1炮||

28.第28次：指示正确||第28次：击中0炮||

29.第29次：指示正确||第29次：击中0炮||

30.第30次：指示正确||第30次：击中0炮||

31.第31次：指示错误，击中0炮||

32.第32次：指示错误，击中0炮||

33.第33次：指示正确||第33次：击中0炮||

34.第34次：指示错误，击中0炮||

35.第35次：指示正确||第35次：击中0炮||

36.第36次：指示正确||第36次：击中0炮||

37.第37次：指示错误，击中0炮||

38.第38次：指示正确||第38次：击中0炮||

39.第39次：指示错误，击中0炮||

40.第40次：指示正确||第40次：击中0炮||

41.第41次：指示正确||第41次：击中1炮||

42.第42次：指示正确||第42次：击中0炮||

43.第43次：指示错误，击中0炮||

44.第44次：指示正确||第44次：击中1炮||

45.第45次：指示正确||第45次：击中0炮||

46.第46次：指示错误，击中0炮||

47.第47次：指示错误，击中0炮||

48.第48次：指示错误，击中0炮||

49.第49次：指示正确||第49次：击中0炮||

50.第50次：指示正确||第50次：击中1炮||

51.

52.out1 =

53.

54.    0.2800

55.

56.

57.out2 =

58.

59.    0.8500

60.

61.

62.out3 =

63.

.    0.2200

65.

66.

67.out4 =

68.

69.    0.0600

一位朋友说要贴出Monte Carlo计算积分的源程序，我就随便做一个简单的吧，复杂的程序完全可以从这个来演化，我想Monte Carlo积分的最大优势就在于高维积分，以及不规则区域， 可以节约很多计算机时。

下面只是演示一个2重积分，可以扩展到20维的只要添加相应的loop项。

被积函数： exp(sqrt(x^2+y^2));

x 上下限：x^2 < sin(y)

y 上下限： y^22010-6-5 15:11:45 上传

下载附件 (42.75 KB) 

1.% MONTE CARLO INT

2.% by caoer

3.clear all

4.N =100000;

5.

6.x = 2*rand(N,1)-1;

7.y = rand(N,1);

8.f = 0;

9.fsq = 0;

10.n = 0;

11.for i=1:N

12.        if x(i)^2 < sin(y(i)) & y(i)^213.        n=n+1;

14.        f=f+exp(sqrt(x(i)^2+y(i)^2)); 

15.        x_plot(n) = x(i);

16.        y_plot(n) = y(i);

17.    end

18.end

19.f = f/N;

20.p=n/N;

21.a = 2*1;

22.I1 = f*a;  %

23.I2 = a*p;  %interesting area

24.I = I1/I2

25.

26.plot(x_plot,y_plot,'o')

复制代码

蒙特卡洛法用于求积分时，与积分重数无关，这点非常重要。虽然四维以下的积分用蒙特卡洛法效率可能不如传统的一些数值积分方法，但是维数高的时候，蒙特卡洛法比传统方法要有效的多，而且实现起来也非常容易。可以说，计算高维积分是蒙特卡洛方法最成功和典型的应用。

基本的蒙特卡洛法具有计算不可重复性的缺点。这里共享采用等序列分布的蒙特卡洛法 等分布序列Monte Carlo积分.doc (92.5 KB, 下载次数: 24) 

2010-6-5 15:09:28 上传

下载次数: 24 

下载积分: 贝壳 -1 

，具有计算可重复性，误差阶比采用基本蒙特卡洛法好的优点。就基本的蒙特卡洛法求积分来说，不管积分重数多少，基本上是计算规模增加100倍，精度提高10倍。