c9_6.m

motecalo仿真实现源代码

mean=sqrt(pi)*(0.5-q((sqrt(2))));
int2=sqrt(pi/2)*(0.5-q(2));
varx=int2-mean*mean;
stdx=sqrt(varx);
alpha=0.1;

nsum=0;
nppseg=1000;
nseg=100;
est=zeros(1,nseg);

for j=1:nseg
    ui=rand(1,nppseg);
    gui=sum(exp(-ui.*ui));
    est(j)=nsum/(j*nppseg);
end

nn=nppseg*(1:nseg);
ub=mean+stdx*qinv(alpha/2)./sqrt(nn);
lb=mean-stdx*qinv(alpha/2)./sqrt(nn);
meanv=mean*ones(1,nseg);
si=1:nseg;
plot(si,est,'k-',si,meanv,'k--',si,ub,'k:',si,lb,'k:')
xlabel('number of segments')
ylabel('estimate of intergal')
legend('estimate value','true value','upperbound','lower bound');