program lynems c Version 4 : containing interface with incubation distribution c modified for AZT treatment according to Brookmeyer & Liao (1990) c Changes to program have been marked : !* c c and c A choice of 4 smoothing functions: Binomial,Epanechnikov,Exponential c and Uniform. October 1990 c A further choice of smoothing: variable. February 1991 c This program estimates HIV intensities from observed AIDS incidences c and an assumed known incubation distribution using the lynems algorithm. c Future AIDS incidences can also be predicted. c It has been written by Lyn Watson, at LaTrobe University, March 1990. c include 'lynems.inc' c lynems.inc contains declaration statements. c Setting parameter defaults epsilon = 1.0e-7 ! Convergence ncon = 0 ! " n = 1 ! Bandwidth (no smoothing) npred = 0 ! No prediction nperyear = 12 ! monthly data nrep = 100 ! A progress report to terminal every 100 iters nsmo = 1 ! Default smoother is Binomial nfix = 1 ! Default smoothing into fixed tail ncond = 1 ! Default conditioning on n_td (old) sim = 0 ! Flag for simulation run hivinit = 50. ! Initial staring lambdas nfixlam = 24 ! No. of fixed lambdas for recent times hivfixed = 50. ! Value of fixed lambdas for recent times nttt = 1 !* Ttt change to incubation dist flag write(*,*) & ' This program estimates HIV intensities from observed AIDS' write(*,*) & ' incidences and an assumed known incubation distribution ' write(*,*) ' using the ems algorithm.' write(*,*) ' Future AIDS incidences can also be predicted.' write(*,*) write(*,*) & ' On the first run, a new set of AIDS data must be used,' write(*,*) ' a new incubation distribution specified and' write(*,*) ' new starting HIVs set.' write(*,*) 160 write(*,*) ' Describe AIDS data, (<30 char),to end type: stop )' read(*,9) heading 9 format(a) if (heading .eq.'stop')stop ' ' CALL namerun ! sets up output files 'runname'.lis (report) ! .cub (incubation probs) and .dat (graph files) if(.not.yesno(' New AIDS data, Y/N ? '))goto 170 CALL aidstime ! Sets length of data array,layout for report ! and reads news AIDS data file. c heading for report c 170 CALL idate(mm,md,my) write (11,79) md,mm,my,heading,aidsfile 79 format(' ! ',3i3,x,a,' file:',a/) iconv = 0 ! Convergence flag diffmax = 0 ! Monitor for diff between old and new lampdas. write(*,29)epsilon,ncon 29 format(' ! Convergence epsilon is ',E10.1,' , ignoring last ', && i4,' estimated HIVs') if(yesno(' Change convergence criteria, Y/N ? '))then write(*,*)' Type in epsilon, then # recent HIVs to be ignored' read(*,*)epsilon,ncon else end if write(11,29)epsilon,ncon if (nfix.eq.1)then if(.not. & yesno(' Smoothing is into fixed tail, change this, Y/N ? ')) & go to 120 else if(.not. & yesno(' Smoothing is indep of tail, change this, Y/N ? ')) & go to 120 endif write(*,*)' Type of endpoint smoothing: 1 into fixed tail' write(*,*)' 2 indep of tail' write(*,*)' Type 1 or 2' read(*,*)nfix 120 if (ncond.eq.1)then if(.not.yesno(' n_td (old) conditioning, change this, Y/N ? ')) & go to 130 else if (ncond.eq.2) then if(.not.yesno(' n_t (new) conditioning, change this, Y/N ? ')) & go to 130 else if (ncond.eq.3) then if(.not. & yesno(' intermediate conditioning, change this, Y/N ? ')) & go to 130 else if(.not.yesno(' varying inter. cond., change this, Y/N ? ')) & go to 130 endif endif endif write(*,*)' Type of conditioning: 1 Old, on n_td' write(*,*)' 2 New, on n_t' write(*,*)' 3 Intermediate' write(*,*)' 4 Varying intermediate' write(*,*)' Type 1, 2, 3 or 4' read(*,*)ncond if ((ncond.eq.3).or.(ncond.eq.4)) then write(*,*) ' Input tau :' read(*,*) tau endif 130 write(*,59)npred 59 format(' Number of prediction intervals are ',i4) if(yesno(' Change # prediction intervals, Y/N ? '))then write(*,39) 39 format( ' Type in # prediction intervals wanted : ',$) read(*,*) npred else endif CALL timevec ! Calculates a time vector for outfile CALL incubate ! Either reads incubation distribution from file ! or calculates it from Weibull with parameters ! set at terminal CALL initlamda ! Sets starting lambdas (ie HIV intensities). CALL window ! Sets weights for smoothing step. CALL expaids ! Calculates expected AIDS counts used in ! denominator of E-step calculation. CALL loglike ! Calculates penalized loglikelihood ! not used although printed in reports. write(*,69)nrep 69 format(' Number of iterations between interim reports is ',i4) if (yesno(' Change this, Y/N ? '))then write(*,109) 109 format(' Number of iterations between interim reports : ',$) read(*,*)nrep else endif ntot = 0 ! Counter for total # iterations for this run. niter = 200 CALL report ! Prints initial interim report to screen & file. 110 write (*,99)niter 99 format(' Number of iterations to try is',i6) if (yesno(' Change this, Y/N ? '))then write(*,89) 89 format(' How many more iterations? 999 to finish ',$) read(*,*)niter else endif if (niter .eq. 999) then iconv = 2 CALL report ! Prints final report to screen & file ! when stopping without convergence. ! Also prints AIDS predictions and std resids. go to 160 ! Try another run else end if do 100 iter = 1,niter ntot = ntot + 1 CALL newlambda ! EM-step. ! Calculates new lambda's (phi's in our paper) if (n.ne.1) CALL smooth ! Smoothes the phi's to form the mu's ! (called new lambda's in paper) CALL converge ! Calculates the difference between old and new ! lambda's then sets the new lambdas into the ! old lambda array. It calculates a new array of ! 'expected AIDS', and the loglikelihood of this ! set of 'expected AIDS'. ! It finally checks convergence via difference ! between old lambda sum and new lambda sum. ! It returns with convergence flag set to one ! if convegence criterion has been met. if (iconv.eq.1)then CALL report ! Prints final report to screen & file of ! estimated HIVs, expected and predicted AIDS ! and standardized residuals. go to 160 ! Return to start for new run else end if nmod = mod(ntot,nrep) ! Checking if interim report need be written. if (nmod.eq.0)CALL report 100 continue go to 110 ! Go ask if more iterations wanted. end LOGICAL FUNCTION yesno(str) character*(*) str character*1 response 100 write(6,10) str 10 format(X,A,X,$) read(5,20) response 20 format(A) if (response.eq.'Y' .or. response.eq.'y') then yesno=.TRUE. else if (response.eq.'N' .or. response.eq.'n') then yesno=.FALSE. else write(6,30) 30 format(X,'Invalid response - Enter Y, y ,N, or n') goto 100 endif return end SUBROUTINE namerun include 'lynems.inc' write (*,79) 79 format(' Write run output filename, <= 8 char without ', & 'extension ',$) read(*,9)runname 9 format(a) do 20 nc = 1,8 if ( runname(nc:nc).eq.' ') go to 30 20 continue 30 nc = nc-1 repfile = runname( :nc)//'.lis' datafile = runname( :nc)//'.dat' incubfile = runname( :nc)//'.cub' open (unit = 11,file = repfile, status = 'new') return end SUBROUTINE AIDStime include 'lynems.inc' write (*,109) 109 format(' AIDS filename ? ',$) read(*,9)aidsfile 9 format(a) open (unit =10,file = aidsfile, status = 'old') write(*,19) 19 format(' Type in #AIDS data points to be used ',$) read(*,*) t write(*,39) 39 format(' Type in #initial zeros for AIDS vector for which HIV ', & 'will be calculated '/' so that time origin is Jan 1979 : ',$) read(*,*) num0 origin = 79 write(11,29)t,num0 29 format(' ! ',i6,' AIDS data points preceeded by ',i4,' zeros') write(*,49)nperyear 49 format( ' ! Number time intervals per year is ',i2) if(yesno( ' Change # time intervals per year, Y/N ? '))then write(*,59) 59 format(' Number time intervals per year is : ',$) read(*,*)nperyear else endif write(11,49)nperyear c Working out lines for report tables x = (t+num0)/nperyear line = nint(x) c Reading new AIDS data c initialising AIDS incidence vector do 10 i = 1,num0 A(i) = 0 10 continue c read (10,*) (A(i),i = num0+1,t + num0) close (unit = 10) ntotaids = 0 do 20 i = num0+1,t+num0 ntotaids = ntotaids + A(i) 20 continue totaids = real(ntotaids) write(*,*)totaids return end c SUBROUTINE timevec include 'lynems.inc' c Working out time vector outdata(1,1) = origin + real(1.0/nperyear) do 9 i = 2,num0+t+npred outdata(1,i) = outdata(1,i-1) + real(1.0/nperyear) 9 continue return end SUBROUTINE incubate c include 'lynems.inc' c if (nttt.eq.0) go to 50 write(*,69)incufile 69 format(' Incubation distribution is not set or is from file ',a) if(yesno(' New incubation distribution, Y/N ? '))go to 20 write(11,69)incufile return c20 if(yesno(' Brookmeyer & Liao''s ttt modification, Y/N ? ')) !* c & go to 40 !* 20 if(yesno(' Is the new prob. dist. to be read from file ? ')) & then i = t + num0 + npred write(*,49)i,i 49 format(' The incubation probability distribution file needs to ' & /' be of length ',i4,' and width ',i4) write(*,59) 59 format(' Incubation pdf filename ? ',$) read(*,9)incufile 9 format(a) CALL readincub nttt = 1 ! Set to use 2 way array for pdf and cdf return else endif 30 write(*,*)' Assuming Weibull, S(t)=exp-(bt)**c,t in yrs,' // & 'enter b & c' read(*,*) b,c write(*,29)b,c 29 format(' b = ',f6.4,' and c = ',f6.4) if(.not.yesno(' Correct, Y/N ? '))go to 30 write(11,19) b,c 19 format(' ! Assuming Weibull, S(t) = exp-(',f6.4,'t)**',f6.4) c c calulating cdf and pdf c if ((ncond.eq.1).or.(ncond.eq.2)) taul=t if ((ncond.eq.3).or.(ncond.eq.4)) taul=tau cf(1) = 0 bpert = b/nperyear do 10 i = 1,taul+num0+npred x = (i*bpert)**c cf(i+1) = 1 - exp(-x) df(i) = cf(i+1) - cf(i) 10 continue nttt = 0 ! Set to indicate indicate Weibull dist if(.not.yesno(' Write incubation distribution to file, Y/N ? ')) & return open (unit = 12,file = incubfile, status = 'new') write(12,*)' ! Cumulative distribution function ' write(*,*)' ! Cumulative distribution function ' CALL printtab(cf,49) ! Prints cdf in a table to screen & file. write(12,*)' ! Probability distribution function ' write(*,*)' ! Probability distribution function ' CALL printtab(df,49) ! Prints pdf in a table to screen & file. close (unit = 12) return 50 write(*,39)b,c 39 format(' Incubation dist currently Weibull (',f6.4,','f6.4,')') if(yesno(' Change incubation distribution, Y/N ? '))go to 20 write(11,19)b,c return c40 CALL tttincub !* c nttt = 1 !* c return !* end c SUBROUTINE readincub include 'lynems.inc' open (unit = 19,file = incufile, status = 'old') do 10 i = 1,t+num0+npred read (19,*) (pdf(i,j),j = 1,t+num0+npred) 10 continue close (unit = 19) write(*,59) 59 format(' Incubation cdf filename ? ',$) read(*,9)incufile 9 format(a) open (unit = 19,file = incufile, status = 'old') do 20 i = 1,t+num0+npred read (19,*) (cdf(i,j),j = 1,t+num0+npred) 20 continue close (unit = 19) return end SUBROUTINE printtab(y,nf) include 'lynems.inc' if(nf.eq.49)then do 40 i = 1,line write (12,49) (y(j),j = (i-1)*nperyear+1,i*nperyear) write (*,49) (y(j),j = (i-1)*nperyear+1,i*nperyear) 40 continue write (12,49) (y(j),j = line*nperyear+1,t+num0) write (*,49) (y(j),j = line*nperyear+1,t+num0) else if (nf.eq.99)then if (nperyear.le.6)then do 20 i = 1,line write (11,29) (y(j),j = (i-1)*nperyear+1,i*nperyear) write (*,29) (y(j),j = (i-1)*nperyear+1,i*nperyear) 20 continue write (11,29) (y(j),j = line*nperyear+1,t+num0) write (*,29) (y(j),j = line*nperyear+1,t+num0) else do 90 i = 1,line write (11,99) (y(j),j = (i-1)*nperyear+1,i*nperyear) write (*,99) (y(j),j = (i-1)*nperyear+1,i*nperyear) 90 continue write (11,99) (y(j),j = line*nperyear+1,t+num0) write (*,99) (y(j),j = line*nperyear+1,t+num0) endif else do 30 i = 1,line write (11,39) (y(j),j = (i-1)*nperyear+1,i*nperyear) write (*,39) (y(j),j = (i-1)*nperyear+1,i*nperyear) 30 continue write (11,39) (y(j),j = line*nperyear+1,t+num0) write (*,39) (y(j),j = line*nperyear+1,t+num0) endif 29 format(6(x,f10.1)) 39 format(12(x,f5.2)) 49 format(12(x,f5.3)) 99 format(12(x,f5.1)) return end c SUBROUTINE initlamda c include 'lynems.inc' character*24 hivfile !** for simulation c write(*,59)hivinit,nfixlam,HIVfixed 59 format(' Starting HIV''s all ',f8.2,' and the last ',i3, & ' set at ',f8.2) if(.not.yesno(' Change starting HIV''s, Y/N ? '))go to 50 write(*,29) 29 format(' Type in # fixed lamdas for recent times : ',$) read(*,*)nfixlam if(nfixlam.le.0)then HIVfixed = 0. go to 30 else write(*,39) 39 format(' Type in value for this fixed lamda : ',$) read(*,*)HIVfixed endif 30 if(yesno(' HIV from input file, Y/N ? '))go to 60 !** write(*,19) 19 format(' Type in starting HIV per time unit (at least 1) ',$) read(*,*)HIVinit 50 write(11,49)nfixlam,HIVfixed 49 format(' ! Last ',i2,' HIVs are fixed at ',f8.2) write(11,9)HIVinit 9 format(' ! Starting HIVs are ',f8.2) if(nfixlam.gt.0) ncon = nfixlam c Initial settings for lamda's do 20 i = 1,t+ num0-nfixlam lamda(i) = HIVinit 20 continue if(nfixlam.le.0)return do 40 i = t+ num0-nfixlam+1,t+num0 lamda(i) = HIVfixed 40 continue return 60 write(*,69) !** 69 format( 'HIV intensity file name ? ' ) !** read(*,79)HIVfile !** 79 format(a) !** open(unit = 15, file = hivfile, status = 'old') !** read(15,*)(lamda(i),i=1,t+num0) !** close(unit = 15) !** return !** end c SUBROUTINE expaids c include 'lynems.inc' c c Calculating expected AIDS from lamda's and incubation distribution c totlam = 0 totsum = 0 c calculate up to t+npred for simulations do 140 i = 1,t+num0+npred totlam = totlam + lamda(i) sum(i) = 0 if (nttt.eq.0) then !* do 30 j = 1,i sum(i) = sum(i) + lamda(j)*df(i-j+1) 30 continue else !* do 31 j = 1,i !* sum(i) = sum(i) + lamda(j)*pdf(i-j+1,j) !* 31 continue !* endif !* 140 continue do 141 i=1,t+num0 totsum = totsum + sum(i) 141 continue return end c SUBROUTINE loglike c c Calculating log likelihood c include 'lynems.inc' c loglkhd = -totsum do 90 i = 1,t+num0 c if (sum(i).gt.0.0) & loglkhd = loglkhd + A(i)*log(sum(i)) 90 continue return end c SUBROUTINE smooth c include 'lynems.inc' c if (n.eq.1) return m = (n-1)/2 c This next do loop sets future mu's equal to final mu for less c biased smoothing if(nfix.eq.1) then do 31 i = 1,m mu(t+num0-nfixlam+i) = hivfixed ! this should enable a ! smooth fit to fixed lambdas 31 continue else do 30 i = 1,m mu(t+num0-nfixlam+i) = mu(t+num0-nfixlam) 30 continue endif if (nsmo.ne.5) then do 10 i = 1,t+num0-nfixlam smu = 0 do 20 j = 1,n c if (((i-m+j-1).gt.0).and.((i-m+j-1).le.(t+num0-nfixlam))) c 1 smu = smu + mu(i-m+j-1)*w(j) if ((i-m+j-1).le.0) smu = smu + mu(1)*w(j) if ((i-m+j-1).gt.0) smu = smu + mu(i-m+j-1)*w(j) 20 continue mu(i) = smu 10 continue else c This is for progressive smoothing, in three steps, for monthly data c starting Jan 79. do 40 i = 1,42 smu = 0. do 45 j = 1,11 if ((i-m+j-1).gt.0) smu = smu + mu(i-m+j-1)*wvar(j,1) ! ** 45 continue mu(i) = smu 40 continue do 50 i = 43,84 smu = 0. do 55 j = 1,11 if ((i-m+j-1).gt.0) smu = smu + mu(i-m+j-1)*wvar(j,2) ! ** 55 continue mu(i) = smu 50 continue do 60 i = 85,t+num0-nfixlam smu = 0. do 65 j = 1,11 if ((i-m+j-1).gt.0) smu = smu + mu(i-m+j-1)*wvar(j,3) ! ** 65 continue mu(i) = smu 60 continue endif return end c SUBROUTINE converge c check convergence include 'lynems.inc' c calculating max difference between individual lamdas for this c iteration and replacing lamdas with new estimates c convergence sum on t + num0 - ncon (ncon usually 24) lamdas c convergence modification by nb 20.4.90 diffmax = 0.0 diff = 0.0 totlamsq = 0.0 do 90 i = 1,t + num0 - ncon ! ncon = nfixlam if nfixlam gt 0 x = abs(mu(i) - lamda(i)) if (x.gt.diffmax)diffmax = x c calculating numerator for new convergence criterion, 2.1.91 diff = diff + x**2.0 c calculating denominator " " " " " totlamsq = totlamsq + lamda(i)**2.0 c replacing old lambdas with new lamda(i) = mu(i) 90 continue 80 CALL expaids ! Updates expected AIDS counts. CALL loglike ! Calculates penalized loglikelihood. if (ntot.lt.100)return ! Inserted 31.8.90 to prevent premature ! convergence conlam = diff/totlamsq if (conlam.ge.epsilon)return iconv = 1 return end SUBROUTINE report c include 'lynems.inc' c if (ntot.ne.0)go to 5 conlam = 0. write(11,109)datafile 109 format(x,' ! Output data from this run is in ',a/) write(*,79) write(11,79) 79 format(/,' ! Iter Loglkhd logdiff HIVs & epsi AIDStot ') hmax = 0 c Finding max HIV 5 hmax = 0. do 6 i = 1,t+num0 if (lamda(i).le.hmax)go to 6 hmax = lamda(i) k = i 6 continue hend = lamda(t+num0-nfixlam) write (*,89) ntot,loglkhd,logdiff,totlam,conlam,totsum write (11,89) ntot,loglkhd,logdiff,totlam,conlam,totsum 89 format (' !',i5,x,f13.4,x,f13.4,x,f11.1,x,E9.2,x,f10.1) if(iconv.eq.0)return if(sim.eq.1)return if(iconv.eq.1)write(11,59)ntot 59 format(' ! Converged after ',i6,' iterations',/) if(iconv.eq.2)write(11,69)ntot 69 format(' ! Stopped after ',i6,' iterations',/) if(nperyear.eq.360)go to 80 ! HIVs and AIDS not printed for daily data c Printing the table of estimated HIVs write(11,*)' ! Estimated HIV incidences' write(*,*)' ! Estimated HIV incidences' CALL printtab(lamda,99) if(yesno(' Write expected AIDS to file, Y/N ? '))then write(11,*)' ! Expected AIDS incidences' write(*,*)' ! Expected AIDS incidences' CALL printtab(sum,99) else endif c Now predict future AIDS incidences if wanted c Writing observed and expected AIDS and estimated HIV to outdata array do 60 i = 1,num0+t outdata(2,i) = A(i) outdata(3,i) = lamda(i) outdata(4,i) = sum(i) 60 continue ncol = 4 if (npred.gt.0)CALL predict if (npred.eq.0)go to 80 do 70 i = num0+t+1,num0+t+npred outdata(2,i) = outdata(4,i) ! future obs as prediction ! from zero future HIVs 70 continue 80 if (yesno(' Goodness of fit, residuals or simulate, Y/N ? ')) & CALL residuals close (unit = 11) c ! Writes to file time,obs AIDS,estimated c ! HIV, expected AIDS predictions and resids if wanted for graphs if(.not.yesno(' Write output datafile, Y/N ? '))return open (unit = 13,file = datafile, status = 'new') write (*,*)num0,t,npred do 10 j = 1,num0+t+npred write(13,99) (outdata(i,j),i=1,ncol) 99 format(f10.3,x,f8.1,x,10(f8.1,x),f7.1) 10 continue close (unit = 13) return end SUBROUTINE predict include 'lynems.inc' c First prediction assuming no new infections extralam = 0. write(11,*)' ! Prediction assuming no new HIV infections ' CALL printpred(extralam) if(nfixlam.ne.0.and.HIVfixed.eq.0.0)return if(nfix.eq.2)go to 60 c Assuming future HIV's are positive fixed value extralam = hivfixed write(11,9)hivfixed 9 format(' ! Prediction assuming future HIV infections are ',f8.1) ncol = ncol + 1 do 11 i = 1,num0+t outdata(ncol,i) = lamda(i) outdata(ncol+1,i) = sum(i) 11 continue ncol = ncol + 1 CALL printpred(extralam) return 60 write(*,89) 89 format(' Choose prediction type: 0 stop, 2 Xtn , ', & ' 3 Hilo '/' Option 3 is not sensible if window is 1 or last ', & ' HIVs are fixed.') read(*,*)nopt if (nopt.eq.0)return if (nopt.eq.2)then extralam = lamda(t+num0) nt = t+num0 write(11,79)nt,extralam 79 format(' ! Prediction assuming HIV infections for each time ', & 'period since',i4,' are ',f8.1) ncol = ncol + 1 do 10 i = 1,num0+t outdata(ncol,i) = lamda(i) outdata(ncol+1,i) = sum(i) 10 continue ncol = ncol + 1 CALL printpred(extralam) else hmax = 0 c Finding max HIV do 5 i = 1,t+num0 if (lamda(i).le.hmax)go to 5 hmax = lamda(i) nmax = i 5 continue c Prediction with HIVs as 2/3 down between max & last value third = (hmax-lamda(t+num0))/3.0 extralam = third + lamda(t+num0) i = nmax do while (lamda(i).gt.extralam) i = i+1 end do nlo = i ncol = ncol + 1 do 20 j = 1,nlo-1 outdata(ncol,j) = lamda(j) outdata(ncol+1,j) = sum(j) 20 continue do 7 j = nlo,t+num0 lamda(j) = extralam outdata(ncol,j) = lamda(j) outdata(ncol+1,j) = sum(j) 7 continue ncol = ncol + 1 totlam = 0 do 8 i = 1,t+num0 totlam = totlam + lamda(i) 8 continue write(11,69)extralam,nlo,totlam 69 format(' ! Prediction assuming ',f8.1,' for HIV for each time ', & 'interval since ',i3/' ! This gives total HIVs as ',f10.1) CALL printpred(extralam) c Prediction of HIVs 1/3 between max & last value extralam = third + extralam i = nmax do while (lamda(i).gt.extralam) i = i+1 end do nmed = i ncol = ncol + 1 do 30 j = 1,nmed-1 outdata(ncol,j) = lamda(j) outdata(ncol+1,j) = sum(j) 30 continue do 17 j = nmed,t+num0 lamda(j) = extralam outdata(ncol,j) = lamda(j) outdata(ncol+1,j) = sum(j) 17 continue ncol = ncol + 1 totlam = 0 do 18 i = 1,t+num0 totlam = totlam + lamda(i) 18 continue write(11,69)extralam,nmed,totlam CALL printpred(extralam) c Predicting HIVs constant from max HIV. This is high extralam = hmax ncol = ncol + 1 do 40 j = 1,nmax-1 outdata(ncol,j) = lamda(j) outdata(ncol+1,j) = sum(j) 40 continue do 27 j = nmax,t+num0 lamda(j) = extralam outdata(ncol,j) = lamda(j) outdata(ncol+1,j) = sum(j) 27 continue ncol = ncol + 1 totlam = 0 do 28 i = 1,t+num0 totlam = totlam + lamda(i) 28 continue write(11,69)extralam,nmax,totlam CALL printpred(extralam) c Some of the lamdas have been changed in the previous prediction c routines and need to be restored in case any other predictions wanted. do 70 i = 1,t + num0 lamda(i) = mu(i) 70 continue endif go to 60 end SUBROUTINE printpred(extralam) include 'lynems.inc' do 6 i = t+num0+1,t+num0+npred lamda(i) = extralam outdata(ncol-1,i) = extralam 6 continue totpred = 0 do 10 j = 1,npred m = t+num0+j sum(m) = 0 if (nttt.eq.0) then !* do 20 i = 1,m sum(m) = sum(m) + lamda(i)*df(m-i+1) 20 continue else !* do 21 i = 1,m !* sum(m) = sum(m) + lamda(i)*pdf(m-i+1,i) !* 21 continue !* endif !* outdata(ncol,m) = sum(m) totpred = totpred + sum(m) 10 continue write(11,25)npred write(*,25)npred 25 format( ' ! Predictions for next ',i6,' time intervals are :'/) x = npred/nperyear nline = nint(x) if(nperyear.gt.6)then do 30 i = 1,nline write(11,29)(sum(j),j=t+num0+1+(i-1)*nperyear,t+num0+i*nperyear) write(*,29)(sum(j),j=t+num0+1+(i-1)*nperyear,t+num0+i*nperyear) 30 continue write (11,29) (sum(j),j = nline*nperyear+t+num0+1,t+num0+npred) write (*,29) (sum(j),j = nline*nperyear+t+num0+1,t+num0+npred) 29 format (' !',12(x,f5.1)) else do 40 i = 1,nline write(11,59)(sum(j),j=t+num0+1+(i-1)*nperyear,t+num0+i*nperyear) write(*,59)(sum(j),j=t+num0+1+(i-1)*nperyear,t+num0+i*nperyear) 40 continue write (11,59) (sum(j),j = nline*nperyear+t+num0+1,t+num0+npred) write (*,59) (sum(j),j = nline*nperyear+t+num0+1,t+num0+npred) 59 format (' !',6(x,f10.0)) endif write(11,39)totpred 39 format(/' ! Total prediction = ',f9.0//) write (*,49)npred,totpred 49 format(' Total prediction for next ',i3,' time intervals is ', & f9.0) return end SUBROUTINE residuals include 'lynems.inc' c if (ntot.eq.0)go to 20 gof = 0 ncol = ncol+1 do 10 i = 1,t+num0 if(sum(i).gt.0) res(i) = (A(i)-sum(i))/sqrt(sum(i)) gof = gof + res(i)**2 outdata(ncol,i) = res(i) 10 continue write(11,19)gof write(*,19)gof 19 format(' ! Goodness of fit = ',f8.2) CALL printtab(res,39) c 20 if (yesno(' Simulate from this run, Y/N ? '))CALL simpois return end SUBROUTINE newlambda c Calculating new mu include 'lynems.inc' if (nttt.eq.0) then !* do 60 i = 1,t+num0-nfixlam sum1 = 0 do 50 j = 1,t-i+num0+1 if (A(i+j-1).ne.0) & sum1 = sum1 + A(i+j-1)*df(j)/sum(i+j-1) 50 continue if (ncond.eq.1) then mu(i) = lamda(i)/cf(t+num0+2-i)*sum1 else if (ncond.eq.2) then mu(i) = lamda(i)*(sum1 + (1 - cf(t+num0+2-i))) else if (ncond.eq.3) then mu(i) = lamda(i)*(sum1 +(cf(tau+num0+2-i)-cf(t+num0+2-i))) mu(i) = mu(i)/cf(tau+num0+2-i) else if ((t-i).gt.(tau-t)) then mu(i)=lamda(i)/cf(t+num0+2-i)*sum1 else mu(i)=lamda(i)*(sum1+(cf(tau-t-num0+2)- & cf(t+num0+2-i))) mu(i)=mu(i)/cf(tau-t-num0+2) endif endif endif 60 continue else !* do 61 i = 1,t+num0-nfixlam !* sum1 = 0 !* do 51 j = 1,t-i+num0+1 !* if (A(i+j-1).ne.0) sum1 = sum1 +A(i+j-1)*pdf(j,i)/sum(i+j-1) !* 51 continue !* if(cdf(t+num0+1-i,i).lt.0.000001)then mu(i) = mu(i-1) else if (ncond.eq.1) then mu(i) = lamda(i)/cdf(t+num0+1-i,i)*sum1 else if (ncond.eq.2) then mu(i) = lamda(i)*(sum1 + (1 - cdf(t+num0+1-i,i))) !* else if (ncond.eq.3) then mu(i) = lamda(i)*(sum1+(cdf(tau+num0+1-i,i)- & cdf(t+num0+1-i,i))) mu(i) = mu(i)/cdf(tau+num0+1-i,i) else if ((t-i).gt.(tau-t)) then mu(i)=lamda(i)/cdf(t+num0+1-i,i)*sum1 else mu(i)=lamda(i)*(sum1+(cdf(tau-t-num0+1,i)- & cdf(t+num0+1-i,i))) mu(i)=mu(i)/cdf(tau-t-num0+1,i) endif endif endif endif 61 continue !* endif !* return end SUBROUTINE window include 'lynems.inc' character*12 smoform(5) real ll,ul,bnd real normconst,g(6) data g/4.0,3.0,2.0,1.0,0.5,0.25/ data smoform/'binomial','Own choice','exponential','uniform', & 'progressive'/ write(*,19) smoform(nsmo),n ! nsmo initially set at 1(=Binomial), c and n as 1(=no smoothing), c in main program 19 format(' ! Smoothing by ',a, ' weights, bandwidth is',i3) if(.not.yesno(' Change, Y/N ? '))then write(11,19) smoform(nsmo),n return else endif write(*,161) 161 format( ' Type of smoothing: 1 Binomial'/ & ' 2 Own choice'/ & ' 3 Exponential'/ & ' 4 Uniform'/ & ' 5 Progressive, only for monthly data'/ & ' Choose 1,2,3 4 or 5') read(*,*)nsmo if (nsmo.eq.5) n=11 if (nsmo.eq.5) goto 80 write(*,162) 162 format( ' Bandwidth? Choose odd 1 to 13 (1 means no smoothing)') read(*,*)n write(11,19) smoform(nsmo),n if (n.eq.1) return ! No smoothing if (nsmo.eq.1)then ! Binomial weights d = n-1 w(1) = 1/(2.0**d) do 10 i = 2,n w(i) = (real(n)-real(i)+1.0)/(real(i)-1)*w(i-1) 10 continue return else if(nsmo-3)2,3,4 endif 2 continue ! Own choice c2 bnd = -real(n/2.0) ! Epanechnikov, parameter = 2 c do 40 i = 1,n c ul = bnd + real(i) c ll = ul - 1.0 c w(i) = 0.75/(-bnd)*(1.0-(ul**2.0+ll**2.0+ul*ll)/(3.0*bnd**2.0)) c40 continue m2 = int((n-1)/2.0) + 1 write(*,*) ' Type in weights in ascending order ' read(*,*)(w(i),i=1,m2) do 40 i = 1,m2-1 w(m2+i) = w(m2-i) 40 continue write(11,*)w return 3 m2 = int((n-1)/2.0) ! Exponential d = n*g(m2)/2.0 normconst = (2.0*(1.0-exp(-d)))**(-1.0) bnd = -real(n)/2.0 do 50 i = 1,m2 ul = bnd + real(i) ll = ul - 1.0 w(i) = normconst*(exp(g(m2)*ul)-exp(g(m2)*ll)) 50 continue w(m2+1) = 2.0*normconst*(1.0-exp(-0.5*g(m2))) do 60 i = n-m2+1,n w(i) = w(n-i+1) 60 continue write(*,*)w return 4 do 70 i = 1,n ! Uniform w(i) = 1.0/real(n) 70 continue write(*,*)w return 80 CALL windowva write(11,19) smoform(nsmo),n return end SUBROUTINE windowva include 'lynems.inc' double precision x(11),wsum do 40 i = 1,11 x(i) = i x(i) = x(i)/12.0 x(i) = x(i) - x(i)**2.0 x(i) = 4.0*x(i) 40 continue wsum1 = 0.0 wsum2 = 0.0 do 20 i = 1,11 wvar(i,1) = x(i)**95.5 ! ad hoc value to mimic var of simvar wsum1 = wsum1 + wvar(i,1) wvar(i,2) = x(i)**5.5 ! " " wsum2 = wsum2 + wvar(i,2) wvar(i,3) = 1./11. ! Uniform smoothing for last 2+ years 20 continue c normalizing the weights do 30 i = 1,11 wvar(i,1) = wvar(i,1)/wsum1 wvar(i,2) = wvar(i,2)/wsum2 30 continue return end SUBROUTINE iterate include 'lynems.inc' iflagx=0 if (iflag.eq.1) then do 228 ix=1,n wy(ix)=w(ix) w(ix)=wx(ix) 228 continue endif 177 continue c re-initialize settings for lambda's do 40 i = 1,t+num0 lamda(i) = orighiv(i) 40 continue do 407 i=t+num0+1,t+num0+npred lamda(i) = lamda(t+num0) 407 continue 60 CALL expaids CALL loglike ntot = 0 iconv = 0 diffmax = 0 do while (iconv.ne.1) ntot = ntot + 1 CALL newlambda if (n.ne.1)CALL smooth CALL converge end do if (iflagx.eq.0) then xl1=loglkhd else logdiff=xl1-loglkhd return endif if (iflag.eq.1) then do 229 ix=1,n w(ix)=wy(ix) 229 continue iflagx=1 goto 177 endif return end