program bcsurvboot
c
c       Box-Cox survival family analysis with percentile method bootstrap CI's.
c
        implicit double precision (a-h,o-z)
        parameter(kmax=21,kbmax=20,nmax=5000,step=25.0d0,isfl=1)
        parameter(maxloop=100,maxboot=1000)
        parameter (tol=1.0d-4,gtol=1.0d-4)
        external pslik, psgrad, marlik, margrd, dbcong, iercd,
     $    erset, du4inf, deqtil, rnopt, rnset

        dimension pars(kmax), t(nmax), del(nmax), z(nmax,kbmax)
        dimension pmin(kmax), pmax(kmax), r(nmax), bta(kbmax)
        dimension rtil(nmax), pold(kmax), rold(nmax)
        dimension rsk(nmax), evcnt(nmax), ilow(nmax), iup(nmax)
        dimension pgues(kmax), pscal(kmax), ipara(7), rpara(7)
        dimension dpar(kmax), parin(kmax)
        dimension epar(maxboot,kmax), ep(maxboot)
        dimension torig(nmax), delorig(nmax), zorig(nmax,kbmax)
        dimension zz(nmax), zsd(kmax)
        dimension qprop(2), quant(2), xlo(2), xhi(2)
        character*20 inpfn

        common /data1/ n, k, kb
        common /data2/ t, del, z
        common /data3/ rsk, evcnt, ilow, iup, ndist
        common /data4/ r
        common /data5/ rtil
        common /data6/ torig, delorig, zorig

        open(1,file='cntlr.crd',status='old')
        open(2,file='bcr.out',status='unknown',access='append')
        read(1,*) inpfn
        read(1,*) imeth
        read(1,*) irecy
        read(1,*) n
        read(1,*) kb
        read(1,*) nboorep
        read(1,*) idum
        read(1,*) big1, rlitl2, big2
        read(1,*) conlev
        if (imeth .eq. 1)
     $   then
           write(2,*) 'Pseudo-likelihood method'
         else
           write(2,*) 'Martingale residual method'
         endif
        write(2,*)
        write(2,*) 'Setup parameters:'
        write(2,*) '  Input file = ', inpfn
        write(2,*) '  Sample size = ', n
        write(2,*) '  Number of covariates = ', kb
        write(2,*) '  Confidence interval level = ', conlev
        write(2,*) '  Number of bootstrap replications = ', nboorep
        write(2,*) '  Bootstrap starting value option = ', irecy
        write(2,*) '  Random number seed = ', idum
        write(2,*) '  Bound on abs(beta) = ', big1
        write(2,*) '  Lower bound on rho = ', rlitl2
        write(2,*) '  Upper bound on rho = ', big2
        write(2,*)

        boot = nboorep
        k = kb + 1
        ibt = 0
        fscal = 1.0d0
        call erset(0,0,0)
        istflg = isfl
        call rnopt(2)
        call rnset(idum)

        read(1,*) (parin(l), l = 1, k)

        do 10  l = 1, kb
          pmin(l) = -big1
          pmax(l) = big1
          pscal(l) = 1.0d0
   10    continue
        pmin(kb+1) = rlitl2
        pmax(kb+1) = big2
        pscal(kb+1) = 1.0d0

        open(3,file=inpfn,status='old')
        do 12  i = 1, n
          read(3,*) (zorig(i,l), l = 1, kb), torig(i), delorig(i)
   12    continue

        do 18  l = 1, kb
          do 14  i = 1, n
            zz(i) = zorig(i,l)
   14      continue
          call simstat(n,zz,zmean,sd)
          zsd(l) = sd
          do 16  i = 1, n
            zorig(i,l) = (zorig(i,l)-zmean)/sd
   16      continue
   18     continue
         zsd(kb+1) = 1.0d0

        iboot = 0
        nboot = 0
        nbgud = 0

   20   continue

        if (iboot .eq. 0)
     $    then
            do 25  i = 1, n
              t(i) = torig(i)
              del(i) = delorig(i)
              do 22  l = 1, kb
                z(i,l) = zorig(i,l)
   22          continue
   25         continue
          else
            nboot = nboot + 1
            call bootgen(n,kb)
         endif

        call sort

        if ((iboot .eq. 0) .or. (irecy .eq. 0))
     $  then
        do 60  l = 1, k
          pars(l) = parin(l)
   60    continue
        else
        do 65  l = 1, k
          pars(l) = dpar(l)
   65    continue
        endif

        do 70  j = 1, kb
          bta(j) = pars(j)
   70    continue
        rho = pars(kb+1)
        call rcomp(bta,rho,iflg)
        if (iflg .eq. 1) goto 110

        nloops = 0

   75   continue

        nloops = nloops + 1

        do 76  l = 1, k
          pold(l) = pars(l)
          pgues(l) = pars(l)
   76    continue
        do 78  i = 1, n
          rold(i) = r(i)
   78    continue

        call du4inf(ipara,rpara)
        rpara(1) = gtol
        if (istflg .eq. 1) rpara(6) = step
        if (imeth .eq. 1)
     $    then
            call dbcong(pslik,psgrad,k,pgues,ibt,pmin,pmax,pscal,fscal,
     $        ipara,rpara,pars,fval)
          else
            call rtilcmp(pars)
            call dbcong(marlik,margrd,k,pgues,ibt,pmin,pmax,pscal,fscal,
     $        ipara,rpara,pars,fval)
          endif

        iflg=0
        ier = iercd()
        if (ier .ne. 0) iflg=1
        rho = pars(kb+1)
        if (rho .ge. big2) iflg=1
        do 80  j = 1, kb
          bta(j) = pars(j)
          if (dabs(pars(j)) .ge. big1) iflg=1
   80    continue
        if (iflg .eq. 1) goto 110
        call rcomp(bta,rho,iflg)
        if (iflg .eq. 1) goto 110

        l = 1
        do 85 while ((l .le. kb) .and. (iflg .eq. 0))
          diff = dabs(pars(l) - pold(l))
          if (diff .ne. 0.0d0) then
            rdiff = 2.0d0*diff/(pold(l)+pars(l))
            if (rdiff .gt. tol) iflg=1
           endif
          l = l + 1
   85    continue
        diff = dabs(pars(kb+1) - pold(kb+1))
        if (diff .gt. tol) iflg=1
        if (iflg .eq. 1) goto 95

        m = 1
        do 90 while ((m .le. ndist) .and. (iflg .eq. 0))
          diff = dabs(r(m) - rold(m))
          if (diff .ne. 0.0d0) then
            rdiff = 2.0d0*diff/(rold(m)+r(m))
            if (rdiff .gt. tol) iflg=1
           endif
          m = m + 1
   90    continue

   95   continue

        if (iflg .eq. 1) then
          if (nloops .lt. maxloop) then
            goto 75
          else
            goto 110
          endif
         endif

        if (iboot .eq. 1) then
          nbgud = nbgud + 1
          do 97  l = 1, k
            epar(nbgud,l)=pars(l)
   97      continue
         endif

  110   continue

        if ((iboot .eq. 0) .and. (iflg .eq. 1)) then
          write(2,*) 'Convergence failure'
          goto 500
         endif

        if (iboot .eq. 0) then
          write(2,*) ' Parameter estimates:'
          do 120  l = 1, k
            dpar(l) = pars(l)
            rpar = pars(l)/zsd(l)
            write(2,115) l, rpar
  115       format(1x,i3,f10.4)
  120      continue
          write(2,*)
          iboot = 1
         endif

        if ((iboot .eq. 1) .and. (nboot .lt. nboorep)) goto 20

        if (nbgud .lt. 2) then
          write(2,*) 'Bootstrap failure'
          goto 500
         endif

        dboo = nbgud
        pconv = dboo/boot
        write(2,*) 'Percent bootstrap reps that converged = ', pconv
        write(2,*)
        write(2,*) 'Percentile method boostrap confidence limits:'
        do 270  l = 1, k
          nqprop = 2
          qprop(1) = (1.0d0-conlev)/2.0d0
          qprop(2) = 1.0d0 - qprop(1)
          do 230  i = 1, nbgud
            ep(i) = epar(i,l)
  230      continue
          call deqtil(nbgud, ep, nqprop, qprop, quant, xlo, xhi, nmiss)
          clo = quant(1)/zsd(l)
          chi = quant(2)/zsd(l)
          write(2,265) l, clo, chi
  265     format(1x,i3,2F10.4)
  270    continue

  500   continue

        write(2,*)
        write(2,*)

        end


        subroutine bootgen(n,kb)
        implicit double precision(a-h,o-z)
        parameter(kbmax=20,nmax=5000,fuzz=2.0d-3)
        external rnund, drnunf

        dimension torig(nmax), delorig(nmax), zorig(nmax,kbmax)
        dimension t(nmax), del(nmax), z(nmax,kbmax), ir(nmax)

        common /data2/ t, del, z
        common /data6/ torig, delorig, zorig

        nn = n
        call rnund(nn,nn,ir)
        do 30  i = 1, n
          t(i) = torig(ir(i))*(1.0d0+fuzz*(drnunf()-0.50d0))
          del(i) = delorig(ir(i))
          do 20  j = 1, kb
            z(i,j) = zorig(ir(i),j)
   20      continue
   30    continue
        return
        end


        subroutine simstat(nn,xx,dmean,sd)
        implicit double precision(a-h,o-z)
        dimension xx(nn)
        dnn = nn
        dmean = 0.0d0
        sumd = 0.0d0
        sumd2 = 0.0d0
        do 130  i = 1, nn
          dmean = dmean + xx(i)
  130    continue
        dmean = dmean/dnn
        do 135  i = 1, nn
          dev = xx(i)-dmean
          sumd = sumd + dev
          sumd2 = sumd2 + (dev**2)
  135    continue
        var = (sumd2-(sumd**2)/dnn)/(dnn-1.0d0)
        sd=dsqrt(var)
        return
        end


        include 'bcsubs.f'