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
    CALL RANDOM_SEED(size = n)
    seed = clock + 37 * (/ (i - 1, i = 1, n) /)
  • 1
    $\begingroup$ 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}$. $\endgroup$ Feb 27 '16 at 21:22
  • $\begingroup$ well that certainly addressed the magnitude of variance issue. I suppose that I am overly concerned with a small fluctuation $\endgroup$
    – Joseph
    Feb 27 '16 at 22:15
  • $\begingroup$ You should check whether the mean indeed converges like $1/\sqrt{N}$. $\endgroup$ Feb 28 '16 at 19:17
  • 1
    $\begingroup$ You might be interested by this web site : pcg-random.org $\endgroup$ Feb 29 '16 at 12:35
  • $\begingroup$ 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. $\endgroup$ Jun 5 '16 at 11:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.