program curerate 
c      bayesian model
c      ming-hui chen, wpi
c      march 4, 1998
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       parameter (iprint=3,maxlag=20,nobs=50001)
       real*8 y(n),d(n),x(4,n)
       real*8 y1(n),d1(n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 start(11),xmin(11),xsec(11),step(11)
       real*8 ynewlo,ysec,reqmin
       real*8 ebeta1(4),elam1,ealpha1
       real*8 beta(4),alam,alpha
       real*8 ebeta(4),sebeta(4),ealam,sealam
       real*8 sbeta(4),s2beta(4)
       real*8 ealpha,sealpha
       real*8 ao
       real*8 betaupp(4),betalow(4)
       real*8 alamupp,alamlow,alphaupp,alphalow
       real*8 aoupp,aolow
       real*8 seqbeta(4,50001),seqlam(50001)
       real*8 seqalpha(50001),seqao(50001)
       real*8 FUZZ,score(50001)
       real*8  ac(maxlag+1),acv(maxlag+1),seac(maxlag+1),
     1         xacmean
       real*8 ahpd(50001)
       integer imean,iseopt
       integer nq1,nq2,nq,nhpd
       integer ITIE,iscore(50001)
       integer nopt,icount
       integer nrep,iseed
       integer nbig(n),nobig(no)
       common  /vecy/y
       common  /vecd/d
       common  /vecx/x
       common  /vecyo/yo
       common  /vecdo/do
       common  /vecxo/xo
       common  /vecao/ao
       common  /vbeta/beta
       common  /valam/alam
       common  /valpha/alpha
       common /vdeltao/deltao
       common /valamo/alamo
       common /vaccept/ip
       common /vnbig/nbig
       common /vnobig/nobig
       external fn1,dacf
c......set deltao and alamo
       api=3.1415927d0
c      deltao=400.0d0
c      alamo=1.0d0
c      deltao=300.0d0
c      alamo=1.0d0
c      deltao=200.0d0
c      alamo=1.0d0
c      deltao=100.0d0
c      alamo=100.0d0
c      deltao=50.0d0
c      alamo=50.0d0
       deltao=400.0d0
       alamo=1.0d0
       ip=0

c......read current data
       open(unit=15,file='comb.dat',status='old')
       do 25 i=1,n
        x(1,i)=1.0d0
        read(15,*) y(i),d(i),y1(i),d1(i),(x(j,i),j=2,4)
        x(2,i)=(x(2,i)-47.0304457d0)/12.9970773d0
        x(3,i)=x(3,i)-1.0d0
 25    continue
       close(15)

c......read historical data
       open(unit=16,file='hist.dat',status='old')
       do 26 i=1,no
        xo(1,i)=1.0d0
        read(16,*) yo(i),do(i),(xo(j,i),j=2,4)
        xo(2,i)=(xo(2,i)-47.0304457d0)/12.9970773d0
        xo(3,i)=xo(3,i)-1.0d0
 26    continue
       close(16)

       do 27 i=1,n
        if (i .eq. 1) then
          ymin=y(1)
        else
         if (ymin .gt. y(i)) ymin=y(i)
        endif
 27    continue
       nmin=0
       do 28 i=1,n
        if ( (y(i) .eq. ymin) .and. (d(i) .eq. 1.0d0) )
     1    then
         nmin=nmin+1
        endif
 28    continue

       do 37 i=1,no
        if (i .eq. 1) then
          yomin=yo(1)
        else
         if (yomin .gt. yo(i)) yomin=yo(i)
        endif
 37    continue
       nomin=0
       do 38 i=1,no
        if ( (yo(i) .eq. yomin) .and. (do(i) .eq. 1.0d0) )
     1    then
         nomin=nomin+1
        endif
 38    continue

c......read the number of Gibbs iterates
       print*,'Please enter number of Gibbs iterates nrep='
       read(*,*) nrep

c......initialize random seed
       iseed=9999999
       
c......obtain mle for new model
c......maximizing beta
       nopt=6
       konvge=5
       kcount=1000
       reqmin=1.0d-10
       start(1)=-1.54843963d0
       start(2)=0.1839883d0
       start(3)=-0.10191943d0
       start(4)=-0.07243823
       start(5)=0.0d0
       start(6)=dlog(1.0d0/1.77326879d0)
       do 32 k=1,6
         step(k)=0.5d0
         xmin(k)=start(k)
         xsec(k)=start(k)
 32    continue
c......call optim.f
       call nelmin(fn1, nopt, start, xmin, ynewlo, reqmin, step,
     1     konvge, kcount, icount, numres, ifault)
     
       do 100 j=1,4
        ebeta1(j)=xmin(j)
 100   continue
       elam1=xmin(5)
       ealpha1=dexp(xmin(6))
c......set starting point
       do 60 j=1,4
        beta(j)=ebeta1(j)
 60    continue
       ao=0.25d0
c      ao=1.0d0
c      ao=0.0d0
c      ao=deltao/(deltao+alamo)
       alam=elam1
       alpha=ealpha1
c......warm up gibbs
       call rnset( iseed )
       icount=0
       do 51 i=1,1000
        call gibbs( iseed )
        if ( (i .eq. 1) .or. (i .eq. icount) ) then
         write(*,*) 'gibbs warm up i=',i
         write(*,980) (beta(j),j=1,4)
 980     format('beta=',4f12.6)
         write(*,981) ao
 981     format('ao=',f10.6)
         write(*,983) alpha
 983     format('alpha=',f12.6)
         write(*,984) alam
 984     format('alam=',f12.6)
         write(*,*) 'nbig(1)=',nbig(1),' nbig(n)=',nbig(n)
         write(*,*) 'nobig(1)=',nobig(1),' nobig(no)=',nobig(no)
         icount=icount+50
        endif
 51    continue
c......initialize cumulative variables
       sao=0.0d0
       s2ao=0.0d0
       do 61 i=1,4
        sbeta(j)=0.0d0
        s2beta(j)=0.0d0
 61    continue
       salam=0.0d0
       s2alam=0.0d0
       salpha=0.0d0
       s2alpha=0.0d0
c......run gibbs
       icount=0
       do 101 i=1,nrep
        call gibbs( iseed )
        if ( (i .eq. 1) .or. (i .eq. icount) ) then
         write(*,*) 'gibbs  i=',i
         write(*,980) (beta(j),j=1,4)
         write(*,981) ao
         write(*,983) alpha
         write(*,984) alam
         write(*,*) 'nbig(1)=',nbig(1),' nbig(n)=',nbig(n)
         write(*,*) 'nobig(1)=',nobig(1),' nobig(no)=',nobig(no)
         icount=icount+50
        endif
        sao=sao+ao
        s2ao=s2ao+ao**2
        seqao(i)=ao
        do 71 j=1,4
         sbeta(j)=sbeta(j)+beta(j)
         s2beta(j)=s2beta(j)+beta(j)**2
         seqbeta(j,i)=beta(j)
 71     continue
        salam=salam+alam
        s2alam=s2alam+alam**2
        seqlam(i)=alam
        salpha=salpha+alpha
        s2alpha=s2alpha+alpha**2
        seqalpha(i)=alpha
 101   continue
       call rnget( iseed )
c......obtain posterior estimates
       eao=sao/real(nrep)
       seao=dsqrt( (s2ao-real(nrep)*eao**2)
     1         / real(nrep-1) )
       do 80 j=1,4
         ebeta(j)=sbeta(j)/real(nrep)
         sebeta(j)=dsqrt( (s2beta(j)-real(nrep)*ebeta(j)**2)
     1         / real(nrep-1) )
 80    continue
       ealpha=salpha/real(nrep)
       sealpha=dsqrt( (s2alpha-real(nrep)*ealpha**2)
     1         / real(nrep-1) )
       ealam=salam/real(nrep)
       sealam=dsqrt( (s2alam-real(nrep)*ealam**2)
     1         / real(nrep-1) )
c......obtain hpd intervals
c......95% HPD intervals for beta
       q1=0.025d0*real(nrep)
       q2=0.975d0*real(nrep)
       nq1=nint(q1)
       nq2=nint(q2)
       nq=nq2-nq1
       FUZZ=0.0d0
       ITIE=1
       imean=1
       iseopt=1
       do 300 jmac=1,4
        write(*,*) 'hpd for beta j=',jmac
        do 301 i=1,nrep
         ahpd(i)=seqbeta(jmac,i)
 301    continue
        write(*,*) 'autocorrelation for beta',jmac
         call dacf(nrep,ahpd,iprint,iseopt,imean,xacmean,maxlag,
     1            acv,ac,seac)
         call DRANKS(nrep,ahpd,FUZZ,ITIE,0,SCORE)
         do 308 i=1,nrep
          iscore(i)=nint(score(i))
 308     continue
         whpd=1000000.0d0
         do 309 j=1,nrep-nq
          do 310 i=1,nrep
           if (iscore(i) .eq. j) pdiff1=ahpd(i)
           if (iscore(i) .eq. j+nq) pdiff2=ahpd(i)
 310      continue
          wb=pdiff2-pdiff1
          if (whpd .gt. wb) then
           whpd=wb
           aupp=pdiff2
           alow=pdiff1
          endif
 309     continue
         betaupp(jmac)=aupp
         betalow(jmac)=alow
 300   continue
c......95% hpd for alpha
        write(*,*) 'hpd for alpha'
        do 401 i=1,nrep
         ahpd(i)=seqalpha(i)
 401    continue
        write(*,*) 'autocorrelation for alpha'
         call dacf(nrep,ahpd,iprint,iseopt,imean,xacmean,maxlag,
     1            acv,ac,seac)
         call DRANKS(nrep,ahpd,FUZZ,ITIE,0,SCORE)
         do 408 i=1,nrep
          iscore(i)=nint(score(i))
 408     continue
         whpd=1000000.0d0
         do 409 j=1,nrep-nq
          do 410 i=1,nrep
           if (iscore(i) .eq. j) pdiff1=ahpd(i)
           if (iscore(i) .eq. j+nq) pdiff2=ahpd(i)
 410      continue
          wb=pdiff2-pdiff1
          if (whpd .gt. wb) then
           whpd=wb
           aupp=pdiff2
           alow=pdiff1
          endif
 409     continue
         alphaupp=aupp
         alphalow=alow
c......95% hpd for alam
        write(*,*) 'hpd for alam'
        do 501 i=1,nrep
         ahpd(i)=seqlam(i)
 501    continue
        write(*,*) 'autocorrelation for alam'
         call dacf(nrep,ahpd,iprint,iseopt,imean,xacmean,maxlag,
     1            acv,ac,seac)
         call DRANKS(nrep,ahpd,FUZZ,ITIE,0,SCORE)
         do 508 i=1,nrep
          iscore(i)=nint(score(i))
 508     continue
         whpd=1000000.0d0
         do 509 j=1,nrep-nq
          do 510 i=1,nrep
           if (iscore(i) .eq. j) pdiff1=ahpd(i)
           if (iscore(i) .eq. j+nq) pdiff2=ahpd(i)
 510      continue
          wb=pdiff2-pdiff1
          if (whpd .gt. wb) then
           whpd=wb
           aupp=pdiff2
           alow=pdiff1
          endif
 509     continue
         alamupp=aupp
         alamlow=alow
        if (ao .gt. 0.0d0) then
c.......95% hpd for ao
        write(*,*) 'hpd for ao'
        do 601 i=1,nrep
         ahpd(i)=seqao(i)
 601    continue
        write(*,*) 'autocorrelation for ao'
c        call dacf(nrep,ahpd,iprint,iseopt,imean,xacmean,maxlag,
c    1            acv,ac,seac)
         call DRANKS(nrep,ahpd,FUZZ,ITIE,0,SCORE)
         do 608 i=1,nrep
          iscore(i)=nint(score(i))
 608     continue
         whpd=1000000.0d0
         do 609 j=1,nrep-nq
          do 610 i=1,nrep
           if (iscore(i) .eq. j) pdiff1=ahpd(i)
           if (iscore(i) .eq. j+nq) pdiff2=ahpd(i)
 610      continue
          wb=pdiff2-pdiff1
          if (whpd .gt. wb) then
           whpd=wb
           aupp=pdiff2
           alow=pdiff1
          endif
 609     continue
         aoupp=aupp
         aolow=alow
         endif


       open(unit=12,file='bayes.out400_1ab',access='sequential',
     1              status='unknown')
       write(12,*) 'cure rate model'
       write(12,*) 'sample size for current study n=',n
       write(12,*) 'sample size for hist. study no=',no
       write(12,*) '------------------------'
       write(12,*) 'MLE for new model'
       write(12,*) '-------------------------------------'
       write(12,*) 'j   ebeta1' 
       do 400 j=1,4
        write(12,998) j,ebeta1(j)
 998    format(1x,I5,f16.6)
 400   continue
       write(12,*) 'elam1=',elam1
       write(12,*) 'ealpha1=',ealpha1
       write(12,*) '-------------------------------------'
       write(12,*) 'number of Gibbs iterations nrep=',nrep
       priorm=deltao/(deltao+alamo)
       priorstd=dsqrt( deltao*alamo/(deltao+alamo+1.0d0) )
       priorstd=priorstd/(deltao+alamo)
       write(12,*) 'deltao=',deltao,' alamo=',alamo
       write(12,*) 'E(ao)=',priorm,' se=',priorstd
       write(12,*) '-------------------------------------'
       write(12,*) 'information needed for setting up a prior'
       write(12,*) 'ymin=',ymin,' nmin=',nmin
       write(12,*) 'yomin=',yomin,' nomin=',nomin
       write(12,*) '-------------------------------------'
       write(12,*) 'posterior estimates'
       write(12,1002) eao,seao,aolow,aoupp
 1002  format('eao=',f10.6,' seao=',f10.6,'(',f10.6,',',
     1         f10.6,')')
       write(12,*) 'acceptance prob for ao=',
     1            real(ip)*10/(1000.0d0+real(nrep))
       write(12,*) '------------------------------------------'
       write(12,*) 'k      ebeta       sebeta         95% HPD Int.'
       do 901 j=1,4
        write(12,1001) j-1, ebeta(j),sebeta(j),betalow(j),betaupp(j)
 1001   format(1x,I4,2f16.8,'(',f12.6,',',f12.6,')')
 901   continue
       write(12,*) '------------------------------------------'
       write(12,*) 'ealpha        sealpha        95% HPD Int.'
       write(12,1003) ealpha,sealpha,alphalow,alphaupp
 1003  format(1x,2f16.8,'(',f12.6,',',f12.6,')')
       write(12,*) '------------------------------------------'
       write(12,*) 'ealam         sealam         95% HPD Int.'
       write(12,1003) ealam,sealam,alamlow,alamupp
       write(12,*) '------------------------------------------'
       close(12)
       stop
       end



       real*8 function fn1(z)
       implicit real*8 (a-h,o-z)
       parameter( n=284 )
       real*8 y(n),d(n),x(4,n)
       real*8 z(6)
       real*8 beta(4),alam,alpha
       common  /vecy/y
       common  /vecd/d
       common  /vecx/x
       do 50 j=1,4
        beta(j)=z(j)
 50    continue
       alam=z(5)
       alpha=dexp(z(6))
       sum=0.0d0
       do 100 i=1,n
        ein=0.0d0
        do 110 j=1,4
         ein=ein+x(j,i)*beta(j)
 110    continue
        ss=dexp(-y(i)**alpha*dexp(alam))
        if (d(i) .eq. 0.0d0) then
         sum=sum-dexp(ein)*(1.0d0-ss)
        endif
        if (d(i) .eq. 1.0d0) then
         sum=sum-dexp(ein)*(1.0d0-ss)
         sum=sum+ein+dlog(alpha)+(alpha-1.0d0)*dlog(y(i))
     1       +alam-y(i)**alpha*dexp(alam)
        endif
 100   continue 
       fn1=-sum
       return
       end


       include 'optim1.f'
       include 'gibbs.f'