! output file has the following structure: ! column 1 -> generation number in units of 100 generations ! columns 2,3: fraction of demes with at least one individual in habitat 1 and 2 at given generation. ! columns 4,5: average allele frequency in habitat 1 and 2 at given generation; averaging performed over all equal-effect loci and over all demes in a given habitat. ! columns 6,7: average heterozygosity in habitat 1 and 2 at given generation; averaging performed over all equal-effect loci and over all demes in a given habitat. ! columns 8,9: average population size of demes in habitat 1 and 2 at given generation; averaging performed over all demes in a given habitat. ! columns 10,11: (average genic variance)/(average genetic variance) in habitat 1 and 2; 1 minus this value indicates the fraction of variance due to LD. program wright_fisher integer(KIND=8) L,K,ndemes,ndemes_hab1,k6,histories,i,j,j1 integer(KIND=8) k5,tsteps,mig,nnew,c3,j3,j6,c1 real(KIND=8) sel1,sel2,migrate,rmax,rho,x,y,pini,t1,t2 real(KIND=8) temp,temp1,temp2,theta,temp3,loglimit real(KIND=8) temp4,temp5,temp6,t2a,t3a,t4a,t1a,x2 integer(KIND=8),allocatable, dimension(:,:,:)::n1_seq,n2_seq integer(KIND=8), allocatable, dimension(:,:)::mig_pool real(KIND=8), allocatable, dimension(:,:)::fit_indivs integer(KIND=8), allocatable, dimension(:):: pop_deme,pop_deme2 real(KIND=8), allocatable, dimension(:):: freq_hab1,pop_hab1 real(KIND=8), allocatable, dimension(:):: gdiv_hab1,gdiv_hab2 real(KIND=8), allocatable, dimension(:):: psurv_hab1,psurv_hab2 real(KIND=8), allocatable, dimension(:):: LDvar_hab1,LDvar_hab2 real(KIND=8), allocatable, dimension(:)::freq_hab2,pop_hab2 real(KIND=8), allocatable, dimension(:)::fit_demes integer(KIND=8) ZBQLPOI K=2000 ! baseline carrying capacity in each deme L=10 ! no. of loci under divergent selection ndemes=100 ! no. of demes in metapopulation rho=0.3d+0 ! fraction of demes in rare habitat ndemes_hab1=nint(ndemes*rho) rmax=40.0d+0/K ! baseline growth rate r0=zeta/K; input zeta sel1=0.2d+0*rmax ! s1=S1*r0; input scaled selective effect S1 per locus in habitat 1 sel2=0.2d+0*rmax ! s2=S2*r0; input scaled selective effect S2 per locus in habitat 2 migrate=rmax*(0.023D+0) ! m=M*r0; input scaled migration rate M pini=0.5d+0 ! initial allele frequency in each habitat at each locus histories=1 ! number of replicates tsteps=300000 ! number of generations simulated theta=1.0D+0/histories loglimit=-log(tiny(double)) allocate(n1_seq(L,K+2000,ndemes),n2_seq(L,2*K+4000,ndemes)) allocate(pop_deme(ndemes),pop_deme2(ndemes)) allocate(fit_indivs(ndemes,K+1000)) allocate(freq_hab1(tsteps),freq_hab2(tsteps),fit_demes(ndemes)) allocate(gdiv_hab1(tsteps),gdiv_hab2(tsteps)) allocate(LDvar_hab1(tsteps),LDvar_hab2(tsteps)) allocate(pop_hab1(tsteps),pop_hab2(tsteps)) allocate(psurv_hab1(tsteps),psurv_hab2(tsteps)) allocate(mig_pool(L,nint(ndemes*K*migrate*100)) ) do i=1,tsteps freq_hab1(i)=0.0d+0 freq_hab2(i)=0.0d+0 gdiv_hab1(i)=0.0d+0 gdiv_hab2(i)=0.0d+0 pop_hab1(i)=0.0d+0 pop_hab2(i)=0.0d+0 ! var_hab1(i)=0.0d+0 ! var_hab2(i)=0.0d+0 LDvar_hab1(i)=0.0d+0 LDvar_hab2(i)=0.0d+0 psurv_hab1(i)=0.0d+0 psurv_hab2(i)=0.0d+0 enddo call cpu_time(t1) call init_random_seed do k6=1,histories print*, k6, "aa" ! create initial condition do j1=1,ndemes_hab1 do i=1,K do j=1,L call random_number(x) if (xnint(ndemes*K*migrate*10000))then print*, "alarm1" endif mig_pool(j3,c3)=n1_seq(j3,pop_deme(i)+1-j,i) enddo c3=c3+1 enddo pop_deme(i)=pop_deme(i)-mig enddo c3=c3-1 if (mod(k5,100)==1)then temp=0.0D+0 do i=1,c3 do j3=1,L temp=temp+mig_pool(j3,i) enddo enddo endif ! decide on number of migrants into demes and their genomes. temp2=(1.0d+0*c3)/ndemes do i=1,ndemes mig=ZBQLPOI(temp2) do j=pop_deme(i)+1,pop_deme(i)+mig call random_number(x2) j3=nint(aint(1+c3*x2)) if (j3>c3)then print*, "severe alarm" endif do j6=1,L n1_seq(j6,j,i)=mig_pool(j6,j3) enddo enddo pop_deme(i)=pop_deme(i)+mig enddo ! selection ! calculate number of delet. mutations carried by each individual in each deme. do i=1,ndemes do j=1,pop_deme(i) fit_indivs(i,j)=0.0d+0 do j3=1,L fit_indivs(i,j)=fit_indivs(i,j)+n1_seq(j3,j,i) enddo enddo enddo ! calculate fitness of each indiv and fitness of deme and new pop. size. do i=1,ndemes_hab1 fit_demes(i)=0.0D+0 if (pop_deme(i)>0)then temp=0.0d+0 do j=1,pop_deme(i) temp1=sel2*fit_indivs(i,j) if (temp10)then temp=0.0d+0 do j=1,pop_deme(i) temp1=sel1*(L-fit_indivs(i,j)) if (temp1fit_indivs(i,j3))then go to 5 else do j6=1,L n2_seq(j6,j,i)=n1_seq(j6,j3,i) temp=temp+(n1_seq(j6,j3,i)/(1.0D+0*L)) enddo endif enddo ! if (pop_deme(i)>0)then ! temp=temp/(2*pop_deme(i)) ! endif ! create new offspring by recombination temp1=0.0D+0 do j=1,pop_deme2(i) do j6=1,L call random_number(y) if (y<0.5D+0)then n1_seq(j6,j,i)=n2_seq(j6,2*j-1,i) else n1_seq(j6,j,i)=n2_seq(j6,2*j,i) endif temp1=temp1+(n1_seq(j6,j,i)/(1.0D+0*L)) enddo enddo ! if (pop_deme(i)>0)then ! temp1=temp1/pop_deme(i) ! endif enddo do i=1,ndemes pop_deme(i)=pop_deme2(i) enddo ! calculate statistics temp3=0.0d+0 temp2=0.0D+0 temp6=0.0d+0 t4a=0.0d+0 t3a=0.0d+0 c1=0 do i=1,ndemes_hab1 temp4=0.0D+0 temp1=0.0D+0 if (pop_deme(i)>0)then c1=c1+1 do j6=1,L temp5=0.0d+0 do j=1,pop_deme(i) temp5=temp5+n1_seq(j6,j,i) enddo temp5=(temp5*1.0d+0)/pop_deme(i) temp1=temp1+temp5 temp4=temp4+(temp5*(1.0D+0-temp5)) enddo temp1=temp1/L temp4=temp4/L t1a=0.0d+0 t2a=0.0D+0 do j=1,pop_deme(i) temp5=0.0d+0 do j6=1,L temp5=temp5+n1_seq(j6,j,i) enddo t1a=t1a+(temp5*1.0D+0) t2a=t2a+(temp5*temp5*1.0d+0) enddo t1a=t1a/pop_deme(i) t2a=t2a/pop_deme(i) t2a=t2a-(t1a*t1a) if (t2a>0.0d+0)then t3a=t3a+((temp4*L)/t2a) endif endif temp3=temp3+temp1 temp6=temp6+temp4 temp2=temp2+pop_deme(i) enddo if (c1>0)then temp3=temp3/c1 temp6=temp6/c1 t4a=t3a/c1 psurv_hab1(k5)=psurv_hab1(k5)+theta endif temp2=temp2/ndemes_hab1 freq_hab1(k5)=freq_hab1(k5)+(theta*temp3) pop_hab1(k5)=pop_hab1(k5)+(theta*temp2) gdiv_hab1(k5)=gdiv_hab1(k5)+(theta*temp6) LDvar_hab1(k5)=LDvar_hab1(k5)+(theta*t4a) ! print*, k5, temp2,temp3,temp6,t4a ! if (temp3<0.000000001D+0)then ! print*, k5, temp3,temp2,"hab1" ! endif temp3=0.0d+0 temp2=0.0D+0 temp6=0.0d+0 t4a=0.0d+0 t3a=0.0d+0 c1=0 do i=ndemes_hab1+1,ndemes temp4=0.0D+0 temp1=0.0D+0 if (pop_deme(i)>0)then c1=c1+1 do j6=1,L temp5=0.0d+0 do j=1,pop_deme(i) temp5=temp5+n1_seq(j6,j,i) enddo temp5=(temp5*1.0d+0)/pop_deme(i) temp1=temp1+temp5 temp4=temp4+(temp5*(1.0D+0-temp5)) enddo temp1=temp1/L temp4=temp4/L t1a=0.0d+0 t2a=0.0D+0 do j=1,pop_deme(i) temp5=0.0d+0 do j6=1,L temp5=temp5+n1_seq(j6,j,i) enddo t1a=t1a+(temp5*1.0D+0) t2a=t2a+(temp5*temp5*1.0d+0) enddo t1a=t1a/pop_deme(i) t2a=t2a/pop_deme(i) t2a=t2a-(t1a*t1a) if (t2a>0.0d+0)then t3a=t3a+((temp4*L)/t2a) endif endif temp3=temp3+temp1 temp6=temp6+temp4 temp2=temp2+pop_deme(i) enddo if (c1>0)then temp3=temp3/c1 temp6=temp6/c1 t4a=t3a/c1 psurv_hab2(k5)=psurv_hab2(k5)+theta endif temp2=temp2/(ndemes-ndemes_hab1) freq_hab2(k5)=freq_hab2(k5)+(theta*temp3) pop_hab2(k5)=pop_hab2(k5)+(theta*temp2) gdiv_hab2(k5)=gdiv_hab2(k5)+(theta*temp6) LDvar_hab2(k5)=LDvar_hab2(k5)+(theta*t4a) ! if (temp3<0.000000001D+0)then ! print*, k5, temp3,temp2,"hab2" ! endif if (mod(k5,5000)==0)then print*,k5 temp=0.0D+0 temp1=0.0d+0 do j6=1,ndemes temp1=temp1+pop_deme(j6) do i=1,pop_deme(j6) do j3=1,L temp=temp+(n1_seq(j3,i,j6)/(1.0D+0*L)) enddo enddo enddo temp=temp/temp1 print*, k5-1,temp, "deme average" endif enddo temp1=0.0D+0 temp2=0.0D+0 do i=1,ndemes_hab1 ! print*, i, pop_deme(i),pop_deme(i+ndemes_hab1) temp1=temp1+pop_deme(i) temp2=temp2+pop_deme(i+ndemes_hab1) enddo print*, k6,temp1,temp2 enddo print*, "b" call cpu_time(t2) print*, t2-t1 19 open(41,file='L10_rho_0.7_S0.2_M0.023_r0.02_1history.txt') do i=1,tsteps write(41,9) i,psurv_hab1(i),psurv_hab2(i),freq_hab1(i), c freq_hab2(i),gdiv_hab1(i),gdiv_hab2(i), c pop_hab1(i)/K,pop_hab2(i)/K,LDvar_hab1(i),LDvar_hab2(i) 9 Format(I8,1X,F15.9,1X,F15.9,1X,F15.9,1X,F15.9,1X,F15.9,1X,F15.9, c 1X,F15.9,1X,F15.9,1X,F15.9,1X,F15.9) enddo close(41) end SUBROUTINE init_random_seed() INTEGER :: z, m, clock INTEGER, DIMENSION(:), ALLOCATABLE :: seed, seedold CALL RANDOM_SEED(size = m) ALLOCATE(seed(m), seedold(m)) CALL SYSTEM_CLOCK(COUNT=clock) seed = clock + 37 * (/ (z - 1, z = 1, m) /) CALL RANDOM_SEED(PUT = seed) CALL RANDOM_SEED(GET = seedold) DEALLOCATE(seed) END SUBROUTINE init_random_seed BLOCK DATA ZBQLBD01 * * Initializes seed array etc. for random number generator. * The values below have themselves been generated using the * NAG generator. * COMMON /ZBQL0001/ ZBQLIX,B,C DOUBLE PRECISION ZBQLIX(43),B,C INTEGER(KIND=8) I DATA (ZBQLIX(I),I=1,43) /8.001441D7,5.5321801D8, +1.69570999D8,2.88589940D8,2.91581871D8,1.03842493D8, +7.9952507D7,3.81202335D8,3.11575334D8,4.02878631D8, +2.49757109D8,1.15192595D8,2.10629619D8,3.99952890D8, +4.12280521D8,1.33873288D8,7.1345525D7,2.23467704D8, +2.82934796D8,9.9756750D7,1.68564303D8,2.86817366D8, +1.14310713D8,3.47045253D8,9.3762426D7 ,1.09670477D8, +3.20029657D8,3.26369301D8,9.441177D6,3.53244738D8, +2.44771580D8,1.59804337D8,2.07319904D8,3.37342907D8, +3.75423178D8,7.0893571D7 ,4.26059785D8,3.95854390D8, +2.0081010D7,5.9250059D7,1.62176640D8,3.20429173D8, +2.63576576D8/ DATA B / 4.294967291D9 / DATA C / 0.0D0 / END ****************************************************************** ****************************************************************** ****************************************************************** SUBROUTINE ZBQLINI(SEED) ****************************************************************** * To initialize the random number generator - either * repeatably or nonrepeatably. Need double precision * variables because integer storage can't handle the * numbers involved ****************************************************************** * ARGUMENTS * ========= * SEED (integer, input). User-input number which generates * elements of the array ZBQLIX, which is subsequently used * in the random number generation algorithm. If SEED=0, * the array is seeded using the system clock if the * FORTRAN implementation allows it. ****************************************************************** * PARAMETERS * ========== * LFLNO (integer). Number of lowest file handle to try when * opening a temporary file to copy the system clock into. * Default is 80 to keep out of the way of any existing * open files (although the program keeps searching till * it finds an available handle). If this causes problems, * (which will only happen if handles 80 through 99 are * already in use), decrease the default value. ****************************************************************** INTEGER(KIND=8) LFLNO PARAMETER (LFLNO=80) ****************************************************************** * VARIABLES * ========= * SEED See above * ZBQLIX Seed array for the random number generator. Defined * in ZBQLBD01 * B,C Used in congruential initialisation of ZBQLIX * SS,MM,} System clock secs, mins, hours and days * HH,DD } * FILNO File handle used for temporary file * INIT Indicates whether generator has already been initialised * INTEGER(KIND=8) SEED,SS,MM,HH,DD,FILNO,I INTEGER(KIND=8) INIT DOUBLE PRECISION ZBQLIX(43),B,C DOUBLE PRECISION TMPVAR1,DSS,DMM,DHH,DDD COMMON /ZBQL0001/ ZBQLIX,B,C SAVE INIT * * Ensure we don't call this more than once in a program * IF (INIT.GE.1) THEN IF(INIT.EQ.1) THEN WRITE(*,1) INIT = 2 ENDIF RETURN ELSE INIT = 1 ENDIF * * If SEED = 0, cat the contents of the clock into a file * and transform to obtain ZQBLIX(1), then use a congr. * algorithm to set remaining elements. Otherwise take * specified value of SEED. * *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> *>>>>>>> NB FOR SYSTEMS WHICH DO NOT SUPPORT THE >>>>>>> *>>>>>>> (NON-STANDARD) 'CALL SYSTEM' COMMAND, >>>>>>> *>>>>>>> THIS WILL NOT WORK, AND THE FIRST CLAUSE >>>>>>> *>>>>>>> OF THE FOLLOWING IF BLOCK SHOULD BE >>>>>>> *>>>>>>> COMMENTED OUT. >>>>>>> *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IF (SEED.EQ.0) THEN *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> *>>>>>>> COMMENT OUT FROM HERE IF YOU DON'T HAVE >>>>>>> *>>>>>>> 'CALL SYSTEM' CAPABILITY ... >>>>>>> *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> CALL SYSTEM(' date +%S%M%H%j > zbql1234.tmp') * * Try all file numbers for LFLNO to 999 * FILNO = LFLNO 10 OPEN(FILNO,FILE='zbql1234.tmp',ERR=11) GOTO 12 11 FILNO = FILNO + 1 IF (FILNO.GT.999) THEN WRITE(*,2) RETURN ENDIF GOTO 10 12 READ(FILNO,'(3(I2),I3)') SS,MM,HH,DD CLOSE(FILNO) CALL SYSTEM('rm zbql1234.tmp') DSS = DINT((DBLE(SS)/6.0D1) * B) DMM = DINT((DBLE(MM)/6.0D1) * B) DHH = DINT((DBLE(HH)/2.4D1) * B) DDD = DINT((DBLE(DD)/3.65D2) * B) TMPVAR1 = DMOD(DSS+DMM+DHH+DDD,B) *<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< *<<<<<<<< ... TO HERE (END OF COMMENTING OUT FOR <<<<<<< *<<<<<<<< USERS WITHOUT 'CALL SYSTEM' CAPABILITY <<<<<<< *<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ELSE TMPVAR1 = DMOD(DBLE(SEED),B) ENDIF ZBQLIX(1) = TMPVAR1 DO 100 I = 2,43 TMPVAR1 = ZBQLIX(I-1)*3.0269D4 TMPVAR1 = DMOD(TMPVAR1,B) ZBQLIX(I) = TMPVAR1 100 CONTINUE 1 FORMAT(//5X,'****WARNING**** You have called routine ZBQLINI ', +'more than',/5X,'once. I''m ignoring any subsequent calls.',//) 2 FORMAT(//5X,'**** ERROR **** In routine ZBQLINI, I couldn''t', +' find an',/5X, +'available file number. To rectify the problem, decrease the ', +'value of',/5X, +'the parameter LFLNO at the start of this routine (in file ', +'randgen.f)',/5X, +'and recompile. Any number less than 100 should work.') END ****************************************************************** FUNCTION ZBQLU01(DUMMY) * * Returns a uniform random number between 0 & 1, using * a Marsaglia-Zaman type subtract-with-borrow generator. * Uses double precision, rather than integer, arithmetic * throughout because MZ's integer constants overflow * 32-bit integer storage (which goes from -2^31 to 2^31). * Ideally, we would explicitly truncate all integer * quantities at each stage to ensure that the double * precision representations do not accumulate approximation * error; however, on some machines the use of DNINT to * accomplish this is *seriously* slow (run-time increased * by a factor of about 3). This double precision version * has been tested against an integer implementation that * uses long integers (non-standard and, again, slow) - * the output was identical up to the 16th decimal place * after 10^10 calls, so we're probably OK ... * DOUBLE PRECISION ZBQLU01,DUMMY,B,C,ZBQLIX(43),X,B2,BINV INTEGER(KIND=8) CURPOS,ID22,ID43 COMMON /ZBQL0001/ ZBQLIX,B,C SAVE /ZBQL0001/ SAVE CURPOS,ID22,ID43 DATA CURPOS,ID22,ID43 /1,22,43/ B2 = B BINV = 1.0D0/B 5 X = ZBQLIX(ID22) - ZBQLIX(ID43) - C IF (X.LT.0.0D0) THEN X = X + B C = 1.0D0 ELSE C = 0.0D0 ENDIF ZBQLIX(ID43) = X * * Update array pointers. Do explicit check for bounds of each to * avoid expense of modular arithmetic. If one of them is 0 the others * won't be * CURPOS = CURPOS - 1 ID22 = ID22 - 1 ID43 = ID43 - 1 IF (CURPOS.EQ.0) THEN CURPOS=43 ELSEIF (ID22.EQ.0) THEN ID22 = 43 ELSEIF (ID43.EQ.0) THEN ID43 = 43 ENDIF * * The integer arithmetic there can yield X=0, which can cause * problems in subsequent routines (e.g. ZBQLEXP). The problem * is simply that X is discrete whereas U is supposed to * be continuous - hence if X is 0, go back and generate another * X and return X/B^2 (etc.), which will be uniform on (0,1/B). * IF (X.LT.BINV) THEN B2 = B2*B GOTO 5 ENDIF ZBQLU01 = X/B2 END ****************************************************************** FUNCTION ZBQLUAB(A,B) * * Returns a random number uniformly distributed on (A,B) * DOUBLE PRECISION A,B,ZBQLU01,ZBQLUAB * * Even if A > B, this will work as B-A will then be -ve * IF (A.NE.B) THEN ZBQLUAB = A + ( (B-A)*ZBQLU01(0.0D0) ) ELSE ZBQLUAB = A WRITE(*,1) ENDIF 1 FORMAT(/5X,'****WARNING**** (function ZBQLUAB) Upper and ', +'lower limits on uniform',/5X,'distribution are identical',/) END ****************************************************************** FUNCTION ZBQLEXP(MU) * * Returns a random number exponentially distributed with * mean MU * DOUBLE PRECISION MU,ZBQLEXP,ZBQLU01 ZBQLEXP = 0.0D0 IF (MU.LT.0.0D0) THEN WRITE(*,1) RETURN ENDIF ZBQLEXP = -DLOG(ZBQLU01(0.0D0))*MU 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLEXP',/) END ****************************************************************** FUNCTION ZBQLNOR(MU,SIGMA) * * Returns a random number Normally distributed with mean * MU and standard deviation |SIGMA|, using the Box-Muller * algorithm * DOUBLE PRECISION THETA,R,ZBQLNOR,ZBQLU01,PI,MU,SIGMA DOUBLE PRECISION SPARE INTEGER(KIND=8) STATUS SAVE STATUS,SPARE,PI DATA STATUS /-1/ IF (STATUS.EQ.-1) PI = 4.0D0*DATAN(1.0D0) IF (STATUS.LE.0) THEN THETA = 2.0D0*PI*ZBQLU01(0.0D0) R = DSQRT( -2.0D0*DLOG(ZBQLU01(0.0D0)) ) ZBQLNOR = (R*DCOS(THETA)) SPARE = (R*DSIN(THETA)) STATUS = 1 ELSE ZBQLNOR = SPARE STATUS = 0 ENDIF ZBQLNOR = MU + (SIGMA*ZBQLNOR) END ****************************************************************** FUNCTION ZBQLBIN(N,P) * * Returns a random number binomially distributed (N,P) * DOUBLE PRECISION P,ZBQLBET1 DOUBLE PRECISION PP,PPP,G,Y,TINY INTEGER(KIND=8) N,ZBQLBIN,ZBQLGEO,IZ,NN TINY = 1.0D-8 ZBQLBIN = 0 IF (.NOT.( (P.GE.0.0D0).AND.(P.LE.1.0D0) )) THEN WRITE(*,1) RETURN ELSEIF(N.LE.0) THEN WRITE(*,1) RETURN ENDIF * * First step: if NP > 10, say, things will be expensive, and * we can get into the right ballpark by guessing a value for * ZBQLBIN (IZ, say), and simulating Y from a Beta distribution * with parameters IZ and NN-IZ+1 (NN starts off equal to N). * If Y is less than PP (which starts off as P) then the IZth order * statistic from NN U(0,1) variates is less than P, and we know * that there are at least IZ successes. In this case we focus on * the remaining (NN-IZ) order statistics and count how many are * less than PP, which is binomial (NN-IZ,(PP-Y)/(1-Y)). * Otherwise, if Y is greater than PP there must be less * than IZ successes, so we can count the number of order statistics * under PP, which is binomial (IZ-1,P/Y). When we've got NN*PP * small enough, we go to the next stage of the algorithm and * generate the final bits directly. * NN = N PP = P 10 IZ = INT(DBLE(NN)*PP) + 1 IF ( (IZ.GT.10).AND.(IZ.LT.NN-10) ) THEN Y = ZBQLBET1(DBLE(IZ),DBLE(NN-IZ+1)) IF (Y.LT.PP) THEN ZBQLBIN = ZBQLBIN + IZ NN = NN - IZ PP = (PP-Y) / (1.0D0-Y) ELSE NN = IZ-1 PP = PP/Y ENDIF GOTO 10 ENDIF * * PP is the probability of the binomial we're currently * simulating from. For the final part, we simulate either number * of failures or number of success, depending which is cheaper. * 20 IF (PP.GT.0.5) THEN PPP = 1.0D0-PP ELSE PPP = PP ENDIF G = 0 IZ = 0 * * ZBQLGEO falls over for miniscule values of PPP, so ignore these * (tiny probability of any successes in this case, anyway) * IF (PPP.GT.TINY) THEN 30 G = G + ZBQLGEO(PPP) IF (G.LE.NN) THEN IZ = IZ + 1 GOTO 30 ENDIF ENDIF IF (PP.GT.0.5) IZ = NN - IZ ZBQLBIN = ZBQLBIN + IZ 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLBIN',/) END ****************************************************************** FUNCTION ZBQLGEO(P) * * Returns a random number geometrically distributed with * parameter P ie. mean 1/P * DOUBLE PRECISION P,ZBQLU01,U,TINY INTEGER(KIND=8) ZBQLGEO TINY = 1.0D-12 ZBQLGEO = 0 IF (.NOT.( (P.GE.0.0D0).AND.(P.LE.1.0D0) )) THEN WRITE(*,1) RETURN ENDIF IF (P.GT.0.9D0) THEN 10 ZBQLGEO = ZBQLGEO + 1 U = ZBQLU01(0.0D0) IF (U.GT.P) GOTO 10 ELSE U = ZBQLU01(0.0D0) * * For tiny P, 1-p will be stored inaccurately and log(1-p) may * be zero. In this case approximate log(1-p) by -p * IF (P.GT.TINY) THEN ZBQLGEO = 1 + INT( DLOG(U)/DLOG(1.0D0-P) ) ELSE ZBQLGEO = 1 + INT(-DLOG(U)/P) ENDIF ENDIF 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLGEO',/) END ****************************************************************** FUNCTION ZBQLPOI(MU) * * Returns a random number Poisson distributed with mean MU * DOUBLE PRECISION ZBQLU01,X,Y,MU,PI DOUBLE PRECISION ZBQLLG,ZBQLGAM,MU1,TMP1,TMP2,T INTEGER(KIND=8) ZBQLPOI,ZBQLBIN,K,INIT SAVE INIT,PI DATA INIT /0/ IF (INIT.EQ.0) THEN PI = 4.0D0*DATAN(1.0D0) INIT = 1 ENDIF ZBQLPOI = 0 IF (MU.LT.0.0D0) THEN WRITE(*,1) RETURN ENDIF * * For small MU, generate exponentials till their sum exceeds 1 * (equivalently, uniforms till their product falls below e^-MU) * IF (MU.LE.1.0D3) THEN MU1 = MU * * For values of MU less than 1000, use order statistics - the Kth * event in a Poisson process of rate MU has a Gamma distribution * with parameters (MU,K); if it's greater than 1 we know that there * are less than K events in (0,1) (and the exact number is binomial) * and otherwise the remaining number is another Poisson. Choose K so * that we'll get pretty close to 1 in the first go but are unlikely * to overshoot it. * 19 IF (MU1.GT.1.0D1) THEN K = INT(MU1-DSQRT(MU1)) Y = ZBQLGAM(DBLE(K),MU1) IF (Y.GT.1.0D0) THEN ZBQLPOI = ZBQLPOI + ZBQLBIN(K-1,(1.0D0/Y)) RETURN ENDIF ZBQLPOI = ZBQLPOI + K MU1 = MU1 * (1.0D0-Y) GOTO 19 ENDIF Y = DEXP(-MU1) X = 1.0D0 20 X = X*ZBQLU01(0.0D0) IF (X.GT.Y) THEN ZBQLPOI = ZBQLPOI + 1 GOTO 20 ENDIF * * For really huge values of MU, use rejection sampling as in * Press et al (1992) - large numbers mean some accuracy may be * lost, but it caps the execution time. * ELSE TMP1 = DSQRT(2.0D0*MU) TMP2 = ZBQLLG(MU+1.0D0)-(MU*DLOG(MU)) 30 Y = DTAN(PI*ZBQLU01(0.0D0)) ZBQLPOI = INT(MU + (TMP1*Y)) IF (ZBQLPOI.LT.0) GOTO 30 X = DBLE(ZBQLPOI) T = (X*DLOG(MU)-ZBQLLG(X+1.0D0)) + TMP2 IF (DABS(T).LT.1.0D2) THEN T = 0.9D0*(1.0D0+(Y*Y))*DEXP(T) IF (ZBQLU01(0.0D0).GT.T) GOTO 30 ELSE T = DLOG(0.9D0*(1.0D0+(Y*Y))) + T IF (DLOG(ZBQLU01(0.0D0)).GT.T) GOTO 30 ENDIF ENDIF 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLPOI',/) END ****************************************************************** FUNCTION ZBQLGAM(G,H) * * Returns a random number with a gamma distribution with mean * G/H and variance G/(H^2). (ie. shape parameter G & scale * parameter H) * DOUBLE PRECISION C,D,R,ZBQLGAM,ZBQLU01,G,H,A,z1,z2,B1,B2,M DOUBLE PRECISION U1,U2,U,V,TEST,X double precision c1,c2,c3,c4,c5,w ZBQLGAM = 0.0D0 IF ( (G.LE.0.0D0).OR.(H.LT.0.0D0) ) THEN WRITE(*,1) RETURN ENDIF IF (G.LT.1.0D0) THEN 889 u=ZBQLU01(0.0d0) v=ZBQLU01(0.0d0) if (u.gt.exp(1.0d0)/(g+exp(1.0d0))) goto 891 ZBQLGAM=((g+exp(1.0d0))*u/exp(1.0d0))**(1.0d0/g) if (v.gt.exp(-ZBQLGAM)) then goto 889 else goto 892 endif 891 ZBQLGAM=-log((g+exp(1.0d0))*(1.0d0-u)/(g*exp(1.0d0))) if (v.gt.ZBQLGAM**(g-1.0)) goto 889 892 ZBQLGAM=ZBQLGAM/h RETURN ELSEIF (G.LT.2.0D0) THEN M = 0.0D0 elseif (g.gt.10.0d0) then c1=g-1.0d0 c2=(g-1.0d0/(6.0d0*g))/c1 c3=2.0d0/c1 c4=c3+2.0d0 c5=1.0d0/sqrt(g) 777 u=ZBQLU01(0.0d0) v=ZBQLU01(0.0d0) if (g.gt.2.50d0) then u=v+c5*(1.0d0-1.860d0*u) endif if (u.le.0.0d0.or.u.ge.1.0d0) goto 777 w=c2*v/u if (c3*u+w+1.0d0/w.le.c4) goto 778 if (c3*log(u)-log(w)+w.ge.1.0d0) goto 777 778 ZBQLGAM=c1*w/h return ELSE M = -(G-2.0D0) ENDIF R = 0.50D0 a = ((g-1.0d0)/exp(1.0d0))**((g-1.0d0)/(r+1.0d0)) C = (R*(M+G)+1.0D0)/(2.0D0*R) D = M*(R+1.0D0)/R z1 = C-DSQRT(C*C-D) * * On some systems (e.g. g77 0.5.24 on Linux 2.4.24), C-DSQRT(C*C) * is not exactly zero - this needs trapping if negative. * IF ((Z1-M.LT.0.0D0).AND.(Z1-M.GT.-1.0D-12)) Z1 = M z2 = C+DSQRT(C*C-D) B1=(z1*(z1-M)**(R*(G-1.0D0)/(R+1.0D0)))*DEXP(-R*(z1-M)/(R+1.0D0)) B2=(z2*(z2-M)**(R*(G-1.0D0)/(R+1.0D0)))*DEXP(-R*(z2-M)/(R+1.0D0)) 50 U1=ZBQLU01(0.0D0) U2=ZBQLU01(0.0D0) U=A*U1 V=B1+(B2-B1)*U2 X=V/(U**R) IF (X.LE.M) GOTO 50 TEST = ((X-M)**((G-1)/(R+1)))*EXP(-(X-M)/(R+1.0D0)) IF (U.LE.TEST) THEN ZBQLGAM = (X-M)/H ELSE GOTO 50 ENDIF 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLGAM',/5X, '(both parameters must be positive)',/) END *************************************************************** FUNCTION ZBQLBET1(NU1,NU2) * * Returns a random number, beta distributed with degrees * of freedom NU1 and NU2. * DOUBLE PRECISION NU1,NU2,ZBQLGAM,ZBQLBET1,ZBQLU01,X1,X2 ZBQLBET1 = 0.0D0 IF ( (NU1.LE.0.0).OR.(NU2.LE.0.0) ) THEN WRITE(*,1) RETURN ENDIF * * If parameters are too small, gamma subroutine tends to return zero * as all the probability goes to the origin and we get rounding * errors, even with double precision. In this case, we use Johnk's * method, suitably scaled to avoid rounding errors as much as possible. * IF ( (NU1.LT.0.9D0).AND.(NU2.LT.0.9D0) ) THEN 10 X1 = ZBQLU01(0.0D0) X2 = ZBQLU01(0.0D0) IF ( (X1**(1.0D0/NU1))+(X2**(1.0D0/NU2)).GT.1.0D0) GOTO 10 X1 = (DLOG(X2)/NU2) - (DLOG(X1)/NU1) ZBQLBET1 = (1.0D0 + DEXP(X1))**(-1) IF (ZBQLBET1.GT.1.0D0) GOTO 10 ELSE X1 = ZBQLGAM(NU1,1.0D0) X2 = ZBQLGAM(NU2,1.0D0) ZBQLBET1 = X1/(X1+X2) ENDIF RETURN 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLBET1',/5X, +'(both degrees of freedom must be positive)',/) END *************************************************************** FUNCTION ZBQLWEI(A,B) * * Returns a random number, Weibull distributed with shape parameter * A and location parameter B, i.e. density is * f(x) = ( A/(B**A) ) * x**(A-1) * EXP( -(x/B)**A ) * DOUBLE PRECISION A,B,ZBQLU01,ZBQLWEI,U ZBQLWEI = 0.0D0 IF ( (A.LE.0.0).OR.(B.LE.0.0) ) THEN WRITE(*,1) RETURN ENDIF U = ZBQLU01(0.0D0) ZBQLWEI = B * ( (-DLOG(U))**(1.0D0/A) ) 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLWEI',/5X, +'(both parameters must be positive)',/) END *************************************************************** FUNCTION ZBQLNB(R,P) * * Returns a pseudo-random number according to a Negative * Binomial distribution with parameters (R,P). NB these are * both DOUBLE - it copes with non-integer R as well. The * form of the distribution is *not* the no. of trials to * the Rth success - see documentation for full spec. * DOUBLE PRECISION R,P,ZBQLGAM,Y INTEGER(KIND=8) ZBQLNB,ZBQLPOI ZBQLNB = 0 IF ( (R.LE.0.0D0).OR.(P.LE.0.0D0).OR.(P.GE.1.0D0) ) THEN WRITE(*,1) RETURN ENDIF Y = ZBQLGAM(R,1.0D0) Y = Y*P/(1.0D0-P) ZBQLNB = ZBQLPOI(Y) 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLNB') END *************************************************************** FUNCTION ZBQLPAR(A,B) * * Returns a random number, Pareto distributed with parameters * A and B. The density is A*(B**A) / (B+X)**(A+1) for X > 0. * (this is slightly nonstandard - see documentation in * randgen.txt). The algorithm is straightforward - it uses the * inverse CDF method. * DOUBLE PRECISION A,B,ZBQLPAR,ZBQLU01,U ZBQLPAR = 0.0D0 IF ( (A.LE.0.0D0).OR.(B.LE.0.0D0) ) THEN WRITE(*,1) RETURN ENDIF U = ZBQLU01(0.0D0) ZBQLPAR = B * (U**(-1.0D0/A)-1.0D0) 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLPAR',/5X, +'(both parameters must be positive)',/) END *************************************************************** FUNCTION ZBQLLG(X) * * Returns log(G(X)) where G is the Gamma function. The algorithm is * that given in Press et al (1992), Section 6.1, although this * version also allows for arguments less than 1. * DOUBLE PRECISION X,Z,Z2,ZBQLLG,PI,RLN2P,C(0:6),TMP,SUM INTEGER(KIND=8) INIT,I SAVE INIT,C,RLN2P,PI DATA INIT /0/ DATA (C(I),I=0,6) / + 1.000000000190015D0,76.18009172947146D0, + -86.50532032941677D0,24.01409824083091D0, + -1.231739572450155D0,0.1208650973866179D-2, + -0.5395239384953D-5/ IF (INIT.EQ.0) THEN PI = 4.0D0*DATAN(1.0D0) RLN2P = 0.5D0*DLOG(2.0D0*PI) INIT = 1 ENDIF * * Compute for x > 1, then use transformation if necessary. Z is * our working argument. * IF (X.GE.1.0D0) THEN Z = X ELSE Z = 2.0D0-X Z2 = 1.0D0-X ENDIF IF (DABS(Z-1.0D0).LT.1.0D-12) THEN ZBQLLG = 0.0D0 RETURN ENDIF TMP = Z + 4.5D0 TMP = ( (Z-0.5D0)*DLOG(TMP) ) - TMP + RLN2P SUM = C(0) DO 50 I=1,6 SUM = SUM + (C(I)/(Z+DBLE(I-1))) 50 CONTINUE ZBQLLG = TMP + DLOG(SUM) * * Transformation required if X<1 * IF (X.LT.1.0D0) THEN TMP = PI*Z2 ZBQLLG = DLOG(TMP/DSIN(TMP)) - ZBQLLG ENDIF END