モンテカルロ法精度 - くもしかつかまない

モンテカルロ法において粒子に対する精度の話になったときに「えー，粒子数を $n$ としたとき精度は $\sqrt{n}$ でしょー」という話を聞いたときから，モンテカルロ法の精度が気になってはいたんですよね．いったい精度が $\sqrt{n}$ とはどういう状態を指すのか．
というわけで今回調べてみました．（モンテカルロ法自体については，itさんの素晴らしい解説があるのでそちらに丸投げで　ﾟωﾟ) つ○）

さて，精度が $\sqrt{n}$ というのは，粒子数が $n$ のモンテカルロ法の解を $a_{n}$ ，真の解を $a$ したとき，

$\left| a_{n} - a \right| \sim \sqrt{n}$

が成り立つことかなー，とか思ってたんですけど，これだと左辺はnに対して単調減少なのに右辺は単調増加でおかしなことになるわなー，ということで，こうか，と．

仮説：モンテカルロ法の精度が $\sqrt{n}$ であるとは，次の式が成り立つことである．
$\left| a_{n} - a \right| \sim \frac{1}{\sqrt{n}}$

つまり，Log10(粒子数)に対して誤差( $\left| a_{n} - a \right|$ )をプロットしたときに，

この曲線に乗るようにプロット点がくるのかなぁと．

さて，実際に粒子数を増やして $\pi$ を順次求めるプログラムを書いてみる．

Program MonteCarlo
  implicit none
  real(8) :: pi, seed, seed2, test, ans,error,fe,re,nr
  real(8),dimension(100000000) :: x,y
  integer :: i,j,k,ntest,nmax,sum,fntest

  pi    = 3.14159265358979d0

  call random_number(x)
  call random_number(y)

  write(*,'(A16,A9,A16,A16,A16)'),"squre of n ratio","ntest","answer","error","error ratio"
  do i = 1,8
     ntest = 10 ** i
     sum   = 0
     error = 0.d0
     do j = 1,ntest
       test = x(j) ** 2 + y(j) ** 2
       if (test < 1.d0) sum = sum + 1
     end do
     ans = 4.d0 * dble(sum) / dble(ntest)
     error = ans - pi
     error = abs(error)
     if(i == 1) then
       fe = error
       re = fe / fe
       fntest = ntest
     else
       re = error / fe
     end if
     nr = dble(ntest) / (fntest)
     nr = sqrt(nr)
     write(*,'(f16.4,I9,f16.8,f16.8,f16.8)'),nr,ntest,ans,error,re
  end do

end Program MonteCarlo