Probleme mit der Einstellung zufälligen Seed

Question

Ich versuche, einen Montecarlo-Algorithmus zu schreiben, um die Interaktion zwischen Agenten in einer Population zu simulieren. Dieser Algorithmus muss bei jeder Iteration zwei Zufallszahlen aufrufen (z. B. 10^9 Iterationen).

Mein Problem hier ist, dass jedes Mal, wenn ich den Samen ändere (um verschiedene Realisierungen zu erhalten), der RNG mir die gleiche Ausgabe gibt (gleiche Montecarlo-Ereignisse). Ich habe verschiedene Möglichkeiten ausprobiert, um es ohne Erfolg zu implementieren.

Ich bin völlig neu in Fortran und kopiere diesen Code von MATLAB. Mache ich etwas falsch in der Art, wie ich diesen Code implementiere?

Unten ist das, was ich versucht:

program Gillespie 

    implicit none 

    integer*8, parameter :: n_max = 10.0**8 ! max. number of iterations 
    integer*8 :: t_ext, I_init, S_init, jump, S_now, I_now, i, u 
    real*8 :: t, N, a0, tau, st, r1, r2 

    real, parameter :: beta = 1000 
    real, parameter :: gammma = 99.98 
    real, parameter :: mu = 0.02 
    real, parameter :: R0 = beta/(gammma+mu) 
    integer :: seed = 11 

    real, dimension(n_max) :: S_new ! susceptible pop. array 
    real, dimension(n_max) :: I_new ! infected pop. array 
    real, dimension(n_max) :: t_new ! time array 
    real, dimension(5) :: events ! events array 

    open(unit=3, file='SIS_output.dat') 

    t = 0      ! initial time 
    N = 40     ! initial population size 

    jump = 1     ! time increment (save in uniform increments) 
    u = 2 

    t_ext = 0     ! extiction time 
    I_init = 2    ! initial infected pop. 
    S_init = N-I_init   ! initial susceptible pop. 

    S_now = S_init 
    I_now = I_init 

    S_new(1) = S_init   ! initialize susceptibles array 
    I_new(1) = I_init   ! initialize infected array 
    t_new(1) = t    ! initialize time array 

    write(3,*) t_new(1), S_new(1), I_new(1) ! write i.c. to array 

    call random_seed(seed) 

    do i=2, n_max 
    call random_number(r1) 
    call random_number(r2) 

    events(1) = mu*N  ! Birth(S) 
    events(2) = mu*S_now ! Death(S) 
    events(3) = mu*I_now ! Death(I) 
    events(4) = (beta*S_now*I_now)/N ! Infection 
    events(5) = gammma*I_now ! Recovery 

    a0 = events(1)+events(2)+events(3)+events(4)+events(5) 

    tau = LOG(1/r1)*(1/a0) ! time increment 
    t = t + tau   ! update time 
    st = r2*a0    ! stochastic time??? 

    ! update the populations 
    if (st .le. events(1)) then 
     S_now = S_now + 1 
    else if (st .gt. events(1) .AND. st .le. 
#(events(1) + events(2))) then 
     S_now = S_now - 1 
    else if (st .gt. (events(1) + events(2)) .AND. st .le. 
#(events(1) + events(2) + events(3))) then 
     I_now = I_now - 1 
    else if (st .gt. (events(1) + events(2) + events(3)) .AND. 
#st .le. (events(1) + events(2) + events(3) + events(4))) then 
     S_now = S_now - 1 
     I_now = I_now + 1 
    else 
     S_now = S_now + 1 
     I_now = I_now - 1 
    end if 

    ! save time in uniform increments 
    if(t .ge. jump) then 
     t_new(u) = floor(t) 
     S_new(u) = S_now 
     I_new(u) = I_now 

     write(3,*) t_new(u), S_new(u), I_new(u) 

     jump = jump+1 
     u = u+1 
    end if 

    ! write(3,*) t_new(i), S_new(i), I_new(i) 

    !N = S_now + I_now ! update population post event 

    if(I_now .le. 0) then ! if extinct, terminate 
     print *, "extinct" 
     goto 2 
    end if 

    end do 

2  end program Gillespie 

ich alle Eingaben zu schätzen wissen. Vielen Dank.

Answer 1

Ihr Anruf

call random_seed(seed) 

ist seltsam. Ich dachte, es sollte nicht ohne ein Schlüsselwortargument erlaubt werden, aber es fragt tatsächlich nach der Größe des zufälligen Seedarrays.

Für eine richtige Art und Weise Samen siehe das Beispiel in

https://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html