2D random number uniformity

Help
1

Hey! I wrote a code to test uniformity of random numbers in 1D and now i need to do the same for 2D and so on until 4, but i’m new to this language and i’m struggling a bit. The code for 1D that i wrote is:

program teste1d
    implicit none

    real(8), dimension(:), allocatable :: v, n 
    real, dimension(:), allocatable :: cell
    real :: c
    integer :: i, k, j
    logical :: cond = .true.

    k = 1000000 !total of random numbers
    j = 20 !number of cells
    
    allocate(v(k)) !vector with the random numbers
    allocate(cell(j)) !how many numbers fall in each cell
    allocate(n(j))
    
    call random_number(v)

    !inicializing the vectors as zeros
    cell = cell*0.0 
    n = n*0.0

    !vector just to plot the result
    do i=1,size(n)
        n(i) = n(i) + i
    end do

    do i=1, size(v) 
        c = (1.0/j)     
        cond = .true.
        do while(cond) 
            if(v(i) < c) then 
                cell(int(j*c)) = cell(int(j*c)) + 1 !cell 1 recieve 1
                cond = .false.
            else
                c = c + (1.0/j)
            end if 
        end do 
    end do 


    !just to plot the result
    open(40, file = 'dados.dat', status='unknown')  
        do i=1,size(n)
            write(40,*)n(i), cell(i)
        end do
    close(40) 
    call system ('gnuplot -p 1d.plt')

end program teste1d

Now, i was trying to do pretty much the same but i’m having difficult to acess the pair of random numbers, my code is avaluating every value instead of the one single pair of two cooordinates. The code is: (In this one i used another RNG):

program teste2d
    use rf_mod, only: rf
    implicit none
    
    real(8), dimension(:,:), allocatable :: v
    real, dimension(:,:), allocatable :: cell
    real :: c 
    integer :: N=5, i, j , k
    logical :: cond = .true.
    
    k = 10

    allocate(v(k,2))
    allocate(cell(N,N))
    cell = cell*0.0 

    do i=1,k
        do j=1, 2
            v(i,j) = rf()
        end do
    end do
    !print *, v

    do i=1, k
        do j=1, 2
            c = 1.0/(N)
            cond = .true.
            do while(cond)
                if(v(i,j) < c) then
                    cell(int(N*c),int(N*c)) = cell(int(N*c),int(N*c)) + 1 !a celula 1 recebe 1
                    cond = .false.
                else
                    c = c + (1.0/N)
                end if 
            end do    
        end do
    end do
    print*, cell, sum(cell)
end program teste2d

module rf_mod
    contains
    DOUBLE PRECISION FUNCTION RF()
        IMPLICIT DOUBLE PRECISION (A-H, O-Z)
        data IX1,IX2,IX3,IX4,IX5,IX6/1500419,1400159,1364,1528,1,3/
        RR1=1.0/FLOAT(IX1)
        RR2=1.0/FLOAT(IX2)
        IX5=MOD(IX5*IX3,IX1)
        IX6=MOD(IX6*IX4,IX2)
        RF=RR1*IX5+RR2*IX6
        IF(RF.GE.1.0)RF=RF-1.0
    END FUNCTION
    end module

The ideia is illustrated in this figure :

image
image881×666 59.3 KB

2 Likes
2

You are progressing at light speed! :smiley:

A couple of tips – you can set an array equal to a scalar value, and that value will be propagated to the whole array, as in

cell = 0.0

In 2-D, you probably want to draw two random numbers at a time, and in 3-D, three at a time. If you are content to use a single RNG, you could define the function to take an argument (nval), and modify the code to return a single real random number for nval = 1, a pair for nval = 2, and a triplet for nval = 3.

module rf_mod
    contains
    FUNCTION RF(n) result (nR)
        IMPLICIT NONE
        integer, parameter :: IX1 = 1500419, IX2 = 1400159, IX3 = 1364, IX4 = 1528
        integer, intent(in) :: n
        integer i, ix5, ix6
        double precision rfv
        double precision, dimension(n) :: nR
        double precision, parameter :: RR1 = 1d0/IX1, RR2 = 1d0/IX2
        data ix5, ix6 /1, 3/
!
        do i = 1, n
           ix5 = MOD(ix5*IX3,IX1)
           ix6 = MOD(ix6*IX4,IX2)
           nR(i) = RR1*IX5+RR2*IX6
        end do
        where (nR >= 1.0d0) nR = nR - 1d0
    END FUNCTION
end module

program tst
use rf_mod
implicit none
double precision :: RFn(10)
integer i, n
!
n = 10
RFn = RF(n)
!
print '(i3,F10.5)',(i,RFn(i),i = 1,n)
end program
1 Like
3

Just a few remarks:

  1. Making the cell array (which is actually used to count the random numbers in your virtual cells) a real array is not a good idea. Why not make it integer? Beside the obvious reason (adding 1 is simpler in integers) there is another one. Compile and the following program:
program main
  implicit none
  real :: x=17000000
  integer :: i, ix=17000000

  do i=1, 20
    x = x + 1
    ix = ix + 1
    print *,ix, x
  enddo
end program main

Due to the constrains of the FP representation, for 4-byte reals (typical) you can successfully add 1 only to numbers up to about 16 million. Compare 2 billion for a 4-byte integer!

  1. Finding the cell for a given number can be as easy as
    call random_number(v)
    do i=1, size(v)
      nc = int(j*v(i))+1  ! nc: 1..j
      cell(nc) = cell(nc)+1
    end do
  1. If using a linear congruential generator (LCG) to build n-tuples (points in n-dimensional unit cube), the points will lie on, at most, (m*n!)**(1/n) hyperplanes (Marsaglia’s theorem) where m is the maximum period of the LCG (the modulo number used) . You must take it into account when planning testing uniformity of random numbers in n-D space.
4

If you want very ‘uniformly’ distributed random number, you may consider quasi-random number. Such as Sobol sequence, etc. Low-discrepancy sequence - Wikipedia
However since they are quasi (so not really random), central limit theorem does not apply anymore. Also due to the ‘uniformity’, they more or less suffer from curse of dimensionality as dimension grows.

1 Like
5

I would like to do that but in this code you send me, the code keeps sending me the message “Symbol ‘nr’ at (1) already has basic type of REAL”, so i don’t know how to make this generator give me an array of random numbers without using a loop.

6

Remove DOUBLE PRECISION from the RF function header: