蒙特卡诺法的matlab实现

原理介绍

简介

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

模拟步骤

产生模拟随机数的matlab命令

1.服从均匀分布

2.服从正态分布

注：适用于正态分布的情况：

3.服从指数分布

4.服从泊松分布

5.指数分布与泊松分布的关系

实战

第一例

clc,clear
count=100000;%进行模拟的次数
a=[20 15 10 5 5 3 3 3 3 1 1 1];%十一个数，代表股东的股份
flag=zeros(1,12);%作为标志数组
c=1:12;
for i=1:count
     b=c(randperm(length(a)));%b就是a中12个数的随机顺序排列,randperm(n)函数随机产生1到n的整数,这里再次在f里索引，为了提高随机性
    for j=1:12
        if(sum(a(b(1:j)))>=50) %判断前j个数是否大于等于50
            break;
        end
    end
    flag(b(j))=flag(b(j))+1; %j停止加一，要追溯到b(j)
end
flag=flag/count;            

结果：

第二例求积分

**题目：**y = x^2 ，y = 12 - x与X轴在第一象限与X轴围成一个曲边三角形。设计一个随机试验，求该图形面积的近似值。

思路：我们可以随机产生很多个点，然后求出点落在所求图形内的概率，然后利用求得的概率乘以矩形区域的面积即可求出面积。

x=unifrnd(0,12,[1,10000000]);%产生服从均匀分布的点
y=unifrnd(0,9,[1,10000000]);
count_point=sum(y<x.^2 & x<=3)+sum(y<12-x & x>=3)%求落入区域的点数
mianji=12*9*count_point/10^7

第三例求解非线性规划

%转化为最小值求解
clc,clear
lb=0;%下界
ub=10;%上界
n=1000;%试验次数
x1=unifrnd(lb,ub, 1,n);%产生服从均匀分布的随机数
x2=unifrnd(lb,ub,1,n);
possible_value=[x1(1) x2(1)];
min=inf;
for i=1:n
     possible_value=[x1(i) x2(i)];
    equ=acc(possible_value);
   
    if equ==1
        f_value=f(possible_value);
        if f_value<min
            min=f_value;
            X_good=possible_value;
        end
    end
end
X_good
result=-min%转化为最大值



%定义函数
function z=f(x)
z=2*x(1)^2+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2);
end
%设定精度值
function accuracy=acc(x)
    if 3*x(1)+x(2)-10<=0.5 && 3*x(1)+x(2)-10>=-0.5 %定义一个接近于0的数
            accuracy=1;
    else
            accuracy=0;
    end 
end

• 本程序重要的一点是添加了一个精度控制

第四例求解整数规划

% 处理整数规划问题
clc,clear
% 初始化随机数发生器
rng('shuffle') % 产生非重复随机数
p0=0;
tic % 计时开始
for i = 1:1e6
    x = randi([0,99],1,5);% 产生范围在[0,99]的整数1*5
    [f,g]=guihua(x);
    if all(g<=0)
        if p0 < f
            x0 = x;
            p0 = f; %记录下当前较好的解
        end
    end
end
p0

toc
% 定义目标函数及约束条件
function [f,g] = guihua(x)
    f = x(1)^2+x(2)^2+3*x(3)^2+4*x(4)^2+2*x(5)^2-8*x(1)-2*x(2)-3*x(3)-x(4)-2*x(5);
    g=[sum(x)-400
        x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-800
        2*x(1)+x(2)+6*x(3)-200
        x(3)+x(4)+5*x(5)-200];
end

• 注意：也可以利用lingo来实现