c9_example5.m

来自「motecalo仿真实现源代码」· M 代码 · 共 17 行

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')