[资源] 【转帖】拟蒙特卡洛方法(Quasi-Monte Carlo)积分实例

%使用Matlab提供的函数求积分,exp(-1/2*x^2)在(0,1)间积分 format long;

syms x

a = sym(1/2);

f = exp(-a*x^2);

ezplot(f)

disp(int(f,-1,1));

fprintf('integral result:%1.18f.\n',double(int(f,0,1)));

%disp(double(int(f,0,1)));
复制代码%使用拟蒙特卡洛方法积分
%得到拟蒙特卡洛序列,即低偏差序列,halton法
%如果有相关的工具箱的话,可以用Matlab里面的haltonset,faureset,sobolset函数实现,

x=halton(10000,2,5577);

n=length(x);

mju=0;

for i=1:n

    mju=mju + exp(-0.5*x(i)^2);

end

mju=mju/n;

fprintf('Quasi-Monte Carlo result:%1.18f.\n',mju);

%disp(mju);

%使用蒙特卡洛方法积分

%得到Uniform序列,

x=random('unif',0,1,10000,1);

n=length(x);

mju=0;

for i=1:n

    mju=mju + exp(-0.5*x(i)^2);

end

mju=mju/n;

fprintf('Monte Carlo result:%1.18f.\n',mju);



%=============生成HALTON序列========================

function result = halton( m,base,seeder )

%生成HALTON序列

% Check inputs

if nargin < 3

seeder = 0;

if nargin < 2

      error('MATLAB:Halton:NotEnoughInputs',...

             'Not enough input arguments. See Halton.');

end

end

res=0;

n=length(base);

for i=1:m

   

    for j=1:n

        element=0;

        temp=seeder+i;

        k=1;

        while temp>0

            element(k)=rem(temp,base(j));

            temp=fix(temp/base(j));

            k=k+1;

        end

        res(i,j)= 0;

        for k=1:length(element)

            res(i,j)=res(i,j)+element(k)/(base(j)^k);

        end

    end

   

end

result=res;