C Metropolis with poisson proposal: PROGRAM METRONET C*******variables: INTEGER IDUM PARAMETER(ITERA=50000) INTEGER K,N,I,W,ROW,COL,ELEM,COUNTY, / X1(100),X2(100),Y(100),COUNT(100), / X1NEW(100),GENERATEDX1(ITERA,100), / GENERATEDX2(ITERA,100),A1(100,100),A2(100,100), / INVA1(100,100),RECORD REAL*8 LAMBDAS(100),GENERATEDLAMBDAS(ITERA,100) REAL*4 GAMDEV,RAN2 C*******interactive input: C PRINT *,'ENTER THE NUMBER OF OBSERVATIONS:' C READ *,K C C PRINT *,'ENTER THE OBSERVATION VALUES:' C READ *, (Y(I),I=1,K) C C C PRINT *,'ENTER THE NUMBER OF O-D PAIRS:' C READ *,N C C PRINT *,'ENTER THE INITIAL GUESS FOR X2:' C READ *, (X2(I),I=1,N-K) C C PRINT *,'ENTER THE INITIAL GUESS FOR X1:' C READ *, (X1(I),I=1,K) C C PRINT *,'ENTER THE ELEMENTS OF A1 ROWWISE:' C READ *,((A1(ROW,COL),COL=1,K),ROW=1,K) C C PRINT *,'ENTER THE ELEMENTS OF A2 ROWWISE:' C READ *,((A2(ROW,COL),COL=1,N-K),ROW=1,K) C C PRINT *,'ENTER THE ELEMENTS OF INVA1 ROWWISE:' C READ *,((INVA1(ROW,COL),COL=1,K),ROW=1,K) C C PRINT *,'ENTER THE NUMBER OF ITERATIONS TO SAVE:' C READ *,RECORD C*******INPUT FROM A FILE OPEN(10,FILE='input',STATUS='OLD') READ(10,*) K READ(10,*) (Y(I),I=1,K) READ(10,*) N READ(10,*) (X2(I),I=1,N-K) READ(10,*) (X1(I),I=1,K) READ(10,*) ((A1(ROW,COL),COL=1,K),ROW=1,K) READ(10,*)((A2(ROW,COL),COL=1,N-K),ROW=1,K) READ(10,*)((INVA1(ROW,COL),COL=1,K),ROW=1,K) READ(10,*) RECORD CLOSE(10) C********MAIN: IDUM=1234567 C*******initial guess for the lambdas: DO I=1,K LAMBDAS(I)=X1(I)+1.D0 END DO DO I=1,N-K LAMBDAS(K+I)=X2(I)+1.D0 END DO C*******initializing COUNT to zero: DO 30 I=1, N-K COUNT(I)=0 30 CONTINUE C*******first iteration: WRITE(*,*)'ITER 1' DO 60 I=1, N-K COUNTY=0 ELEM=0 CALL SIMULAX2(X2,A1,A2,INVA1,Y,LAMBDAS, / I,ELEM,COUNTY,K,N,IDUM) X2(I)=ELEM COUNT(I)=COUNT(I)+COUNTY 60 CONTINUE CALL SIMULAX1(A1,A2,INVA1,X2,Y,X1NEW,K,N) DO 90 I=1,K LAMBDAS(I)=DBLE(GAMDEV(X1NEW(I)+1.D0,IDUM)) 90 CONTINUE DO 100 J=K+1,N LAMBDAS(J)=DBLE(GAMDEV(X2(J-K)+1.D0,IDUM)) 100 CONTINUE C*******all the other iterations before starting recording DO 110 W=2, ITERA-RECORD DO 140 I=1, N-K COUNTY=0 ELEM=0 CALL SIMULAX2(X2,A1,A2,INVA1,Y,LAMBDAS, / I,ELEM,COUNTY,K,N,IDUM) X2(I)=ELEM COUNT(I)=COUNT(I)+COUNTY 140 CONTINUE CALL SIMULAX1(A1,A2,INVA1,X2,Y,X1NEW,K,N) DO 170 I=1, K LAMBDAS(I)=DBLE(GAMDEV(X1NEW(I)+1.D0,IDUM)) 170 CONTINUE DO 180 J=K+1,N LAMBDAS(J)=DBLE(GAMDEV(X2(J-K)+1.D0,IDUM)) 180 CONTINUE 110 CONTINUE C*********all the iterations to be recorded DO 1110 W=ITERA-RECORD+1,ITERA DO 1140 I=1, N-K COUNTY=0 ELEM=0 CALL SIMULAX2(X2,A1,A2,INVA1,Y,LAMBDAS, / I,ELEM,COUNTY,K,N,IDUM) X2(I)=ELEM COUNT(I)=COUNT(I)+COUNTY 1140 CONTINUE DO 1150 I=1, N-K GENERATEDX2(W-ITERA+RECORD,I)=X2(I) 1150 CONTINUE CALL SIMULAX1(A1,A2,INVA1,X2,Y,X1NEW,K,N) DO 1160 I=1, K GENERATEDX1(W-ITERA+RECORD,I)=X1NEW(I) 1160 CONTINUE DO 1170 I=1, K LAMBDAS(I)=DBLE(GAMDEV(X1NEW(I)+1.D0, IDUM)) 1170 CONTINUE DO 1180 J=K+1,N LAMBDAS(J)=DBLE(GAMDEV(X2(J-K)+1.D0,IDUM)) 1180 CONTINUE DO 1190 J=1,N GENERATEDLAMBDAS(W-ITERA+RECORD,J)=LAMBDAS(J) 1190 CONTINUE 1110 CONTINUE OPEN(9,FILE='output.poiss',STATUS='UNKNOWN') WRITE(9,99) 'LAMBDAS' 99 FORMAT(// 1X, A //) DO 999 I=1,RECORD WRITE(9,9999) I,(GENERATEDLAMBDAS(I,J),J=1,N) 9999 FORMAT(1X, I6, 64E10.5) 999 CONTINUE WRITE(9,88) 'X1X2' 88 FORMAT(// 1X, A //) DO 888 I=1, RECORD WRITE(9,8888) I,(GENERATEDX1(I,J),J=1,K),(GENERATEDX2(I,J), / J=1,(N-K)) 8888 FORMAT(1X, I6, 64I6) 888 CONTINUE WRITE(9,77) 'ACCEPTED X2' 77 FORMAT(// 1X, A //) DO 777 I=1, N-K WRITE(9,7777) COUNT(I) 7777 FORMAT(1X, 44I6) 777 CONTINUE STOP END C***********GAMMA RANDOM GENERATION******************************************* real*4 function gamdev(a,ix) c c Generates a random gamma variable with shape a>0 and scale=1 c ix : random seed c requires uniform random generator, ran2 c integer ix real*4 aa,ea,u0,u1,c1,c2,c3,c4,c5,w,ran2 real*8 a aa=sngl(a) if ( aa .eq. 1.0 ) then gamdev = -alog( ran2( ix ) ) return endif if ( aa .lt. 1.0 ) then ea = 2.7182818 / ( aa + 2.7182818 ) 11 u0 = ran2( ix ) u1 = ran2( ix ) if ( u0 .le. ea ) then gamdev = ( u0 / ea ) ** (1.0/aa) if ( u1 .gt. exp( -gamdev ) ) then goto 11 else return endif else gamdev = -alog( (1.0-u0)/(ea*aa) ) if( u1 .gt. gamdev**(aa-1.0) ) then go to 11 else return endif endif else c1 = aa - 1.0 c2 = ( aa- 1.0/(6.0*aa) ) / c1 c3 = 2.0 / c1 c4 = c3 + 2.0 c5 = 1.0 / sqrt(aa) 12 u1 = ran2( ix ) u2 = ran2( ix ) u1 = u2 + c5 * ( 1.0-1.86 * u1 ) if ( abs( u1- 0.5 ) .ge. 0.5 ) goto 12 w = c2 * u2 / u1 if ( ( c3*alog( u1 ) - alog( w ) + w) .ge. 1.0 ) then goto 12 else gamdev = c1 * w endif endif return end C**********POISSON GENERATION************************************************** SUBROUTINE RPOISSON(PHI,IX,ISEED) c c generates x ~ Poisson(phi) c INTEGER IX,ISEED REAL*8 PHI,Z REAL*4 RAN2 Z=0.D0 IX=-1 DO WHILE(Z.LT.PHI) Z=Z-DLOG(1.D0-RAN2(ISEED)) IX=IX+1 ENDDO RETURN END C**********RANDOM NUMBER GENERATION******************************************** REAL*4 FUNCTION RAN2(ISEED) c c ran2 from 2nd Edn of Numerical Recipes c integer iseed,im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv real*4 am,eps,rnmx parameter(im1=2147483563,im2=2147483399,am=1./im1,imm1=im1-1, & ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=12211, & ir2=3791,ntab=32,ndiv=1+imm1/ntab,eps=1.2e-7,rnmx=1.-eps) integer idum2,j,k,iv(ntab),iy c SAVE IV,IY,IDUM2 c data idum2/123456789/,iv/ntab*0/,iy/0/ if (iseed.le.0) then iseed=max(-iseed,1) idum2=iseed do j=ntab+8,1,-1 k=iseed/iq1 iseed=ia1*(iseed-k*iq1)-k*ir1 if (iseed.lt.0) iseed=iseed+im1 if (j.le.ntab) iv(j)=iseed enddo iy=iv(1) endif k=iseed/iq1 iseed=ia1*(iseed-k*iq1)-k*ir1 if (iseed.lt.0) iseed=iseed+im1 k=idum2/iq2 idum2=ia2*(idum2-k*iq2)-k*ir2 if (idum2.lt.0) idum2=idum2+im2 j=1+iy/ndiv iy=iv(j)-idum2 iv(j)=iseed if (iy.lt.1) iy=iy+imm1 ran2=min(am*iy,rnmx) RETURN END c C*******MAT-VEC MULTIPLICATION************************************************* SUBROUTINE MATVEC(MAT1,VEC2,PROD,LIM,M,N) INTEGER LIM,M,N,I,K,MAT1(LIM,LIM),VEC2(LIM), / PROD(LIM),SUM DO 600 I=1, M SUM=0 DO 610 K=1, N SUM=SUM + MAT1(I,K)*VEC2(K) 610 CONTINUE PROD(I)=SUM 600 CONTINUE END C*******SIMULA.X1************************************************************** SUBROUTINE SIMULAX1(A1,A2,INVA1,X2,Y,X1NEW,K,N) INTEGER K,N,X2(100),Y(100),I,A1(100,100),X1NEW(100), / A2(100,100),INVA1(100,100),A2X2(100),A1X1(100),LIMIT DO 4000 I=1,K A2X2(I)=0 4000 CONTINUE LIMIT=100 C WRITE(*,*)(X2(I),I=1,N-K) CALL MATVEC(A2,X2,A2X2, LIMIT,K,N-K) C WRITE(*,*)(A2X2(I),I=1,K) DO 5000 I=1,K A1X1(I)=Y(I)-A2X2(I) 5000 CONTINUE C WRITE(*,*)(A1X1(I),I=1,K) LIMIT=100 CALL MATVEC(INVA1,A1X1,X1NEW,LIMIT,K,K) C WRITE(*,*)(X1NEW(J),J=1,K) END C*******POISSON DISTR********************************************************** REAL*8 FUNCTION POISS(LAMBDA,N) INTEGER N REAL*8 LAMBDA real*4 gammln IF(N .GE.0)THEN POISS=dexp(-lambda+n*dlog(lambda)-gammln(n+1.0)) ELSE POISS=0 END IF END C*******SIMULAX2*************************************************************** SUBROUTINE SIMULAX2(X2,A1,A2,INVA1,Y,LAMBDAS, / I,ELEM,COUNTY,K,N,IDUM) INTEGER LIMSUP,IDUM,COUNTY,K,N,RESID(100), / X2(100),Y(100),I,J,ELEM,OLDELEM,A2X2(100),PROP,LIMIT, / A1(100,100),A2(100,100),INVA1(100,100),MINI REAL*8 LAMBDAS(100),NUM,DEN,FLY,PSI,VOILA REAL*4 RAN2 OLDELEM=X2(I) X2(I)=0 LIMIT=100 CALL MATVEC(A2,X2,A2X2,LIMIT,K,N-K) DO 999 J=1,100 IF(A2(J,I).EQ.1)THEN RESID(J)=Y(J)-A2X2(J) ELSE RESID(J)=10000 END IF 999 CONTINUE MINI=0 CALL MINIMO(K,RESID,MINI) LIMSUP=MINI IF(LIMSUP.NE.0)THEN COUNTY=0 PROP=0 CALL RPOISSON(LAMBDAS(K+I),PROP,IDUM) FLY=DBLE(RAN2(IDUM)) VOILA=0 CALL PSI(A1,A2,INVA1,Y,PROP,I,LAMBDAS,X2,K,N, / VOILA) NUM=VOILA CALL PSI(A1,A2,INVA1,Y,OLDELEM,I,LAMBDAS,X2,K,N, / VOILA) DEN=VOILA IF(FLY .LE. NUM/DEN)THEN ELEM=PROP COUNTY=1 ELSE ELEM=OLDELEM COUNTY=0 END IF ELSE ELEM=0 COUNTY=1 END IF END C*******PSI FOR METROPOLIS***************************************************** SUBROUTINE PSI(A1,A2,INVA1,Y,EL,I,LAMBDAS,X2,K,N, / VOILA) INTEGER K,N,X2(100),EL,I,Y(100),J,X1NEW(100),A1(100,100), / A2(100),INVA1(100),A2X2(100),A1X1(100) REAL *8 LAMBDAS(100),POISS,PROD,VOILA PROD=1 X2(I)=EL DO 414 J=1,K X1NEW(J)=0 414 CONTINUE CALL SIMULAX1(A1,A2,INVA1,X2,Y,X1NEW,K,N) DO 515 J=1,K PROD=PROD*POISS(LAMBDAS(J),X1NEW(J)) 515 CONTINUE VOILA=PROD C WRITE(*,*)'END PSI' END C*******FIND MINIMUM*********************************************************** SUBROUTINE MINIMO(N,RA,MINI) INTEGER RA(100),N,TEMP,I I=1 TEMP=RA(1) DO 386 I=2,N IF(RA(I).LT.TEMP)THEN TEMP=RA(I) ELSE TEMP=TEMP+0 END IF 386 CONTINUE MINI=TEMP END c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ REAL*4 FUNCTION GAMMLN(XX) REAL*4 XX REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, * -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ X=XX-ONE TMP=X+FPF TMP=(X+HALF)*LOG(TMP)-TMP SER=ONE DO J=1,6 X=X+ONE SER=SER+COF(J)/X ENDDO GAMMLN=TMP+LOG(STP*SER) RETURN END c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++