本帖转载于Matlab技术论坛，原帖参见http://www.matlabsky.com/thread-673-1-1.html

%使用拟蒙特卡洛方法积分
%得到拟蒙特卡洛序列,即低偏差序列,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;