昨天,呃,不,是今天凌晨,我问一位大牛有没有什么既简单又强势的算法,他说了蒙特卡洛。今天查了些资料,见识到它的强大与及简洁。参考了这篇文章,和维基百科。
现在用MATLAB实现蒙特卡洛方法的几个应用。蒙特卡洛算法MCA
1.计算圆周率
status = 10;
for i=1:6 %6次模拟
point=status.^i; %模拟的随机点数
RandData=rand(2,point); % 根据随机点数,产生随机的(x,y)散点
x = RandData(1,:);
y = RandData(2,:);
inner=find(x.^2+y.^2 <= 1); % find返回符合条件(在圆内)的点的索引值,构成数组
Outcome(i)=length(inner)/length(RandData); % 计算圆内的点占所有点的比例,即为pi/4
my_pi(i) = Outcome(i)*4;
end
my_pi
输出如下:
可以看到,随着随机点数的增加,算得的pi精度越来越高。(注意这里,我将脚本文件命名为pi.m,系统会有冲突警告,所以最好改下文件名。)
2.计算积分
这里以计算y=x^2在[0,1]上的定积分为例。
status = 10;
for i=1:4 %4次模拟
point=status.^i; %模拟的随机点数
RandData=rand(2,point); % 根据随机点数,产生随机的(x,y)散点
x = RandData(1,:);
y = RandData(2,:);
%scatter(x,y)
Below=find(x.^2>y); % y=x^2曲线下的点(的下标)
Outcome(i)=length(Below)/length(RandData);
subplot(2,2,i)
scatter(x(Below),y(Below),'r')
end
axis([0,1,0,1])
Outcome
输出如下:
3.计算测度
这里参考一篇论文,讲解关于测度的MATLAB实现。
相关代码如下:
clear all
n = input('please input the dimensionn n of variable=:');
m=0;
for i=1:n
x=rand(1);
y=rand(1);
z=rand(1);
if x^2+sin(y)<=z
if x-z+exp(y)<=1
m=m+1;
end
end
end
ration=m/n
运行结果如下:
可见结果运行还是很稳定的。
4.计算非线性规划的优化解
代码如下:
clear all
L=10000; %要求最小值,设置足够大的上限
M = input('set the number of dots M=')
for N=1:M
% 将不等关系巧妙地转化为乘积关系
x2=-9+13*rand(1);
x1=((x2)^2+9)*rand(1);
x3=3*(x1)+(x2)-5+(12-(x1)*(x2)-3*(x1)-(x2))*rand(1);
if 12-(x1)*(x2)-3*x1-x2>=0
f=3*(x1)^2-4*(x2)^2+5*(x1)*(x3);
if f<L
L=f;
end
end
end
L
运行输出结果:
并且可行域和目标函数的复杂性不对MC方法造成影响。
最后,看到知乎上有句话总结了蒙特卡洛方法——尽量找好的,但不保证是最好的。多么值得思考的一句话啊。