c9_example5.m

motecalo仿真实现源代码

M=10;
N=10000;
u=rand(N,M);
uu=1./(1+u.*u);
data=zeros(N,M);


data(1,:)=4*uu(1,:);
for j=2:N
    data(j,:)=4*sum(uu(1:j,:))/j;
end
est=data(N,:);
est1=sum(est)/M
plot(data,'k')
xlabel('number of trial')
ylabel('estimate of pi')