# Generating Random Numbers in Fortran for Metropolis method

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

• You compute the sum of random values, which in general diverges like $\sqrt{N}$. You probably want to compute the mean of the random values, which should converge like $1/\sqrt{N}$. Feb 27 '16 at 21:22
• well that certainly addressed the magnitude of variance issue. I suppose that I am overly concerned with a small fluctuation Feb 27 '16 at 22:15
• You should check whether the mean indeed converges like $1/\sqrt{N}$. Feb 28 '16 at 19:17
• You might be interested by this web site : pcg-random.org Feb 29 '16 at 12:35
• You may also write the EAx, EAy and EAz to a file and plot a histogram, this will give you an idea of the distribution. An example in Python: "import numpy as np ; import matplotlib.pyplot as plt ; data = np.loadtxt('fort.50').reshape((-1,3)) ; plt.hist(data); plt.show();". Further hint: to debug RNG problems, make the seed an input variable. Jun 5 '16 at 11:35