I am having a really hard time getting any kind of reliable / consistent result from my Metropolis code. I have torn it apart and am now examining just the randomness in my random number generator.
I decided I would just run many iterations over a spread from -1 to 1 and it should give me an average of 0. I was completely surprised by the lack of precision. My code does give results about zero, but, from -200 to 200 with almost an equal distribution throughout the range.
I'm trying to find a way to bring this closer together. I have tried more runs (knowing 10,000 isn't really that many) but even up to a billion runs, the averages can range just as much if not more.
Any suggestions on how to fix this randomness?
I have been using:
program test_rand real*8 x, y, z, EAx, EAy, EAz, Tx, Ty, Tz integer j, i CALL init_random_seed Do i = 1, 10000 EAx = 0.D0 EAy = 0.D0 EAz = 0.D0 Tx = 0.D0 Ty = 0.D0 Tz = 0.D0 Do j=1,10000 CALL RANDOM_NUMBER(x) CALL RANDOM_NUMBER(y) CALL RANDOM_NUMBER(z) EAx = 2.0D0 * x - 1.0D0 EAy = 2.0D0 * y - 1.0D0 EAz = 2.0D0 * z - 1.0D0 Tx = EAx + Tx Ty = EAy + Ty Tz = EAz + Tz end do write(50,*) i, Tx, Ty, Tz end do end program test_rand ! initialize a random seed from the system clock at every run (fortran 95 code) subroutine init_random_seed() INTEGER :: i, n, clock INTEGER, DIMENSION(:), ALLOCATABLE :: seed CALL RANDOM_SEED(size = n) ALLOCATE(seed(n)) CALL SYSTEM_CLOCK(COUNT=clock) seed = clock + 37 * (/ (i - 1, i = 1, n) /) CALL RANDOM_SEED(PUT = seed) DEALLOCATE(seed) end