C*****FILE RPVAX07.F C JOB SPECIFIC SOURCE FILE FOR PROGRAM MTG C*****BROWNIAN MOTION PARAMETERS OF TRAJECTORY WITH KNOWN CENSUS ERROR STD C.....ESTIMATE DENNIS MODEL PVA PARAMETERS FROM POPULATION TRAJECTORY C ASSUMING UNDERLYING DENNIS-LIKE MODEL WITH LOG-NORMAL C OBSERVATION ERROR PRESENT ON THE CENSUS OBSERVATIONS WITH KNOWN C LOG SPACE STD DEV SUPPLIED FOR EACH CENSUS AS PART OF THE DATA, C AND WITH PRIOR INCLUDING COMPONENTS FOR STATE AT EACH TIME STEP C AS WELL AS FOR STATE TRANSITIONS C.....NOTATION C NDT NUMBER OF TIME STEPS C NC NUMBER OF CENSUSES C EN POPULATION SIZE C AN LOG POPULATION SIZE C X CENSUS POINT ESTIMATE LOGNORMAL DISTRIBUTED AND UNBIASED C ABOUT AN IN THE LOG SPACE C XS ORIGINAL SPACE ERROR STD OF THE POINT ESTIMATE C AX LOG OF CENSUS POINT ESTIMATE UNBIASED IN THE LOG SPACE C AXS LOG SPACE ERROR STD OF THE (LOG) POINT ESTIMATE C ER MEAN EXPONENTIAL GROWTH RATE C SR STD DEV OF NORMAL DISTRIBUTION OF EXPONENTIAL GROWTH RATE C ASR LOG OF STD DEV OF NORMAL DISTRIBUTION OF EXPONENTIAL GROWTH RATE C ERBAR MEAN FOR NORMAL PRIOR ON ER C ERSTD STD DEV FOR NORMAL PRIOR ON ER C SRBR MEAN FOR NORMAL PRIOR ON SR (IGNORING TRUNCATION, SR>0) C SRST STD DEV OF NORMAL PRIOR ON SR (IGNORING TRUNCATION, SR>0) C ENB MEAN OF NORMAL PRIOR DISTRIBUTION OF POPULATION SIZE C (IGNORING TRUNCATION EN>0) C ENS STD OF NORMAL PRIOR DISTRIBUTION OF POPULATION SIZE C (IGNORING TRUNCATION EN>0) C.....PROCESS MODEL (STATE EQUATION) C AN(T+1)~N([AN(T)+ER],[SR*SR]) C.....OBSERVATION MODEL C AX(T)~N(AN(T),(AXS(T)**2) C IN OTHER WORDS C AX=AN+EPS WHERE EPS~N(0,AXS**2) C.....PRIOR COMPONENTS C ER~N(ERBAR,ERSTD**2) C SR~N(SRBR,SRST**2) (SUBJECT TO SR>0) C EN(T)~N(ENB(T),ENS(T)**2) (SUBJECT TO EN>0) C AND PROCESS MODEL PRODUCT SEQUENCE C.....JOB SET UP C DATA FILE OF NDT RECORDS AND 4 FIELDS C WHERE NDT IS THE NUMBER OF TIME STEPS STARTING FROM THE FIRST C CENSUS AND ENDING WITH THE LAST CENSUS, AND THE FIELDS ARE C (1) LOG CENSUS POINT ESTIMATE (IN LOG SPACE) C CODED AS -99 FOR NO CENSUS THAT TIME STEP C (2) LOG SPACE ERROR STD DEV OF POINT ESTIMATE C CODED AS -99 FOR NO CENSUS THAT TIME STEP C (3) MODE OF NORMAL PRIOR ON POPULATION SIZE IN C THAT TIME STEP (ORIGINAL SPACE) C (4) STD DEV OF NORMAL PRIOR ON POPULATION SIZE IN C THAT TIME STEP (ORIGINAL SPACE), IGNORING THE TRUNCATION AT 0 C.....JOB FILE SUPPLIES C SPV(1) ERBAR MEAN FOR NORMAL PRIOR ON ER C SPV(2) ERSTD STD DEV FOR NORMAL PRIOR ON ER C SPV(3) SRBR MEAN FOR NORMAL PRIOR ON SR (IGNORING TRUNCATION) C SPV(4) SRST STD FOR NORMAL PRIOR ON SR (IGNORING TRUNCATION) C.....PARAMETERS C [1] POPULATION SIZE IN TIME STEP 1 C ... C [NDT] POPULATION SIZE IN TIME STEP NDT C [NDT+1] MEAN OF EXPONENTIAL GROWTH RATE, ER C [NDT+2] STD OF EXPONENTIAL GROWTH RATE, SR C***** C VERSION OF APR 26, 2011 C***** SUBROUTINE MTGIN(LUNS,LUNF,SPV,NDT,MDT,DTV,CFNAM,FWV,IWV) C*****JOB-SPECIFIC DEDICATED SUBROUTINE THAT INITIALIZES FIXED PARAMETERS C FOR PROGRAM MTG C.....GIVEN C LUNS DEVICE NUMBER FOR SCREEN, PRESUMED OPEN C LUNF DEVICE NUMBER FOR OUTPUT FILE, PRESUMED OPEN C SPV VECTOR OF FIXED PARAMETER VALUES, READ FROM JOB FILE C NDT NUMBER OF OBSERVATIONS IN DATA, SPECIFIED IN JOB FILE C MDT NUMBER OF VARIABLES PER OBSERVATION IN DATA, SPECIFIED C IN JOB FILE C DTV FLOATING POINT ARRAY OF DATA ELEMENTS, IN SINGLE SUBSCRIPT C STORAGE, WHERE ELEMENT IJ=(I-1)*M+J IS THE J-TH VARIABLE C FROM OBSERVATION I, AS READ FROM THE DATA FILE SPECIFIED C IN THE JOB FILE C.....SELF-IDENTIFIES BY WRITING NAME OF THE FILE OF JOB-SPECIFIC C SUBROUTINES, AS HARD CODED IN THE FIRST USER-MODIFIABLE C EXECUTABLE STATEMENT, TO LUNS AND LUNF, AND C.....RETURNS C CFNAM CHARACTER*8 NAME, SET IN THE FIRST USER-MODIFIABLE C EXECUTABLE STATEMENT, THAT MUST MATCH THE JOB SPECIFIC C SUBROUTINE FILE NAME SPECIFIED IN THE JOB FILE, AND MUST C MATCH THE FILE NAME OF THE ACTUAL FILE (THIS FILE) CONTAINING C THE SOURCE CODE FOR MTGIN AND MTGPP C SPV POSSIBLY REARRANGED, TRANSFORMED, SUBSET OR AUGMENTED; WILL C BE COMMUNICATED TO SUBROUTINE MTGPP IN THIS FORM C DTV POSSIBLY REARRANGED, TRANSFORMED, SUBSET OR AUGMENTED; WILL C BE COMMUNICATED TO SUBROUTINE MTGPP IN THIS FORM C FWV VECTOR AVAILABLE FOR FLOATING POINT WORKSPACE AND COMMUNICATING C ADDITIONAL FLOATING POINT VALUES TO SUBROUTINE MTGPP C IWV VECTOR AVAILABLE FOR INTEGER WORKSPACE AND COMMUNICATING C INTEGER VALUES TO SUBROUTINE MTGPP C***** CHARACTER*8 CFNAM DIMENSION DTV(1),SPV(1),FWV(1),IWV(1) C========BEGIN BLOCK OF JOB-SPECIFIC CODE=============================== C.....SELF-IDENTIFY WITH CHARACTER STRING OF UP TO 8 CHARACTERS, IN C SINGLE QUOTES (MANDATORY; MUST MATCH NAME OF FILE OF JOB CODE). C DO NOT CHANGE THE NAME OF THE VARIABLE CFNAM. C..... CFNAM='rpvax07 ' C.....RECEIVE C NDT NUMBER OF TIME STEPS IN DATA C DTV ELEMENT IJ=(I-1)*4+J IS DATA RECORD I WHERE FOR C J=1: LOG CENSUS POINT ESTIMATE (IN LOG SPACE) C CODED AS -99 FOR NO CENSUS THAT TIME STEP C J=2: LOG SPACE ERROR STD DEV OF POINT ESTIMATE C CODED AS -99 FOR NO CENSUS THAT TIME STEP C J=3: MODE OF NORMAL PRIOR ON POPULATION SIZE IN C THAT TIME STEP (ORIGINAL SPACE) C J=4: ORIGINAL SPACE STD DEV OF NORMAL PRIOR ON POPULATION SIZE C IN THAT TIME STEP, IGNORING THE TRUNCATION AT 0 C C.....LOOP THROUGH DATA TIME STEPS C AND CODE IWV 1-NDT AS C 1 FOR CENSUS PRESENT C 0 FOR MISSING C..... IJ=0 DO 1040 I=1,NDT IJ=IJ+1 AX=DTV(IJ) IF (AX+98.0) 1010,1020,1020 1010 IWV(I)=0 GO TO 1030 1020 IWV(I)=1 1030 IJ=IJ+3 1040 CONTINUE C.....LOAD NDT-1 IN SPV(5) SPV(5)=FLOAT(NDT-1) C========END THIS BLOCK OF JOB SPECIFIC CODE============================ C.....OPERATION COMPLETE 9999 RETURN END C----- SUBROUTINE MTGPP(NDT,MDT,DTV,NN,PRV,SPV,FWV,IWV,PPP,IFLG) C*****LOG PROPORTIONAL POSTERIOR PROBABILITY FOR PROGRAM MTG C.....GIVEN C NDT NUMBER OF OBSERVATIONS IN DATA, SPECIFIED IN THE JOB FILE C MDT NUMBER OF VARIABLES PER OBSERVATION IN DATA, SPECIFIED C IN THE JOB FILE C DTV FLOATING POINT ARRAY OF DATA ELEMENTS, IN SINGLE SUBSCRIPT C STORAGE, WHERE ELEMENT IJ=(I-1)*M+J IS THE J-TH VARIABLE C FROM OBSERVATION I, AS READ FROM THE DATA FILE SPECIFIED C IN THE JOB FILE (AND POSSIBLY MODIFIED BY MTGIN) C NN NUMBER OF UNKNOWN PARAMETERS FOR INFERENCE, SPECIFIED C IN THE JOB FILE C PRV VECTOR CURRENT TRIAL VALUES OF THE NN PARAMETERS C SPV VECTOR OF FIXED PARAMETER VALUES, AS SPECIFIED IN THE C JOB FILE AND POSSIBLY MODIFIED IN MTGIN C FWV VECTOR OF FLOATING WORKSPACE AND POSSIBLE FLOATING POINT C VALUES COMMUNICATED FROM SUBROUTINE MTGIN C IWV VECTOR OF INTEGER WORKSPACE AND POSSIBLE INTEGER VALUES C COMMUNICATED FROM SUBROUTINE MTGIN C.....RETURNS C PPP LOG PROPORTIONAL POSTERIOR PROBABILITY C IFLG FLAG: 0 SIGNALS ZERO PROBABILITY, 1 OTHERWISE C***** DIMENSION DTV(1),PRV(1),SPV(1),FWV(1),IWV(1) C========BEGIN BLOCK OF JOB-SPECIFIC CODE=============================== C.....UNLOAD GLOBAL PRIOR SPECIFICATIONS C SPV(1) ERBAR MEAN FOR NORMAL PRIOR ON ER C SPV(2) ERSTD STD DEV FOR NORMAL PRIOR ON ER C SPV(3) SRBR MEAN FOR NORMAL PRIOR ON SR C SPV(4) SRST STD FOR NORMAL PRIOR ON SR C SPV(5) FNM NUMBER OF STATE TRANSITIONS C..... ERBAR=SPV(1) ERSTD=SPV(2) SRBR=SPV(3) SRST=SPV(4) FNM=SPV(5) C.....UNLOAD DYNAMIC PARAMETERS ER=PRV(NDT+1) SR=PRV(NDT+2) IF (SR) 1010,1010,1020 1010 IFLG=0 GO TO 9999 1020 ASR=ALOG(SR) C.....PRIOR COMPONENT ON DYNAMIC PARAMETERS PPP=-((ER-ERBAR)/ERSTD)**2-((SR-SRBR)/SRST)**2 C.....FIRST DATA RECORD AND ASSOCIATED STATE DESCRIPTION IJ=0 I=1 ISW=IWV(I) EN=PRV(I) IF (EN) 1010,1010,1030 1030 IJ=IJ+1 AX=DTV(IJ) IJ=IJ+1 AXS=DTV(IJ) IJ=IJ+1 ENB=DTV(IJ) IJ=IJ+1 ENS=DTV(IJ) C.....PRIOR CONTRIBUTION OF STATE PPP=PPP-((EN-ENB)/ENS)**2 AN=ALOG(EN) C.....LIKELIHOOD CONTRIBUTION OF OBSERVATION IF (ISW) 9999,1050,1040 1040 PPP=PPP-((AX-AN)/AXS)**2 C.....ADVANCE STATE AND CACHE 1050 ANC=AN+ER C.....LOOP THROUGH REMAINING DATA RECORDS AND ASSOCIATED STATE DESCRIPTION DO 1090 I=2,NDT ISW=IWV(I) EN=PRV(I) IF (EN) 1010,1010,1060 1060 IJ=IJ+1 AX=DTV(IJ) IJ=IJ+1 AXS=DTV(IJ) IJ=IJ+1 ENB=DTV(IJ) IJ=IJ+1 ENS=DTV(IJ) C.....PRIOR CONTRIBUTION OF STATE PPP=PPP-((EN-ENB)/ENS)**2 AN=ALOG(EN) C.....PRIOR CONTRIBUTION OF STATE TRANSITION IGNORING DENOMINATOR TERM PPP=PPP-((AN-ANC)/SR)**2 C.....LIKELIHOOD CONTRIBUTION OF OBSERVATION IF (ISW) 9999,1080,1070 1070 PPP=PPP-((AX-AN)/AXS)**2 C.....ADVANCE STATE AND CACHE 1080 ANC=AN+ER 1090 CONTINUE C.....SCALE NORMAL EXPONENT (NUMERATOR) TERMS PPP=PPP/2.0 C.....INCORPORATE DENOMINATOR TERM FOR THE TRANSITIONS PPP=PPP-FNM*ASR C.....POSTERIOR PROBABILITY COMPLETE IFLG=1 C========END THIS BLOCK OF JOB SPECIFIC CODE============================ C.....OPERATION COMPLETE 9999 RETURN END C***** C END OF FILE C*****