beta.f000644 000423 000017 00000014135 06470175400 012637 0ustar00hughesusers000000 032510 C The following routines are appropriate if beta(i,j), the conditional C log-odds between stations i and j (see equation 1 in Hughes et al (1998)) C is LINEAR in the parameters b. See the README file for further details. C C The following routines are called from various areas in nhmm code and C must be replicated for any model for beta(i,j). C C Only the routines fbetaij, lbeta, dbetadb, d2bdb2 and setup need to be C rewritten. real*8 function fbetaij(b,s,i,j) C returns beta_ij as a function of b and distance between stations include 'nhmm.inc' include 'dist.inc' integer s,i,j real*8 b(*) c fbetaij = b(1) + b(2)*d(s,i,j) return end logical function lbeta(b) C Check to see if conditional independence model holds (i.e. beta_ij = 0 C for all i and j). In this case, that means b(1)=b(2)=0. real*8 b(*) c lbeta = .FALSE. if (b(1).eq.0.0.and.b(2).eq.0.0) lbeta=.TRUE. return end real*8 function dbetadb(b,s,k,i,j) c Compute first derivative of beta(i,j) as a function of b(k) include 'nhmm.inc' include 'dist.inc' integer s,i,j,k real*8 b(*) c if (i.eq.j) then dbetadb = 0.0 elseif (k.eq.1) then dbetadb = 1.0 elseif (k.eq.2) then dbetadb = d(s,i,j) endif return end real*8 function d2bdb2(b,s,k,h,i,j) C compute second derivative of beta(i,j) as a function of b(k) integer s,i,j,k,h real*8 b(*) c d2bdb2 = 0.0 return end subroutine setup(raiopt) include 'nhmm.inc' include 'dist.inc' include 'dims.inc' include 'dhash.inc' c If needed, this gets called when nhmm is started, just after parameters are c read. In this case, the distance scale is transformed to achieve c approximate isotropy. It is more efficeint to do this once rather than c each tine fbetaij is called. integer s,i,j,raiopt real*8 phi(MAXM),e(MAXM) data phi/7*0.0/, * e/7*1.0/ c if (raiopt.eq.3) then do 10 s = 1,m do 10 i = 1,nsta-1 do 10 j = i+1,nsta d(s,i,j) = log(dis(i,j)*sqrt(cos(1.570796-dir(i,j)+phi(s))**2 + * sin(1.570796-dir(i,j)+phi(s))**2/e(s))/100.) d(s,j,i) = d(s,i,j) 10 continue endif c do 1035 j=0,MAXD1 table(j) = -1 1035 continue ierrd = 0 return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC real*8 function gbetai(b,s,i,nsta,nb,ib) C less efficient version of fbetai, but doesn't use up room in sumd integer s,i,j,nsta,ib(*),nb real*8 b(*),fbetaij c gbetai = 0.0 do 10 j = 1,i-1 10 if (ib(j).eq.1) gbetai = gbetai + fbetaij(b,s,i,j) do 20 j = i+1,nsta 20 if (ib(j).eq.1) gbetai = gbetai + fbetaij(b,s,i,j) return end subroutine fbetas(b,s,nsta,nb,ib,sumb) include 'nhmm.inc' include 'dhash.inc' C returns the vector sumb where sumb(i) = \sum_j beta_ij r_j integer s,i,nsta,ib(*),nb integer loc,lookup real*8 b(*),sumb(*) c loc = lookup(ib) do 10 i = 1,nsta sumb(i) = 0.0 do 10 k = 1,nb 10 sumb(i) = sumb(i) + b(k)*sumd(k,s,i,loc) return end subroutine fbeta(ldb,b,m,nsta,nb,ib,sumb) include 'nhmm.inc' include 'dhash.inc' C returns the matrix sumb where sumb(i,s) = \sum_j beta_sij r_j integer ldb,i,nsta,ib(*),nb,m real*8 b(ldb,*),sumb(MAXSTA,MAXM) c loc = lookup(ib) do 10 s=1,m do 10 i = 1,nsta sumb(i,s) = 0.0 do 10 k = 1,nb 10 sumb(i,s) = sumb(i,s) + b(k,s)*sumd(k,s,i,loc) return end subroutine dbeta(b,b0,s,nsta,nb,ib,db) include 'nhmm.inc' C returns the vector db where db(i) = \sum_j (beta_ij - beta0_ij) r_j C this could be done with just calls to fbetai, but hopefully this is C a little faster integer s,i,j,nsta,ib(*),nb real*8 b(*),b0(*),db(*),dbetadb,dbet(MAXB) c do 10 k = 1,nb 10 dbet(k) = b(k)-b0(k) do 100 i = 1,nsta db(i) = 0.0 do 20 j = i+1,nsta if (ib(j).eq.1) then do 30 k = 1,nb 30 db(i) = db(i) + dbet(k)*dbetadb(b,s,k,i,j) endif 20 continue 100 continue return end subroutine dbdb(ldb,b,nb,nsta,s,ib,db) include 'nhmm.inc' include 'dhash.inc' integer ldb,s,i,ib(*),nsta,nb integer loc,lookup real*8 b(*),db(ldb,*) c loc = lookup(ib) do 10 k = 1,nb do 10 i = 1,nsta db(k,i) = sumd(k,s,i,loc) 10 continue return end integer function lookup(ib) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'dist.inc' include 'dhash.inc' integer l,imis,i,j,k,ib(*) integer loc,nincr,nmax integer s real*8 b(MAXB,MAXM),dbetadb data nmax/MAXD/ c call icode(nsta,ib,l,imis) loc = mod(l,nmax) nincr = nmax-2 - mod(l,(nmax-2)) c 175 if (table(loc).eq.l) then c write(*,*) 'Found ',l goto 160 elseif (table(loc).eq.-1) then ierrd=ierrd+1 table(loc) = l c write(*,*) 'New ',l do 165 i=1,nsta do 120 k = 1,nb do 120 s = 1,m 120 sumd(k,s,i,loc)=0.0 do 135 j=1,i-1 if (ib(j).eq.1) then do 130 k=1,nb do 130 s = 1,m c Note: value of b is immaterial in linear version 130 sumd(k,s,i,loc)=sumd(k,s,i,loc)+ * dbetadb(b(1,s),s,k,i,j) endif 135 continue do 145 j=i+1,nsta if (ib(j).eq.1) then do 140 k=1,nb do 140 s = 1,m 140 sumd(k,s,i,loc)=sumd(k,s,i,loc)+ * dbetadb(b(1,s),s,k,i,j) endif 145 continue 165 continue goto 160 else c write(*,*) 'Collision' if (ierrd.ge. MAXD) stop 'Error! Increase maxd' loc = mod(loc+nincr,nmax) goto 175 endif 160 lookup = loc return end cmlib.f000644 000423 000017 00000052055 06242170130 012777 0ustar00hughesusers000000 031010 ccccccc Routines from CMLIB ccccccccccccccc double precision function dasum(n,dx,incx) c c takes the sum of the absolute values. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dtemp integer i,incx,m,mp1,n,nincx c dasum = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dtemp = dtemp + dabs(dx(i)) 10 continue dasum = dtemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dabs(dx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) 50 continue 60 dasum = dtemp return end subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dy(*),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end double precision function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end SUBROUTINE DGECO(A,LDA,N,IPVT,RCOND,Z) INTEGER LDA,N,IPVT(*) DOUBLE PRECISION A(LDA,*),Z(*) DOUBLE PRECISION RCOND C C DGECO FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION C AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, DGEFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW DGECO BY DGESL. C TO COMPUTE INVERSE(A)*C , FOLLOW DGECO BY DGESL. C TO COMPUTE DETERMINANT(A) , FOLLOW DGECO BY DGEDI. C TO COMPUTE INVERSE(A) , FOLLOW DGECO BY DGEDI. C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND DOUBLE PRECISION C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z DOUBLE PRECISION(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK DGEFA C BLAS DAXPY,DDOT,DSCAL,DASUM C FORTRAN DABS,DMAX1,DSIGN C C INTERNAL VARIABLES C DOUBLE PRECISION DDOT,EK,T,WK,WKM DOUBLE PRECISION ANORM,S,DASUM,SM,YNORM INTEGER INFO,J,K,KB,KP1,L C C C COMPUTE 1-NORM OF A C ANORM = 0.0D0 DO 10 J = 1, N ANORM = DMAX1(ANORM,DASUM(N,A(1,J),1)) 10 CONTINUE C C FACTOR C CALL DGEFA(A,LDA,N,IPVT,INFO) C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E . C TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE C TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID C OVERFLOW. C C SOLVE TRANS(U)*W = E C EK = 1.0D0 DO 20 J = 1, N Z(J) = 0.0D0 20 CONTINUE DO 100 K = 1, N IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) IF (DABS(EK-Z(K)) .LE. DABS(A(K,K))) GO TO 30 S = DABS(A(K,K))/DABS(EK-Z(K)) CALL DSCAL(N,S,Z,1) EK = S*EK 30 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = DABS(WK) SM = DABS(WKM) IF (A(K,K) .EQ. 0.0D0) GO TO 40 WK = WK/A(K,K) WKM = WKM/A(K,K) GO TO 50 40 CONTINUE WK = 1.0D0 WKM = 1.0D0 50 CONTINUE KP1 = K + 1 IF (KP1 .GT. N) GO TO 90 DO 60 J = KP1, N SM = SM + DABS(Z(J)+WKM*A(K,J)) Z(J) = Z(J) + WK*A(K,J) S = S + DABS(Z(J)) 60 CONTINUE IF (S .GE. SM) GO TO 80 T = WKM - WK WK = WKM DO 70 J = KP1, N Z(J) = Z(J) + T*A(K,J) 70 CONTINUE 80 CONTINUE 90 CONTINUE Z(K) = WK 100 CONTINUE S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) C C SOLVE TRANS(L)*Y = W C DO 120 KB = 1, N K = N + 1 - KB IF (K .LT. N) Z(K) = Z(K) + DDOT(N-K,A(K+1,K),1,Z(K+1),1) IF (DABS(Z(K)) .LE. 1.0D0) GO TO 110 S = 1.0D0/DABS(Z(K)) CALL DSCAL(N,S,Z,1) 110 CONTINUE L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T 120 CONTINUE S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) C YNORM = 1.0D0 C C SOLVE L*V = Y C DO 140 K = 1, N L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T IF (K .LT. N) CALL DAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) IF (DABS(Z(K)) .LE. 1.0D0) GO TO 130 S = 1.0D0/DABS(Z(K)) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM 130 CONTINUE 140 CONTINUE S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE U*Z = V C DO 160 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 150 S = DABS(A(K,K))/DABS(Z(K)) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM 150 CONTINUE IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K) IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 T = -Z(K) CALL DAXPY(K-1,T,A(1,K),1,Z(1),1) 160 CONTINUE C MAKE ZNORM = 1.0 S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 RETURN END SUBROUTINE DGEDI(A,LDA,N,IPVT,DET,WORK,JOB) INTEGER LDA,N,IPVT(*),JOB DOUBLE PRECISION A(LDA,*),DET(2),WORK(*) C C DGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX C USING THE FACTORS COMPUTED BY DGECO OR DGEFA. C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGECO OR DGEFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM DGECO OR DGEFA. C C WORK DOUBLE PRECISION(N) C WORK VECTOR. CONTENTS DESTROYED. C C JOB INTEGER C = 11 BOTH DETERMINANT AND INVERSE. C = 01 INVERSE ONLY. C = 10 DETERMINANT ONLY. C C ON RETURN C C A INVERSE OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE UNCHANGED. C C DET DOUBLE PRECISION(2) C DETERMINANT OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE NOT REFERENCED. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. DABS(DET(1)) .LT. 10.0 C OR DET(1) .EQ. 0.0 . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED. C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY C AND IF DGECO HAS SET RCOND .GT. 0.0 OR DGEFA HAS SET C INFO .EQ. 0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS DAXPY,DSCAL,DSWAP C FORTRAN DABS,MOD C C INTERNAL VARIABLES C DOUBLE PRECISION T DOUBLE PRECISION TEN INTEGER I,J,K,KB,KP1,L,NM1 C C C COMPUTE DETERMINANT C IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = 1.0D0 DET(2) = 0.0D0 TEN = 10.0D0 DO 50 I = 1, N IF (IPVT(I) .NE. I) DET(1) = -DET(1) DET(1) = A(I,I)*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0D0) GO TO 60 10 IF (DABS(DET(1)) .GE. 1.0D0) GO TO 20 DET(1) = TEN*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 10 20 CONTINUE 30 IF (DABS(DET(1)) .LT. TEN) GO TO 40 DET(1) = DET(1)/TEN DET(2) = DET(2) + 1.0D0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE C C COMPUTE INVERSE(U) C IF (MOD(JOB,10) .EQ. 0) GO TO 150 DO 100 K = 1, N A(K,K) = 1.0D0/A(K,K) T = -A(K,K) CALL DSCAL(K-1,T,A(1,K),1) KP1 = K + 1 IF (N .LT. KP1) GO TO 90 DO 80 J = KP1, N T = A(K,J) A(K,J) = 0.0D0 CALL DAXPY(K,T,A(1,K),1,A(1,J),1) 80 CONTINUE 90 CONTINUE 100 CONTINUE C C FORM INVERSE(U)*INVERSE(L) C NM1 = N - 1 IF (NM1 .LT. 1) GO TO 140 DO 130 KB = 1, NM1 K = N - KB KP1 = K + 1 DO 110 I = KP1, N WORK(I) = A(I,K) A(I,K) = 0.0D0 110 CONTINUE DO 120 J = KP1, N T = WORK(J) CALL DAXPY(N,T,A(1,J),1,A(1,K),1) 120 CONTINUE L = IPVT(K) IF (L .NE. K) CALL DSWAP(N,A(1,K),1,A(1,L),1) 130 CONTINUE 140 CONTINUE 150 CONTINUE RETURN END SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO) INTEGER LDA,N,IPVT(*),INFO DOUBLE PRECISION A(LDA,*) C C DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION. C C DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) . C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE C INDICATION OF SINGULARITY. C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS DAXPY,DSCAL,IDAMAX C C INTERNAL VARIABLES C DOUBLE PRECISION T INTEGER IDAMAX,J,K,KP1,L,NM1 C C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C INFO = 0 NM1 = N - 1 IF (NM1 .LT. 1) GO TO 70 DO 60 K = 1, NM1 KP1 = K + 1 C C FIND L = PIVOT INDEX C L = IDAMAX(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (A(L,K) .EQ. 0.0D0) GO TO 40 C C INTERCHANGE IF NECESSARY C IF (L .EQ. K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE C C COMPUTE MULTIPLIERS C T = -1.0D0/A(K,K) CALL DSCAL(N-K,T,A(K+1,K),1) C C ROW ELIMINATION WITH COLUMN INDEXING C DO 30 J = KP1, N T = A(L,J) IF (L .EQ. K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (A(N,N) .EQ. 0.0D0) INFO = N RETURN END SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(*),JOB DOUBLE PRECISION A(LDA,*),B(*) C C DGESL SOLVES THE DOUBLE PRECISION SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY DGECO OR DGEFA. C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGECO OR DGEFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM DGECO OR DGEFA. C C B DOUBLE PRECISION(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE TRANS(A)*X = B WHERE C TRANS(A) IS THE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0 C OR DGEFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL DGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL DGESL(A,LDA,N,IPVT,C(1,J),0) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS DAXPY,DDOT C C INTERNAL VARIABLES C DOUBLE PRECISION DDOT,T INTEGER K,KB,L,NM1 C NM1 = N - 1 IF (JOB .NE. 0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (NM1 .LT. 1) GO TO 30 DO 20 K = 1, NM1 L = IPVT(K) T = B(L) IF (L .EQ. K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1, N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL DAXPY(K-1,T,A(1,K),1,B(1),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K = 1, N T = DDOT(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/A(K,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (NM1 .LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1) L = IPVT(K) IF (L .EQ. K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c double precision da,dx(*) integer i,incx,m,mp1,n,nincx c if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end subroutine dswap (n,dx,incx,dy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end integer function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dmax integer i,incx,ix,n c idamax = 0 if( n .lt. 1 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end condsim.f000644 000423 000017 00000003237 06242170130 013353 0ustar00hughesusers000000 031036 subroutine condsim(nsim,iseed) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'data.inc' include 'params.inc' include 'dims.inc' include 'nsave.inc' c integer nsim,iseed integer ijk,t,s,isample,iten,lseq(10),year,day real*8 delta(MAXM),a(MAXM,MAXM),prob(MAXM) data iten/10/ c do 900 ijk = 1,nsim t = 0 do 140 year=1,nyears t = t+1 call steady(m,gam,natm,mu,isigma,det,adata(1,t),delta) s = isample(m,delta) call seqgen(iseed,pipars(1,1,s),beta(1,s), * s,nsta,nb,ista(1,s),iten,lseq) write(12,'(i2, 31 i2)') s,(ista(i,s),i=1,nsta) do 140 day = 2,nobs t = t+1 call makeA(m,gam,natm,mu,isigma,det,adata(1,t),a) do 160 j=1,m 160 prob(j) = a(s,j) s = isample(m,prob) call seqgen(iseed,pipars(1,1,s),beta(1,s), * s,nsta,nb,ista(1,s),iten,lseq) write(12,'(i2, 31 i2)') s,(ista(i,s),i=1,nsta) 140 continue write(12,*) write(*,*) 'Finished simulation ',ijk 900 continue return end integer function isample(n,prob) c random sample an integer from 1 to n integer iseed real*8 prob(n),u,ran1,sum c c write(*,*) 'isample',n c write(*,*) (prob(i),i=1,10) iseed = 1 1 sum = 0 u = ran1(iseed) c write(*,*) u do 10 i=1,n sum = sum + prob(i) if (u.le.sum) then isample=i goto 99 endif 10 continue write(*,*) 'Error in isample' write(*,*) n,u,sum write(*,*) 'Trying again' goto 1 99 return end data.inc000644 000423 000017 00000000206 06242170131 013150 0ustar00hughesusers000000 031506 integer nyears,nobs integer rdata(2,MAXOBS) real*8 adata(MAXATM,MAXOBS) common /data/ adata,rdata,nyears,nobs dhash.inc000644 000423 000017 00000000174 06242170131 013333 0ustar00hughesusers000000 031507 integer ierrd,table(0:MAXD1) real*8 sumd(MAXB,MAXM,MAXSTA,0:MAXD1) common /dhash/ sumd,table,ierrd dims.inc000644 000423 000017 00000000100 06242170131 013167 0ustar00hughesusers000000 030753 integer m,nsta,natm,nb common /dims/ m,nsta,natm,nb dist.inc000644 000423 000017 00000000164 06242170131 013211 0ustar00hughesusers000000 030754 real*8 dis(MAXSTA,MAXSTA),dir(MAXSTA,MAXSTA) real*8 d(MAXM,MAXSTA,MAXSTA) common /dist/ d,dir,dis dldp.f000644 000423 000017 00000013331 06470172240 012646 0ustar00hughesusers000000 033260 subroutine dlset1(t,gam,mu,isigma,det,sum,eQ,denom,delta,ia) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer s0,s1,t integer ipvt(MAXM) real*8 sum(MAXM),eQ(MAXM,MAXM),denom(MAXM) real*8 atmf real*8 delta(MAXM),a(MAXM,MAXM),ia(MAXM,MAXM) real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM) real*8 isigma(MAXATM,MAXATM) real*8 det,work(MAXM) c call steady(m,gam,natm,mu,isigma,det,adata(1,t),delta) call makeA(m,gam,natm,mu,isigma,det,adata(1,t),a) call augment(MAXM,m,a) do 202 s0=1,m sum(s0) = 0.0 do 200 s1=1,m eQ(s0,s1) = atmf(natm,adata(1,t),s0,s1,mu,isigma,det) sum(s0) = sum(s0) + gam(s0,s1)*eQ(s0,s1) 200 continue denom(s0) = sum(s0)*sum(s0) 202 continue call invsig(MAXM,m,a,det,ia,work,ipvt) return end subroutine dldgm1(t,ifix,gam,mu,isigma,det,sum,eQ,denom,delta, * ia,Sgm,n) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer s0,t,ifix(5),n real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM) real*8 isigma(MAXATM,MAXATM),det real*8 dQ,dQdmu real*8 sum(MAXM),eQ(MAXM,MAXM),Sgm(MAXM,*) real*8 denom(MAXM),ap1(MAXM),dp1(MAXM) real*8 delta(MAXM),ia(MAXM,MAXM) c n=0 if (ifix(1).ne.0) then do 320 i = 1,m do 320 j = 1,m n = n+1 call dAdgam(i,j,MAXM,m,sum,denom,gam,eQ,ap1) call ddelta(i,MAXM,m,ia,ap1,delta,dp1) do 320 s0 = 1,m Sgm(s0,n) = dp1(s0)/delta(s0) 320 continue endif if (ifix(4).ne.0) then do 325 i=1,m do 325 j=1,m call dAdmu(i,j,MAXM,m,sum,denom,gam,eQ,ap1) do 325 k =1,natm n = n+1 dQ = dQdmu(t,i,j,k,isigma,mu) call ddelta(i,MAXM,m,ia,ap1,delta,dp1) do 325 s0 = 1,m Sgm(s0,n) = dQ*dp1(s0)/delta(s0) 325 continue endif return end subroutine dlset2(t,gam,mu,isigma,det,sum,eQ,denom,a) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer s0,s1,t real*8 sum(MAXM),eQ(MAXM,MAXM),denom(MAXM) real*8 atmf real*8 a(MAXM,MAXM) real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM) real*8 isigma(MAXATM,MAXATM) real*8 det c call makeA(m,gam,natm,mu,isigma,det,adata(1,t),a) do 312 s1=1,m sum(s1) = 0.0 do 310 s0=1,m eQ(s1,s0) = atmf(natm,adata(1,t),s1,s0,mu,isigma,det) sum(s1) = sum(s1) + gam(s1,s0)*eQ(s1,s0) 310 continue denom(s1) = sum(s1)*sum(s1) 312 continue return end subroutine dldgm2(t,ifix,gam,mu,isigma,det,sum,eQ,denom,a,Sgm,n) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer s0,t,ifix(5),n real*8 sum(MAXM),eQ(MAXM,MAXM),Sgm(MAXM,*) real*8 denom(MAXM),ap1(MAXM),a(MAXM,MAXM) real*8 dQ,dQdmu real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM) real*8 isigma(MAXATM,MAXATM),det c n = 0 if (ifix(1).ne.0) then do 315 i=1,m do 315 j=1,m n=n+1 call dAdgam(i,j,MAXM,m,sum,denom,gam,eQ,ap1) do 315 s0 = 1,m Sgm(s0,n) = ap1(s0)/a(i,s0) 315 continue endif if (ifix(4).ne.0) then do 325 i=1,m do 325 j=1,m call dAdmu(i,j,MAXM,m,sum,denom,gam,eQ,ap1) do 325 k =1,natm n=n+1 dQ = dQdmu(t,i,j,k,isigma,mu) do 325 s0 = 1,m Sgm(s0,n) = ap1(s0)*dQ/a(i,s0) 325 continue endif return end subroutine dldpb(t,expect,Spb,npb) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'fb.inc' include 'params.inc' include 'dims.inc' include 'data.inc' include 'options.inc' c integer s,t,npb integer i,ib1(MAXSTA),ib2(MAXSTA) real*8 Spb(*),expect(*) real*8 sum real*8 newbeta(MAXB),dbetadb real*8 zero data zero/0.0d0/ npb = m*(nsta+nb) call vbit(ib1,rdata(1,t),nsta) call vbit(ib2,rdata(2,t),nsta) if (raiopt.eq.1) then c independence model i1 = 0 do 1000 s = 1,m do 140 i = 1,nsta i1 = i1 + 1 if (ib2(i).eq.0) then if (ib1(i).eq.0) then Spb(i1) = -1.0/(1.0-pipars(2,i,s)) else Spb(i1) = 1.0/pipars(2,i,s) endif else Spb(i1) = 0.0 endif 140 continue 1000 continue elseif (raiopt.eq.3) then c autologistic model i1 = 0 do 2000 s = 1,m do 145 k = 1,nb 145 newbeta(k) = beta(k,s) do 150 i = 1,nsta i1 = i1+1 Spb(i1) = ib1(i) - expect(i1) 150 continue do 170 k = 1,nb i1 = i1+1 sum = 0.0 do 160 i = 1,nsta-1 do 160 j = i+1,nsta if (ib1(i).eq.1.and.ib1(j).eq.1) * sum = sum + dbetadb(newbeta,s,k,i,j) 160 continue Spb(i1) = sum - expect(i1) 170 continue 2000 continue endif return end emalg.f000644 000423 000017 00000026272 06470174733 013024 0ustar00hughesusers000000 033030 subroutine emalg(lsets,lsetr) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' include 'norm.inc' include 'options.inc' include 'npsetup.inc' include 'weights.inc' include 'fb.inc' logical lsets,lsetr integer ns,nr real*8 xs(MAXNS),xr(MAXNR) real*8 diff,tol,l2norm,pctdif real*8 rlik,loglik,oldlik character*80 fintmd c data tol/.0000000001/ c ifails = 4 ifailr = 4 call gptox(ifix,m,nsta,natm,nb,gam,mu,pipars,beta,xs,ns,xr,nr) if (lsets) call npsets(nclins,ncnlns,leniws,lenws,as,bls,bus, * lsets) if (lsetr) call npsetr(nclinr,ncnlnr,leniwr,lenwr,ar,blr,bur, * lsetr) call e04uef('Cold Start') call e04uef('Derivative Level = 3') c call e04uef('Derivative Level = 0') c call e04uef('Verify Level = 3') call e04uef('Verify Level = -1') call e04uef('Major Print Level = 0') call e04uef('Hessian = Yes') call getnam(fintmd) c calculate forward-backward probs (E step) call forbac() oldlik = loglik() call weights(oldlik) write(*,*) 'Current Likelihood = ',oldlik call prtpar(15) write(15,*) 'Current Likelihood = ',oldlik close(15) 1 continue open(15,file=fintmd,status="OLD",access="APPEND") c maximize psi_s (M step) if (ifix(1).eq.1.or.ifix(4).eq.1) then c if (ifails.eq.0) then c call e04uef('Warm Start') c else c call e04uef('Cold Start') c endif call stepps(xs,ns,nclins,ncnlns,leniws,lenws,as,bls,bus,ifails) call e04uef('Warm Start') if (ifails.eq.-1) call gtox(ifix,m,natm,gam,mu,xs,ns) endif c maximize psi_r (M step) if (ifix(2).eq.1.or.ifix(3).eq.1) then call steppr(xr,nr,blr,bur,ifailr) c if (ifailr.eq.-1) call ptox(ifix,m,nsta,nb,pipars,beta,xr,nr) endif c diff = l2norm(xs,xr,gam,mu,pipars,beta,ifix) call xtogp(ifix,m,nsta,natm,nb,gam,mu,pipars,beta,xs,xr) call rehash() if (raiopt.eq.3) call norman(m,nsta,nb,pipars,beta,lnorm) call prtpar(15) call forbac() rlik = loglik() call weights(rlik) pctdif = dabs((rlik-oldlik)/oldlik) c write(*,*) 'Parameter change = ',diff write(*,*) 'Likelihood % change = ', pctdif write(*,*) 'Current Likelihood = ',rlik write(15,*) 'Likelihood % change = ', pctdif write(15,*) 'Current Likelihood = ',rlik close(15) oldlik=rlik if (pctdif.gt.tol) goto 1 write(*,*) 'Parameter change on final iteration = ',diff return end subroutine stepps(x,n,nclin,ncnln,leniw,lenw,a,bl,bu,ifail) include 'nhmm.inc' include 'params.inc' include 'dims.inc' integer n,nclin,ncnln integer nrowa,nrowj,nrowr integer iter,leniw,lenw integer istates(MAXPAR),iw(MAXIW),ifail,iuser(1) integer mode,nstate real*8 x(*),objf,objf1,grad(MAXNS) real*8 a(MAXLIN,MAXNS),bl(MAXPAR),bu(MAXPAR) real*8 c(MAXNLN),cjac(MAXNLN,MAXNS),clambs(MAXPAR) real*8 rs(MAXNS,MAXNS) real*8 w(MAXW),user(1) common /nps/ clambs,rs,istates external confun,objfuns data nrowa /MAXLIN/ data nrowj /MAXNLN/ data nrowr /MAXNS/ c mode = 2 nstate = 0 ifail = -1 call objfuns(mode,n,x,objf1,grad,nstate,iuser,user) c write(*,*) objf1 write(*,*) 'Gradients of hidden model params at start of M-step' write(*,*) (grad(i),i=1,n) call e04ucf(n,nclin,ncnln,nrowa,nrowj,nrowr, * a,bl,bu,confun,objfuns,iter, * istates,c,cjac,clambs,objf,grad,rs,x, * iw,leniw,w,lenw,iuser,user,ifail) write(*,*) 'M-step error code for hidden model params ',ifail if (ifail.eq.2.or.ifail.eq.3.or.ifail.eq.7.or.ifail.eq.9) stop if (ifail.ne.0.and.objf1.lt.objf) then c write(*,*) objf ifail = -1 endif return end subroutine weights(rlik) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'params.inc' include 'data.inc' include 'fb.inc' include 'weights.inc' integer s,s2,t,year,day real*8 rlik,scale,scale2 real*8 pi(MAXM),a(MAXM,MAXM) c t = 0 do 100 year = 1,nyears t = t+1 scale = exp(forw(m+1,t)+back(m+1,t)-rlik) c call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pi) c call makeA(m,gam,natm,mu,isigma,det,adata(1,t),a) do 10 s = 1,m 10 c(s,t) = forw(s,t)*back(s,t)*scale do 100 day = 2,nobs t = t+1 scale = exp(forw(m+1,t)+back(m+1,t)-rlik) scale2 = exp(forw(m+1,t-1)+back(m+1,t)-rlik) call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pi) call makeA(m,gam,natm,mu,isigma,det,adata(1,t),a) do 100 s = 1,m c(s,t) = forw(s,t)*back(s,t)*scale do 100 s2 = 1,m c2(s,s2,t) = forw(s,t-1)*a(s,s2)*pi(s2)*back(s2,t)*scale2 100 continue return end subroutine steppr(x,n,bl,bu,ifail) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'options.inc' include 'dims.inc' include 'params.inc' c integer n,ifail integer nrowr,mode,s,info real*8 x(*),objf,grad(MAXIR) real*8 d(MAXIR),alpha,newx(MAXIR),objf1 real*8 bl(*),bu(*),tol real*8 rr(MAXIR,MAXIR) external objfr1,objfr3 data nrowr /MAXIR/ data tol/1D-14/ c write(*,*) 'Gradients of output model params at start of M-step' goto (100,200,300) raiopt C Independence model 100 continue do 110 s = 1,m ind=nsta*(s-1)+1 call objfr0(nsta,s,x(ind),grad,bl(ind),bu(ind)) write(*,*) 'state ',s write(*,*) (grad(i),i=1,nsta) 110 continue ifail=0 goto 999 C Saturated model not implemented 200 continue stop 'Saturated model not implemented' C Autologistic model 300 continue ifail = 0 nr = ifix(2)*nsta + ifix(3)*nb do 3050 s=1,m mode = 0 ind=nr*(s-1)+1 if (iest.le.2) then call objfr1(mode,nr,s,nrowr,x(ind),objf,grad,rr) elseif (iest.eq.3) then call objfr3(mode,nr,s,nrowr,x(ind),objf,grad,rr) endif write(*,*) 'state ',s c write(*,*) objf c write(*,*) (x(ind+i-1),i=1,nr) write(*,*) (grad(i),i=1,nr) c do 1 i = 1,nr c write(*,*) (rr(i,j),j=1,nr) c1 continue call delta(nr,nrowr,grad,rr,d,info) if (info.gt.0) then write(*,*) 'Unable to make progress' goto 3050 endif c write(*,*) (d(i),i=1,nr) alpha = 2.0 3005 alpha = alpha/2.0 if (alpha.lt.tol) then ifail=-1 goto 3050 endif do 3010 i = 1,nr 3010 newx(i) = x(ind+i-1) + d(i)*alpha mode = 2 if (iest.le.2) then call objfr1(mode,nr,s,nrowr,newx,objf1,grad,rr) elseif (iest.eq.3) then call objfr3(mode,nr,s,nrowr,newx,objf1,grad,rr) endif write(*,*) objf1 if (objf1.gt.objf) goto 3005 call copyx(newx,x(ind),nr) c write(*,*) (x(ind+i-1),i=1,nr) 3050 continue 999 return end subroutine delta(nr,nrowr,grad,rr,d,info) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' integer nr,nrowr integer ipvt(MAXIR),info,job real*8 grad(*),d(*),rr(nrowr,*),det(2),work(MAXIR) real*8 a(MAXIR,MAXIR) data job/01/,zero/0.0D0/ c n = nr do 1 i = 1,n do 1 j = 1,n 1 a(i,j) = rr(i,j) call DGEFA(a,nrowr,n,ipvt,info) c write(*,*) info if (info.gt.0) then n = nr - nb do 2 i = 1,n do 2 j = 1,n 2 a(i,j) = rr(i,j) call DGEFA(a,nrowr,n,ipvt,info) endif if (info.gt.0) write(*,*) 'Error inverting hessian' call DGEDI(a,nrowr,n,ipvt,det,work,job) c do 5 i = 1,nr 5 d(i) = zero do 10 i = 1,n do 10 j = 1,n 10 d(i) = d(i) - grad(j)*a(i,j) c return end subroutine npsets(nclin,ncnln,leniw,lenw,a,bl,bu,lsets) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' logical lsets integer n,nclin,ncnln,leniw,lenw,ncladd integer nclin1,nclin2,i,j,k,i1,i2 real*8 a(MAXLIN,MAXNS), bl(*),bu(*) real*8 zero,one,BIGBND,eps data zero,one,BIGBND,eps/0.0,1.0,1.0D15,1.0D-7/ c n = 0 nclin = 0 nclin1 = 0 nclin2 = 0 ncnln = 0 if (ifix(1).ne.0) then nclin1= m nclin = m c Linear constraints - rows of gam must sum to 1 do 140 i=1,nclin do 140 j = (i-1)*m+1,(i-1)*m+m 140 a(i,j) = one c Elements of gam must be between 0 and 1 do 150 i=n+1,n+m*m bl(i) = zero + eps 150 bu(i) = one - eps n = n + m*m endif if (ifix(4).ne.0) then c mu_i sum to 0 nclin2 = m*natm do 145 i=1,m do 145 j=1,natm i1 = nclin + (i-1)*natm + j do 145 k=1,m i2 = n + (i-1)*m*natm + (k-1)*natm + j 145 a(i1,i2) = one nclin = nclin + nclin2 do 151 i=n+1,n+natm*m*m bl(i) = -BIGBND 151 bu(i) = BIGBND n = n + natm*m*m endif c Linear constraints on gamma do 155 i=n+1, n+nclin1 bl(i) = one 155 bu(i) = one c Linear constraints on mu do 156 i=n+nclin1+1, n+nclin1+nclin2 bl(i) = zero 156 bu(i) = zero c Additional user defined linear constraints read(11,*) ncladd do 158 i=1,ncladd 158 read(11,*) (a(nclin+i,j),j=1,n),bl(n+nclin+i),bu(n+nclin+i) nclin = nclin + ncladd leniw = 3*n + nclin +2*ncnln lenw = 2*n**2 + n*nclin + 2*n*ncnln + 20*n + 11*nclin + 21*ncnln lsets = .FALSE. return end subroutine npsetr(nclin,ncnln,leniw,lenw,a,bl,bu,lsetr) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' logical lsetr integer n,nclin,ncnln,leniw,lenw integer i real*8 a(1,MAXNR), bl(*),bu(*) real*8 BIGBND data PIBND,BIGBND/17.0,1.0D15/ c n = 0 nclin = 0 ncnln = 0 do 180 i = 1,m c if (ifix(2).ne.0) then do 171 j=1,nsta n = n+1 bl(n) = -PIBND 171 bu(n) = PIBND endif c if (ifix(3).ne.0) then do 175 j=1,nb n = n+1 bl(n) = -BIGBND bu(n) = BIGBND 175 continue endif c 180 continue leniw = 3*n + nclin +2*ncnln lenw = 2*n**2 + n*nclin + 2*n*ncnln + 20*n + 11*nclin + 21*ncnln lsetr = .FALSE. return end c====================================================================== c CONFUN is required by e04ucf c---------------------------------------------------------------------- subroutine CONFUN(mode,ncnln,n,nrowj,needc,x,c,cjac,nstate) Implicit double precision (A-H,O-Z) integer mode,ncnln,n,nrowj,nstate,needc(*) real*8 x(n),c(*),cjac(nrowj,*) c return end subroutine getnam(fintmd) character*80 fintmd character*9 ints data ints/"123456789"/ c fintmd = "nhmm000.out" do 15 i = 1,9 fintmd(7:7) = ints(i:i) open(15,file=fintmd,status="NEW",err=15) goto 16 15 continue 16 write(*,*) "Intermediate results will be stored in ",fintmd return end fb.inc000644 000423 000017 00000000117 06242170132 012630 0ustar00hughesusers000000 031416 real*8 forw(MAXM1,MAXOBS),back(MAXM1,MAXOBS) common /fb/ forw,back forbac.f000644 000423 000017 00000007260 06302403124 013153 0ustar00hughesusers000000 031550 subroutine forbac() Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'params.inc' include 'data.inc' include 'fb.inc' include 'options.inc' integer i,j,t,year,day real*8 zero,one,scale,tscale real*8 pisum(MAXM),delta(MAXM),a(MAXM,MAXM) data zero,one/0.0,1.0/ c c Compute the forward-backward probabilities for the NHMM. forw and back c are the forward and backward probabilities. For day t, forw(1:m,t) and c back(1:m,t) are the scaled f/b probabilities corresponding to that day. c forw(m+1,t) and back(m+1,t) are the (log)scale factors. Thus, the true c forward probabilities for day t are forw(1:m,t)*exp(forw(m+1,t)), and c similarly for the backward probabilities. c t=0 c Take advantage of the fact that fortran arrays are stored in column c major order by (effectively) passing a pointer to rdata and adata in c sumpi, steady and makeA do 102 year=1,nyears c Forward probabilities t=t+1 call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pisum) call steady(m,gam,natm,mu,isigma,det,adata(1,t),delta) if (year.eq.1) then forw(m+1,t) = zero else forw(m+1,t) = forw(m+1,t-1) + log(scale) endif scale = zero do 101 i = 1,m forw(i,t) = delta(i)*pisum(i) scale = scale + forw(i,t) 101 continue do 103 day=2,nobs t=t+1 call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pisum) call makeA(m,gam,natm,mu,isigma,det,adata(1,t),a) tscale = zero do 104 i=1,m forw(i,t) = zero do 105 j=1,m forw(i,t) = forw(i,t) + forw(j,t-1)*a(j,i) 105 continue forw(i,t) = forw(i,t)*pisum(i)/scale tscale = tscale + forw(i,t) 104 continue forw(m+1,t) = forw(m+1,t-1) + log(scale) scale = tscale if (scale.eq.zero) write(*,*) * 'Warning. Scale = 0 at year=',year,' day = ',day 103 continue 102 continue c Backward probabilities (only compute if needed) t = t+1 do 212 year = nyears, 1, -1 t=t-1 if (year.eq.nyears) then do 201 j=1,m back(j,t) = one 201 continue back(m+1,t) = 0 else call sumpi(rdata(1,t+1),m,nsta,nb,pipars,beta,pisum) call steady(m,gam,natm,mu,isigma,det,adata(1,t+1),delta) tscale = zero do 202 i = 1,m back(i,t) = zero do 222 j=1,m 222 back(i,t) = back(i,t)+delta(j)*pisum(j)*back(j,t+1) back(i,t) = back(i,t)/scale tscale = tscale + back(i,t) 202 continue back(m+1,t) = back(m+1,t+1) + log(scale) scale = tscale endif scale = m do 203 day = nobs-1, 1, -1 t=t-1 call sumpi(rdata(1,t+1),m,nsta,nb,pipars,beta,pisum) call makeA(m,gam,natm,mu,isigma,det,adata(1,t+1),a) tscale = zero do 205 i=1,m back(i,t) = zero do 204 j=1,m 204 back(i,t) = back(i,t)+a(i,j)*pisum(j)*back(j,t+1) back(i,t) = back(i,t)/scale tscale = tscale + back(i,t) 205 continue back(m+1,t) = back(m+1,t+1) + log(scale) scale = tscale if (scale.eq.zero) write(*,*) * 'Warning. Scale = 0 at year=',year,' day = ',day 203 continue 212 continue return end c====================================================================== c This computes the likelihood for the data c---------------------------------------------------------------------- real*8 function loglik() Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'fb.inc' include 'data.inc' include 'dims.inc' c integer day data zero/0.0d0/ c loglik = zero day = nyears*nobs do 40 i=1,m 40 loglik = loglik + forw(i,day) loglik=log(loglik)+forw(m+1,day) return end hash.inc000644 000423 000017 00000000146 06242170132 013167 0ustar00hughesusers000000 031417 integer hvec(2,0:MAXPI1),ierrh real*8 pi(0:MAXPI1,MAXM) common /HASH/ hvec,pi,ierrh makea.f000644 000423 000017 00000006647 06302400314 013005 0ustar00hughesusers000000 030760 c====================================================================== c This subroutine makes the matrix A(x) c---------------------------------------------------------------------- subroutine makeA(m,gam,natm,mu,isigma,det,x,a) include 'nhmm.inc' Implicit double precision (A-H,O-Z) integer m,natm,i,j real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM) real*8 isigma(MAXATM,MAXATM),det,x(*) real*8 a(MAXM,MAXM),sumrow real*8 atmf c write(*,*) 'makeA' do 40 i=1,m sumrow = 0.0 do 25 j=1,m a(i,j) = atmf(natm,x,i,j,mu,isigma,det)*gam(i,j) sumrow = sumrow + a(i,j) 25 continue do 30 j=1,m a(i,j) = a(i,j)/sumrow 30 continue c write(*,*) (a(i,j),j=1,m) 40 continue return end c====================================================================== c This routine computes P(atmospheric data | S(t),S(t-1)). c---------------------------------------------------------------------- real*8 function atmf(natm,x,i,j,mu,isigma,det) include 'nhmm.inc' Implicit double precision (A-H,O-Z) real*8 x(*), mu(MAXM,MAXM,MAXATM), isigma(MAXATM,MAXATM) real*8 det,qf c real*8 pi integer i,j,k,l,natm c data pi/3.1415926535/ c c write(*,*) 'atmf' c Multivariate normal qf=0.0 do 10 k = 1,natm do 10 l = 1,natm 10 qf = qf + isigma(k,l)*(x(k)-mu(i,j,k))*(x(l)-mu(i,j,l)) c atmf = (1/dsqrt(2*pi))**natm * (1/dsqrt(det)) * dexp(-.5*qf) atmf = dexp(-.5*qf) c write(*,*) i,j,atmf return end c====================================================================== c This routine inverts sigma and gets its determinant. LINPACKD routines c DGEFA and DGEDI are used. c---------------------------------------------------------------------- subroutine invsig(lda,n,sigma,det,isigma,work,ipvt) Implicit double precision (A-H,O-Z) integer lda,n,ipvt(*) real*8 det,sigma(lda,*),isigma(lda,*) real*8 work(*),d(2) integer job,info,i,j data job/11/ do 20 i=1,n do 20 j=1,n 20 isigma(i,j) = sigma(i,j) call DGEFA(isigma,lda,n,ipvt,info) c write(*,*) 'After DGEFA' call DGEDI(isigma,lda,n,ipvt,d,work,job) c write(*,*) 'After DGEDI' if (d(1).lt.0.5) then det = 0.0 else det = d(1)*10.0**d(2) endif return end c====================================================================== c This computes the steady state probabilites. c To solve the linear equations the cmlib LINPACKD subroutines c DGECO and DGESL are used. c---------------------------------------------------------------------- subroutine steady(m,gam,natm,mu,isigma,det,x,delta) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer m,natm integer i,j integer lda,ipvt(MAXM),job real*8 gam(MAXM,MAXM),delta(*) real*8 x(*) real*8 mu(MAXM,MAXM,MAXATM),isigma(MAXATM,MAXATM),det real*8 zero,one real*8 a(MAXM,MAXM),z(MAXM),rcond,t(MAXM,MAXM) c external DGESL,DGECO data zero,one/0.0d0,1.0d0/ data lda,job/MAXM,0/ c write(*,*) 'steady' call makeA(m,gam,natm,mu,isigma,det,x,t) c write(*,6001) ((t(i,j),j=1,m),i=1,m) c6001 format(" A = ", 2f8.5) do 20 i=1,m do 10 j=1,m a(i,j)=t(j,i) 10 continue 20 continue do 30 i=1,m delta(i)=zero a(i,i)=a(i,i)-one a(m,i)=one 30 continue delta(m)=one call DGECO(a,lda,m,ipvt,rcond,z) call DGESL(a,lda,m,ipvt,delta,job) c write(*,*) (delta(i),i=1,m) return end nhmm.f000644 000423 000017 00000074734 06435141157 012703 0ustar00hughesusers000000 031540 c====================================================================== c nhmm.m4 c====================================================================== c Version 13 m-state nsta-station model with dependence between stations c c General program for estimating paramters of m-state, nsta-station c nonnstationary hidden markov model. The nonstationary refers to the c transition matrix which is assumed to depend on atmospheric data. The c EM algorithm with modification as described in Rai and Matthews c (Biometrics 49,587-591, 1993) is used. The user may specify a c model for between station dependence. c Missing values in the rainfall data are allowed. c Also, the user can enter arbitrary parameter constraints. c c Estimation can be by exact maximum likelihood or maximum Pseudo- c Likelihood c c NAME: nhmm.for c c====================================================================== c INPUT: c c Input consists of 2 or 3 files, depending on the options chosen. These c are the c c PARAMETER CONTROL FILE - sets user-defined options, selects model, c specifies initial parameter values c c RAW DATA FILE - contains the (daily) rain and atmospheric data. c c DISTANCE FILE - contains the matrix of distances between rain c stations. Used only if the spatial model c is selected for rainfall. c c Each of these files is described further below. c c---------------------------------------------------------------------- c PARAMETER CONTROL FILE. c Name: supplied by the user interactively c Function: sets user-defined options, selects model, specifies c initial parameter values c Format: all input is free format. c c Record Description c 1 number of climate states (m), number of stations (nsta), number c of years of data (nyear), number of time periods per year (ndays), c number of atmospheric variables (natm), atmospheric model c option (atmopt), rain model option (raiopt) c Note: 1 <= dbeg < dend <= ndays c Note: Years are treated as independent sequences. Thus, one c could have 10 years of Jan. data. However, if the data c do not form independent sequences (e.g. 5 complete c years of data), then a "year" may be arbitrarily long. c For example, one "year" of 1825 days. The only constraint c is that nyears*ndays < MAXOBS c Note: atmopt = 1 is Bayes model c = 2 is Logistic model (not implemented) c Note: raiopt = 1 is Independence model c = 3 is autologistic model (a spatial model) c c 2 the data file name (datfil) c c 3 the format of the data in datfil (dfmt); this should read 1 day c per line and the first `ndays' lines should be the first year, c etc. Each line gives the rain indicators and the atmospheric data. c See format of RAW DATA FILE below. c c------------------------------------------------------------ c ENTER THIS ONLY IF raiopt = 3 c 4 The number of parameters in the spatial model for each state (nb) c and the name of file containing distances/directions between c stations (disfil) c------------------------------------------------------------ c c next An indicator variable which determines whether the next set of c parameters (gam) should be treated as fixed c or estimated. Values are 0 = fixed, 1 = estimate. c c next m (gam) c each line should contain the m initial transition probabilities c from state i (i=1,m) to states 1 - m c c next An indicator variable which determines whether the next set of c parameters (pipars) should be treated as fixed or estimated. c Values are 0 = fixed, 1 = estimate. c c next nsta (p) c each line should contain m values, one for each state. There c will be nsta such lines, one for each station. The value given c is the initial value for p_is = 1/(1+exp(-alpha_is)) where i is c the station number and s is the state. alpha is the parameter in c equation 4.1 in Hughes (1993). Under the independence model, p_is c is the probability of rain in state s (s=1,m) for station i c (i=1,nsta). c C ----------------------------------------------------------------------- C ENTER THESE ONLY IF raiopt = 3 C next An indicator variable which determines whether the next set of C parameters (beta) should be treated as fixed or estimated. C Values are 0 = fixed, 1 = estimate. C C next nb (beta) C these lines should contain the m*nb inital values for beta C ----------------------------------------------------------------------- C ENTER THESE ONLY IF natm > 0 c next An indicator variable which determines whether the next set of c parameters (mu) should be treated as fixed c or estimated. Values are 0 = fixed, 1 = estimate. c c next m*m (mu or b) c if model is Bayes (mu) c each line should contain the natm initial means for the c atmospheric variables when S(t-1) = i, S(t) = j (i,j = 1,m) c with j varying more quickly. c c next natm (sigma) c the (symmetric) variance-covariance matrix for the atmospheric c variables C ----------------------------------------------------------------------- c c next beginning of time period to analyse (dbeg), c end of time period to analyse (dend) c c next estimation option (iest) c Note: iest = 1 is exact maximum likelihood c = 2 is mcml c = 2 is mcml + mpsl c if iest >= 2, also enter, on same line, a negative random number c seed and the relative variability desired in the c monte carlo estimate of the norm c C ----------------------------------------------------------------------- C ENTER THESE ONLY IF iest >= 2 c c next random number seed (iseed), relative error for estimating the c norm (eps=se(norm)/norm). c c next the m log(normalizing constants) corresponding to the initial c values of alpha and beta. If beta = 0, these are ignored and c the program computes the norms anyway, so just enter 0. c c next m bcd random numbers from each of the autologistic distributions c corresponding to the initial values of alpha and beta. If c beta = 0 these are ignored, so just enter 0. C ----------------------------------------------------------------------- c c next ncladd - number of additional linear constraints on the c gamma and mu parameters (besides the usual) c c ----------------------------------------------------------------------- c ENTER THESE ONLY IF ncladd > 0 c c next a(i,*), b - each row is a set of n coefficients that form c a constraint on the parameters. The last element of the c row is the constant that the constraint sums to. c That is, ax = b. c c next ncladd - number of additional linear constraints on the c pipars and beta parameters (NOT IMPLEMENTED AT PRESENT) c c next a(i,*), b - each row is a set of n coefficients that form c a constraint on the parameters. The last element of the c row is the constant that the constraint sums to. c That is, ax = b. (NOT IMPLEMENTED AT PRESENT) c------------------------------------------------------------------------ c---------------------------------------------------------------------- c RAW DATA FILE c c Name: Given in PARAMETER CONTROL FILE c Function: supplies the data c Format: This file should contain ndays*nyears lines. Each line c represents one day and contains the nsta rainfall c indicators (0 = dry, 1 = wet, 9 = missing) and the c natm atmospheric vars for that day. Each line is read c with the format specified by dfmt (see PARAMETER CONTROL FILE) c------------------------------------------------------------------------ c---------------------------------------------------------------------- c DISTANCE FILE c c Name: Given in PARAMETER CONTROL FILE c Function: supplies the distance and direction matrix for the spatial c model for c rainfall. IGNORED UNLESS SPATIAL MODEL (RAIOPT=3) IS USED. c Format: There are 2*nsta lines and each line contains nsta numbers. c The first nsta lines give distances; the second nsta give c directions (in degrees). c The j'th number in line i is the distance from station j to i. c====================================================================== c c OUTPUT: c gam (stationary transition probs for the hidden mc) c pipars(parameters used to compute pi matrix) c mu (mean of atmospheric variable k for climate transition c j -> i) c beta (spatial parameters) c sigma (var-covar matrix of atmospheric vars) c====================================================================== c c REMARKS: c c The libraries CMLIB and NPSOL are accessed by this program. See c the Makefile for more information. c c The file nhmm.inc contains all key dimensions of the program. To c expand program to handle more states or atmospheric vars, change c the dimensions there and recompile. Increasing the number of c stations over 32 would require reprogramming. c c The name of the input data file can be up to 80 characters long, at c present. To change, just change the declaration c character*80 datfil c and read formats to the desired value and recompile. c c If all the parameters are estimated the estimates c will not be unique. Typically, sigma is known apriori and fixed. c c To estimate a standard hmm set atmopt = 1, natm = 0, c and do not enter mu and sigma c c There is a problem if many rain station obervations are missing c on a given day. Since the program sums over all possible missing c rain states there are 2**missing probabilities to sum up. This c will create a disaster if missing is large. c c HISTORY: c date who? what? c 10/93 Jim Hughes This program is a revision of the est14 c program used for my dissertation. The c major change is that estimation is now c based on the EM algorithm instead of c direct maximization of the likelihood. c c================================================================= c c A few important variables: c c m = number of states (max = MAXM) c nsta = number of stations (max = MAXSTA <= 31) c nyears = number of years of data c ndays = number of time units per year in data file c dbeg = day of year that marks beginning of season c dend = day of year that marks end of season c nobs = dend-dbeg+1 = number of days to analyze (max nyears*ndays = MAXOBS) c natm = number of atmospheric variables (max = MAXATM) c atmopt = 1 use Bayes model; 2 use Logistic model c gam(i,j) = prob(S(t)=j|S(t-1)=i) c pipars(1,i,j) = exp(pipars(2,i,j))/(1+exp(pipars(2,i,j))) c pipars(2,i,j) = prob(rain at station i|S(t)=j) (for independence model) c beta(i,j) = dependence parameter i for state j c mu(i,j,k) = E(atm var(k)|S(t)=j,S(t-1)=i) c sigma(i,j) = Cov(atm var(i), atm var(j)) c rdata(*,t) = Rain state on day t (t = (year-1)*nobs + day) c rdata is two 4 byte integers. The first gives a binary c coded decimal number which indicates which stations had rain. c The second is a binary coded decimal number which c indicates which stations have missing values for rainfall. c adata(k,t) = Atmospheric var k on day t (t = (year-1)*nobs + day) c Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'data.inc' include 'params.inc' include 'options.inc' include 'dims.inc' include 'fb.inc' include 'norm.inc' include 'dist.inc' logical lsets,lsetr integer iopt,ndays,dbeg,dend,iseed2 integer ipvt(MAXATM) real*8 work(MAXATM) real*8 loglik,rlik character*80 dfmt,datfil,disfil,fname write(*,3003) 3003 format( *'***************************************************************'/ *'| Hidden Markov Model Parameter Estimation |'/ *'| via the EM algorithm |'/ *'***************************************************************') write(*,*) 'Enter name of parameter control file.' read(*,501) datfil write(*,*) 'Reading parameters ...' call getpar(m,nsta,nyears,ndays,natm,atmopt,raiopt,nb, * datfil,dfmt,disfil,ifix,dis,dir, * gam,pipars,beta,mu,sigma,dbeg,dend,iest,iseed, * neps) c call prtmdl(m,nsta,nyears,ndays,natm,dbeg,dend,atmopt,raiopt, * iest,datfil,dfmt,disfil,ifix) call prtpar(6) c Do some setup and parameter transformations if (natm.gt.0) call invsig(MAXATM,natm,sigma,det,isigma,work,ipvt) call setup(raiopt) call rehash() if (raiopt.eq.3) call norman(m,nsta,nb,pipars,beta,lnorm) c Read the data write(*,*) write(*,*) 'Reading data ... ' call getdat(nyears,ndays,nsta,natm,dbeg,dend, * rdata,adata,dfmt,datfil) nobs = dend-dbeg+1 call sumdat(nyears,nobs,nsta,natm,rdata,adata) c if (ifix(1).eq.0.and.ifix(2).eq.0.and. * ifix(3).eq.0.and.ifix(4).eq.0) then write(*,*) 'All parameter values fixed' stop endif lsets = .TRUE. lsetr = .TRUE. c 100 write(*,*) 133 write(*,*) write(*,*) ' (1) Run EM algorithm from current parameters' write(*,*) ' (2) Evaluate likelihood at current parameters' write(*,*) ' (3) Compute covariance matrix of parameters' write(*,*) ' (4) Print current parameters to the screen' write(*,*) ' (5) Print current parameters to a file' write(*,*) ' (6) Run Viterbi algorithm for current parameters' write(*,*) ' (7) Generate conditional simulations from current * model' write(*,*) ' (9) Exit' read(*,*,err=133) iopt c goto (1000,2000,3000,4000,5000,6000,7000,8000,9000) iopt 1000 call emalg(lsets,lsetr) goto 133 2000 call forbac() rlik = loglik() write(*,*) 'Likelihood = ',rlik goto 133 3000 call forbac() call variance(lsets,lsetr) goto 133 4000 call prtpar(6) goto 133 5000 write(*,*) 'Enter file name ' read(*,501) fname open(12,file=fname,iostat=ierr,status="NEW") call prtpar(12) close(12) goto 133 6000 call viterbi() goto 133 7000 continue write(*,*) 'How many simulations?' read(*,*) nsim write(*,*) 'Enter file name for simulation(s) ' read(*,501) fname open(12,file=fname,iostat=ierr,status="NEW") 5 write(*,*) 'Enter a neg. integer to seed random number generator.' read(*,*) iseed2 if (iseed2.ge.0) goto 5 call condsim(nsim,iseed2) goto 133 8000 continue goto 133 9000 continue stop 501 format(a80) end c====================================================================== c This prints current parameter values c---------------------------------------------------------------------- subroutine prtpar(iunit) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' include 'options.inc' include 'norm.inc' include 'nsave.inc' integer i,j,k,lsamp(MAXM) character*5 par1(2),par2(2) data par1/'gam',' a '/, par2/' Mu ',' b '/ write(iunit,102) par1(atmopt) 102 format(/' ',a5) do 60 i=1,m write(iunit,6040)(gam(i,j),j=1,m) 6040 format( 6g12.4) 60 continue write(iunit,104) 104 format(/' Pipars') do 70 i=1,nsta write(iunit,6040)(pipars(2,i,j),j=1,m) 70 continue if (raiopt.eq.3) then write(iunit,105) 105 format(/' Beta') do 80 i = 1,nb 80 write(iunit,6040) (beta(i,j),j=1,m) endif write(iunit,106) par2(atmopt) 106 format(/' ',a5/' i j ') do 90 i=1,m do 90 j=1,m write(iunit,6070) i,j,(mu(i,j,k),k=1,natm) 6070 format(2i4, 6g12.4) 90 continue write(iunit,107) 107 format(/' Sigma') do 100 i=1,natm write(iunit,6040) (sigma(i,j),j=1,natm) 100 continue if (raiopt.eq.3) then do 200 i = 1,m 200 call icode(nsta,ista(1,i),lsamp(i),k) write(iunit,108) 108 format(/'Current value of the log(normalizing constants)') write(iunit,6050) (lnorm(i),i=1,m) 6050 format( 6e15.8) write(iunit,109) 109 format(/'A random no. from each state (use to restart EM)') write(iunit,6041) (lsamp(i),i=1,m) 6041 format(6(i10,x)) endif return end c====================================================================== c This converts gam and pipars to x c---------------------------------------------------------------------- subroutine gptox(ifix,m,nsta,natm,nb,gam,mu,pipars,beta, * xs,ns,xr,nr) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer m,nsta,natm,ns,nr,ifix(5),nb real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM),pipars(2,MAXSTA,MAXM) real*8 beta(MAXB,MAXM) real*8 xs(*),xr(*) call gtox(ifix,m,natm,gam,mu,xs,ns) call ptox(ifix,m,nsta,nb,pipars,beta,xr,nr) return end subroutine gtox(ifix,m,natm,gam,mu,xs,ns) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer m,natm,ns,ifix(5) integer i,j,k real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM) real*8 xs(*) c write(*,*) 'gtox' ns=0 if (ifix(1).ne.0) then do 10 i=1,m do 10 j=1,m ns = ns + 1 xs(ns) = gam(i,j) 10 continue endif C if (ifix(4).ne.0) then do 30 i=1,m do 30 j=1,m do 30 k=1,natm ns = ns + 1 xs(ns) = mu(i,j,k) 30 continue endif return end subroutine ptox(ifix,m,nsta,nb,pipars,beta,xr,nr) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer m,nsta,nr,ifix(5),nb integer i,j real*8 pipars(2,MAXSTA,MAXM),beta(MAXB,MAXM) real*8 xr(*) c write(*,*) 'ptox' nr=0 C do 30 i=1,m if (ifix(2).ne.0) then do 20 j=1,nsta nr = nr + 1 xr(nr) = pipars(1,j,i) 20 continue endif C if (ifix(3).ne.0) then do 25 j=1,nb nr = nr + 1 xr(nr) = beta(j,i) 25 continue endif c 30 continue return end c====================================================================== c This converts x to gam and pipars c---------------------------------------------------------------------- subroutine xtogp(ifix,m,nsta,natm,nb,gam,mu,pipars,beta,xs,xr) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer m,nsta,natm,ifix(5),nb real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM),pipars(2,MAXSTA,MAXM) real*8 beta(MAXB,MAXM) real*8 xs(*),xr(*) c write(*,*) 'xtogp' call xtog(ifix,m,natm,gam,mu,xs) call xtop(ifix,m,nsta,nb,pipars,beta,xr) return end subroutine xtop(ifix,m,nsta,nb,pipars,beta,xr) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer m,nsta,ifix(5),nb integer i,j,ind real*8 pipars(2,MAXSTA,MAXM) real*8 beta(MAXB,MAXM) real*8 xr(*) C ind = 0 do 30 i = 1,m if (ifix(2).ne.0) then do 20 j = 1,nsta ind=ind+1 pipars(2,j,i)=1.0/(1.0+exp(-xr(ind))) pipars(1,j,i)=xr(ind) 20 continue endif C if (ifix(3).ne.0) then do 25 j = 1,nb ind=ind+1 beta(j,i)=xr(ind) 25 continue endif c 30 continue return end subroutine xtog(ifix,m,natm,gam,mu,x) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer m,natm,ifix(5) integer i,j,ind,k real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM) real*8 x(*) c ind=0 if (ifix(1).ne.0) then do 10 i = 1,m do 10 j = 1,m ind=ind+1 gam(i,j) = x(ind) 10 continue endif C if (ifix(4).ne.0) then do 30 i=1,m do 30 j=1,m do 30 k=1,natm ind=ind+1 mu(i,j,k) = x(ind) 30 continue endif return end c====================================================================== c This gets initial parameter estimates c---------------------------------------------------------------------- subroutine getpar(m,nsta,nyears,ndays,natm,atmopt,raiopt,nb, * datfil,dfmt,disfil,ifix,dis,dir, * gam,pipars,beta,mu,sigma,dbeg,dend,iest,iseed, * neps) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'norm.inc' include 'nsave.inc' c integer m,nsta,nyears,ndays,natm,dbeg,dend,atmopt,nb integer raiopt,iseed integer i,j,ierr,k integer ifix(5),iest,lsamp(MAXM) real*8 gam(MAXM,*),pipars(2,MAXSTA,*),beta(MAXB,*),neps real*8 mu(MAXM,MAXM,MAXATM),sigma(MAXATM,MAXATM) real*8 dis(MAXSTA,MAXSTA),dir(MAXSTA,MAXSTA) character*80 dfmt,datfil,disfil c open(11,file=datfil,iostat=ierr,status="OLD") if (ierr.ne.0) stop 'Error opening parameter control file' read(11,*) m, nsta, nyears, ndays, natm, atmopt, raiopt if (m.gt. MAXM) stop 'Error. Too many states. max = MAXM ' if (nsta.gt. MAXSTA) * stop 'Error. Too many stations. max = MAXSTA ' if (natm.gt. MAXATM) stop 'Error. Too many atm vars. max = MAXATM' if (atmopt.ne.1) stop 'Error. Only works with Bayes model.' if (raiopt.eq.2) stop 'Error. Saturated model not implemented.' read(11,'(a80)') datfil read(11,'(a80)') dfmt if (raiopt.eq.3) then read(11,'(i5)') nb read(11,'(a80)') disfil open(12,file=disfil,status='OLD',iostat=ierr) if (ierr.ne.0) stop 'Error opening distance file' do 15 i=1,nsta 15 read(12,*)(dis(i,j),j=1,nsta) do 25 i=1,nsta 25 read(12,*)(dir(i,j),j=1,nsta) endif c read(11,*) ifix(1) do 10 i=1,m read(11,*)(gam(i,j),j=1,m) 10 continue read(11,*) ifix(2) do 20 i=1,nsta read(11,*)(pipars(2,i,j),j=1,m) do 120 j=1,m pipars(1,i,j) = dlog(pipars(2,i,j)/(1.0-pipars(2,i,j))) 120 continue 20 continue if (raiopt.eq.3) then read(11,*) ifix(3) do 28 i = 1,nb 28 read(11,*) (beta(i,j),j=1,m) else ifix(3) = 0 endif if (natm.gt.0) then read(11,*) ifix(4) do 30 i=1,m do 30 j=1,m read(11,*)(mu(i,j,k),k=1,natm) 30 continue read(11,*) ifix(5) if (ifix(5).ne.0) stop 'Error. Sigma cannot be estimated' do 40 i=1,natm 40 read(11,*)(sigma(i,j),j=1,natm) else ifix(4) = 0 ifix(5) = 0 endif read(11,*)dbeg,dend if (nyears*(dbeg-dend+1).gt.MAXOBS) * stop 'ERROR. too many observations' read(11,*) iest if (iest.ge.2) then read(11,*) iseed,neps read(11,*) (norm0(j),j=1,m) read(11,*) (lsamp(j),j=1,m) do 410 j = 1,m lnorm(j) = norm0(j) do 400 i = 1,nsta pipar0(1,i,j) = pipars(1,i,j) pipar0(2,i,j) = pipars(2,i,j) 400 continue do 405 i = 1,nb 405 beta0(i,j) = beta(i,j) call vbit(ista(1,j),lsamp(j),nsta) 410 continue endif return end subroutine icode(nsta,ivec,idat,imis) c c Function to code the rain data into two integers. c The lowest nsta bits of the first integer are cleared (0) or set (1) c according to the rain state at station i (i=1,nsta). The lowest c nsta bits in the second integer are cleared or set according to whether c or not the rain data for station i (i=1,nsta) is missing (0 = not c missing 1 = missing). Since integers are 32 bits typically, there c can be at most 31 rain stations with this scheme (don't want to overwrite c the sign bit). c Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer nsta,ivec(*),idat,imis,i,i0 c idat = 0 imis = 0 do 10 i=1,nsta i0 = i-1 call setbit(i0,imis,ivec(i).eq.9) call setbit(i0,idat,ivec(i).eq.1) 10 continue return end subroutine rehash() Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'hash.inc' do 1 i=0, MAXPI1 hvec(1,i) = -1 hvec(2,i) = -1 1 continue ierrh = 0 return end subroutine prtmdl(m,nsta,nyears,ndays,natm,dbeg,dend, * atmopt,raiopt,iest, * datfil,dfmt,disfil,ifix) Implicit double precision (A-H,O-Z) integer m,nsta,nyears,ndays,natm,dbeg,dend,atmopt,raiopt,iest integer ifix(5) character*80 datfil,dfmt,disfil character*5 par1(2),par2(2) character*8 atmmdl(2) character*11 raimdl(3) character*4 estmdl(3) data par1/'gam',' a '/, par2/' Mu ',' b '/ data atmmdl/' Bayes ','Logistic'/ data raimdl/'Independent',' Saturated ',' Spatial '/ data estmdl/'MLE ','MCML','MPLE'/ c write(6,101) m,nsta,nyears,ndays,natm,dbeg,dend, * atmmdl(atmopt), * raimdl(raiopt),estmdl(iest), * datfil,dfmt if (raiopt.eq.3) write(6,107) disfil 101 format(/' Model: ',i2,' states, ',i2,' stations',/, * ' ',i3,' years of data with ',i3,' obs/year',/, * ' ',i2,' atmospheric variables',/, * ' Season: day ',i3,' to day ',i3,/, * ' Model for state transitions: ',a8,/, * ' Model for rain states: ',a11,/, * ' Estimation method: ',a4,//, * ' Data file: ',/,a80,/, * ' Format: ',/,a80) 107 format(/' Distance file: ',/,a80) c write(6,102) 102 format(//' Initial values'/) if (ifix(1).eq.0) write(6,1021) par1(atmopt) 1021 format('Note: ',a5,' is fixed') if (ifix(2).eq.0) write(6,*) 'Note: Pipars is fixed' if (raiopt.eq.3.and.ifix(3).eq.0) write(6,*) 'Note: Beta is fixed' if (ifix(4).eq.0) write(6,1022) par2(atmopt) 1022 format('Note: ',a5,' is fixed') if (atmopt.eq.1.and.ifix(5).eq.0)write(6,*) 'Note: Sigma is fixed' return end subroutine getdat(nyears,ndays,nsta,natm,dbeg,dend, * rdata,adata,dfmt,datfil) Implicit double precision (A-H,O-Z) include 'nhmm.inc' logical first integer nyears,year,ndays,day,dbeg,dend,nsta,natm integer rdata(2,MAXOBS),itemp(MAXSTA) real*8 adata(MAXATM,MAXOBS) character*60 dfmt,datfil c open(10,file=datfil,status="OLD") t=0 do 30 year=1,nyears call skip(10,dbeg-1) do 25 day=dbeg,dend t=t+1 if (dfmt.eq.'*') then read(10,*,end=99) (itemp(i),i=1,nsta),(adata(i,t),i=1,natm) else read(10,dfmt,end=99) (itemp(i),i=1,nsta),(adata(i,t),i=1,natm) endif call icode(nsta,itemp,idat,imis) rdata(1,t)=idat rdata(2,t)=imis c if (t.eq.1) then write(*,*) ' First line is ...' write(*,*) (itemp(i),i=1,nsta),(adata(i,1),i=1,natm) first = .FALSE. endif 25 continue call skip(10,ndays-dend) 30 continue 99 return end subroutine skip(iunit,nrec) integer iunit,nrec,i c do 10 i=1,nrec read(iunit,*) 10 continue return end subroutine sumdat(nyears,nobs,nsta,natm,rdata,adata) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer t,nyears,nobs,nsta,natm,year,day,i0,i integer rdata(2,MAXOBS) integer status(3,MAXSTA) integer bit real*8 adata(MAXATM,MAXOBS) real*8 s1(MAXATM),s2(MAXATM),amin(MAXATM),amax(MAXATM) c c Summarize the rain and atmospheric data in the following manner c c Rain data c station c 1 2 3 4 5 6 .... nsta c dry(0) c Status wet(1) c miss(9) c c Atmospheric data c variable c name(1) name(2) name(3) .... name(natm) c mean c std dev c median c min c max c t=0 do 5 i=1,natm s1(i) = 0.0 s2(i) = 0.0 amin(i) = adata(i,1) amax(i) = adata(i,1) 5 continue do 100 year = 1,nyears do 100 day = 1,nobs t = t+1 do 10 i=1,nsta i0 = i-1 if (bit(i0,rdata(1,t)).eq.1) then status(2,i) = status(2,i) + 1 elseif (bit(i0,rdata(2,t)).eq.1) then status(3,i) = status(3,i) + 1 else status(1,i) = status(1,i) + 1 endif 10 continue do 20 i = 1,natm s1(i) = s1(i) + adata(i,t) s2(i) = s2(i) + adata(i,t)*adata(i,t) if (adata(i,t).lt.amin(i)) amin(i) = adata(i,t) if (adata(i,t).gt.amax(i)) amax(i) = adata(i,t) 20 continue 100 continue do 30 i=1,natm s1(i) = s1(i)/t s2(i) = (s2(i) - t*s1(i)*s1(i))/(t-1) 30 continue write(*,*) write(*,*) "Rainfall Summary" write(*,*) "Station Dry Wet Miss" write(*,"(i2,5x,3i5)") (i,(status(j,i),j=1,3),i=1,nsta) write(*,*) write(*,*) "Atmospheric Data Summary" write(*,*) " Mean SD Min Max" do 40 i=1,natm 40 write(*,"(i2,1x,4g10.4)") i,s1(i),sqrt(s2(i)),amin(i),amax(i) return end real*8 function l2norm(xs,xr,gam,mu,pipars,beta,ifix) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' c integer ind,i,j,k,ifix(5) real*8 xs(MAXNS),xr(MAXNR) real*8 gam(MAXM,*),pipars(2,MAXSTA,*),beta(MAXB,*) real*8 mu(MAXM,MAXM,MAXATM) c l2norm=0.0 ind = 0 if (ifix(1).ne.0) then do 10 i = 1,m do 10 j = 1,m ind=ind+1 l2norm = l2norm + ((gam(i,j) - xs(ind)))**2 10 continue endif C if (ifix(4).ne.0) then do 30 i=1,m do 30 j=1,m do 30 k=1,natm ind=ind+1 l2norm = l2norm + ((mu(i,j,k) - xs(ind)))**2 30 continue endif C ind = 0 do 40 i = 1,m if (ifix(2).ne.0) then do 20 j = 1,nsta ind=ind+1 l2norm = l2norm + ((pipars(1,j,i)-xr(ind)))**2 20 continue endif C if (ifix(3).ne.0) then do 25 j = 1,nb ind=ind+1 l2norm = l2norm + ((beta(j,i)-xr(ind)))**2 25 continue endif 40 continue l2norm = sqrt(l2norm) return end integer function bit(i0,inum) integer i0,inum c return 1 if the i0 bit of inum is set, 0 otherwise if (btest(inum,i0)) then bit=1 else bit=0 endif return end subroutine vbit(ib,inum,n) integer ib(*),inum,n c do 10 i=1,n i0 = i-1 if (btest(inum,i0)) then ib(i) = 1 else ib(i) = 0 endif 10 continue return end subroutine setbit(i0,inum,test) integer i0,inum logical test c set the i0 bit of inum if test is .TRUE. otherwise clear it. if (test) then inum = ibset(inum,i0) else inum = ibclr(inum,i0) endif return end nhmm.inc000644 000423 000017 00000001536 06471075507 013224 0ustar00hughesusers000000 031407 parameter(MAXM=7, MAXATM=5, MAXSTA=31, MAXB=5, MAXOBS=6000) parameter(MAXSTA2=MAXSTA*(MAXSTA-1)/2) parameter(MAXMA=MAXM*(MAXATM+1)) c MAXPI and MAXD must be prime numbers; used in hashing routines parameter(MAXPI=3697, MAXD=14983) parameter(MAXM1=MAXM+1, MAXPI1=MAXPI-1, MAXD1=MAXD-1) parameter(MAXLIN=MAXM*(MAXATM+1)+10,MAXNLN=1) parameter(MAXNS=MAXM*MAXM*(MAXATM+1),MAXNR=MAXM*(MAXSTA+MAXB)) parameter(MAXN=max(MAXNS,MAXNR)) parameter(MAXPAR=MAXN+MAXLIN+MAXNLN) c see NAG/NPSOL documentation for MAXIW and MAXW parameter(MAXIW=3*MAXN+MAXLIN+2*MAXNLN) parameter(MAXW=20*MAXN*(MAXN+1)+11*MAXLIN) parameter(MAXP2=MAXNR+1+1) c MAXSEQ = 10000*log(MAXSTA) parameter(MAXSEQ=34340) parameter(MAXTOT=MAXNS+MAXNR) parameter(MAXIS=MAXM*(MAXATM+1),MAXIR=MAXSTA+MAXB) norm.inc000644 000423 000017 00000000063 06242170133 013222 0ustar00hughesusers000000 030755 real*8 lnorm(MAXM) common /norm/ lnorm npsetup.inc000644 000423 000017 00000000536 06242170133 013753 0ustar00hughesusers000000 030756 integer nclins,ncnlns integer leniws,lenws integer nclinr,ncnlnr integer leniwr,lenwr real*8 as(MAXLIN,MAXNS),bls(MAXPAR),bus(MAXPAR) real*8 ar(1,MAXNR),blr(MAXP2),bur(MAXP2) common /npsetup/ nclins,ncnlns,nclinr,ncnlnr,leniws,lenws, * leniwr,lenwr,as,bls,bus,ar,blr,bur nsave.inc000644 000423 000017 00000000215 06242170134 013365 0ustar00hughesusers000000 030757 integer ista(MAXSTA,MAXM) real*8 pipar0(2,MAXSTA,MAXM),norm0(MAXM),beta0(MAXB,MAXM) common /nsave/ pipar0,norm0,beta0,ista objfun.f000644 000423 000017 00000057076 06470172152 013222 0ustar00hughesusers000000 031420 subroutine OBJFUNS(mode,n,x,objf,objgr,nstate,iuser,user) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'fb.inc' include 'params.inc' include 'dims.inc' include 'data.inc' include 'weights.inc' c integer mode,n,nstate integer t,i,ii,year,day,s0,s1,nchk integer istate,iuser(*) real*8 objf,x(*),objgr(*),user(*) real*8 delta(MAXM),a(MAXM,MAXM) real*8 sum(MAXM),eQ(MAXM,MAXM),denom(MAXM) real*8 Sgm(MAXM,MAXNS) real*8 newgam(MAXM,MAXM),newmu(MAXM,MAXM,MAXATM) c c write(*,*) (x(i),i=1,n) call newgm(m,natm,x,newgam,newmu) c t = 0 objf=0.0 do 130 i = 1,n 130 objgr(i) = 0.0 c loop over time do 300 year = 1,nyears t=t+1 c day = 1 call dlset1(t,newgam,newmu,isigma,det,sum,eQ,denom,delta,a) call dldgm1(t,ifix,newgam,newmu,isigma,det,sum,eQ,denom, * delta,a,Sgm,nchk) do 220 s0=1,m objf = objf - c(s0,t)*log(delta(s0)) do 220 i=1,n objgr(i) = objgr(i) - c(s0,t)*Sgm(s0,i) 220 continue c day > 1 do 300 day = 2,nobs t=t+1 call dlset2(t,newgam,newmu,isigma,det,sum,eQ,denom,a) do 305 s0 = 1,m do 305 s1 = 1,m objf = objf - c2(s0,s1,t)*log(a(s0,s1)) 305 continue call dldgm2(t,ifix,newgam,newmu,isigma,det,sum,eQ,denom, * a,Sgm,nchk) do 310 i=1,n ii = istate(i) do 310 s = 1,m objgr(i) = objgr(i) - c2(ii,s,t)*Sgm(s,i) 310 continue 300 continue c write(*,*) objf,(objgr(i),i=1,n) return end integer function istate(i) include 'dims.inc' integer i c if (i.le.m*m) then istate = (i-1)/m + 1 else istate = (i-m*m-1)/(natm*m) + 1 endif return end integer function istate2(i) include 'dims.inc' integer i c istate2 = (i-1)/(nsta+nb) + 1 return end real*8 function dQdmu(t,i,j,k,isigma,mu) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'data.inc' integer i,j,k,l,t real*8 isigma(MAXATM,MAXATM),mu(MAXM,MAXM,MAXATM) dQdmu = 0.0 do 10 l=1,natm 10 dQdmu = dQdmu + isigma(k,l)*(adata(l,t)-mu(i,j,l)) return end subroutine objfr0(n,s,x,dx,bl,bu) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'fb.inc' include 'dims.inc' include 'data.inc' include 'weights.inc' c integer n,s integer ib1(MAXSTA),ib2(MAXSTA),year,day,i,t real*8 x(n),dx(n),bl(n),bu(n) real*8 p,zero real*8 num(MAXSTA),denom(MAXSTA) data zero/0.0D0/ c do 100 i=1,n num(i) = zero denom(i) = zero dx(i) = zero 100 continue t = 0 do 110 year = 1,nyears do 110 day = 1,nobs t = t+1 call vbit(ib1,rdata(1,t),nsta) call vbit(ib2,rdata(2,t),nsta) do 110 i=1,n p = 1.0/(1.0+exp(-x(i))) if (ib2(i).eq.1) then num(i) = num(i) + c(s,t)*p denom(i) = denom(i) + c(s,t)*(1-p) else if (ib1(i).eq.1) then num(i) = num(i) + c(s,t) dx(i) = dx(i) + c(s,t)/p else denom(i) = denom(i) + c(s,t) dx(i) = dx(i) - c(s,t)/(1-p) endif endif 110 continue do 120 i=1,n x(i) = log(num(i)/denom(i)) if (x(i).gt.bu(i)) x(i)=bu(i) if (x(i).lt.bl(i)) x(i)=bl(i) 120 continue return end subroutine OBJFR1(mode,n,s,nrowr,x,objf,objgr,rr) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'options.inc' include 'dims.inc' include 'data.inc' include 'weights.inc' include 'nsave.inc' integer mode,n,s,nrowr real*8 objf,x(*),objgr(*),rr(nrowr,*) c integer t,year,day integer i,ind,ib1(MAXSTA) real*8 p,dpda(MAXSTA),dpdb(MAXB),b(MAXSTA) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 exp(MAXSTA),expb(MAXB),lnorm real*8 sum,wt real*8 zero real*8 newpi(MAXSTA),newbeta(MAXB) real*8 dbetadb data zero/0.0d0/ c call newpb(nsta,nb,s,ifix,x,pipars,beta,newpi,newbeta) call rehash() c objf = zero if (mode.eq.0) then do 50 i = 1,nsta dpda(i) = zero 50 continue do 60 i = 1,nb dpdb(i) = zero 60 continue endif c if (iest.eq.1) then call deriv1(mode,nsta,nb,s,newpi,newbeta,lnorm,exp,expb, * dp2da2,dp2dab,dp2db2) elseif (iest.eq.2) then call deriv2(mode,nsta,nb,s,ista(1,s),iseed,neps, * newpi,newbeta,lnorm,exp,expb, * dp2da2,dp2dab,dp2db2,pipar0,beta0,norm0) endif c wt = zero t = 0 c loop over time do 200 year = 1,nyears do 200 day = 1,nobs t=t+1 if (rdata(2,t).gt.0) stop 'Missing data not allowed' call vbit(ib1,rdata(1,t),nsta) call fbetas(newbeta,s,nsta,nb,ib1,b) p = zero do 140 i = 1,nsta if (ib1(i).eq.1) p = p + newpi(i)+b(i)/2.0 140 continue objf = objf - c(s,t)*(p - lnorm) if (mode.eq.0) then wt = wt + c(s,t) do 150 i = 1,nsta dpda(i) = dpda(i) - c(s,t)*(ib1(i) - exp(i)) 150 continue do 170 k = 1,nb sum = 0.0 do 160 i = 1,nsta-1 do 160 j = i+1,nsta if (ib1(i).eq.1.and.ib1(j).eq.1) * sum = sum + dbetadb(newbeta,s,k,i,j) 160 continue dpdb(k) = dpdb(k) - c(s,t)*(sum-expb(k)) 170 continue endif 200 continue c if (mode.eq.0) then ind = 0 if (ifix(2).eq.1) then do 260 i = 1,nsta ind=ind+1 objgr(ind) = dpda(i) ind1 = ind-1 do 260 j = i,nsta ind1 = ind1+1 rr(ind,ind1) = wt*dp2da2(i,j) rr(ind1,ind) = rr(ind,ind1) 260 continue endif if (ifix(3).eq.1) then do 280 k = 1,nb ind=ind+1 objgr(ind) = dpdb(k) ind1 = ind-1 do 270 l = k,nb ind1 = ind1+1 rr(ind,ind1) = wt*dp2db2(k,l) rr(ind1,ind) = rr(ind,ind1) 270 continue if (ifix(2).eq.1) then ind1 = 0 do 275 i = 1,nsta ind1 = ind1+1 rr(ind,ind1) = wt*dp2dab(k,i) rr(ind1,ind) = rr(ind,ind1) 275 continue endif 280 continue endif endif c 999 continue c write(*,*) (x(i),i=1,n) c write(*,*) mode,objf c write(*,*) (objgr(i),i=1,n) return end subroutine begder(nsta,nb,expect,expectb,dp2da2,dp2dab,dp2db2) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer nsta,nb real*8 expect(*),expectb(*) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 zero data zero/0.0d0/ c do 100 i = 1,nsta expect(i) = zero do 105 j = 1,nsta 105 dp2da2(i,j) = zero do 100 k = 1,nb 100 dp2dab(k,i) = zero do 110 k = 1,nb expectb(k) = zero do 110 l = 1,nb 110 dp2db2(k,l) = zero return end subroutine endder(nsta,nb,denom,expect,expectb, * dp2da2,dp2dab,dp2db2) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer nsta,nb real*8 expect(*),denom,expectb(*) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) c do 160 i = 1,nsta expect(i) = expect(i)/denom do 160 j = 1,i dp2da2(i,j) = dp2da2(i,j)/denom - expect(i)*expect(j) dp2da2(j,i) = dp2da2(i,j) 160 continue do 180 k = 1,nb expectb(k) = expectb(k)/denom do 170 l = 1,k dp2db2(k,l)= (dp2db2(k,l)/denom)-expectb(k)*expectb(l) dp2db2(l,k) = dp2db2(k,l) 170 continue do 175 i = 1,nsta dp2dab(k,i) = (dp2dab(k,i)/denom)- expect(i)*expectb(k) 175 continue 180 continue return end subroutine moments(nsta,nb,s,beta,p,ib,expect,expectb, * dp2da2,dp2dab,dp2db2) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer nsta,nb,ib(*),s real*8 p,expect(*),expectb(*),beta(*) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 sum(MAXB),sum2(MAXB,MAXB),zero real*8 dbetadb,d2bdb2 data zero/0.0d0/ c do 125 k = 1,nb sum(k) = zero do 125 h = k,nb sum2(k,h) = zero 125 continue c do 135 i = 1,nsta if (ib(i).eq.1) then expect(i) = expect(i) + p dp2da2(i,i) = dp2da2(i,i) + p do 140 j = i+1,nsta if (ib(j).eq.1) then dp2da2(i,j) = dp2da2(i,j) + p dp2da2(j,i) = dp2da2(i,j) do 145 k = 1,nb sum(k) = sum(k) + dbetadb(beta,s,k,i,j) do 145 h = k,nb sum2(k,h) = sum2(k,h) + d2bdb2(beta,s,k,h,i,j) 145 continue endif 140 continue endif 135 continue c do 155 k = 1,nb expectb(k) = expectb(k)+sum(k)*p do 150 h = k,nb dp2db2(k,h) = dp2db2(k,h)+(sum(k)*sum(h)+sum2(k,h))*p dp2db2(h,k) = dp2db2(k,h) 150 continue do 155 i =1,nsta if (ib(i).eq.1) then dp2dab(k,i) = dp2dab(k,i)+sum(k)*p endif 155 continue c return end subroutine moment2(nsta,nb,s,beta,p,expect,expectb, * dp2da2,dp2dab,dp2db2) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer nsta,nb,s integer i,j,ii,jj,k,h real*8 p(*),expect(*),expectb(*) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 dbetadb,zero data zero/0.0d0/ c do 210 i=1,nsta expect(i) = p(i) dp2da2(i,i) = p(i)*(1.0-p(i)) 210 continue do 220 k = 1,nb expectb(k) = zero do 225 i = 1,nsta dp2dab(k,i) = zero do 230 j = 1,nsta expectb(k)=expectb(k)+p(i)*p(j)*dbetadb(beta,s,k,i,j)/2. dp2dab(k,i)=dp2dab(k,i)+p(j)*dbetadb(beta,s,k,i,j) 230 continue dp2dab(k,i)=dp2dab(k,i)*p(i)*(1.0-p(i)) 225 continue do 235 h = k,nb dp2db2(k,h) = zero do 235 i = 1,nsta do 235 j = 1,nsta do 235 ii = 1,nsta do 235 jj = 1,nsta if (i.ne.j.and.ii.ne.jj) then if (i.eq.ii) then if (j.eq.jj) then dp2db2(k,h) = dp2db2(k,h)+p(i)*p(j)*(1.0-p(i)*p(j))* * dbetadb(beta,s,k,i,j)*dbetadb(beta,s,h,ii,jj)/4. else dp2db2(k,h) = dp2db2(k,h)+p(i)*p(j)*p(jj)*(1.0-p(i))* * dbetadb(beta,s,k,i,j)*dbetadb(beta,s,h,ii,jj)/4. endif endif if (i.eq.jj) then if (j.eq.ii) then dp2db2(k,h) = dp2db2(k,h)+p(i)*p(j)*(1.0-p(i)*p(j))* * dbetadb(beta,s,k,i,j)*dbetadb(beta,s,h,ii,jj)/4. else dp2db2(k,h) = dp2db2(k,h)+p(i)*p(j)*p(ii)*(1.0-p(i))* * dbetadb(beta,s,k,i,j)*dbetadb(beta,s,h,ii,jj)/4. endif endif if (j.eq.jj) then if (i.ne.ii) then dp2db2(k,h) = dp2db2(k,h)+p(i)*p(j)*p(ii)*(1.0-p(j))* * dbetadb(beta,s,k,i,j)*dbetadb(beta,s,h,ii,jj)/4. endif endif if (j.eq.ii) then if (i.ne.jj) then dp2db2(k,h) = dp2db2(k,h)+p(i)*p(j)*p(jj)*(1.0-p(j))* * dbetadb(beta,s,k,i,j)*dbetadb(beta,s,h,ii,jj)/4. endif endif endif 235 continue dp2db2(h,k) = dp2db2(k,h) 220 continue return end subroutine deriv1(mode,nsta,nb,s,pipars,beta,lnorm, * expect,expectb,dp2da2,dp2dab,dp2db2) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer nsta,s,mode,nb integer ib(MAXSTA),i,l real*8 pipars(*),beta(*),lnorm,expect(*),expectb(*) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 zero,p,b(MAXSTA),denom data zero/0.0d0/ c denom = zero if (mode.eq.0) call begder(nsta,nb,expect,expectb, * dp2da2,dp2dab,dp2db2) c do 120 l = 0, 2**nsta-1 call vbit(ib,l,nsta) call fbetas(beta,s,nsta,nb,ib,b) p = zero do 130 i = 1,nsta if (ib(i).eq.1) p = p + pipars(i)+b(i)/2.0 130 continue p = exp(p) denom = denom + p if (mode.eq.0) call moments(nsta,nb,s,beta,p,ib,expect,expectb, * dp2da2,dp2dab,dp2db2) 120 continue c if (mode.eq.0) call endder(nsta,nb,denom,expect,expectb, * dp2da2,dp2dab,dp2db2) lnorm = log(denom) return end subroutine deriv2(mode,nsta,nb,s,ista,iseed,neps, * pipars,beta,lnorm,expect,expectb, * dp2da2,dp2dab,dp2db2,pipar0,beta0,norm0) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer nsta,s,mode,nb,iseed,ista(*) integer i,imode real*8 pipars(*),beta(*),expect(*),expectb(*),lnorm real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 pipar0(2,MAXSTA,MAXM),beta0(MAXB,MAXM),norm0(MAXM) real*8 sum0 real*8 zero,p(MAXSTA) real*8 err,dnorm,step(MAXNR),neps,ep real*8 ran1 logical lbeta data zero/0.0d0/ c if (lbeta(beta)) then c Compute norm analytically by assuming independence; lnorm = zero if (mode.eq.0) call begder(nsta,nb,expect,expectb, * dp2da2,dp2dab,dp2db2) do 200 i=1,nsta ep = exp(pipars(i)) p(i) = ep/(1.0+ep) lnorm = lnorm + log(1.0+ep) ista(i) = 0 if (ran1(iseed).le.p(i)) ista(i)=1 200 continue if (mode.eq.0) call moment2(nsta,nb,s,beta,p,expect,expectb, * dp2da2,dp2dab,dp2db2) else c Estimate change in norm from previous value via Monte Carlo lnorm = norm0(s) call stepsiz(nsta,nb,s,pipars,beta,pipar0,beta0,step,sum0) if (sum0.eq.0.0.and.mode.ne.0) return n = min0(int(10000*log(float(nsta))),MAXSEQ) i1 = 1 i0 = 0 write(*,*) 'Estimating normalizing constant for state ',s 300 continue if (i1.eq.1) then imode = mode else imode = 2 endif call delnrm(imode,nsta,nb,n,iseed,s,pipar0,beta0,step,i0, * expect,expectb,dp2da2,dp2dab,dp2db2, * ista,dnorm,err) write(*,*) i1,dnorm,err if (err.le.neps) then i1 = i1 - 1 i0 = i0 + 1 lnorm = lnorm + log(dnorm) if (i1.eq.0) goto 310 else i0 = 2*i0 i1 = 2*i1 do 305 i = 1,nsta+nb 305 step(i) = step(i)/2.0 endif goto 300 c 310 continue endif return end subroutine OBJFR3(mode,n,s,nrowr,x,objf,objgr,rr) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' include 'data.inc' include 'weights.inc' integer mode,n,nrowr,s real*8 objf,x(*),objgr(*),rr(nrowr,*) integer i,ind,ind1,t,year,day,ib1(MAXSTA) real*8 p,ex real*8 dpda(MAXSTA),dpdb(MAXB) real*8 dp2da2(MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 zero real*8 newpi(MAXSTA),newbeta(MAXB) real*8 v,term,b(MAXSTA),db(MAXB,MAXSTA) data zero/0.0d0/ c call newpb(nsta,nb,s,ifix,x,pipars,beta,newpi,newbeta) call rehash() c objf = zero if (mode.eq.0) then do 50 i = 1,nsta dp2da2(i) = zero dpda(i) = zero do 50 k = 1,nb dp2dab(k,i) = zero 50 continue do 60 i = 1,nb dpdb(i) = zero do 60 j = 1,nb dp2db2(i,j) = zero 60 continue endif c loop over time t = 0 do 200 year = 1,nyears do 200 day = 1,nobs t=t+1 if (rdata(2,t).gt.0) stop 'Missing data not allowed' call vbit(ib1,rdata(1,t),nsta) call fbetas(newbeta,s,nsta,nb,ib1,b) call dbdb(MAXB,newbeta,nb,nsta,s,ib1,db) c loop over stations do 200 i=1,nsta ex = exp(newpi(i)+b(i)) p = ex/(1+ex) if (ib1(i).eq.1) then objf = objf - c(s,t)*log(p) else objf = objf - c(s,t)*log(1.0-p) endif if (mode.eq.0) then term = c(s,t)*(ib1(i)-p) v = c(s,t)*p*(1-p) dpda(i) = dpda(i) - term dp2da2(i) = dp2da2(i) + v do 205 k = 1,nb dpdb(k) = dpdb(k) - db(k,i)*term dp2dab(k,i) = dp2dab(k,i) + db(k,i)*v do 205 l = k,nb dp2db2(k,l) = dp2db2(k,l) + db(k,i)*db(l,i)*v 205 continue endif 200 continue c if (mode.eq.0) then ind=0 if (ifix(2).eq.1) then do 260 i = 1,nsta ind=ind+1 objgr(ind) = dpda(i) rr(ind,ind) = dp2da2(i) ind1 = ind do 260 j = i+1,nsta ind1 = ind1+1 rr(ind,ind1) = zero rr(ind1,ind) = zero 260 continue endif if (ifix(3).eq.1) then do 280 k = 1,nb ind=ind+1 objgr(ind) = dpdb(k) ind1 = ind-1 do 270 l = k,nb ind1 = ind1+1 rr(ind,ind1) = dp2db2(k,l) rr(ind1,ind) = dp2db2(k,l) 270 continue if (ifix(2).eq.1) then ind1 = 0 do 275 i = 1,nsta ind1 = ind1+1 rr(ind,ind1) = dp2dab(k,i) rr(ind1,ind) = dp2dab(k,i) 275 continue endif 280 continue endif endif c 999 continue c write(*,*) (x(i),i=1,n) c write(*,*) mode,objf c write(*,*) (objgr(i),i=1,n) return end subroutine missum1(m,nsta,nb,ib1,ib2,pipars,beta,s,sum1, * sumb,sum2) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer m,nsta,ib1(*),ib2(*),s integer i,j,iii(MAXSTA),nmis real*8 pipars(*),beta(*) real*8 sumb(*),sum1(*),sum2,sum real*8 zero,p,b(MAXSTA),dbetadb data zero/0.0d0/ c call whatmiss(nsta,ib2,iii,nmis) if (nmis.eq.0) goto 999 sum2 = zero do 185 i = 1,nsta 185 sum1(i) = zero do 190 i = 1,nb 190 sumb(i) = zero do 195 l=0,2**nmis-1 call ifill(nmis,l,iii,ib1) call fbetas(beta,s,nsta,nb,ib1,b) p = zero do 198 i=1,nsta if (ib1(i).eq.1) p = p + pipars(i) + b(i)/2.0 198 continue p = exp(p) sum2 = sum2 + p do 197 k = 1,nb sum = zero do 196 i = 1,nsta do 196 j = i+1,nsta if (ib1(i).eq.1.and.ib1(j).eq.1) * sum = sum + dbetadb(beta,s,k,i,j) 196 continue sumb(k) = sumb(k) + p*sum 197 continue do 199 i=1,nmis 199 sum1(iii(i)) = sum1(iii(i)) + p*ib1(iii(i)) 195 continue 999 return end subroutine augment(lda,m,a) real*8 a(lda,m),temp c do 10 i = 1,m-1 a(i,i) = 1.0 - a(i,i) do 10 j = i+1,m temp = a(i,j) a(i,j) = -a(j,i) a(j,i) = -temp 10 continue do 20 j = 1,m a(m,j) = -1.0 20 continue return end subroutine dAdgam(i,s1,lda,m,sum,denom,gam,eQ,aprime) integer i,j,s1,m,lda real*8 sum(m),denom(m),gam(lda,m),eQ(lda,m),aprime(m) c derivative of A wrt gamma(i,s1) c return i'th row of aprime only do 30 j = 1,s1-1 30 aprime(j) = -gam(i,j)*eQ(i,j)*eQ(i,s1)/denom(i) aprime(s1) = eQ(i,s1)*(sum(i)-gam(i,s1)*eQ(i,s1))/denom(i) do 35 j = s1+1,m 35 aprime(j) = -gam(i,j)*eQ(i,j)*eQ(i,s1)/denom(i) return end subroutine dAdmu(i,s1,lda,m,sum,denom,gam,eQ,aprime) integer i,j,s1,m,lda real*8 sum(m),denom(m),gam(lda,m),eQ(lda,m),aprime(m) c derivative of A wrt mu(i,s1,*) c return i'th row of aprime only do 30 j = 1,s1-1 30 aprime(j) = -gam(i,j)*gam(i,s1)*eQ(i,j)*eQ(i,s1)/denom(i) aprime(s1) = gam(i,s1)*eQ(i,s1)*(sum(i)-gam(i,s1)*eQ(i,s1)) * /denom(i) do 35 j = s1+1,m 35 aprime(j) = -gam(i,j)*gam(i,s1)*eQ(i,j)*eQ(i,s1)/denom(i) return end subroutine ddelta(i,lda,m,ia,ap,delta,dp) integer i,lda,m integer j real*8 ia(lda,*),ap(*),delta(*),dp(*) c do 10 s = 1,m dp(s) = 0.0 do 15 j = 1,m-1 15 dp(s) = dp(s) + ia(s,j)*ap(j) dp(s) = dp(s) * delta(i) 10 continue return end subroutine weight1(t,rlik,c1) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'fb.inc' integer t integer s1 real*8 rlik,c1(*) real*8 scale1 c c c1 = P(W_t1 | Y) c scale1 = exp(forw(m+1,t)+back(m+1,t)-rlik) do 312 s1=1,m c1(s1) = forw(s1,t)*back(s1,t)*scale1 312 continue return end subroutine weight2(t,rlik,c2) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer t integer s0,s1 real*8 rlik,c2(MAXM,*),pi(MAXM),a(MAXM,MAXM) real*8 scale2 c c c2 = p(W_t1, W_{t1-1} | Y) c scale2 = exp(forw(m+1,t-1)+back(m+1,t)-rlik) call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pi) call makeA(m,gam,natm,mu,isigma,det,adata(1,t),a) do 310 s1=1,m do 310 s0=1,m c2(s1,s0) = forw(s1,t-1)*a(s1,s0)*pi(s0)*back(s0,t)*scale2 310 continue return end subroutine newpb(nsta,nb,s,ifix,x,pipars,beta,newpi,newbeta) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer nsta,nb,s integer i,ifix(5) real*8 pipars(2,MAXSTA,MAXM),beta(MAXB,MAXM) real*8 x(*),newpi(*),newbeta(*) c ind=0 if (ifix(2).eq.0) then do 5 i =1,nsta 5 newpi(i) = pipars(1,i,s) else do 10 i = 1,nsta ind=ind+1 10 newpi(i) = x(ind) endif c if (ifix(3).eq.0) then do 15 i =1,nb 15 newbeta(i) = beta(i,s) else do 20 i = 1,nb ind=ind+1 20 newbeta(i) = x(ind) endif c return end subroutine newgm(m,natm,x,newgam,newmu) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' integer m,natm integer i,j,k real*8 x(*),newgam(MAXM,MAXM),newmu(MAXM,MAXM,MAXATM) c ind=0 if (ifix(1).eq.0) then do 5 i = 1,m do 5 j = 1,m 5 newgam(i,j) = gam(i,j) else do 10 i = 1,m do 10 j = 1,m ind=ind+1 10 newgam(i,j) = x(ind) endif c if (ifix(4).eq.0) then do 15 i =1,m do 15 j =1,m do 15 k =1,natm 15 newmu(i,j,k) = mu(i,j,k) else do 20 i=1,m do 20 j=1,m do 20 k=1,natm ind=ind+1 20 newmu(i,j,k) = x(ind) endif c return end options.inc000644 000423 000017 00000000161 06242170134 013730 0ustar00hughesusers000000 031510 integer atmopt,raiopt,iest,iseed real*8 neps common /options/ neps,atmopt,raiopt,iest,iseed params.inc000644 000423 000017 00000000402 06242170134 013517 0ustar00hughesusers000000 031511 real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM) real*8 pipars(2,MAXSTA,MAXM),beta(MAXB,MAXM) real*8 sigma(MAXATM,MAXATM),isigma(MAXATM,MAXATM) real*8 det integer ifix(5) common /params/ gam,mu,pipars,beta,sigma,isigma,det,ifix sumpi.f000644 000423 000017 00000024256 06242170135 013066 0ustar00hughesusers000000 031460 subroutine sumpi(l,m,nsta,nb,pipars,beta,pisum) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'hash.inc' include 'options.inc' c c Compute the probability of this rain state (i.e. this l(1) and l(2)) c Use a hashing scheme so that pisum is only computed once for any c given l(1) and l(2) - then stored for future access. Make sure c MAXPI is greater than the maximum number of rain patterns and that c MAXPI is prime. MAXPI1 = MAXPI + 1 c integer i,l(2),nmax,m,nsta,nb integer loc,nincr real*8 pipars(2,MAXSTA,MAXM),beta(MAXB,MAXM) real*8 pisum(*) data nmax/MAXPI/ c loc = mod(l(1),nmax) nincr = nmax-2 - mod(l(1),(nmax-2)) c write(*,*) l(1),l(2),loc,nincr c 5 if (hvec(1,loc).eq.l(1) .and. hvec(2,loc).eq. l(2)) then do 10 i=1,m 10 pisum(i) = pi(loc,i) c write(*,*) 'Found! pisum = ',(pisum(i),i=1,m) goto 900 elseif (hvec(1,loc).eq.-1) then ierrh = ierrh + 1 hvec(1,loc) = l(1) hvec(2,loc) = l(2) c goto (1,2,3) raiopt c Independence model 1 call indep(l,m,nsta,nb,pipars,beta,pisum) goto 100 c 2 stop 'Saturated model not implemented' c 3 call spatial(l,m,nsta,nb,pipars,beta,pisum) c 100 do 15 i=1,m 15 pi(loc,i) = pisum(i) goto 900 else if (ierrh .ge. MAXPI) stop 'Error! Increase maxpi' loc = mod(loc+nincr,nmax) endif goto 5 900 return end subroutine indep(l,m,nsta,nb,pipars,beta,pisum) Implicit double precision (A-H,O-Z) include 'nhmm.inc' c integer l(2),m,nsta,i,j,i0,bit,nb real*8 pisum(*),one real*8 pipars(2,MAXSTA,MAXM),beta(MAXB,MAXM) data one/1.0d0/ c c Independence model do 10 j=1,m pisum(j)=one 10 continue do 50 i=1,nsta i0 = i-1 if (bit(i0,l(2)).eq.0) then c if not missing if (bit(i0,l(1)).eq.1) then do 30 j=1,m pisum(j)=pisum(j)*pipars(2,i,j) 30 continue else do 40 j=1,m pisum(j)=pisum(j)*(1-pipars(2,i,j)) 40 continue endif c otherwise pisum = pisum*1.0 endif 50 continue return end subroutine spatial(l,m,nsta,nb,pipars,beta,p) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer l(2),m,nsta,i,j,ltp,ind(MAXSTA),nmis,imis integer ib1(MAXSTA),ib2(MAXSTA),nb real*8 p(*),temp(MAXM) real*8 pipars(2,MAXSTA,MAXM),beta(MAXB,MAXM) c if (l(2).eq.0) then c no missing data call autolog(l(1),m,nsta,nb,pipars,beta,p) else c sum over missing data call vbit(ib1,l(1),nsta) call vbit(ib2,l(2),nsta) call whatmiss(nsta,ib2,ind,nmis) do 105 i =1,m 105 p(i) = 0.0 do 110 i = 0,2**nmis-1 call ifill(nmis,i,ind,ib1) call icode(nsta,ib1,ltp,imis) call autolog(ltp,m,nsta,nb,pipars,beta,temp) do 155 j=1,m 155 p(j) = p(j) + temp(j) 110 continue endif return end subroutine whatmiss(nsta,ib,ind,nmis) integer nsta,nmis,ib(*),ind(*),i c nmis=0 do 100 i=1,nsta if (ib(i).eq.1) then nmis=nmis+1 ind(nmis) = i endif 100 continue return end subroutine ifill(nmis,i,ind,ib1) integer nmis,i,ind(*),ib1(*),bit c do 115 k =1,nmis 115 ib1(ind(k)) = bit(k-1,i) return end subroutine autolog(l,m,nsta,nb,pipars,beta,p) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'norm.inc' integer l,m,nsta,s,nb real*8 p(*),b(MAXSTA,MAXM) real*8 pipars(2,MAXSTA,MAXM),beta(MAXB,MAXM) integer i,ib(MAXSTA) c call vbit(ib,l,nsta) call fbeta(MAXB,beta,m,nsta,nb,ib,b) do 120 s=1,m p(s) = 0.0 do 130 i=1,nsta if (ib(i).eq.1) p(s) = p(s) + pipars(1,i,s) + b(i,s)/2.0 130 continue p(s) = exp(p(s) - lnorm(s)) 120 continue return end subroutine norman(m,nsta,nb,pipars,beta,lnorm) Implicit double precision (A-H,O-Z) include 'nhmm.inc' c c program for computing the norm for the autologistic model c integer m,nsta,nb real*8 pipars(2,MAXSTA,MAXM),beta(MAXB,MAXM),lnorm(MAXM) real*8 pips(MAXSTA),bets(MAXB) integer i,k,s c do 100 s = 1,m do 10 i = 1,nsta 10 pips(i) = pipars(1,i,s) do 20 k = 1,nb 20 bets(k) = beta(k,s) call norms(nsta,nb,s,pips,bets,lnorm(s)) 100 continue return end subroutine norms(nsta,nb,s,pipars,beta,lnorm) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'options.inc' include 'nsave.inc' c c program for computing the norm for the autologistic model c integer nsta,nb,s integer mode real*8 pipars(*),beta(*),lnorm real*8 expect(MAXSTA),expectb(MAXB) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) c mode = 2 if (iest.eq.1) then c Compute norm analytically (non-independent case) by summing over states call deriv1(mode,nsta,nb,s,pipars,beta,lnorm, * expect,expectb,dp2da2,dp2dab,dp2db2) c Compute norm by Monte Carlo elseif (iest.ge.2) then call deriv2(mode,nsta,nb,s,ista(1,s),iseed,neps, * pipars,beta,lnorm,expect,expectb, * dp2da2,dp2dab,dp2db2,pipar0,beta0,norm0) call copy0(nsta,nb,s,pipars,beta,lnorm,pipar0,beta0,norm0) endif return end subroutine stepsiz(nsta,nb,s,pipars,beta,pipar0,beta0,step,sum0) include 'nhmm.inc' integer nsta,nb,s real*8 pipars(*),beta(*),step(*) real*8 pipar0(2,MAXSTA,MAXM),beta0(MAXB,MAXM),sum0 real*8 zero data zero/0.0d0/ c ind = 0 sum0 = zero do 320 i=1,nsta ind = ind+1 step(ind) = pipars(i) - pipar0(1,i,s) 320 sum0 = sum0 + abs(step(ind)) do 330 i=1,nb ind = ind+1 step(ind) = beta(i) - beta0(i,s) 330 sum0 = sum0 + abs(step(ind)) return end subroutine copy0(nsta,nb,s,pipars,beta,lnorm,pipar0,beta0,norm0) Implicit double precision (A-H,O-Z) include 'nhmm.inc' integer nsta,nb,s,i real*8 pipars(*),beta(*),lnorm real*8 pipar0(2,MAXSTA,MAXM),beta0(MAXB,MAXM),norm0(MAXM) c do 350 i = 1,nsta pipar0(1,i,s)=pipars(i) 350 pipar0(2,i,s)=1.0/(1.0 + exp(-pipars(i))) do 360 i = 1,nb 360 beta0(i,s) = beta(i) norm0(s) = lnorm return end subroutine delnrm(mode,nsta,nb,n,iseed,s,pipar0,beta0,step,i0, * expect,expectb,dp2da2,dp2dab,dp2db2, * ista,dnorm,err) Implicit double precision (A-H,O-Z) include 'nhmm.inc' c c program for computing the change in the norm for the autologistic model c integer mode,nsta,nb,n,iseed,s,ista(*),i0 integer ib(MAXSTA),lseq(MAXSEQ) real*8 pipar0(2,MAXSTA,MAXM),beta0(MAXB,MAXM),step(*),dnorm,err real*8 expect(*),expectb(*) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 pipars(2,MAXSTA),beta(MAXB),beta1(MAXB) real*8 var,dpi(MAXSTA),db(MAXSTA),p,zero data zero/0.D0/ c dnorm = zero if (mode.eq.0) call begder(nsta,nb,expect,expectb, * dp2da2,dp2dab,dp2db2) var = zero do 305 i = 1,nsta pipars(1,i) = pipar0(1,i,s) + i0*step(i) pipars(2,i) = 1.0/(1.0+exp(-pipars(1,i))) 305 dpi(i) = step(i) do 310 i = 1,nb beta(i) = beta0(i,s) + i0*step(nsta+i) 310 beta1(i) = beta(i) + step(nsta+i) call seqgen(iseed,pipars,beta,s,nsta,nb,ista,n,lseq) do 320 ijk = 1,n call vbit(ib,lseq(ijk),nsta) call dbeta(beta1,beta,s,nsta,nb,ib,db) p = zero do 335 i=1,nsta if (ib(i).eq.1) p = p + dpi(i) + db(i) 335 continue p = exp(p) dnorm = dnorm + p var = var + p*p if (mode.eq.0) call moments(nsta,nb,s,step(nsta+1), * p,ib,expect,expectb, * dp2da2,dp2dab,dp2db2) 320 continue if (mode.eq.0) call endder(nsta,nb,dnorm,expect,expectb, * dp2da2,dp2dab,dp2db2) dnorm = dnorm/n var = sqrt(((var-n*dnorm**2)/n)/n) err = var/dnorm c write(*,*) ' N del.norm se(del.norm) rel. se.' c write(*,*) n,dnorm,var,err return end subroutine seqgen(iseed,pipars,beta,s,nsta,nb,ista,n,lseq) Implicit double precision (A-H,O-Z) include 'nhmm.inc' c c Generate a random sequence of numbers from the distribution given c by pipars and beta c integer n, s, nsta, lseq(*), iseed, ista(*) integer i,k,imis real*8 ran1, rnum, gbetai real*8 pipars(2,MAXSTA),beta(MAXB) logical lbeta data imis/0/ c c if (lbeta(beta)) then c Independence - simple bernoulli sampling do 20 i = 1,n do 10 k = 1, nsta rnum = ran1(iseed) ista(k) = 0 if (rnum.le.pipars(2,k)) ista(k)=1 10 continue call icode(nsta,ista,lseq(i),imis) 20 continue else c Autologistic - Gibb's sampler; first value is in ista do 120 i = 1,n do 110 k = 1, nsta rnum = ran1(iseed) ista(k)=1 xp = exp(pipars(1,k) + gbetai(beta,s,k,nsta,nb,ista)) if (rnum.gt.xp/(1.0 + xp)) ista(k)=0 110 continue call icode(nsta,ista,lseq(i),imis) 120 continue endif return end real*8 function ran1(idum) dimension r(97) parameter (m1=259200,ia1=7141,ic1=54773,rm1=1./m1) parameter (m2=134456,ia2=8121,ic2=28411,rm2=1./m2) parameter (m3=243000,ia3=4561,ic3=51349) data iff/0/ c if (idum.lt.0 .or. iff.eq.0) then iff=1 ix1=mod(ic1-idum,m1) ix1=mod(ia1*ix1+ic1,m1) ix2=mod(ix1,m2) ix1=mod(ia1*ix1+ic1,m1) ix3=mod(ix1,m3) do 11 j=1,97 ix1=mod(ia1*ix1+ic1,m1) ix2=mod(ia2*ix2+ic2,m2) r(j)=(float(ix1)+float(ix2)*rm2)*rm1 11 continue idum=1 endif ix1=mod(ia1*ix1+ic1,m1) ix2=mod(ia2*ix2+ic2,m2) ix3=mod(ia3*ix3+ic3,m3) j=1+(97*ix3)/m3 if (j.gt.97 .or. j.lt.1) pause ran1=r(j) r(j)=(float(ix1)+float(ix2)*rm2)*rm1 return end variance.f000644 000423 000017 00000073272 06467643122 013532 0ustar00hughesusers000000 031520 c variance4.f 10/20/95 Subroutine variance(lsets,lsetr) c c From Louis (1981) we have c c Obs. data info = E(complete data info) - E(complete data score^2) c c Redundent parameters are then projected out of the Obs. data info c (see subroutine project) and the result is inverted to give c parameter variances c Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' include 'npsetup.inc' c MAXWK = MAXTOT*(MAXTOT-MAXLIN) PARAMETER(MAXWK=171072) logical lsets,lsetr integer ns,nr,npar,ncon,n,ierr,job,lag integer iout(MAXTOT) real*8 temp(MAXTOT,MAXTOT), info(MAXTOT,MAXTOT) real*8 xs(MAXNS),xr(MAXNR),work(MAXWK) real*8 a(MAXLIN,MAXTOT),d(2) real*8 rlik,loglik real*8 eps data eps/.1/,job/11/ write(*,*) 'Enter maximum lag for computation of observed info' read(*,*) lag call gptox(ifix,m,nsta,natm,nb,gam,mu,pipars,beta,xs,ns,xr,nr) if (lsets) call npsets(nclins,ncnlns,leniws,lenws,as,bls,bus, * lsets) if (lsetr) call npsetr(nclinr,ncnlnr,leniwr,lenwr,ar,blr,bur, * lsetr) c c WRITE(*,*) ns,nr c Assume all ifix(1) - ifix(4) = 1 npar = ns+nr c rlik = loglik() call weights(rlik) write(*,*) 'Expected info' call einfo(temp) c open(12,file="variance.out") write(*,*) 'Observed info' call oinfo(temp,eps,lag) do 20 i = 1,npar 20 write(12,'(10(e14.8,x))') (temp(i,j),j=1,npar) close(12) c c find constrained observed information matrix c ncon = nclins do 25 i = 1,ncon do 23 j = 1,ns 23 a(i,j) = as(i,j) do 24 j = ns+1, ns+nr 24 a(i,j) = 0.0 25 continue do 30 i = 1,m 30 iout(i) = i*m do 35 i = 1,m do 35 j = 1,natm 35 iout(m+(i-1)*natm+j) = m*m + i*m*natm - natm + j write(*,*) 'Project' write(*,*) npar,ncon write(*,*) (iout(i),i=1,ncon) write(*,*) (a(1,i),i=1,npar) call project(npar,ncon,MAXLIN,a,iout,MAXTOT,temp,info,work) c open(13,file="variance2.out") do 19 i = 1,npar-ncon 19 write(13,'(10(e14.8,x))') (info(i,j),j=1,npar-ncon) close(13) c c invert to get variance and log determinant c n = npar - ncon write(*,*) "Inverting..." call DGEFA(info,MAXTOT,n,iout,ierr) call DGEDI(info,MAXTOT,n,iout,d,work,job) write(*,*) "log determinant = ",log(d(1)) + d(2)*log(10.0) write(*,*) "standard errors" write(*,201) (sqrt(info(i,i)),i=1,n) 201 format(8e10.3) return end subroutine einfo(info) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'options.inc' include 'dims.inc' include 'data.inc' include 'weights.inc' include 'nsave.inc' integer i,j,t,s,i1,i2,year,day integer ib1(MAXSTA),ib2(MAXSTA) real*8 info(MAXTOT,*) real*8 sum(MAXM),eQ(MAXM,MAXM) real*8 denom(MAXM) real*8 hss(MAXM) real*8 Sgm(MAXM,MAXNS) real*8 delta(MAXM),a(MAXM,MAXM) real*8 newpi(MAXSTA),newbeta(MAXB) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 exp(MAXSTA),expb(MAXB),lnorm real*8 second(MAXIR,MAXIR,MAXM) mode = 0 ns = m*m*(natm+1) nr = m*(nsta+nb) npar = ns+nr do 2 i = 1,npar do 2 j = 1,npar info(i,j) = 0.0 2 continue if (raiopt.eq.3) then do 10 s = 1,m do 12 i = 1,nsta 12 newpi(i) = pipars(1,i,s) do 14 i = 1,nb 14 newbeta(i) = beta(i,s) if (iest.eq.1) then call deriv1(mode,nsta,nb,s,newpi,newbeta,lnorm,exp,expb, * dp2da2,dp2dab,dp2db2) elseif (iest.eq.2) then call deriv2(mode,nsta,nb,s,ista(1,s),iseed,neps, * newpi,newbeta,lnorm,exp,expb, * dp2da2,dp2dab,dp2db2,pipar0,beta0,norm0) endif do 16 i=1,nsta do 15 j = 1,nsta 15 second(i,j,s) = dp2da2(i,j) do 16 j = 1,nb second(i,nsta+j,s) = dp2dab(j,i) 16 second(nsta+j,i,s) = dp2dab(j,i) do 18 i = 1,nb do 18 j = 1,nb 18 second(nsta+i,nsta+j,s) = dp2db2(i,j) 10 continue endif c t = 0 do 1000 year = 1,nyears write(*,*) year do 1000 day = 1,nobs t=t+1 c gamma and mu if (day.eq.1) then call dlset1(t,gam,mu,isigma,det,sum,eQ,denom,delta,a) c a is really the inverse of the augmented version of a call dldgm1(t,ifix,gam,mu,isigma,det,sum,eQ,denom, * delta,a,Sgm,nchk) do 310 i1 = 1,ns ii = istate(i1) do 310 i2 = i1,ns jj = istate(i2) call d2ldgm1(t,ii,jj,i1,i2,gam,mu,isigma,det,sum,eQ,denom, * a,delta,Sgm,hss) do 317 s = 1,m info(i1,i2) = info(i1,i2) - c(s,t)* * (hss(s)/delta(s)-Sgm(s,i1)*Sgm(s,i2)) info(i2,i1) = info(i1,i2) 317 continue 310 continue else call dlset2(t,gam,mu,isigma,det,sum,eQ,denom,a) call dldgm2(t,ifix,gam,mu,isigma,det,sum,eQ,denom, * a,Sgm,nchk) do 320 i1=1,ns ii = istate(i1) do 320 i2 = i1,ns if (istate(i2).ne.ii) goto 320 * call d2ldgm2(t,ii,i1,i2,gam,mu,isigma,det,sum,eQ,denom,hss) do 318 s = 1,m 318 info(i1,i2) = info(i1,i2) - c2(ii,s,t)* * (hss(s)/a(ii,s)-Sgm(s,i1)*Sgm(s,i2)) info(i2,i1) = info(i1,i2) 320 continue endif c c pipars and beta call vbit(ib1,rdata(1,t),nsta) if (raiopt.eq.1) then c independence model i1 = ns do 140 s = 1,m do 140 i = 1,nsta i1 = i1 + 1 if (ib2(i).eq.0) then if (ib1(i).eq.0) then info(i1,i1) = info(i1,i1)+c(s,t)/(1-pipars(2,i,s))**2 else info(i1,i1) = info(i1,i1)+c(s,t)/pipars(2,i,s)**2 endif endif 140 continue elseif (raiopt.eq.3) then c autologistic model i1 = ns do 400 s = 1,m do 180 i = 1,nsta i1 = i1+1 c dpi dpi do 185 j = i,nsta i2 = i1 + j - i info(i1,i2)=info(i1,i2)-c(s,t)*second(i,j,s) info(i2,i1) = info(i1,i2) 185 continue c dpi dbeta do 190 j = 1,nb i2 = m*m + (s-1)*nsta + j info(i1,i2) = info(i1,i2)-c(s,t)*second(i,nsta+j,s) info(i2,i1) = info(i1,i2) 190 continue 180 continue c dbeta dbeta do 200 i = 1,nb i1 = i1+1 do 200 j = i,nb i2 = i1 + j - i info(i1,i2) = info(i1,i2)-c(s,t)*second(nsta+i,nsta+j,s) info(i2,i1) = info(i1,i2) 200 continue 400 continue endif 1000 continue c At this point info contains the expected complete data information return end subroutine d2ldgm1(t,ii,jj,i1,i2,gam,mu,isigma,det,sum,eQ,denom, * ia,delta,Sgm,hss) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer s0,s1,t,ii,jj,i1,i2,j,k,l,h real*8 sum(MAXM),eQ(MAXM,MAXM),Sgm(MAXM,*) real*8 hss(MAXM),temp(MAXM),ap1(MAXM),ap2(MAXM) real*8 denom(MAXM) real*8 dQ,dQk,dQh,dQdmu real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM),delta(MAXM) real*8 isigma(MAXATM,MAXATM),det,ia(MAXM,MAXM) c do 10 i = 1,m 10 hss(i) = 0.0 if (i1.le.m*m) then j = i1 - (ii-1)*m if (i2.le.m*m) then c gamma gamma l = i2 - (jj-1)*m call dAdgam(ii,j,MAXM,m,sum,denom,gam,eQ,ap1) call dAdgam(jj,l,MAXM,m,sum,denom,gam,eQ,ap2) if (ii.eq.jj) call d2Adg2(ii,j,l,MAXM,m,sum,denom,gam,eQ,hss) else c gamma mu l = (i2 - m*m - (jj-1)*m*natm - 1)/natm + 1 k = mod(i2 - m*m - (jj-1)*m*natm - 1,natm) + 1 call dAdgam(ii,j,MAXM,m,sum,denom,gam,eQ,ap1) call dAdmu(jj,l,MAXM,m,sum,denom,gam,eQ,ap2) dQ = dQdmu(t,jj,l,k,isigma,mu) if (ii.eq.jj) call d2Adgm(ii,j,l,MAXM,m,sum,denom,gam,eQ, * dQ,hss) endif else c mu mu j = (i1 - m*m - (ii-1)*m*natm - 1)/natm + 1 h = mod(i1 - m*m - (ii-1)*m*natm - 1,natm) + 1 l = (i2 - m*m - (jj-1)*m*natm - 1)/natm + 1 k = mod(i2 - m*m - (jj-1)*m*natm - 1,natm) + 1 call dAdmu(ii,j,MAXM,m,sum,denom,gam,eQ,ap1) call dAdmu(jj,l,MAXM,m,sum,denom,gam,eQ,ap2) dQ = -isigma(k,h) dQk = dQdmu(t,jj,l,k,isigma,mu) dQh = dQdmu(t,ii,j,h,isigma,mu) if (ii.eq.jj) call d2Adm2(ii,j,l,MAXM,m,sum,denom,gam,eQ,dQ, * dQk,dQh,hss) endif c do 100 s0 = 1,m-1 100 temp(s0) = hss(s0)*delta(ii)+ap1(s0)*Sgm(ii,i2)*delta(ii) * +ap2(s0)*Sgm(jj,i1)*delta(jj) do 110 s1 = 1,m do 110 s0 = 1,m-1 hss(s1) = ia(s1,s0)*temp(s0) 110 continue return end subroutine d2ldgm2(t,ii,i1,i2,gam,mu,isigma,det,sum,eQ,denom,hss) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer t,ii,i1,i2 real*8 sum(MAXM),eQ(MAXM,MAXM),hss(MAXM) real*8 denom(MAXM) real*8 dQ,dQk,dQh,dQdmu real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM) real*8 isigma(MAXATM,MAXATM),det c if (i1.le.m*m) then j = i1 - (ii-1)*m if (i2.le.m*m) then c gamma gamma l = i2 - (ii-1)*m call d2Adg2(ii,j,l,MAXM,m,sum,denom,gam,eQ,hss) else c gamma mu l = (i2 - m*m - (ii-1)*m*natm - 1)/natm + 1 k = mod(i2 - m*m - (ii-1)*m*natm - 1,natm) + 1 dQ = dQdmu(t,ii,l,k,isigma,mu) call d2Adgm(ii,j,l,MAXM,m,sum,denom,gam,eQ,dQ,hss) endif else c mu mu j = (i1 - m*m - (ii-1)*m*natm - 1)/natm + 1 h = mod(i1 - m*m - (ii-1)*m*natm - 1,natm) + 1 l = (i2 - m*m - (ii-1)*m*natm - 1)/natm + 1 k = mod(i2 - m*m - (ii-1)*m*natm - 1,natm) + 1 dQ = -isigma(k,h) dQk = dQdmu(t,ii,l,k,isigma,mu) dQh = dQdmu(t,ii,j,h,isigma,mu) call d2Adm2(ii,j,l,MAXM,m,sum,denom,gam,eQ,dQ,dQk,dQh,hss) endif return end subroutine d2Adg2(i,j,l,lda,m,sum,denom,gam,eQ,hss) Implicit double precision (A-H,O-Z) integer i,j,l,m,lda,k real*8 sum(*),denom(*),gam(lda,*),eQ(lda,*),hss(m),ddd c find derivative of A(i,*) with respect to gamma(i,j) and gamma(i,l) c return i'th row of hss only ddd = sum(i)*denom(i) if (j.eq.l) then do 20 k=1,j-1 20 hss(k) = 2*gam(i,k)*eQ(i,k)*eQ(i,j)*eQ(i,l)/ddd hss(j) = -2*eQ(i,j)*eQ(i,j)*(sum(i)-gam(i,j)*eQ(i,j))/ddd do 25 k=j+1,m 25 hss(k) = 2*gam(i,k)*eQ(i,k)*eQ(i,j)*eQ(i,l)/ddd else do 30 k = 1,m if (k.eq.j) then hss(k) = -eQ(i,k)*eQ(i,l)*(sum(i)-2*gam(i,k)*eQ(i,k))/ddd elseif (k.eq.l) then hss(k) = -eQ(i,k)*eQ(i,j)*(sum(i)-2*gam(i,k)*eQ(i,k))/ddd else hss(k) = 2*gam(i,k)*eQ(i,k)*eQ(i,j)*eQ(i,l)/ddd endif 30 continue endif return end subroutine d2Adm2(i,j,l,lda,m,sum,denom,gam,eQ,dQ,dQk,dQh,hss) Implicit double precision (A-H,O-Z) integer i,j,l,m,lda,k real*8 sum(*),denom(*),gam(lda,*),eQ(lda,*),hss(m),ddd real*8 dQ,dQk,dQh,term c find derivative of A(i,*) wrt mu(i,j,k) and mu(i,l,h) c return i'th row of hss only ddd = sum(i)*denom(i) if (j.eq.l) then do 20 k=1,j-1 term = gam(i,k)*gam(i,j)*eQ(i,k)*eQ(i,j) 20 hss(k) = term*(2*gam(i,j)*eQ(i,j)*dQh*dQk * -sum(i)*(dQ+dQh*dQk))/ddd hss(j) = gam(i,j)*eQ(i,j)*(sum(i)-gam(i,j)*eQ(i,j))* * (sum(i)-2*gam(i,j)*eQ(i,j))*dQh*dQk/ddd + * gam(i,j)*eQ(i,j)*dQ* * (sum(i)-gam(i,j)*eQ(i,j))/denom(i) do 25 k=j+1,m term = gam(i,k)*gam(i,j)*eQ(i,k)*eQ(i,j) 25 hss(k) = term*(2*gam(i,j)*eQ(i,j)*dQk*dQh * -sum(i)*(dQ+dQh*dQk))/ddd else do 30 k=1,m if (k.eq.j) then hss(k) = gam(i,k)*eQ(i,k)*gam(i,l)*eQ(i,l)*dQk*dQh* * (sum(i)-2*gam(i,k)*eQ(i,k))/ddd elseif (k.eq.l) then hss(k) = gam(i,k)*eQ(i,k)*gam(i,j)*eQ(i,j)*dQk*dQh* * (sum(i)-2*gam(i,k)*eQ(i,k))/ddd else hss(k) = 2*gam(i,k)*eQ(i,k)*gam(i,l)*eQ(i,l)*dQk*dQh* * gam(i,j)*eQ(i,j)/ddd endif 30 continue endif return end subroutine d2Adgm(i,j,l,lda,m,sum,denom,gam,eQ,dQ,hss) Implicit double precision (A-H,O-Z) integer i,j,l,m,lda,k real*8 sum(*),denom(*),gam(lda,*),eQ(lda,*),hss(m),ddd,dQ c c find derivative of A(i,*) with respect to gamma(i,j) and mu(i,l,.) c return i'th row of hss only ddd = sum(i)*denom(i) if (j.eq.l) then do 20 k=1,j-1 20 hss(k) = -gam(i,k)*eQ(i,k)*eQ(i,j)*dQ* * (sum(i)-2*gam(i,j)*eQ(i,j))/ddd hss(j) = eQ(i,j)*(sum(i)-2*gam(i,j)*eQ(i,j))*dQ* * (sum(i)-gam(i,j)*eQ(i,j))/ddd do 25 k=j+1,m 25 hss(k) = -gam(i,k)*eQ(i,k)*eQ(i,j)*dQ* * (sum(i)-2*gam(i,j)*eQ(i,j))/ddd else do 30 k = 1,m if (k.eq.j) then hss(k) = -gam(i,l)*eQ(i,k)*eQ(i,l)*dQ* * (sum(i)-2*gam(i,k)*eQ(i,k))/ddd elseif (k.eq.l) then hss(k) = -gam(i,j)*eQ(i,k)*eQ(i,j)*dQ* * (sum(i)-2*gam(i,k)*eQ(i,k))/ddd else hss(k) = 2*gam(i,k)*gam(i,l)*eQ(i,k)*eQ(i,j)*dQ* * eQ(i,l)/ddd endif 30 continue endif return end subroutine project(n,p,nrowa,a,iout,nrowr,r,rz,z) c c Computation of the projected hessian. In the case of a linearly c constrained mle problem, this is the observed information. c c The problem is to minimize -log l(x) c c subject to Ax = b c c The vector of parameters x has length n. Assume that the (p) c constraints are all of the simple form c c x_i + x_{i+1} + ... + x_{i+k} = b c c i.e. the rows of A are all 1's and 0's. The point is to eliminate c the redundent parameters. c c Reorder the x's so they are in the order x = (x_B, x_S) where x_B c are the p redundent parameters which you wish to eliminate. Then c reorder the columns of A in a similar manner. Now we can write c c A = B | S c pxn pxp px(n-p) c c and define c c Z = ( -B^{-1} S ) c ( I ) c c where I is the (n-p)x(n-p) identity matrix. So Z is n x (n-p). c c Then the projected hessian (= observed information) is c c Z^{T} r Z c c where r is the hessian (with rows and columns reordered in the c appropriate manner). The inverse of this is an estimate of c the covariance matrix. c c n = number of parameters c p = number of constraints c nrowa = row dimension of a c a = constraint matrix of dimension p x n c iout = vector of length n. On input the first p elements give the c indicies (from 1 to n) of the parameters to be eliminated. c nrowr = row dimension of r c r = hessian c rz = projected hessian c z = temporary storage vector of length at least n*(n-p) c integer n,nrowa,nrowr,p,iout(*),job,ipvt(1000) real*8 a(nrowa,*),r(nrowr,*),rz(nrowr,*) real*8 z(n,*),d(2),work(1000) data job/11/ c icnt = p do 10 i=1,n do 20 j=1,p if (i.eq.iout(j)) goto 10 20 continue icnt = icnt + 1 iout(icnt) = i 10 continue c write(*,*) (iout(i),i=1,n) c do 15 i=1,p c15 write(*,*) (a(i,j),j=1,n) c compute z = B^{-1} S do 50 i=1,p do 51 j = 1,p z(i,j) = a(i,iout(j)) 51 continue c write(*,*) (z(i,j),j=1,p) 50 continue call DGEFA(z,n,p,ipvt,info) call DGEDI(z,n,p,ipvt,d,work,job) do 55 i=1,p c write(*,*) (z(i,j),j=1,p) do 58 j = 1,n-p work(j) = 0 do 57 k=1,p 57 work(j) = work(j) - z(i,k)*a(k,iout(p+j)) 58 continue do 59 j = 1,n-p 59 z(i,j) = work(j) c write(*,*) (z(i,j),j=1,n-p) 55 continue do 60 i = p+1,n do 60 j = 1,n-p z(i,j) = 0.0 if (i-p.eq.j) z(i,j) = 1.0 60 continue c do 333 i=1,n c333 write(*,"(12f5.1)") (z(i,j),j=1,n-p) c compute projected hessian = Z^T r Z do 70 i = 1,n-p do 70 j = 1,n-p rz(i,j) = 0.0 do 80 k = 1,n do 80 l = 1,n c if (i.eq.1.and.j.eq.1) c * write(*,*) l,k,z(l,i),z(k,j),r(iout(l),iout(k)) 80 rz(i,j) = rz(i,j) + z(l,i)*r(iout(l),iout(k))*z(k,j) 70 continue c Done return end subroutine oinfo(info,eps,lag) c On input, info contains the expected complete data information c On output, info contains the (approx) observed information Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' include 'data.inc' include 'options.inc' include 'weights.inc' include 'nsave.inc' integer t1,t2,ngm,npb,year,day1,day2,s,lag real*8 info(MAXTOT,*) real*8 rlik,loglik,dmatmin real*8 k2(MAXM,MAXM),k3a(MAXM,MAXM,MAXM),k3b(MAXM,MAXM,MAXM) real*8 k4(MAXM,MAXM,MAXM,MAXM) real*8 Sgm1(MAXM,MAXNS),Spb1(MAXNR) real*8 Sgm2(MAXM,MAXNS),Spb2(MAXNR) real*8 Bnd real*8 delta(MAXM),a(MAXM,MAXM),b(MAXM,MAXM) real*8 sum(MAXM),eQ(MAXM,MAXM),denom(MAXM) real*8 newpi(MAXSTA),newbeta(MAXB) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 exp(MAXSTA),expb(MAXB),lnorm real*8 expect(MAXNR) rlik = loglik() if (raiopt.eq.3) then i1 = 0 do 10 s = 1,m do 12 i = 1,nsta 12 newpi(i) = pipars(1,i,s) do 14 i = 1,nb 14 newbeta(i) = beta(i,s) if (iest.eq.1) then call deriv1(mode,nsta,nb,s,newpi,newbeta,lnorm,exp,expb, * dp2da2,dp2dab,dp2db2) elseif (iest.eq.2) then call deriv2(mode,nsta,nb,s,ista(1,s),iseed,neps, * newpi,newbeta,lnorm,exp,expb, * dp2da2,dp2dab,dp2db2,pipar0,beta0,norm0) endif do 16 i = 1,nsta i1 = i1+1 16 expect(i1) = exp(i) do 18 i = 1,nb i1 = i1+1 18 expect(i1) = expb(i) 10 continue endif c t1 = 0 do 1001 year = 1,nyears write(*,*) year t1=t1+1 call dlset1(t1,gam,mu,isigma,det,sum,eQ,denom,delta,a) call dldgm1(t1,ifix,gam,mu,isigma,det,sum,eQ,denom, * delta,a,Sgm1,ngm) call dldpb(t1,expect,Spb1,npb) call SSt1(ngm,npb,c(1,t1),Sgm1,Spb1,delta,info) do 1000 day1 = 2,nobs t1=t1+1 call dlset2(t1,gam,mu,isigma,det,sum,eQ,denom,a) call dldgm2(t1,ifix,gam,mu,isigma,det,sum,eQ,denom, * a,Sgm1,ngm) call dldpb(t1,expect,Spb1,npb) call SSt2(ngm,npb,c(1,t1),c2(1,1,t1),Sgm1,Spb1,a,info) t2 = t1 c Bnd = 1.0 do 2000 day2 = day1-1,1,-1 t2 = t2 - 1 if (day2.lt.day1-lag) goto 1000 call weight3(day2,t2,t1,k2,k3a,k3b,k4,rlik,b) c Bnd = Bnd * (1-m*dmatmin(MAXM,m,m,b)) c if (Bnd.lt.eps) goto 1000 if (day2.eq.1) then call dlset1(t2,gam,mu,isigma,det,sum,eQ,denom,delta,a) call dldgm1(t2,ifix,gam,mu,isigma,det,sum,eQ,denom, * delta,a,Sgm2,ngm) else call dlset2(t2,gam,mu,isigma,det,sum,eQ,denom,a) call dldgm2(t2,ifix,gam,mu,isigma,det,sum,eQ,denom, * a,Sgm2,ngm) endif call dldpb(t2,expect,Spb2,npb) call SSe(day2,ngm,npb,k2,k3a,k3b,k4,Sgm1,Sgm2, * Spb1,Spb2,info) c write(*,*) t1,t2,info(1,1)+info(2,2)-2*info(1,2) 2000 continue 1000 continue 1001 continue return end subroutine SSt1(ngm,npb,c1,Sgm,Spb,delta,info) c Compute E(SS') - E(S)E(S') for day 1 of season Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' integer ngm,npb,i,k,l,i1,i2,ii,jj real*8 c1(*),delta(*) real*8 d1(MAXNS),d2(MAXNR) real*8 Sgm(MAXM,*),Spb(*) real*8 info(MAXTOT,*) do 100 i1 = 1,ngm d1(i1) = 0.0 do 105 l = 1,m 105 d1(i1) = d1(i1) + c1(l)*Sgm(l,i1) do 100 i2 = 1,i1 do 110 l = 1,m 110 info(i1,i2) = info(i1,i2) - c1(l)*Sgm(l,i1)*Sgm(l,i2) info(i1,i2) = info(i1,i2) + d1(i1)*d1(i2) info(i2,i1) = info(i1,i2) 100 continue c write(*,*) (d1(i),i=1,ngm) c do 200 i1 = 1, npb ii = ngm + i1 i = istate2(i1) d2(i1) = c1(i)*Spb(i1) do 200 i2 = 1,i1 jj = ngm + i2 k = istate2(i2) if (i.eq.k) info(ii,jj)=info(ii,jj)-c1(i)*Spb(i1)*Spb(i2) info(ii,jj) = info(ii,jj) + d2(i1)*d2(i2) info(jj,ii) = info(ii,jj) 200 continue c write(*,*) (d2(i),i=1,npb) c do 300 i1 = 1, ngm do 300 i2 = 1, npb jj = ngm + i2 k = istate2(i2) info(i1,jj) = info(i1,jj) - c1(k)*Sgm(k,i1)*Spb(i2) info(i1,jj) = info(i1,jj) + d1(i1)*d2(i2) info(jj,i1) = info(i1,jj) 300 continue return end subroutine SSt2(ngm,npb,c1,c2,Sgm,Spb,a,info) c Compute E(SS') - E(S)E(S') for day 2 - end of season Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' integer ngm,npb,i,k,l,i1,i2,ii,jj,istate real*8 c1(*),c2(MAXM,*),a(MAXM,*) real*8 d1(MAXNS),d2(MAXNR) real*8 Sgm(MAXM,*),Spb(*) real*8 info(MAXTOT,*) do 100 i1 = 1,ngm i = istate(i1) d1(i1) = 0.0 do 105 l = 1,m 105 d1(i1) = d1(i1) + c2(i,l)*Sgm(l,i1) do 100 i2 = 1,i1 k = istate(i2) if (i.eq.k) then do 110 l = 1,m 110 info(i1,i2) = info(i1,i2) - c2(i,l)*Sgm(l,i1)*Sgm(l,i2) endif info(i1,i2) = info(i1,i2) + d1(i1)*d1(i2) info(i2,i1) = info(i1,i2) 100 continue c write(*,*) (d1(i),i=1,ngm) c do 200 i1 = 1, npb ii = ngm + i1 i = istate2(i1) d2(i1) = c1(i)*Spb(i1) do 200 i2 = 1,i1 jj = ngm + i2 k = istate2(i2) if (i.eq.k) info(ii,jj)=info(ii,jj)-c1(i)*Spb(i1)*Spb(i2) info(ii,jj) = info(ii,jj) + d2(i1)*d2(i2) info(jj,ii) = info(ii,jj) 200 continue c write(*,*) (d2(i),i=1,npb) c do 300 i1 = 1, ngm i = istate(i1) do 300 i2 = 1, npb jj = ngm +i2 k = istate2(i2) info(i1,jj) = info(i1,jj) - c2(i,k)*Sgm(k,i1)*Spb(i2) info(i1,jj) = info(i1,jj) + d1(i1)*d2(i2) info(jj,i1) = info(i1,jj) 300 continue return end real*8 function dmatmin(lda,m,n,a) integer lda,m,n real*8 a(lda,*),small c small = a(1,1) do 10 i=1,m do 10 j=1,n if (a(i,j).lt.small) small = a(i,j) 10 continue dmatmin = small return end subroutine weight3(day,t2,t1,k1,k2a,k2b,k3,rlik,bmat) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer t1,t2,day integer s0,s1,s2,s3 real*8 rlik,k1(MAXM,*),k2a(MAXM,MAXM,*),k2b(MAXM,MAXM,*) real*8 k3(MAXM,MAXM,MAXM,*) real*8 bmat(MAXM,*) real*8 pw1(MAXM),pw2(MAXM),pww1(MAXM,MAXM),pww2(MAXM,MAXM) real*8 a(MAXM,MAXM),a1(MAXM,MAXM) real*8 temp1(MAXM) real*8 pi(MAXM) real*8 scale,scale2,scale3,scale4 c c k1 = P(W_t2, W_t1 | Y) - P(W_t1 | Y)P(W_t2 | Y) c k2a = P(W_{t2-1}, W_t2, W_t1 | Y) - P(W_{t2-1}, W_t2 | Y)P(W_t1 | Y) c k2b = P(W_t2, W_{t1-1}, W_t1 | Y) - P(W_{t1-1}, W_t1 | Y)P(W_t2 | Y) c k3 = P(W_{t2-1},W_t2,W_{t1-1},W_t1 | Y) - c P(W_{t2-1},W_t2 | Y)P(W_{t1-1},W_t1 | Y) do 10 s0 = 1,m do 10 s1 = 1,m k1(s0,s1) = 0.0 do 10 s2 = 1,m k2a(s2,s1,s0) = 0.0 k2b(s2,s1,s0) = 0.0 do 10 s3 = 1,m k3(s3,s2,s1,s0) = 0.0 10 continue c scale = exp(forw(m+1,t2)+back(m+1,t2)-rlik) call makeA(m,gam,natm,mu,isigma,det,adata(1,t2+1),a1) do 200 s0=1,m pw2(s0) = forw(s0,t2)*back(s0,t2)*scale do 200 s1=1,m k1(s0,s1) = forw(s0,t2)*a1(s0,s1) C the following only gets used if t1 = t2+1 k2b(s0,s0,s1) = k1(s0,s1) 200 continue call sumpi(rdata(1,t2),m,nsta,nb,pipars,beta,pi) call makeA(m,gam,natm,mu,isigma,det,adata(1,t2),a) if (day.gt.1) then scale2 = exp(forw(m+1,t2-1)+back(m+1,t2)-rlik) do 220 s1=1,m do 220 s0=1,m pww2(s1,s0)=forw(s1,t2-1)*a(s1,s0)*pi(s0)*back(s0,t2)*scale2 bmat(s1,s0) = a(s1,s0)*pi(s0)*back(s0,t2)/back(s0,t2-1) do 220 s2 = 1,m k2a(s1,s0,s2) = forw(s1,t2-1)*a(s1,s0)*pi(s0)*a1(s0,s2) C the following only gets used if t1 = t2+1 k3(s1,s0,s0,s2) = k2a(s1,s0,s2) 220 continue endif do 300 t = t2+1,t1-1 call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pi) call makeA(m,gam,natm,mu,isigma,det,adata(1,t+1),a1) do 320 s0=1,m do 310 s1=1,m temp1(s1) = 0 do 310 k = 1,m 310 temp1(s1) = temp1(s1) + k1(s0,k)*pi(k)*a1(k,s1) do 320 s1=1,m if (t .eq. t1-1) then do 325 k = 1,m k2b(s0,s1,k) = k1(s0,s1)*pi(s1)*a1(s1,k) 325 continue endif k1(s0,s1) = temp1(s1) 320 continue do 340 s1=1,m do 340 s0=1,m do 350 s2=1,m temp1(s2) = 0 do 350 k = 1,m temp1(s2) = temp1(s2) + k2a(s1,s0,k)*pi(k)*a1(k,s2) 350 continue do 340 s2=1,m if (t .eq. t1-1) then do 345 k = 1,m k3(s1,s0,s2,k) = k2a(s1,s0,s2)*pi(s2)*a1(s2,k) 345 continue endif k2a(s1,s0,s2) = temp1(s2) 340 continue 300 continue scale = exp(forw(m+1,t1)+back(m+1,t1)-rlik) scale2 = exp(forw(m+1,t1-1)+back(m+1,t1)-rlik) scale3 = exp(forw(m+1,t2)+back(m+1,t1)-rlik) if (day.gt.1) scale4 = exp(forw(m+1,t2-1)+back(m+1,t1)-rlik) call sumpi(rdata(1,t1),m,nsta,nb,pipars,beta,pi) call makeA(m,gam,natm,mu,isigma,det,adata(1,t1),a) do 400 s0=1,m pw1(s0) = forw(s0,t1)*back(s0,t1)*scale do 400 s1=1,m k1(s0,s1) = k1(s0,s1)*pi(s1)*back(s1,t1)*scale3 400 continue do 420 s1=1,m do 420 s0=1,m pww1(s1,s0)=forw(s1,t1-1)*a(s1,s0)*pi(s0)*back(s0,t1)*scale2 do 420 s2 = 1,m k2b(s1,s0,s2) = k2b(s1,s0,s2)*pi(s2)*back(s2,t1)*scale3 if (day.eq.1) goto 420 k2a(s1,s0,s2) = k2a(s1,s0,s2)*pi(s2)*back(s2,t1)*scale4 do 410 s3 = 1,m 410 k3(s1,s0,s2,s3)=k3(s1,s0,s2,s3)*pi(s3)*back(s3,t1)*scale4 420 continue c do 500 s0 = 1,m do 500 s1 = 1,m k1(s0,s1) = k1(s0,s1) - pw2(s0)*pw1(s1) do 500 s2 = 1,m k2b(s0,s1,s2) = k2b(s0,s1,s2) - pw2(s0)*pww1(s1,s2) if (day.eq.1) goto 500 k2a(s0,s1,s2) = k2a(s0,s1,s2) - pww2(s0,s1)*pw1(s2) do 510 s3 = 1,m 510 k3(s0,s1,s2,s3) = k3(s0,s1,s2,s3) - pww2(s0,s1)*pww1(s2,s3) 500 continue return end subroutine SSe(day,ngm,npb,k1,k2a,k2b,k3,Sgm1,Sgm2, * Spb1,Spb2,info) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' integer ngm,npb,istate,istate2,s0,s1,day real*8 k1(MAXM,*),k2a(MAXM,MAXM,*),k2b(MAXM,MAXM,*) real*8 k3(MAXM,MAXM,MAXM,*) real*8 Sgm1(MAXM,*),Sgm2(MAXM,*),Spb1(*),Spb2(*) real*8 info(MAXTOT,*) do 100 i1 = 1,ngm i = istate(i1) do 100 i2 = 1,i1 k = istate(i2) do 110 s0 = 1,m do 110 s1 = 1,m if (day.eq.1) then info(i1,i2) = info(i1,i2) - * k2b(s1,i,s0)*Sgm1(s0,i1)*Sgm2(s1,i2)- * k2b(s1,k,s0)*Sgm1(s0,i2)*Sgm2(s1,i1) else info(i1,i2) = info(i1,i2) - * k3(k,s1,i,s0)*Sgm1(s0,i1)*Sgm2(s1,i2)- * k3(i,s1,k,s0)*Sgm1(s0,i2)*Sgm2(s1,i1) endif 110 continue info(i2,i1) = info(i1,i2) 100 continue c do 200 i1 = 1, npb ii = ngm + i1 i = istate2(i1) do 200 i2 = 1,i1 jj = ngm + i2 k = istate2(i2) c write(*,*) day,ii,jj,i1,i2,i,k,info(ii,jj) info(ii,jj)=info(ii,jj)-k1(i,k)*Spb2(i1)*Spb1(i2) * -k1(k,i)*Spb2(i2)*Spb1(i1) info(jj,ii) = info(ii,jj) 200 continue c do 300 i1 = 1, ngm i = istate(i1) do 300 i2 = 1, npb jj = ngm +i2 k = istate2(i2) do 310 s0 = 1,m if (day.eq.1) then info(i1,jj)=info(i1,jj)-k1(s0,k)*Sgm2(s0,i1)*Spb1(i2) * -k2b(k,i,s0)*Sgm1(s0,i1)*Spb2(i2) else info(i1,jj)=info(i1,jj)-k2a(i,s0,k)*Sgm2(s0,i1)*Spb1(i2) * -k2b(k,i,s0)*Sgm1(s0,i1)*Spb2(i2) endif 310 continue info(jj,i1) = info(i1,jj) 300 continue return end viterbi.f000644 000423 000017 00000011666 06470140122 013373 0ustar00hughesusers000000 031470 subroutine viterbi() Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'data.inc' include 'dims.inc' include 'params.inc' include 'fb.inc' c integer year,day,s integer Q(MAXM,MAXOBS) integer argmax,jmax,nnn(MAXM,MAXM) real*8 f(MAXM),g(MAXM) real*8 a(MAXM,MAXM),pisum(MAXM),arow(MAXM) real*8 xbar(MAXM,MAXM,MAXATM),xss(MAXM,MAXM,MAXATM,MAXATM) real*8 sighat(MAXATM,MAXATM),gamsum real*8 arg,temp(MAXATM),newmu(MAXM,MAXM,MAXATM),newgam(MAXM,MAXM) character*80 outfil write(*,*) 'Enter name of output file for state sequence.' read(*,'(a80)') outfil open(4,file=outfil) do 40 i=1,m do 40 j=1,m nnn(i,j)=0 do 40 k=1,natm xbar(i,j,k)=0.0 do 40 l=1,natm 40 xss(i,j,k,l)=0.0 c do 60 year=1,nyears t = year*nobs call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pisum) call copyx(pisum,g,m) call scale(g,m) c now do it for t = dend-1 ... dbeg do 50 day=nobs-1,1,-1 t = t-1 call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pisum) call makeA(m,gam,natm,mu,isigma,det,adata(1,t+1),a) c do 55 s=1,m do 53 i = 1,m 53 arow(i) = a(s,i) jmax = argmax(arow,g,m) Q(s,day+1) = jmax f(s) = pisum(s)*arow(jmax)*g(jmax) c write(*,*) s,jmax,f(s) 55 continue c call scale(f,m) call copyx(f,g,m) 50 continue c and finish up call steady(m,gam,natm,mu,isigma,det,adata(1,t),delta) jmax = argmax(delta,f,m) Q(jmax,1) = jmax write(4,'(365i2)') (Q(jmax,day),day=1,nobs) c c accumulate sums of squares for pooled estimate of sigma c t = (year-1)*nobs + 1 do 62 day=2,nobs t = t+1 i = Q(jmax,day-1) j = Q(jmax,day) nnn(i,j) = nnn(i,j)+1 do 62 k=1,natm xbar(i,j,k) = xbar(i,j,k) + adata(k,t) do 62 l=1,natm 62 xss(i,j,k,l) = xss(i,j,k,l) + adata(k,t)*adata(l,t) 60 continue c c compute pooled estimate of sigma c do 70 i=1,m do 70 j=1,m do 70 k=1,natm 70 if (nnn(i,j).gt.0) xbar(i,j,k) = xbar(i,j,k)/nnn(i,j) df = nyears*nobs-m*m do 71 k=1,natm do 71 l=1,natm sighat(k,l)=0.0 do 72 i=1,m do 72 j=1,m c write(*,*) k,l,i,j,xss(i,j,k,l),xbar(i,j,k),xbar(i,j,l) 72 sighat(k,l) = sighat(k,l) + * (xss(i,j,k,l) - nnn(i,j)*xbar(i,j,k)*xbar(i,j,l)) 71 sighat(k,l) = sighat(k,l)/df c compute rescaled estimates of the parameters c do 370 i=1,natm c370 write(*,*) (isigma(i,j),j=1,natm) do 80 i=1,m gamsum=0.0 do 85 j=1,m do 90 k=1,natm temp(k) = 0.0 do 93 l=1,natm 93 temp(k) = temp(k) + mu(i,j,l)*isigma(l,k) 90 continue c write(*,*) (temp(l),l=1,natm) arg = 0.0 do 95 k=1,natm newmu(i,j,k) = 0.0 do 96 l=1,natm 96 newmu(i,j,k) = newmu(i,j,k) + temp(l)*sighat(l,k) arg = arg + temp(k)*(mu(i,j,k) - newmu(i,j,k)) 95 continue newgam(i,j) = gam(i,j) * dexp(-.5*arg) gamsum = gamsum + newgam(i,j) 85 continue c write(*,*) (newgam(i,j),j=1,m) do 87 j=1,m 87 newgam(i,j)=newgam(i,j)/gamsum 80 continue c print parameters - Bayes form write(6,202) 202 format(//' Rescaled parameters'//' Gamma') do 160 i=1,m write(6,6040)(newgam(i,j),j=1,m) 6040 format(6g12.4) 160 continue write(6,204) 204 format(/' Pipars') do 170 i=1,nsta write(6,6050)(pipars(2,i,j),j=1,m) 6050 format(6g12.4) 170 continue write(6,205) 205 format(/' Beta') do 180 i=1,nb write(6,6050)(beta(i,j),j=1,m) 180 continue write(6,206) 206 format(/' Mu'/' row col of gamma') do 190 i=1,m do 190 j=1,m write(6,6070) i,j,(newmu(i,j,k),k=1,natm) 6070 format(2i4,5g12.4) 190 continue write(6,207) 207 format(/' Sigma') do 200 i=1,natm write(6,6080) (sighat(i,j),j=1,natm) 6080 format(5g12.4) 200 continue c stop end subroutine copyx(x1,x2,n) Implicit double precision (A-H,O-Z) real*8 x1(*),x2(*) do 10 i=1,n 10 x2(i)=x1(i) return end integer function argmax(a,f,m) Implicit double precision (A-H,O-Z) real*8 a(*),f(*),rmax,rcand c c write(*,*) 'argmax' c write(*,*) (a(i),i=1,m) c write(*,*) (f(i),i=1,m) argmax = 1 rmax = a(1)*f(1) do 10 i = 2,m rcand = a(i)*f(i) if (rcand.gt.rmax) then rmax = rcand argmax=i endif 10 continue c write(*,*) argmax return end subroutine scale(f,m) integer m,s real*8 ts,f(*) c ts = 0.0 do 57 s=1,m 57 ts = ts+f(s) ts=ts/dble(m) do 58 s=1,m 58 f(s) = f(s)/ts return end weights.inc000644 000423 000017 00000000115 06242170135 013713 0ustar00hughesusers000000 031406 real*8 c(MAXM,MAXOBS),c2(MAXM,MAXM,MAXOBS) common /weight/ c,c2 variance.f000644 000423 000017 00000073360 06704463370 013526 0ustar00hughesusers000000 101170 c variance4.f 10/20/95 Subroutine variance(lsets,lsetr) c c From Louis (1981) we have c c Obs. data info = E(complete data info) - E(complete data score^2) c c Redundent parameters are then projected out of the Obs. data info c (see subroutine project) and the result is inverted to give c parameter variances c Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' include 'npsetup.inc' c MAXWK = MAXTOT*(MAXTOT-MAXLIN) PARAMETER(MAXWK=171072) logical lsets,lsetr integer ns,nr,npar,ncon,n,ierr,job,lag integer iout(MAXTOT) real*8 temp(MAXTOT,MAXTOT), info(MAXTOT,MAXTOT) real*8 xs(MAXNS),xr(MAXNR),work(MAXWK) real*8 a(MAXLIN,MAXTOT),d(2) real*8 rlik,loglik real*8 eps data eps/.1/,job/11/ write(*,*) 'Enter maximum lag for computation of observed info' read(*,*) lag call gptox(ifix,m,nsta,natm,nb,gam,mu,pipars,beta,xs,ns,xr,nr) if (lsets) call npsets(nclins,ncnlns,leniws,lenws,as,bls,bus, * lsets) if (lsetr) call npsetr(nclinr,ncnlnr,leniwr,lenwr,ar,blr,bur, * lsetr) c c WRITE(*,*) ns,nr c Assume all ifix(1) - ifix(4) = 1 npar = ns+nr c rlik = loglik() call weights(rlik) write(*,*) 'Expected info' call einfo(temp) open(12,file="variance.out") write(*,*) 'Observed info' call oinfo(temp,eps,lag) do 20 i = 1,npar 20 write(12,'(10(e14.8,x))') (temp(i,j),j=1,npar) close(12) c c find constrained observed information matrix c ncon = nclins do 25 i = 1,ncon do 23 j = 1,ns 23 a(i,j) = as(i,j) do 24 j = ns+1, ns+nr 24 a(i,j) = 0.0 25 continue do 30 i = 1,m 30 iout(i) = i*m do 35 i = 1,m do 35 j = 1,natm 35 iout(m+(i-1)*natm+j) = m*m + i*m*natm - natm + j write(*,*) 'Project' write(*,*) npar,ncon write(*,*) (iout(i),i=1,ncon) write(*,*) (a(1,i),i=1,npar) call project(npar,ncon,MAXLIN,a,iout,MAXTOT,temp,info,work) open(13,file="variance2.out") do 19 i = 1,npar-ncon 19 write(13,'(10(e14.8,x))') (info(i,j),j=1,npar-ncon) close(13) c c invert to get variance and log determinant c n = npar - ncon write(*,*) "Inverting..." call DGEFA(info,MAXTOT,n,iout,ierr) call DGEDI(info,MAXTOT,n,iout,d,work,job) write(*,*) "log determinant = ",log(d(1)) + d(2)*log(10.0) write(*,*) "standard errors" write(*,201) (sqrt(info(i,i)),i=1,n) 201 format(8e10.3) return end subroutine einfo(info) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'options.inc' include 'dims.inc' include 'data.inc' include 'weights.inc' include 'nsave.inc' integer i,j,t,s,i1,i2,year,day,ns,nr,npar integer ib1(MAXSTA),ib2(MAXSTA) real*8 info(MAXTOT,*) real*8 sum(MAXM),eQ(MAXM,MAXM) real*8 denom(MAXM) real*8 hss(MAXM) real*8 Sgm(MAXM,MAXNS) real*8 delta(MAXM),a(MAXM,MAXM) real*8 newpi(MAXSTA),newbeta(MAXB) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 exp(MAXSTA),expb(MAXB),lnorm real*8 second(MAXIR,MAXIR,MAXM) mode = 0 ns = m*m*(natm+1) nr = m*(nsta+nb) npar = ns+nr do 2 i = 1,npar do 2 j = 1,npar info(i,j) = 0.0 2 continue if (raiopt.eq.3) then do 10 s = 1,m do 12 i = 1,nsta 12 newpi(i) = pipars(1,i,s) do 14 i = 1,nb 14 newbeta(i) = beta(i,s) if (iest.eq.1) then call deriv1(mode,nsta,nb,s,newpi,newbeta,lnorm,exp,expb, * dp2da2,dp2dab,dp2db2) elseif (iest.eq.2) then call deriv2(mode,nsta,nb,s,ista(1,s),iseed,neps, * newpi,newbeta,lnorm,exp,expb, * dp2da2,dp2dab,dp2db2,pipar0,beta0,norm0) endif do 16 i=1,nsta do 15 j = 1,nsta 15 second(i,j,s) = dp2da2(i,j) do 16 j = 1,nb second(i,nsta+j,s) = dp2dab(j,i) 16 second(nsta+j,i,s) = dp2dab(j,i) do 18 i = 1,nb do 18 j = 1,nb 18 second(nsta+i,nsta+j,s) = dp2db2(i,j) 10 continue endif c t = 0 do 1000 year = 1,nyears write(*,*) year do 1000 day = 1,nobs t=t+1 c gamma and mu if (day.eq.1) then call dlset1(t,gam,mu,isigma,det,sum,eQ,denom,delta,a) c a is really the inverse of the augmented version of a call dldgm1(t,ifix,gam,mu,isigma,det,sum,eQ,denom, * delta,a,Sgm,nchk) do 310 i1 = 1,ns ii = istate(i1) do 310 i2 = i1,ns jj = istate(i2) call d2ldgm1(t,ii,jj,i1,i2,gam,mu,isigma,det,sum,eQ,denom, * a,delta,Sgm,hss) do 317 s = 1,m info(i1,i2) = info(i1,i2) - c(s,t)* * (hss(s)/delta(s)-Sgm(s,i1)*Sgm(s,i2)) info(i2,i1) = info(i1,i2) 317 continue 310 continue else call dlset2(t,gam,mu,isigma,det,sum,eQ,denom,a) call dldgm2(t,ifix,gam,mu,isigma,det,sum,eQ,denom, * a,Sgm,nchk) do 320 i1=1,ns ii = istate(i1) do 320 i2 = i1,ns if (istate(i2).ne.ii) goto 320 * call d2ldgm2(t,ii,i1,i2,gam,mu,isigma,det,sum,eQ,denom,hss) do 318 s = 1,m 318 info(i1,i2) = info(i1,i2) - c2(ii,s,t)* * (hss(s)/a(ii,s)-Sgm(s,i1)*Sgm(s,i2)) info(i2,i1) = info(i1,i2) 320 continue endif c c pipars and beta call vbit(ib1,rdata(1,t),nsta) if (raiopt.eq.1) then c independence model i1 = ns do 140 s = 1,m do 140 i = 1,nsta i1 = i1 + 1 if (ib2(i).eq.0) then if (ib1(i).eq.0) then info(i1,i1) = info(i1,i1)+c(s,t)/(1-pipars(2,i,s))**2 else info(i1,i1) = info(i1,i1)+c(s,t)/pipars(2,i,s)**2 endif endif 140 continue elseif (raiopt.eq.3) then c autologistic model i1 = ns do 400 s = 1,m do 180 i = 1,nsta i1 = i1+1 c dpi dpi do 185 j = i,nsta i2 = i1 + j - i info(i1,i2)=info(i1,i2)+c(s,t)*second(i,j,s) info(i2,i1) = info(i1,i2) 185 continue c dpi dbeta do 190 j = 1,nb i2 = ns + s*nsta + (s-1)*nb + j info(i1,i2) = info(i1,i2)+c(s,t)*second(i,nsta+j,s) info(i2,i1) = info(i1,i2) 190 continue 180 continue c dbeta dbeta do 200 i = 1,nb i1 = i1+1 do 200 j = i,nb i2 = i1 + j - i info(i1,i2) = info(i1,i2)+c(s,t)*second(nsta+i,nsta+j,s) info(i2,i1) = info(i1,i2) 200 continue 400 continue endif 1000 continue c At this point info contains the expected complete data information return end subroutine d2ldgm1(t,ii,jj,i1,i2,gam,mu,isigma,det,sum,eQ,denom, * ia,delta,Sgm,hss) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer s0,s1,t,ii,jj,i1,i2,j,k,l,h real*8 sum(MAXM),eQ(MAXM,MAXM),Sgm(MAXM,*) real*8 hss(MAXM),temp(MAXM),ap1(MAXM),ap2(MAXM) real*8 denom(MAXM) real*8 dQ,dQk,dQh,dQdmu real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM),delta(MAXM) real*8 isigma(MAXATM,MAXATM),det,ia(MAXM,MAXM) c do 10 i = 1,m 10 hss(i) = 0.0 if (i1.le.m*m) then j = i1 - (ii-1)*m if (i2.le.m*m) then c gamma gamma l = i2 - (jj-1)*m call dAdgam(ii,j,MAXM,m,sum,denom,gam,eQ,ap1) call dAdgam(jj,l,MAXM,m,sum,denom,gam,eQ,ap2) if (ii.eq.jj) call d2Adg2(ii,j,l,MAXM,m,sum,denom,gam,eQ,hss) else c gamma mu l = (i2 - m*m - (jj-1)*m*natm - 1)/natm + 1 k = mod(i2 - m*m - (jj-1)*m*natm - 1,natm) + 1 call dAdgam(ii,j,MAXM,m,sum,denom,gam,eQ,ap1) call dAdmu(jj,l,MAXM,m,sum,denom,gam,eQ,ap2) dQ = dQdmu(t,jj,l,k,isigma,mu) if (ii.eq.jj) call d2Adgm(ii,j,l,MAXM,m,sum,denom,gam,eQ, * dQ,hss) endif else c mu mu j = (i1 - m*m - (ii-1)*m*natm - 1)/natm + 1 h = mod(i1 - m*m - (ii-1)*m*natm - 1,natm) + 1 l = (i2 - m*m - (jj-1)*m*natm - 1)/natm + 1 k = mod(i2 - m*m - (jj-1)*m*natm - 1,natm) + 1 call dAdmu(ii,j,MAXM,m,sum,denom,gam,eQ,ap1) call dAdmu(jj,l,MAXM,m,sum,denom,gam,eQ,ap2) dQ = -isigma(k,h) dQk = dQdmu(t,jj,l,k,isigma,mu) dQh = dQdmu(t,ii,j,h,isigma,mu) if (ii.eq.jj) call d2Adm2(ii,j,l,MAXM,m,sum,denom,gam,eQ,dQ, * dQk,dQh,hss) endif c do 100 s0 = 1,m-1 100 temp(s0) = hss(s0)*delta(ii)+ap1(s0)*Sgm(ii,i2)*delta(ii) * +ap2(s0)*Sgm(jj,i1)*delta(jj) do 110 s1 = 1,m do 110 s0 = 1,m-1 hss(s1) = ia(s1,s0)*temp(s0) 110 continue return end subroutine d2ldgm2(t,ii,i1,i2,gam,mu,isigma,det,sum,eQ,denom,hss) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer t,ii,i1,i2 real*8 sum(MAXM),eQ(MAXM,MAXM),hss(MAXM) real*8 denom(MAXM) real*8 dQ,dQk,dQh,dQdmu real*8 gam(MAXM,MAXM),mu(MAXM,MAXM,MAXATM) real*8 isigma(MAXATM,MAXATM),det c if (i1.le.m*m) then j = i1 - (ii-1)*m if (i2.le.m*m) then c gamma gamma l = i2 - (ii-1)*m call d2Adg2(ii,j,l,MAXM,m,sum,denom,gam,eQ,hss) else c gamma mu l = (i2 - m*m - (ii-1)*m*natm - 1)/natm + 1 k = mod(i2 - m*m - (ii-1)*m*natm - 1,natm) + 1 dQ = dQdmu(t,ii,l,k,isigma,mu) call d2Adgm(ii,j,l,MAXM,m,sum,denom,gam,eQ,dQ,hss) endif else c mu mu j = (i1 - m*m - (ii-1)*m*natm - 1)/natm + 1 h = mod(i1 - m*m - (ii-1)*m*natm - 1,natm) + 1 l = (i2 - m*m - (ii-1)*m*natm - 1)/natm + 1 k = mod(i2 - m*m - (ii-1)*m*natm - 1,natm) + 1 dQ = -isigma(k,h) dQk = dQdmu(t,ii,l,k,isigma,mu) dQh = dQdmu(t,ii,j,h,isigma,mu) call d2Adm2(ii,j,l,MAXM,m,sum,denom,gam,eQ,dQ,dQk,dQh,hss) endif return end subroutine d2Adg2(i,j,l,lda,m,sum,denom,gam,eQ,hss) Implicit double precision (A-H,O-Z) integer i,j,l,m,lda,k real*8 sum(*),denom(*),gam(lda,*),eQ(lda,*),hss(m),ddd c find derivative of A(i,*) with respect to gamma(i,j) and gamma(i,l) c return i'th row of hss only ddd = sum(i)*denom(i) if (j.eq.l) then do 20 k=1,j-1 20 hss(k) = 2*gam(i,k)*eQ(i,k)*eQ(i,j)*eQ(i,l)/ddd hss(j) = -2*eQ(i,j)*eQ(i,j)*(sum(i)-gam(i,j)*eQ(i,j))/ddd do 25 k=j+1,m 25 hss(k) = 2*gam(i,k)*eQ(i,k)*eQ(i,j)*eQ(i,l)/ddd else do 30 k = 1,m if (k.eq.j) then hss(k) = -eQ(i,k)*eQ(i,l)*(sum(i)-2*gam(i,k)*eQ(i,k))/ddd elseif (k.eq.l) then hss(k) = -eQ(i,k)*eQ(i,j)*(sum(i)-2*gam(i,k)*eQ(i,k))/ddd else hss(k) = 2*gam(i,k)*eQ(i,k)*eQ(i,j)*eQ(i,l)/ddd endif 30 continue endif return end subroutine d2Adm2(i,j,l,lda,m,sum,denom,gam,eQ,dQ,dQk,dQh,hss) Implicit double precision (A-H,O-Z) integer i,j,l,m,lda,k real*8 sum(*),denom(*),gam(lda,*),eQ(lda,*),hss(m),ddd real*8 dQ,dQk,dQh,term c find derivative of A(i,*) wrt mu(i,j,k) and mu(i,l,h) c return i'th row of hss only ddd = sum(i)*denom(i) if (j.eq.l) then do 20 k=1,j-1 term = gam(i,k)*gam(i,j)*eQ(i,k)*eQ(i,j) 20 hss(k) = term*(2*gam(i,j)*eQ(i,j)*dQh*dQk * -sum(i)*(dQ+dQh*dQk))/ddd hss(j) = gam(i,j)*eQ(i,j)*(sum(i)-gam(i,j)*eQ(i,j))* * (sum(i)-2*gam(i,j)*eQ(i,j))*dQh*dQk/ddd + * gam(i,j)*eQ(i,j)*dQ* * (sum(i)-gam(i,j)*eQ(i,j))/denom(i) do 25 k=j+1,m term = gam(i,k)*gam(i,j)*eQ(i,k)*eQ(i,j) 25 hss(k) = term*(2*gam(i,j)*eQ(i,j)*dQk*dQh * -sum(i)*(dQ+dQh*dQk))/ddd else do 30 k=1,m if (k.eq.j) then hss(k) = gam(i,k)*eQ(i,k)*gam(i,l)*eQ(i,l)*dQk*dQh* * (sum(i)-2*gam(i,k)*eQ(i,k))/ddd elseif (k.eq.l) then hss(k) = gam(i,k)*eQ(i,k)*gam(i,j)*eQ(i,j)*dQk*dQh* * (sum(i)-2*gam(i,k)*eQ(i,k))/ddd else hss(k) = 2*gam(i,k)*eQ(i,k)*gam(i,l)*eQ(i,l)*dQk*dQh* * gam(i,j)*eQ(i,j)/ddd endif 30 continue endif return end subroutine d2Adgm(i,j,l,lda,m,sum,denom,gam,eQ,dQ,hss) Implicit double precision (A-H,O-Z) integer i,j,l,m,lda,k real*8 sum(*),denom(*),gam(lda,*),eQ(lda,*),hss(m),ddd,dQ c c find derivative of A(i,*) with respect to gamma(i,j) and mu(i,l,.) c return i'th row of hss only ddd = sum(i)*denom(i) if (j.eq.l) then do 20 k=1,j-1 20 hss(k) = -gam(i,k)*eQ(i,k)*eQ(i,j)*dQ* * (sum(i)-2*gam(i,j)*eQ(i,j))/ddd hss(j) = eQ(i,j)*(sum(i)-2*gam(i,j)*eQ(i,j))*dQ* * (sum(i)-gam(i,j)*eQ(i,j))/ddd do 25 k=j+1,m 25 hss(k) = -gam(i,k)*eQ(i,k)*eQ(i,j)*dQ* * (sum(i)-2*gam(i,j)*eQ(i,j))/ddd else do 30 k = 1,m if (k.eq.j) then hss(k) = -gam(i,l)*eQ(i,k)*eQ(i,l)*dQ* * (sum(i)-2*gam(i,k)*eQ(i,k))/ddd elseif (k.eq.l) then hss(k) = -gam(i,j)*eQ(i,k)*eQ(i,j)*dQ* * (sum(i)-2*gam(i,k)*eQ(i,k))/ddd else hss(k) = 2*gam(i,k)*gam(i,l)*eQ(i,k)*eQ(i,j)*dQ* * eQ(i,l)/ddd endif 30 continue endif return end subroutine project(n,p,nrowa,a,iout,nrowr,r,rz,z) c c Computation of the projected hessian. In the case of a linearly c constrained mle problem, this is the observed information. c c The problem is to minimize -log l(x) c c subject to Ax = b c c The vector of parameters x has length n. Assume that the (p) c constraints are all of the simple form c c x_i + x_{i+1} + ... + x_{i+k} = b c c i.e. the rows of A are all 1's and 0's. The point is to eliminate c the redundent parameters. c c Reorder the x's so they are in the order x = (x_B, x_S) where x_B c are the p redundent parameters which you wish to eliminate. Then c reorder the columns of A in a similar manner. Now we can write c c A = B | S c pxn pxp px(n-p) c c and define c c Z = ( -B^{-1} S ) c ( I ) c c where I is the (n-p)x(n-p) identity matrix. So Z is n x (n-p). c c Then the projected hessian (= observed information) is c c Z^{T} r Z c c where r is the hessian (with rows and columns reordered in the c appropriate manner). The inverse of this is an estimate of c the covariance matrix. c c n = number of parameters c p = number of constraints c nrowa = row dimension of a c a = constraint matrix of dimension p x n c iout = vector of length n. On input the first p elements give the c indicies (from 1 to n) of the parameters to be eliminated. c nrowr = row dimension of r c r = hessian c rz = projected hessian c z = temporary storage vector of length at least n*(n-p) c integer n,nrowa,nrowr,p,iout(*),job,ipvt(1000) integer icnt, i, j, k, l real*8 a(nrowa,*),r(nrowr,*),rz(nrowr,*) real*8 z(n,*),d(2),work(1000) data job/11/ c icnt = p do 10 i=1,n do 20 j=1,p if (i.eq.iout(j)) goto 10 20 continue icnt = icnt + 1 iout(icnt) = i 10 continue c write(*,*) (iout(i),i=1,n) c do 15 i=1,p c15 write(*,*) (a(i,j),j=1,n) c compute z = B^{-1} S do 50 i=1,p do 51 j = 1,p z(i,j) = a(i,iout(j)) 51 continue c write(*,*) (z(i,j),j=1,p) 50 continue call DGEFA(z,n,p,ipvt,info) call DGEDI(z,n,p,ipvt,d,work,job) do 55 i=1,p c write(*,*) (z(i,j),j=1,p) do 58 j = 1,n-p work(j) = 0 do 57 k=1,p 57 work(j) = work(j) - z(i,k)*a(k,iout(p+j)) 58 continue do 59 j = 1,n-p 59 z(i,j) = work(j) c write(*,*) (z(i,j),j=1,n-p) 55 continue do 60 i = p+1,n do 60 j = 1,n-p z(i,j) = 0.0 if (i-p.eq.j) z(i,j) = 1.0 60 continue c do 333 i=1,n c333 write(*,"(12f5.1)") (z(i,j),j=1,n-p) c compute projected hessian = Z^T r Z do 70 i = 1,n-p do 70 j = 1,n-p rz(i,j) = 0.0 do 80 k = 1,n do 80 l = 1,n c if (i.eq.1.and.j.eq.1) c * write(*,*) l,k,z(l,i),z(k,j),r(iout(l),iout(k)) 80 rz(i,j) = rz(i,j) + z(l,i)*r(iout(l),iout(k))*z(k,j) 70 continue c Done return end subroutine oinfo(info,eps,lag) c On input, info contains the expected complete data information c On output, info contains the (approx) observed information Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' include 'data.inc' include 'options.inc' include 'weights.inc' include 'nsave.inc' integer t1,t2,ngm,npb,year,day1,day2,s,lag real*8 info(MAXTOT,*),eps real*8 rlik,loglik,dmatmin real*8 k2(MAXM,MAXM),k3a(MAXM,MAXM,MAXM),k3b(MAXM,MAXM,MAXM) real*8 k4(MAXM,MAXM,MAXM,MAXM) real*8 Sgm1(MAXM,MAXNS),Spb1(MAXNR) real*8 Sgm2(MAXM,MAXNS),Spb2(MAXNR) real*8 Bnd real*8 delta(MAXM),a(MAXM,MAXM),b(MAXM,MAXM) real*8 sum(MAXM),eQ(MAXM,MAXM),denom(MAXM) real*8 newpi(MAXSTA),newbeta(MAXB) real*8 dp2da2(MAXSTA,MAXSTA),dp2dab(MAXB,MAXSTA),dp2db2(MAXB,MAXB) real*8 exp(MAXSTA),expb(MAXB),lnorm real*8 expect(MAXNR) rlik = loglik() if (raiopt.eq.3) then i1 = 0 do 10 s = 1,m do 12 i = 1,nsta 12 newpi(i) = pipars(1,i,s) do 14 i = 1,nb 14 newbeta(i) = beta(i,s) if (iest.eq.1) then call deriv1(mode,nsta,nb,s,newpi,newbeta,lnorm,exp,expb, * dp2da2,dp2dab,dp2db2) elseif (iest.eq.2) then call deriv2(mode,nsta,nb,s,ista(1,s),iseed,neps, * newpi,newbeta,lnorm,exp,expb, * dp2da2,dp2dab,dp2db2,pipar0,beta0,norm0) endif do 16 i = 1,nsta i1 = i1+1 16 expect(i1) = exp(i) do 18 i = 1,nb i1 = i1+1 18 expect(i1) = expb(i) 10 continue endif c t1 = 0 do 1001 year = 1,nyears write(*,*) year t1=t1+1 call dlset1(t1,gam,mu,isigma,det,sum,eQ,denom,delta,a) call dldgm1(t1,ifix,gam,mu,isigma,det,sum,eQ,denom, * delta,a,Sgm1,ngm) call dldpb(t1,expect,Spb1,npb) call SSt1(ngm,npb,c(1,t1),Sgm1,Spb1,delta,info) do 1000 day1 = 2,nobs t1=t1+1 call dlset2(t1,gam,mu,isigma,det,sum,eQ,denom,a) call dldgm2(t1,ifix,gam,mu,isigma,det,sum,eQ,denom, * a,Sgm1,ngm) call dldpb(t1,expect,Spb1,npb) call SSt2(ngm,npb,c(1,t1),c2(1,1,t1),Sgm1,Spb1,a,info) t2 = t1 c Bnd = 1.0 do 2000 day2 = day1-1,1,-1 t2 = t2 - 1 if (day2.lt.day1-lag) goto 1000 call weight3(day2,t2,t1,k2,k3a,k3b,k4,rlik,b) c Bnd = Bnd * (1-m*dmatmin(MAXM,m,m,b)) c if (Bnd.lt.eps) goto 1000 if (day2.eq.1) then call dlset1(t2,gam,mu,isigma,det,sum,eQ,denom,delta,a) call dldgm1(t2,ifix,gam,mu,isigma,det,sum,eQ,denom, * delta,a,Sgm2,ngm) else call dlset2(t2,gam,mu,isigma,det,sum,eQ,denom,a) call dldgm2(t2,ifix,gam,mu,isigma,det,sum,eQ,denom, * a,Sgm2,ngm) endif call dldpb(t2,expect,Spb2,npb) call SSe(day2,ngm,npb,k2,k3a,k3b,k4,Sgm1,Sgm2, * Spb1,Spb2,info) c write(*,*) t1,t2,info(1,1)+info(2,2)-2*info(1,2) 2000 continue 1000 continue 1001 continue return end subroutine SSt1(ngm,npb,c1,Sgm,Spb,delta,info) c Compute E(SS') - E(S)E(S') for day 1 of season Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' integer ngm,npb,i,k,l,i1,i2,ii,jj real*8 c1(*),delta(*) real*8 d1(MAXNS),d2(MAXNR) real*8 Sgm(MAXM,*),Spb(*) real*8 info(MAXTOT,*) do 100 i1 = 1,ngm d1(i1) = 0.0 do 105 l = 1,m 105 d1(i1) = d1(i1) + c1(l)*Sgm(l,i1) do 100 i2 = 1,i1 do 110 l = 1,m 110 info(i1,i2) = info(i1,i2) - c1(l)*Sgm(l,i1)*Sgm(l,i2) info(i1,i2) = info(i1,i2) + d1(i1)*d1(i2) info(i2,i1) = info(i1,i2) 100 continue c write(*,*) (d1(i),i=1,ngm) c do 200 i1 = 1, npb ii = ngm + i1 i = istate2(i1) d2(i1) = c1(i)*Spb(i1) do 200 i2 = 1,i1 jj = ngm + i2 k = istate2(i2) if (i.eq.k) info(ii,jj)=info(ii,jj)-c1(i)*Spb(i1)*Spb(i2) info(ii,jj) = info(ii,jj) + d2(i1)*d2(i2) info(jj,ii) = info(ii,jj) 200 continue c write(*,*) (d2(i),i=1,npb) c do 300 i1 = 1, ngm do 300 i2 = 1, npb jj = ngm + i2 k = istate2(i2) info(i1,jj) = info(i1,jj) - c1(k)*Sgm(k,i1)*Spb(i2) info(i1,jj) = info(i1,jj) + d1(i1)*d2(i2) info(jj,i1) = info(i1,jj) 300 continue return end subroutine SSt2(ngm,npb,c1,c2,Sgm,Spb,a,info) c Compute E(SS') - E(S)E(S') for day 2 - end of season Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' integer ngm,npb,i,k,l,i1,i2,ii,jj,istate real*8 c1(*),c2(MAXM,*),a(MAXM,*) real*8 d1(MAXNS),d2(MAXNR) real*8 Sgm(MAXM,*),Spb(*) real*8 info(MAXTOT,*) do 100 i1 = 1,ngm i = istate(i1) d1(i1) = 0.0 do 105 l = 1,m 105 d1(i1) = d1(i1) + c2(i,l)*Sgm(l,i1) do 100 i2 = 1,i1 k = istate(i2) if (i.eq.k) then do 110 l = 1,m 110 info(i1,i2) = info(i1,i2) - c2(i,l)*Sgm(l,i1)*Sgm(l,i2) endif info(i1,i2) = info(i1,i2) + d1(i1)*d1(i2) info(i2,i1) = info(i1,i2) 100 continue c write(*,*) (d1(i),i=1,ngm) c do 200 i1 = 1, npb ii = ngm + i1 i = istate2(i1) d2(i1) = c1(i)*Spb(i1) do 200 i2 = 1,i1 jj = ngm + i2 k = istate2(i2) if (i.eq.k) info(ii,jj)=info(ii,jj)-c1(i)*Spb(i1)*Spb(i2) info(ii,jj) = info(ii,jj) + d2(i1)*d2(i2) info(jj,ii) = info(ii,jj) 200 continue c write(*,*) (d2(i),i=1,npb) c do 300 i1 = 1, ngm i = istate(i1) do 300 i2 = 1, npb jj = ngm +i2 k = istate2(i2) info(i1,jj) = info(i1,jj) - c2(i,k)*Sgm(k,i1)*Spb(i2) info(i1,jj) = info(i1,jj) + d1(i1)*d2(i2) info(jj,i1) = info(i1,jj) 300 continue return end real*8 function dmatmin(lda,m,n,a) integer lda,m,n real*8 a(lda,*),small c small = a(1,1) do 10 i=1,m do 10 j=1,n if (a(i,j).lt.small) small = a(i,j) 10 continue dmatmin = small return end subroutine weight3(day,t2,t1,k1,k2a,k2b,k3,rlik,bmat) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'params.inc' include 'dims.inc' include 'data.inc' include 'fb.inc' integer t1,t2,day integer s0,s1,s2,s3 real*8 rlik,k1(MAXM,*),k2a(MAXM,MAXM,*),k2b(MAXM,MAXM,*) real*8 k3(MAXM,MAXM,MAXM,*) real*8 bmat(MAXM,*) real*8 pw1(MAXM),pw2(MAXM),pww1(MAXM,MAXM),pww2(MAXM,MAXM) real*8 a(MAXM,MAXM),a1(MAXM,MAXM) real*8 temp1(MAXM) real*8 pi(MAXM) real*8 scale,scale2,scale3,scale4 c c k1 = P(W_t2, W_t1 | Y) - P(W_t1 | Y)P(W_t2 | Y) c k2a = P(W_{t2-1}, W_t2, W_t1 | Y) - P(W_{t2-1}, W_t2 | Y)P(W_t1 | Y) c k2b = P(W_t2, W_{t1-1}, W_t1 | Y) - P(W_{t1-1}, W_t1 | Y)P(W_t2 | Y) c k3 = P(W_{t2-1},W_t2,W_{t1-1},W_t1 | Y) - c P(W_{t2-1},W_t2 | Y)P(W_{t1-1},W_t1 | Y) do 10 s0 = 1,m do 10 s1 = 1,m k1(s0,s1) = 0.0 do 10 s2 = 1,m k2a(s2,s1,s0) = 0.0 k2b(s2,s1,s0) = 0.0 do 10 s3 = 1,m k3(s3,s2,s1,s0) = 0.0 10 continue c scale = exp(forw(m+1,t2)+back(m+1,t2)-rlik) call makeA(m,gam,natm,mu,isigma,det,adata(1,t2+1),a1) do 200 s0=1,m pw2(s0) = forw(s0,t2)*back(s0,t2)*scale do 200 s1=1,m k1(s0,s1) = forw(s0,t2)*a1(s0,s1) C the following only gets used if t1 = t2+1 k2b(s0,s0,s1) = k1(s0,s1) 200 continue call sumpi(rdata(1,t2),m,nsta,nb,pipars,beta,pi) call makeA(m,gam,natm,mu,isigma,det,adata(1,t2),a) if (day.gt.1) then scale2 = exp(forw(m+1,t2-1)+back(m+1,t2)-rlik) do 220 s1=1,m do 220 s0=1,m pww2(s1,s0)=forw(s1,t2-1)*a(s1,s0)*pi(s0)*back(s0,t2)*scale2 bmat(s1,s0) = a(s1,s0)*pi(s0)*back(s0,t2)/back(s0,t2-1) do 220 s2 = 1,m k2a(s1,s0,s2) = forw(s1,t2-1)*a(s1,s0)*pi(s0)*a1(s0,s2) C the following only gets used if t1 = t2+1 k3(s1,s0,s0,s2) = k2a(s1,s0,s2) 220 continue endif do 300 t = t2+1,t1-1 call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pi) call makeA(m,gam,natm,mu,isigma,det,adata(1,t+1),a1) do 320 s0=1,m do 310 s1=1,m temp1(s1) = 0 do 310 k = 1,m 310 temp1(s1) = temp1(s1) + k1(s0,k)*pi(k)*a1(k,s1) do 320 s1=1,m if (t .eq. t1-1) then do 325 k = 1,m k2b(s0,s1,k) = k1(s0,s1)*pi(s1)*a1(s1,k) 325 continue endif k1(s0,s1) = temp1(s1) 320 continue do 340 s1=1,m do 340 s0=1,m do 350 s2=1,m temp1(s2) = 0 do 350 k = 1,m temp1(s2) = temp1(s2) + k2a(s1,s0,k)*pi(k)*a1(k,s2) 350 continue do 340 s2=1,m if (t .eq. t1-1) then do 345 k = 1,m k3(s1,s0,s2,k) = k2a(s1,s0,s2)*pi(s2)*a1(s2,k) 345 continue endif k2a(s1,s0,s2) = temp1(s2) 340 continue 300 continue scale = exp(forw(m+1,t1)+back(m+1,t1)-rlik) scale2 = exp(forw(m+1,t1-1)+back(m+1,t1)-rlik) scale3 = exp(forw(m+1,t2)+back(m+1,t1)-rlik) if (day.gt.1) scale4 = exp(forw(m+1,t2-1)+back(m+1,t1)-rlik) call sumpi(rdata(1,t1),m,nsta,nb,pipars,beta,pi) call makeA(m,gam,natm,mu,isigma,det,adata(1,t1),a) do 400 s0=1,m pw1(s0) = forw(s0,t1)*back(s0,t1)*scale do 400 s1=1,m k1(s0,s1) = k1(s0,s1)*pi(s1)*back(s1,t1)*scale3 400 continue do 420 s1=1,m do 420 s0=1,m pww1(s1,s0)=forw(s1,t1-1)*a(s1,s0)*pi(s0)*back(s0,t1)*scale2 do 420 s2 = 1,m k2b(s1,s0,s2) = k2b(s1,s0,s2)*pi(s2)*back(s2,t1)*scale3 if (day.eq.1) goto 420 k2a(s1,s0,s2) = k2a(s1,s0,s2)*pi(s2)*back(s2,t1)*scale4 do 410 s3 = 1,m 410 k3(s1,s0,s2,s3)=k3(s1,s0,s2,s3)*pi(s3)*back(s3,t1)*scale4 420 continue c do 500 s0 = 1,m do 500 s1 = 1,m k1(s0,s1) = k1(s0,s1) - pw2(s0)*pw1(s1) do 500 s2 = 1,m k2b(s0,s1,s2) = k2b(s0,s1,s2) - pw2(s0)*pww1(s1,s2) if (day.eq.1) goto 500 k2a(s0,s1,s2) = k2a(s0,s1,s2) - pww2(s0,s1)*pw1(s2) do 510 s3 = 1,m 510 k3(s0,s1,s2,s3) = k3(s0,s1,s2,s3) - pww2(s0,s1)*pww1(s2,s3) 500 continue return end subroutine SSe(day,ngm,npb,k1,k2a,k2b,k3,Sgm1,Sgm2, * Spb1,Spb2,info) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' integer ngm,npb,istate,istate2,s0,s1,day real*8 k1(MAXM,*),k2a(MAXM,MAXM,*),k2b(MAXM,MAXM,*) real*8 k3(MAXM,MAXM,MAXM,*) real*8 Sgm1(MAXM,*),Sgm2(MAXM,*),Spb1(*),Spb2(*) real*8 info(MAXTOT,*) do 100 i1 = 1,ngm i = istate(i1) do 100 i2 = 1,i1 k = istate(i2) do 110 s0 = 1,m do 110 s1 = 1,m if (day.eq.1) then info(i1,i2) = info(i1,i2) - * k2b(s1,i,s0)*Sgm1(s0,i1)*Sgm2(s1,i2)- * k2b(s1,k,s0)*Sgm1(s0,i2)*Sgm2(s1,i1) else info(i1,i2) = info(i1,i2) - * k3(k,s1,i,s0)*Sgm1(s0,i1)*Sgm2(s1,i2)- * k3(i,s1,k,s0)*Sgm1(s0,i2)*Sgm2(s1,i1) endif 110 continue info(i2,i1) = info(i1,i2) 100 continue c do 200 i1 = 1, npb ii = ngm + i1 i = istate2(i1) do 200 i2 = 1,i1 jj = ngm + i2 k = istate2(i2) c write(*,*) day,ii,jj,i1,i2,i,k,info(ii,jj) info(ii,jj)=info(ii,jj)-k1(i,k)*Spb2(i1)*Spb1(i2) * -k1(k,i)*Spb2(i2)*Spb1(i1) info(jj,ii) = info(ii,jj) 200 continue c do 300 i1 = 1, ngm i = istate(i1) do 300 i2 = 1, npb jj = ngm +i2 k = istate2(i2) do 310 s0 = 1,m if (day.eq.1) then info(i1,jj)=info(i1,jj)-k1(s0,k)*Sgm2(s0,i1)*Spb1(i2) * -k2b(k,i,s0)*Sgm1(s0,i1)*Spb2(i2) else info(i1,jj)=info(i1,jj)-k2a(i,s0,k)*Sgm2(s0,i1)*Spb1(i2) * -k2b(k,i,s0)*Sgm1(s0,i1)*Spb2(i2) endif 310 continue info(jj,i1) = info(i1,jj) 300 continue return end viterbi.f000644 000423 000017 00000011666 06704463377 013407 0ustar00hughesusers000000 101230 subroutine viterbi() Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'data.inc' include 'dims.inc' include 'params.inc' include 'fb.inc' c integer year,day,s integer Q(MAXM,MAXOBS) integer argmax,jmax,nnn(MAXM,MAXM) real*8 f(MAXM),g(MAXM) real*8 a(MAXM,MAXM),pisum(MAXM),arow(MAXM) real*8 xbar(MAXM,MAXM,MAXATM),xss(MAXM,MAXM,MAXATM,MAXATM) real*8 sighat(MAXATM,MAXATM),gamsum real*8 arg,temp(MAXATM),newmu(MAXM,MAXM,MAXATM),newgam(MAXM,MAXM) character*80 outfil write(*,*) 'Enter name of output file for state sequence.' read(*,'(a80)') outfil open(4,file=outfil) do 40 i=1,m do 40 j=1,m nnn(i,j)=0 do 40 k=1,natm xbar(i,j,k)=0.0 do 40 l=1,natm 40 xss(i,j,k,l)=0.0 c do 60 year=1,nyears t = year*nobs call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pisum) call copyx(pisum,g,m) call scale(g,m) c now do it for t = dend-1 ... dbeg do 50 day=nobs-1,1,-1 t = t-1 call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pisum) call makeA(m,gam,natm,mu,isigma,det,adata(1,t+1),a) c do 55 s=1,m do 53 i = 1,m 53 arow(i) = a(s,i) jmax = argmax(arow,g,m) Q(s,day+1) = jmax f(s) = pisum(s)*arow(jmax)*g(jmax) c write(*,*) s,jmax,f(s) 55 continue c call scale(f,m) call copyx(f,g,m) 50 continue c and finish up call steady(m,gam,natm,mu,isigma,det,adata(1,t),delta) jmax = argmax(delta,f,m) Q(jmax,1) = jmax write(4,'(365i2)') (Q(jmax,day),day=1,nobs) c c accumulate sums of squares for pooled estimate of sigma c t = (year-1)*nobs + 1 do 62 day=2,nobs t = t+1 i = Q(jmax,day-1) j = Q(jmax,day) nnn(i,j) = nnn(i,j)+1 do 62 k=1,natm xbar(i,j,k) = xbar(i,j,k) + adata(k,t) do 62 l=1,natm 62 xss(i,j,k,l) = xss(i,j,k,l) + adata(k,t)*adata(l,t) 60 continue c c compute pooled estimate of sigma c do 70 i=1,m do 70 j=1,m do 70 k=1,natm 70 if (nnn(i,j).gt.0) xbar(i,j,k) = xbar(i,j,k)/nnn(i,j) df = nyears*nobs-m*m do 71 k=1,natm do 71 l=1,natm sighat(k,l)=0.0 do 72 i=1,m do 72 j=1,m c write(*,*) k,l,i,j,xss(i,j,k,l),xbar(i,j,k),xbar(i,j,l) 72 sighat(k,l) = sighat(k,l) + * (xss(i,j,k,l) - nnn(i,j)*xbar(i,j,k)*xbar(i,j,l)) 71 sighat(k,l) = sighat(k,l)/df c compute rescaled estimates of the parameters c do 370 i=1,natm c370 write(*,*) (isigma(i,j),j=1,natm) do 80 i=1,m gamsum=0.0 do 85 j=1,m do 90 k=1,natm temp(k) = 0.0 do 93 l=1,natm 93 temp(k) = temp(k) + mu(i,j,l)*isigma(l,k) 90 continue c write(*,*) (temp(l),l=1,natm) arg = 0.0 do 95 k=1,natm newmu(i,j,k) = 0.0 do 96 l=1,natm 96 newmu(i,j,k) = newmu(i,j,k) + temp(l)*sighat(l,k) arg = arg + temp(k)*(mu(i,j,k) - newmu(i,j,k)) 95 continue newgam(i,j) = gam(i,j) * dexp(-.5*arg) gamsum = gamsum + newgam(i,j) 85 continue c write(*,*) (newgam(i,j),j=1,m) do 87 j=1,m 87 newgam(i,j)=newgam(i,j)/gamsum 80 continue c print parameters - Bayes form write(6,202) 202 format(//' Rescaled parameters'//' Gamma') do 160 i=1,m write(6,6040)(newgam(i,j),j=1,m) 6040 format(6g12.4) 160 continue write(6,204) 204 format(/' Pipars') do 170 i=1,nsta write(6,6050)(pipars(2,i,j),j=1,m) 6050 format(6g12.4) 170 continue write(6,205) 205 format(/' Beta') do 180 i=1,nb write(6,6050)(beta(i,j),j=1,m) 180 continue write(6,206) 206 format(/' Mu'/' row col of gamma') do 190 i=1,m do 190 j=1,m write(6,6070) i,j,(newmu(i,j,k),k=1,natm) 6070 format(2i4,5g12.4) 190 continue write(6,207) 207 format(/' Sigma') do 200 i=1,natm write(6,6080) (sighat(i,j),j=1,natm) 6080 format(5g12.4) 200 continue c stop end subroutine copyx(x1,x2,n) Implicit double precision (A-H,O-Z) real*8 x1(*),x2(*) do 10 i=1,n 10 x2(i)=x1(i) return end integer function argmax(a,f,m) Implicit double precision (A-H,O-Z) real*8 a(*),f(*),rmax,rcand c c write(*,*) 'argmax' c write(*,*) (a(i),i=1,m) c write(*,*) (f(i),i=1,m) argmax = 1 rmax = a(1)*f(1) do 10 i = 2,m rcand = a(i)*f(i) if (rcand.gt.rmax) then rmax = rcand argmax=i endif 10 continue c write(*,*) argmax return end subroutine scale(f,m) integer m,s real*8 ts,f(*) c ts = 0.0 do 57 s=1,m 57 ts = ts+f(s) ts=ts/dble(m) do 58 s=1,m 58 f(s) = f(s)/ts return end 40 continue 300 continue scale = exp(forw(m+1,t1)+back(m+1,t1)-rlik) scale2 = exp(forw(m+1,t1-1)+back(m+1,t1)-rlik) scale3 = exp(forw(m+1,t2)+back(m+1,t1)-rlik) if (day.gt.1) scale4 = exp(forw(m+1,t2-1)+back(m+1,t1)-rlik) call sumpi(rdata(1,t1),m,nsta,nb,pipars,beta,pi) call makeA(m,gam,natm,mu,isigma,det,adata(1,t1),a) do 400 s0=1,m pw1(s0) = forw(s0,t1)*back(s0,t1)*scale do 400 s1=1,m k1(s0,s1) = k1(s0,s1)*pi(s1)*back(s1,t1)*scale3 400 continue do 420 s1=1,m do 420 s0=1,m pww1(s1,s0)=forw(s1,t1-1)*a(s1,s0)*pi(s0)*back(s0,t1)*scale2 do 420 s2 = 1,m k2b(s1,s0,s2) = k2b(s1,s0,s2)*pi(s2)*back(s2,t1)*scale3 if (day.eq.1) goto 420 k2a(s1,s0,s2) = k2a(s1,s0,s2)*pi(s2)*back(s2,t1)*scale4 do 410 s3 = 1,m 410 k3(s1,s0,s2,s3)=k3(s1,s0,s2,s3)*pi(s3)*back(s3,t1)*scale4 420 continue c do 500 s0 = 1,m do 500 s1 = 1,m k1(s0,s1) = k1(s0,s1) - pw2(s0)*pw1(s1) do 500 s2 = 1,m k2b(s0,s1,s2) = k2b(s0,s1,s2) - pw2(s0)*pww1(s1,s2) if (day.eq.1) goto 500 k2a(s0,s1,s2) = k2a(s0,s1,s2) - pww2(s0,s1)*pw1(s2) do 510 s3 = 1,m 510 k3(s0,s1,s2,s3) = k3(s0,s1,s2,s3) - pww2(s0,s1)*pww1(s2,s3) 500 continue return end subroutine SSe(day,ngm,npb,k1,k2a,k2b,k3,Sgm1,Sgm2, * Spb1,Spb2,info) Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'dims.inc' integer ngm,npb,istate,istate2,s0,s1,day real*8 k1(MAXM,*),k2a(MAXM,MAXM,*),k2b(MAXM,MAXM,*) real*8 k3(MAXM,MAXM,MAXM,*) real*8 Sgm1(MAXM,*),Sgm2(MAXM,*),Spb1(*),Spb2(*) real*8 info(MAXTOT,*) do 100 i1 = 1,ngm i = istate(i1) do 100 i2 = 1,i1 k = istate(i2) do 110 s0 = 1,m do 110 s1 = 1,m if (day.eq.1) then info(i1,i2) = info(i1,i2) - * k2b(s1,i,s0)*Sgm1(s0,i1)*Sgm2(s1,i2)- * k2b(s1,k,s0)*Sgm1(s0,i2)*Sgm2(s1,i1) else info(i1,i2) = info(i1,i2) - * k3(k,s1,i,s0)*Sgm1(s0,i1)*Sgm2(s1,i2)- * k3(i,s1,k,s0)*Sgm1(s0,i2)*Sgm2(s1,i1) endif 110 continue info(i2,i1) = info(i1,i2) 100 continue c do 200 i1 = 1, npb ii = ngm + i1 i = istate2(i1) do 200 i2 = 1,i1 jj = ngm + i2 k = istate2(i2) c write(*,*) day,ii,jj,i1,i2,i,k,info(ii,jj) info(ii,jj)=info(ii,jj)-k1(i,k)*Spb2(i1)*Spb1(i2) * -k1(k,i)*Spb2(i2)*Spb1(i1) info(jj,ii) = info(ii,jj) 200 continue c do 300 i1 = 1, ngm i = istate(i1) do 300 i2 = 1, npb jj = ngm +i2 k = istate2(i2) do 310 s0 = 1,m if (day.eq.1) then info(i1,jj)=info(i1,jj)-k1(s0,k)*Sgm2(s0,i1)*Spb1(i2) * -k2b(k,i,s0)*Sgm1(s0,i1)*Spb2(i2) else info(i1,jj)=info(i1,jj)-k2a(i,s0,k)*Sgm2(s0,i1)*Spb1(i2) * -k2b(k,i,s0)*Sgm1(s0,i1)*Spb2(i2) endif 310 continue info(jj,i1) = info(i1,jj) 300 continue return end viterbi.f000644 000423 000017 00000011666 06704463377 013407 0ustar00hughesusers000000 101230 subroutine viterbi() Implicit double precision (A-H,O-Z) include 'nhmm.inc' include 'data.inc' include 'dims.inc' include 'params.inc' include 'fb.inc' c integer year,day,s integer Q(MAXM,MAXOBS) integer argmax,jmax,nnn(MAXM,MAXM) real*8 f(MAXM),g(MAXM) real*8 a(MAXM,MAXM),pisum(MAXM),arow(MAXM) real*8 xbar(MAXM,MAXM,MAXATM),xss(MAXM,MAXM,MAXATM,MAXATM) real*8 sighat(MAXATM,MAXATM),gamsum real*8 arg,temp(MAXATM),newmu(MAXM,MAXM,MAXATM),newgam(MAXM,MAXM) character*80 outfil write(*,*) 'Enter name of output file for state sequence.' read(*,'(a80)') outfil open(4,file=outfil) do 40 i=1,m do 40 j=1,m nnn(i,j)=0 do 40 k=1,natm xbar(i,j,k)=0.0 do 40 l=1,natm 40 xss(i,j,k,l)=0.0 c do 60 year=1,nyears t = year*nobs call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pisum) call copyx(pisum,g,m) call scale(g,m) c now do it for t = dend-1 ... dbeg do 50 day=nobs-1,1,-1 t = t-1 call sumpi(rdata(1,t),m,nsta,nb,pipars,beta,pisum) call makeA(m,gam,natm,mu,isigma,det,adata(1,t+1),a) c do 55 s=1,m do 53 i = 1,m 53 arow(i) = a(s,i) jmax = argmax(arow,g,m) Q(s,day+1) = jmax f(s) = pisum(s)*arow(jmax)*g(jmax) c write(*,*) s,jmax,f(s) 55 continue c call scale(f,m) call copyx(f,g,m) 50 continue c and finish up call steady(m,gam,natm,mu,isigma,det,adata(1,t),delta) jmax = argmax(delta,f,m) Q(jmax,1) = jmax write(4,'(365i2)') (Q(jmax,day),day=1,nobs) c c accumulate sums of squares for pooled estimate of sigma c t = (year-1)*nobs + 1 do 62 day=2,nobs t = t+1 i = Q(jmax,day-1) j = Q(jmax,day) nnn(i,j) = nnn(i,j)+1 do 62 k=1,natm xbar(i,j,k) = xbar(i,j,k) + adata(k,t) do 62 l=1,natm 62 xss(i,j,k,l) = xss(i,j,k,l) + adata(k,t)*adata(l,t) 60 continue c c compute pooled estimate of sigma c do 70 i=1,m do 70 j=1,m do 70 k=1,natm 70 if (nnn(i,j).gt.0) xbar(i,j,k) = xbar(i,j,k)/nnn(i,j) df = nyears*nobs-m*m do 71 k=1,natm do 71 l=1,natm sighat(k,l)=0.0 do 72 i=1,m do 72 j=1,m c write(*,*) k,l,i,j,xss(i,j,k,l),xbar(i,j,k),xbar(i,j,l) 72 sighat(k,l) = sighat(k,l) + * (xss(i,j,k,l) - nnn(i,j)*xbar(i,j,k)*xbar(i,j,l)) 71 sighat(k,l) = sighat(k,l)/df c compute rescaled estimates of the parameters c do 370 i=1,natm c370 write(*,*) (isigma(i,j),j=1,natm) do 80 i=1,m gamsum=0.0 do 85 j=1,m do 90 k=1,natm temp(k) = 0.0 do 93 l=1,natm 93 temp(k) = temp(k) + mu(i,j,l)*isigma(l,k) 90 continue c write(*,*) (temp(l),l=1,natm) arg = 0.0 do 95 k=1,natm newmu(i,j,k) = 0.0 do 96 l=1,natm 96 newmu(i,j,k) = newmu(i,j,k) + temp(l)*sighat(l,k) arg = arg + temp(k)*(mu(i,j,k) - newmu(i,j,k)) 95 continue newgam(i,j) = gam(i,j) * dexp(-.5*arg) gamsum = gamsum + newgam(i,j) 85 continue c write(*,*) (newgam(i,j),j=1,m) do 87 j=1,m 87 newgam(i,j)=newgam(i,j)/gamsum 80 continue c print parameters - Bayes form write(6,202) 202 format(//' Rescaled parameters'//' Gamma') do 160 i=1,m write(6,6040)(newgam(i,j),j=1,m) 6040 format(6g12.4) 160 continue write(6,204) 204 format(/' Pipars') do 170 i=1,nsta write(6,6050)(pipars(2,i,j),j=1,m) 6050 format(6g12.4) 170 continue write(6,205) 205 format(/' Beta') do 180 i=1,nb write(6,6050)(beta(i,j),j=1,m) 180 continue write(6,206) 206 format(/' Mu'/' row col of gamma') do 190 i=1,m do 190 j=1,m write(6,6070) i,j,(newmu(i,j,k),k=1,natm) 6070 format(2i4,5g12.4) 190 continue write(6,207) 207 format(/' Sigma') do 200 i=1,natm write(6,6080) (sighat(i,j),j=1,natm) 6080 format(5g12.4) 200 continue c stop end subroutine copyx(x1,x2,n) Implicit double precision (A-H,O-Z) real*8 x1(*),x2(*) do 10 i=1,n 10 x2(i)=x1(i) return end integer function argmax(a,f,m) Implicit double precision (A-H,O-Z) real*8 a(*),f(*),rmax,rcand c c write(*,*) 'argmax' c write(*,*) (a(i),i=1,m) c write(*,*) (f(i),i=1,m) argmax = 1 rmax = a(1)*f(1) do 10 i = 2