       program partionweightestimation 
c      ming-hui chen
c      wpi, 7/17/99
       implicit real*8 (a-h,o-z)
       parameter (n=183)
       real*8  y(n),x(8,n)
       real*8  beta(7)
       real*8  sbeta(7),s2beta(7)
       real*8  ebeta(7),sebeta(7)
       real*8  sbetap(7),s2betap(7)
       real*8  ebetap(7),sebetap(7)
       real*8  FUZZ,score(100001)
       real*8  seqbeta(7,100001)
       real*8  seqb(100001)
       real*8  api
       real*8  amuc(7)
       real*8  q1,q2,pb1,pb2,wb,whpd,aupp,alow
       real*8  betaup(4),betalow(4)
       real   tarray1(2), tarray2(2)
       real   tarray3(2)
       real*8  seqlike(100001)
       real*8  b(10),a(10),p(10),w(10)
       real*8  sb(10),sp(10)
       real*8  b6(10),a6(10),p6(10),w6(10)
       real*8  sb6(10),sp6(10)
       integer nq1,nq2,nq,nhpd
       integer iseed,idum
       integer ITIE,iscore(100001)
       integer nobs,nrep,npr
       integer idumd(7)
       COMMON /vecy/y
       COMMON /vecx/x
       COMMON /vecbeta/beta
       COMMON /vecjmac/jmac
       external DLINDS
       external DLFTDS,DLFDDS
       api=3.1415927d0
c
c
c......set true values
       iseed=8888888

c......read the current data
       open(unit=16,file='current.data',status='old')
       do 100 i=1,183
        x(1,i)=1.0d0
        read(16,*) (x(j,i),j=2,7),y(i)
        x(2,i)=(x(2,i)-334.6658566d0)/109.3475625d0
        x(3,i)=(x(3,i)-34.6403402d0)/7.6790501d0
        x(4,i)=(x(4,i)-0.5078979d0)/0.5002416d0
        x(5,i)=(x(5,i)-0.9137303d0)/0.2809327d0
        x(6,i)=(x(6,i)-0.8907104d0)/0.3128582d0
        x(7,i)=(x(7,i)-0.4262295d0)/0.4958847d0
 100   continue
c      close(16)
c......read the number of Gibbs iterates
       print*,'Please enter number of Gibbs iterates nrep='
       read(*,*) nrep
c
c......calculate asymptotic variance for current data
       write(*,*) 'asymptotic variance for current data'
       amuc(1)=-4.4397d0
       amuc(2)=-1.5502d0
       amuc(3)=0.1190d0
       amuc(4)=-0.0427d0
       amuc(5)=0.0500d0
       amuc(6)=-0.1829d0
       amuc(7)=0.1288d0
c


       call etime(tarray1)
c......initialize
       call rnset( iseed )
c......initial values of beta(2),b
       do 44 j=1,7
        beta(j)=amuc(j)
        if (j .gt. 5) beta(j)=0.0d0
 44    continue
       ao=0.5d0
       bo=0.0d0
c......warm up gibbs
       icount=0
       ipao=0
       ipbo=0
       do 200 i=1,500
        call gibbslog(iseed)
        if ( (i .eq. 1) .or. (i .eq. icount) ) then
         write(*,*) 'warm up (log) i=',i
         write(*,888) (beta(j),j=1,7)
 888     format('beta=',1x,7f16.6)
         icount=icount+100
        endif
 200   continue
c......initialize the cumulative variables
       do 211 j=1,7
         sbeta(j)=0.0d0
         s2beta(j)=0.0d0
 211   continue
       icount=0
       sLn=0.0d0
       s2Ln=0.0d0
       sr=0.0d0
       s2r=0.0d0
       do 210 i=1,nrep
        do 216 j=1,10
         jmac=j
         call gibbslog(iseed)
 216    continue
        if ( (i .eq. 1) .or. (i .eq. icount) ) then
         write(*,*) 'Gibbs (log) i=',i
         write(*,888) (beta(j),j=1,7)
         icount=icount+20
        endif
        do 212 j=1,5
         seqbeta(j,i)=beta(j)
         sbeta(j)=sbeta(j)+beta(j)
         s2beta(j)=s2beta(j)+beta(j)**2
 212    continue
        slog=0.0d0
        do 217 i1=1,n
         ein=0.0d0
         do 218 j=1,5
          ein=ein+beta(j)*x(j,i1)
 218     continue
         if (ein .gt. -70.0d0) then
          ein1=dexp(ein)
         else
          ein1=0.0d0
         endif
         if (y(i1) .eq. 1.0d0) then
          slog=slog+dlog(1.0d0-dexp(-ein1))
         endif
         if (y(i1) .eq. 0.0d0) then
          slog=slog-ein1
         endif
         sum1=0.0d0
         if (ein .gt. 0.0d0) then
           sum1=sum1+ein*(y(i1)-1.0d0)
           if (ein .le. 50.0d0) then
            sum1=sum1-dlog(1.0d0+dexp(-ein))
           endif
         else
           sum1=sum1+ein*y(i1)
           if (ein .gt. -50.0d0) then
            sum1=sum1-dlog(1.0d0+dexp(ein))
           endif
         endif
         slog=slog-sum1
 217    continue
        sLn=sLn-slog
        s2Ln=s2Ln+slog**2
        seqlike(i)=dexp(-slog)
        if (-slog .gt. -80.0d0) then
         sr=sr+dexp(-slog)
         s2r=s2r+dexp(-slog)**2
        endif
 210   continue
       er=sr/real(nrep)
       ser=dsqrt( (s2r-real(nrep)*er**2)
     1       / real(nrep*(nrep-1)) )


       q1=0.005d0*real(nrep)
       q2=0.995d0*real(nrep)
       nq1=nint(q1)
       nq2=nint(q2)
       FUZZ=0.0d0
       ITIE=1
       call DRANKS(nrep,seqlike,FUZZ,ITIE,0,SCORE)
       do 1541 i=1,nrep
        iscore(i)=nint(score(i))
 1541  continue
       ert=0.0d0 
       s2=0.0d0
       num=0
       do 1542 i=1,nrep
        if (iscore(i) .gt. nq2) goto 1542
        if (iscore(i) .le. nq1) goto 1542 
        num=num+1
        ert=ert+seqlike(i)
        s2=s2+seqlike(i)**2
 1542  continue
       ert=ert/real(num)
       sert=dsqrt( (s2-real(num)*ert**2)
     1       / real(num*(num-1)) )
c......obtain the estimates
c......for beta
       do 220 i=1,5
        ebeta(i)=sbeta(i)/real(nrep)
        sebeta(i)=dsqrt( (s2beta(i)-real(nrep)*ebeta(i)**2)
     1                    /real(nrep-1) )
 220   continue
c......to obtian 95% hpd interval of beta
       q1=0.025d0*real(nrep)
       q2=0.975d0*real(nrep)
       nq1=nint(q1)
       nq2=nint(q2)
       nq=nq2-nq1
       do 300 k=1,5
       FUZZ=0.0d0
       ITIE=1
       do 540 j=1,nrep
        seqb(j)=seqbeta(k,j)
 540   continue
       call DRANKS(nrep,seqb,FUZZ,ITIE,0,SCORE)
       do 541 i=1,nrep
        iscore(i)=nint(score(i))
 541   continue
       whpd=1000000.0d0
       do 641 j=1,nrep-nq
        do 542 i=1,nrep
         if (iscore(i) .eq. j) pdiff1=seqb(i)
         if (iscore(i) .eq. j+nq) pdiff2=seqb(i)
 542    continue
        wb=pdiff2-pdiff1
        if (whpd .gt. wb) then
         whpd=wb
         aupp=pdiff2
         alow=pdiff1
        endif
 641   continue
       betaup(k)=aupp
       betalow(k)=alow
 300   continue




       call etime(tarray2)
c
       a(1)=0.75d0
       a(2)=1.5d0
       a(3)=2.5d0
       a(4)=3.5d0
       a6(1)=0.75d0
       a6(2)=1.00d0
       a6(3)=1.25d0
       a6(4)=1.50d0
       a6(5)=2.0d0
       a6(6)=2.5d0
       a6(7)=3.0d0
       a6(8)=3.50d0
       a6(9)=4.0d0

c
       sb(1)=0.0d0
       sb(2)=0.0d0
       sb(3)=0.0d0
       sb(4)=0.0d0
       sb(5)=0.0d0
       sb6(1)=0.0d0
       sb6(2)=0.0d0
       sb6(3)=0.0d0
       sb6(4)=0.0d0
       sb6(5)=0.0d0
       sb6(6)=0.0d0
       sb6(7)=0.0d0
       sb6(8)=0.0d0
       sb6(9)=0.0d0
       sb6(10)=0.0d0
       do 410 i=1,nrep
        slog=seqlike(i)
c
        if (slog .le. a(1)) then
         sb(1)=sb(1)+seqlike(i)**2 
        endif
        if ((slog .gt. a(1)) .and.
     1      (slog .le. a(2)) )  then
         sb(2)=sb(2)+seqlike(i)**2 
        endif
        if ((slog .gt. a(2)) .and.
     1      (slog .le. a(3)) )  then
         sb(3)=sb(3)+seqlike(i)**2 
        endif
        if ((slog .gt. a(3)) .and.
     1      (slog .le. a(4)) )  then
         sb(4)=sb(4)+seqlike(i)**2 
        endif
        if (slog .gt. a(4)) then
         sb(5)=sb(5)+seqlike(i)**2 
        endif
c
        if (slog .le. a6(1)) then
         sb6(1)=sb6(1)+seqlike(i)**2
        endif
        if ((slog .gt. a6(1)) .and.
     1      (slog .le. a6(2)) )  then
         sb6(2)=sb6(2)+seqlike(i)**2
        endif
        if ((slog .gt. a6(2)) .and.
     1      (slog .le. a6(3)) )  then
         sb6(3)=sb6(3)+seqlike(i)**2
        endif
        if ((slog .gt. a6(3)) .and.
     1      (slog .le. a6(4)) )  then
         sb6(4)=sb6(4)+seqlike(i)**2
        endif
        if ((slog .gt. a6(4)) .and.
     1      (slog .le. a6(5)) )  then
         sb6(5)=sb6(5)+seqlike(i)**2
        endif
        if ((slog .gt. a6(5)) .and.
     1      (slog .le. a6(6)) )  then
         sb6(6)=sb6(6)+seqlike(i)**2
        endif
        if ((slog .gt. a6(6)) .and.
     1      (slog .le. a6(7)) )  then
         sb6(7)=sb6(7)+seqlike(i)**2
        endif
        if ((slog .gt. a6(7)) .and.
     1      (slog .le. a6(8)) )  then
         sb6(8)=sb6(8)+seqlike(i)**2
        endif
        if ((slog .gt. a6(8)) .and.
     1      (slog .le. a6(9)) )  then
         sb6(9)=sb6(9)+seqlike(i)**2
        endif
        if (slog .gt. a6(9)) then
         sb6(10)=sb6(10)+seqlike(i)**2
        endif

 410   continue
       do 411 j=1,10
        if (j .le. 5) then
         b(j)=sb(j)/real(nrep)
        endif
        b6(j)=sb6(j)/real(nrep)
 411   continue
        
       do 1200 i=1,500
        call gibbs(iseed)
        if ( (i .eq. 1) .or. (i .eq. icount) ) then
         write(*,*) 'warm up (logit) i=',i
         write(*,888) (beta(j),j=1,7)
         icount=icount+100
        endif
 1200  continue
c......initialize the cumulative variables
       icount=0
       sp(1)=0.0d0
       sp(2)=0.0d0
       sp(3)=0.0d0
       sp(4)=0.0d0
       sp(5)=0.0d0
       sp6(1)=0.0d0
       sp6(2)=0.0d0
       sp6(3)=0.0d0
       sp6(4)=0.0d0
       sp6(5)=0.0d0
       sp6(6)=0.0d0
       sp6(7)=0.0d0
       sp6(8)=0.0d0
       sp6(9)=0.0d0
       sp6(10)=0.0d0
c      do 1210 i=1,50000
c      do 1210 i=1,5000
       do 1210 i=1,nrep
         do 1222 j=1,10
          call gibbs(iseed)
 1222    continue 
        if ( (i .eq. 1) .or. (i .eq. icount) ) then
         write(*,*) 'Gibbs (logit) i=',i
         write(*,888) (beta(j),j=1,7)
         icount=icount+50
        endif
        slog=0.0d0
        do 1217 i1=1,n
         ein=0.0d0
         do 1218 j=1,7
          ein=ein+beta(j)*x(j,i1)
 1218    continue
         if (ein .gt. -70.0d0) then
          ein1=dexp(ein)
         else
          ein1=0.0d0
         endif
         if (y(i1) .eq. 1.0d0) then
          slog=slog+dlog(1.0d0-dexp(-ein1))
         endif
         if (y(i1) .eq. 0.0d0) then
          slog=slog-ein1
         endif
         sum1=0.0d0
         if (ein .gt. 0.0d0) then
           sum1=sum1+ein*(y(i1)-1.0d0)
           if (ein .le. 50.0d0) then
            sum1=sum1-dlog(1.0d0+dexp(-ein))
           endif
         else
           sum1=sum1+ein*y(i1)
           if (ein .gt. -50.0d0) then
            sum1=sum1-dlog(1.0d0+dexp(ein))
           endif
         endif
         slog=slog-sum1
 1217   continue
        ein=dexp(-slog)
        if (ein .le. a(1)) then
         sp(1)=sp(1)+1.0d0
        endif
        if ((ein .gt. a(1)) .and.
     1      (ein .le. a(2)) )  then
         sp(2)=sp(2)+1.0d0
        endif
        if ((ein .gt. a(2)) .and.
     1      (ein .le. a(3)) )  then
         sp(3)=sp(3)+1.0d0
        endif
        if ((ein .gt. a(3)) .and.
     1      (ein .le. a(4)) )  then
         sp(4)=sp(4)+1.0d0
        endif
        if (ein .gt. a(4)) then
         sp(5)=sp(5)+1.0d0
        endif
c
        if (ein .le. a6(1)) then
         sp6(1)=sp6(1)+1.0d0
        endif
        if ((ein .gt. a6(1)) .and.
     1      (ein .le. a6(2)) )  then
         sp6(2)=sp6(2)+1.0d0
        endif
        if ((ein .gt. a6(2)) .and.
     1      (ein .le. a6(3)) )  then
         sp6(3)=sp6(3)+1.0d0
        endif
        if ((ein .gt. a6(3)) .and.
     1      (ein .le. a6(4)) )  then
         sp6(4)=sp6(4)+1.0d0
        endif
        if ((ein .gt. a6(4)) .and.
     1      (ein .le. a6(5)) )  then
         sp6(5)=sp6(5)+1.0d0
        endif
        if ((ein .gt. a6(5)) .and.
     1      (ein .le. a6(6)) )  then
         sp6(6)=sp6(6)+1.0d0
        endif
        if ((ein .gt. a6(6)) .and.
     1      (ein .le. a6(7)) )  then
         sp6(7)=sp6(7)+1.0d0
        endif
        if ((ein .gt. a6(7)) .and.
     1      (ein .le. a6(8)) )  then
         sp6(8)=sp6(8)+1.0d0
        endif
        if ((ein .gt. a6(8)) .and.
     1      (ein .le. a6(9)) )  then
         sp6(9)=sp6(9)+1.0d0
        endif
        if (ein .gt. a6(9)) then
         sp6(10)=sp6(10)+1.0d0
        endif
 1210  continue
       ss=0.0d0
       ss6=0.0d0
       do 1250 j=1,10
        if (j .le. 5) then
c        p(j)=sp(j)/5000.0d0
         p(j)=sp(j)/real(nrep)
         ss=ss+p(j)**2/b(j)
        endif
c       p6(j)=sp6(j)/5000.0d0
        p6(j)=sp6(j)/real(nrep)
        ss6=ss6+p6(j)**2/b6(j)
 1250  continue
       do 1251 j=1,10
        if (j .le. 5) then
         w(j)=(p(j)/b(j))/ss
        endif
        w6(j)=(p6(j)/b6(j))/ss6
 1251  continue

       swr=0.0d0
       sw6r=0.0d0
       do 1260 i=1,nrep
        slog=seqlike(i)
        if (slog .le. a(1)) then
         ng=1 
        endif
        if ((slog .gt. a(1)) .and.
     1      (slog .le. a(2)) )  then
         ng=2 
        endif
        if ((slog .gt. a(2)) .and.
     1      (slog .le. a(3)) )  then
         ng=3 
        endif
        if ((slog .gt. a(3)) .and.
     1      (slog .le. a(4)) )  then
         ng=4 
        endif
        if (slog .gt. a(4)) then
         ng=5 
        endif
c
        if (slog .le. a6(1)) then
         ng6=1
        endif
        if ((slog .gt. a6(1)) .and.
     1      (slog .le. a6(2)) )  then
         ng6=2
        endif
        if ((slog .gt. a6(2)) .and.
     1      (slog .le. a6(3)) )  then
         ng6=3
        endif
        if ((slog .gt. a6(3)) .and.
     1      (slog .le. a6(4)) )  then
         ng6=4
        endif
        if ((slog .gt. a6(4)) .and.
     1      (slog .le. a6(5)) )  then
         ng6=5
        endif
        if ((slog .gt. a6(5)) .and.
     1      (slog .le. a6(6)) )  then
         ng6=6
        endif
        if ((slog .gt. a6(6)) .and.
     1      (slog .le. a6(7)) )  then
         ng6=7
        endif
        if ((slog .gt. a6(7)) .and.
     1      (slog .le. a6(8)) )  then
         ng6=8
        endif
        if ((slog .gt. a6(8)) .and.
     1      (slog .le. a6(9)) )  then
         ng6=9
        endif
        if (slog .gt. a6(9)) then
         ng6=10
        endif
c
        swr=swr+seqlike(i)*w(ng)
        sw6r=sw6r+seqlike(i)*w6(ng6)
 1260  continue
       ewr=swr/real(nrep)
       varw=(1.0d0/ss - ewr**2)
       ew6r=sw6r/real(nrep)
       varw6=(1.0d0/ss6 - ew6r**2)

       

       call etime(tarray3)
       open(unit=12,file='log.out50000a',access='sequential',
     1              status='unknown')
       write(12,*) 'loglog regression model'
       write(12,*) 'x1---T4COUNT'
       write(12,*) 'x2---AGE'
       write(12,*) 'x3---TREAT'
       write(12,*) 'Gibbs iterates for post and prior nrep=',nrep
       write(12,*) 'MLE of beta for current data'
       write(12,1996) (amuc(j),j=1,7)
 1996  format(1x, 7f16.6)
       write(12,*) '-------------------------------------------'
       write(12,*) 'posterior estimate of beta for post full model'
       write(12,*) 'index        ebeta          sebeta  95% hpd'
       do 900 i=1,5
         write(12,1998) i-1, ebeta(i), sebeta(i),
     1          betalow(i),betaup(i)
 1998    format(1x,I5, 4f16.6)
 900   continue
       write(12,*) 'runtime 1=', (tarray1(i), i=1,2)
       write(12,*) 'runtime 2=', (tarray2(i), i=1,2)
       write(12,*) 'runtime 3=', (tarray3(i), i=1,2)
       write(12,*) 'er=',er,' ser=',ser 
       write(12,*) 'var=', (ser**2)*real(nrep) 
       write(12,*) 'ert=',ert,' sert=',sert 
       write(12,*) 'var=', (sert**2)*real(num),' num=',num 
       write(12,*) 'ewr=',ewr,' sewr=', dsqrt(varw/real(nrep))
       write(12,*) 'varw=', varw 
       write(12,*) 'ew6r=',ew6r,' sew6r=', dsqrt(varw6/real(nrep))
       write(12,*) 'varw6=', varw6 
       write(12,*) 'eLn=',eLn,' seLn=',seLn 
       write(12,*) 'b6     p6     b     p'
       do 910 j=1,10
        if (j .le. 5) then
         write(12,920) j, b6(j),p6(j),b(j),p(j)
 920     format(I6,4f16.6)
        else
         write(12,930) j, b6(j),p6(j)
 930     format(I6,2f16.6)
        endif
 910   continue
       write(12,*) ' j,    a6'   
       do 940 j=1,10
        write(12,950) j,a6(j)
 950    format(I6,1x,f16.6)
 940   continue
       close(12)
       open(unit=14,file='log50000a.dat',access='sequential',
     1              status='unknown')
       do 970 i=1,nrep
                slog=seqlike(i)
        if (slog .le. a(1)) then
         ng=1
        endif
        if ((slog .gt. a(1)) .and.
     1      (slog .le. a(2)) )  then
         ng=2
        endif
        if ((slog .gt. a(2)) .and.
     1      (slog .le. a(3)) )  then
         ng=3
        endif
        if ((slog .gt. a(3)) .and.
     1      (slog .le. a(4)) )  then
         ng=4
        endif
        if (slog .gt. a(4)) then
         ng=5
        endif
c
        if (slog .le. a6(1)) then
         ng6=1
        endif
        if ((slog .gt. a6(1)) .and.
     1      (slog .le. a6(2)) )  then
         ng6=2
        endif
        if ((slog .gt. a6(2)) .and.
     1      (slog .le. a6(3)) )  then
         ng6=3
        endif
        if ((slog .gt. a6(3)) .and.
     1      (slog .le. a6(4)) )  then
         ng6=4
        endif
        if ((slog .gt. a6(4)) .and.
     1      (slog .le. a6(5)) )  then
         ng6=5
        endif
        if ((slog .gt. a6(5)) .and.
     1      (slog .le. a6(6)) )  then
         ng6=6
        endif
        if ((slog .gt. a6(6)) .and.
     1      (slog .le. a6(7)) )  then
         ng6=7
        endif
        if ((slog .gt. a6(7)) .and.
     1      (slog .le. a6(8)) )  then
         ng6=8
        endif
        if ((slog .gt. a6(8)) .and.
     1      (slog .le. a6(9)) )  then
         ng6=9
        endif
        if (slog .gt. a6(9)) then
         ng6=10
        endif
         write(14,*) seqlike(i),w(ng)*seqlike(i),w6(ng6)*seqlike(i)
 970   continue
       close(14)

       stop
       end

c......including Gibbs sampling codes
       include 'gibbslog.f'
       include 'gibbslogit.f'
