I'm running a Fortran code which performs a stochastic simulation of a marked Poisson cluster process. In practice, event properties (eg. time of occurrences) are generated by inversion method, i.e. by random sampling of the cumulative distribution function. Because of the Poissonian randomness, I expect each generated sequence to be different, but this is not the case. I guess the reason is that the seed for the pseudorandom number generator is the same at each simulation. I do not know Fortran, so I have no idea how to solve this issue. Here is the part of the code involved with the pseudorandom number generator, any idea?
subroutine pseud0(r)
c generation of pseudo-random numbers
c data ir/584287/
data ir/574289/
ir=ir*48828125
if(ir) 10,20,20
10 ir=(ir+2147483647)+1
20 r=float(ir)*0.4656613e-9
return
end
subroutine pseudo(random)
c wichmann+hill (1982) Appl. Statist 31
data ix,iy,iz /1992,1111,1151/
ix=171*mod(ix,177)-2*(ix/177)
iy=172*mod(iy,176)-35*(iy/176)
iz=170*mod(iz,178)-63*(iz/178)
if (ix.lt.0) ix=ix+30269
if (iy.lt.0) iy=iy+30307
if (iz.lt.0) iz=iz+30323
random=mod(float(ix)/30269.0+float(iy)/30307.0+
& float(iz)/30323.0,1.0)
return
end
First, I would review the modern literature for PRNG and pick a modern implementation. Second, I would rewrite the code in modern Fortran.
You need to follow @francescalus advice and have a method for updating the seed. Without attempting to modernizing your code, here is one method for the pseud0
prng
subroutine init0(i)
integer, intent(in) :: i
common /myseed0/iseed
iseed = i
end subroutine init0
subroutine pseud0(r)
common /myseed0/ir
ir = ir * 48828125
if (ir) 10,20,20
10 ir = (ir+2147483647)+1
20 r = ir*0.4656613e-9
end subroutine pseud0
program foo
integer i
real r1
call init0(574289) ! Original seed
do i = 1, 10
call pseud0(r1)
print *, r1
end do
print *
call init0(289574) ! New seed
do i = 1, 10
call pseud0(r1)
print *, r1
end do
print *
end program foo