Search code examples
randomfortransimulationrandom-seedstochastic-process

How to change seed number in Fortran stochastic simulator code


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

Solution

  • 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