       subroutine gibbs( iseed )
c      ming-hui chen, wpi
c      march 6, 1998
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       integer 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
c......generate nbig
       call Gnbig( iseed )
c      write(*,*) 'after Gnbig'
c......generate nobig
       call Gnobig( iseed )
c      write(*,*) 'after Gnobig'
c......generate beta
       call Gbeta( iseed )
c      write(*,*) 'after Gbeta'
c......generate alam
       call Galam( iseed )
c      write(*,*) 'after Galam'
c......generate alpha
       call Galpha( iseed )
c      write(*,*) 'after Galpha'
c......generate ao
       call Gao( iseed )
c      write(*,*) 'after ao'
       return
       end

       subroutine Gnbig( iseed )
c      ming-hui chen, wpi
c      march 6, 1998
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       integer 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 DRNUNF
       do 100 i=1,n
c.......START of Generation of a Poisson RV
        ein=0.0d0
        do 101 j=1,4
         ein=ein+x(j,i)*beta(j)
 101    continue
        theta=dexp(ein)
        ss=y(i)**alpha*dexp(alam)
        amu=theta*dexp(-ss)
        np=0
        sum=0.0d0
 33     call rnset( iseed )
         u=DRNUNF()
        call rnget( iseed )
        rne=-dlog(u)
        sum=sum+rne
        if (sum .lt. amu) then
         np=np+1
         goto 33
        endif
        nbig(i)=np+nint(d(i))
c       write(*,*) 'i=',i,' nbig=',nbig(i)
c.......END of Generation of a Poisson RV
 100   continue
c      write(*,*) 'nbig(1)=',nbig(1),' nbig(n)=',nbig(n)
       return
       end

       subroutine Gnobig( iseed )
c      ming-hui chen, wpi
c      march 6, 1998
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       integer 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 DRNUNF
       do 100 i=1,no
c.......START of Generation of a Poisson RV
        ein=0.0d0
        do 101 j=1,4
         ein=ein+xo(j,i)*beta(j)
 101    continue
        theta=dexp(ein)
        ss=yo(i)**alpha*dexp(alam)
        amu=ao*theta*dexp(-ss)
        np=0
        sum=0.0d0
 33     call rnset( iseed )
         u=DRNUNF()
        call rnget( iseed )
        rne=-dlog(u)
        sum=sum+rne
        if (sum .lt. amu) then
         np=np+1
         goto 33
        endif
        nobig(i)=np
c.......END of Generation of a Poisson RV
c       write(*,*) 'i=',i,' nobig=',nobig(i)
 100   continue
c      write(*,*) 'nobig(1)=',nobig(1),' nobig(n)=',nobig(no)
       return
       end

       subroutine Gbeta( iseed )
c      ming-hui chen, wpi
c      march 6, 1998
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
c      use Gilks log-concave sampling
C *************************** MODEL PARAMETERS *****************************
C -------------------------- Commonly used variables ---------------------
       PARAMETER(maxtries=25,im=2,NS=10,np=7,in=1000)
       real*8  emax, step
C ------------------------ Working variables ------------------------------
       real*8  val,a(im+in),ha(im+in),hpa(im+in),
     1         xlb,xub,rwv(6*NS+15+in)
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       integer nbig(n),nobig(no)
       integer ifault,iwv(NS+7+in),ifcount
       integer iseed,idum
       logical lb,ub
       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
       common /dummy/idum
       external eval
       emax=60.0d0
       DO 80 j=1,4
        idum=j
        step=1.0d0
        a(1)=beta(j)-step
        a(2)=beta(j)+step
20      CALL eval(a(1), ha(1), hpa(1))
        CALL eval(a(2), ha(2), hpa(2))
c        write(*,*) 'a(1)=',a(1),' ha(1)=',ha(1)
c        write(*,*) ' hpa(1)=',hpa(1)
c ------------------ Get starting points for Gilks -----------------------
         IF (hpa(2) .gt. 0.0) THEN
c        Both points to the left of the mode; push a(2) to the right
30          a(2)=a(2)+step
            CALL EVAL(a(2), ha(2), hpa(2))
c        If a(2) still to the left of the mode repeat, else continue
            IF (hpa(2).gt.0.0) THEN
               GOTO 30
            ELSE
               a(1)=a(2)-step
            ENDIF
         ELSE
            IF (hpa(1) .le. 0.0) THEN
c        Both points to the right of the mode
c        Insert new point to the left of a(1)
               a(2)=a(1)
               ha(2)=ha(1)
               hpa(2)=hpa(1)
40             a(1)=a(1)-step
               CALL EVAL(a(1), ha(1), hpa(1))
               IF (hpa(1) .le. 0.0) THEN
                  GOTO 40
c              ELSE
c                 a(2)=a(1)+step
               ENDIF
             ENDIF
         ENDIF
c -----------------------------------------------------------------------
         ifcount=0
50       CALL INITIAL(ns, im, emax, a, ha, hpa, lb, xlb, ub, xub,
     +                ifault, iwv, rwv, eval)
         IF((ifault.eq.3) .or. (ifault.eq.4)) GOTO 20
C
         CALL SAMPLE(iwv,rwv,eval,val,ifault,iseed)
         IF(ifault.eq.5) THEN
c          PRINT*, 'ifault =', ifault, ' trying again'
           a(1)=a(1)-step
           a(2)=a(2)+step
           ifcount=ifcount+1
           IF(ifcount.le.maxtries) THEN
             GOTO 50
           ELSE
c             PRINT*, 'Too many tries resetting starting values'
              a(1)=0.0
              a(2)=0.0
              step=1.2d0*step
              GOTO 20
           ENDIF
         ENDIF
         IF(ifault.eq.7) THEN
           a(1)=a(1)-step
           a(2)=a(2)+step
           ifcount=ifcount+1
           IF(ifcount.le.maxtries) THEN
             GOTO 50
           ELSE
c            PRINT*, 'Too many tries resetting starting values'
             a(1)=0.0
             a(2)=0.0
             step=1.2d0*step
             GOTO 20
           ENDIF
         ENDIF
         IF (IFAULT .NE. 0) THEN
c           PRINT*, ' The sampling status is ', IFAULT
         ENDIF
C
         beta(j)=val
C
 80    continue
       return
       end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                             SUBROUTINE EVAL(x, hx, hpx)                   C
C Arguments:                                                                C
C       z                     Input.           The point of evaluation      C
C       hz                    Output           hz=h(z)                      C
C       hpz                   Output           hpz=h'(z)                    C
C                                                                           C
C       Subroutine EVAL(z,hz,hpz) evaluates the log-likelihood function     C
C       at a predetermined point Z (input), such that hz=h(z) (output).     C
C       The only information that is COMMON to it and the outside world     C
C       is J, the element of the vector of parameters BETA, with respect    C
C       of which the derivative h'(.) is taken.                             C
C       This derivative is also evaluated at a, i.e., hpz=h'(z) (output).   C
C       EVAL(.) is declared EXTERNAL in Gilks's adaptive-rejection method   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       subroutine eval(z1,ahz,ahpz)
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       real*8 z1,ahz,ahpz
       integer 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
       common /dummy/idum
       beta(idum)=z1
c......calculate log likelihood and its derivative
       sum1=0.0d0
       sum2=0.0d0
       do 300 i=1,n
        ein=0.0d0
        do 301 j=1,4
         ein=ein+x(j,i)*beta(j)
 301    continue
        sum1=sum1+ein*d(i)
        sum2=sum2+x(idum,i)*d(i)
        ss=y(i)**alpha*dexp(alam)
        ein1=ein
        if (ss .le. 80.0d0) then
         ein1=ein1+dlog(1.0d0-dexp(-ss))
        endif
        if (ein1 .gt. -80.0d0) then
          sum1=sum1-dexp(ein1)
          sum2=sum2-x(idum,i)*dexp(ein1)
        endif
 300   continue
c......from historical data
       sum3=0.0d0
       sum4=0.0d0
       if (ao .eq. 0.0d0) goto 322
       do 310 i=1,no
        ein=0.0d0
        do 311 j=1,4
         ein=ein+xo(j,i)*beta(j)
 311    continue
        sum3=sum3+ein*do(i)
        sum4=sum4+xo(idum,i)*do(i)
        ss=yo(i)**alpha*dexp(alam)
        ein1=ein
        if (ss .le. 80.0d0) then
         ein1=ein1+dlog(1.0d0-dexp(-ss))
        endif
        if (ein1 .gt. -80.0d0) then
          sum3=sum3-dexp(ein1)
          sum4=sum4-xo(idum,i)*dexp(ein1)
        endif
 310   continue
c322   ein1=-0.0001*beta(idum)**2/2.0d0
c      ein2=-0.0001*beta(idum)
 322   sum3=sum3*ao
       sum4=sum4*ao
       ahz=sum1+sum3
       ahpz=sum2+sum4
       return
       end

       include 'gilks1.f'

       subroutine Galam( iseed )
c      ming-hui chen, wpi
c      march 6, 1998
       implicit real*8 (a-h,o-z)
c      use Gilks log-concave sampling
C *************************** MODEL PARAMETERS *****************************
C -------------------------- Commonly used variables ---------------------
       PARAMETER(maxtries=25,im=2,NS=10,np=7,in=1000)
       real*8  emax, step
C ------------------------ Working variables ------------------------------
       real*8  val,a(im+in),ha(im+in),hpa(im+in),
     1         xlb,xub,rwv(6*NS+15+in)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       integer ifault,iwv(NS+7+in),ifcount
       integer 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 evala1
       emax=60.0d0
        step=5.0d0
        a(1)=alam-step
        a(2)=alam+step
20      CALL evala1(a(1), ha(1), hpa(1))
        CALL evala1(a(2), ha(2), hpa(2))
c        write(*,*) 'a(1)=',a(1),' ha(1)=',ha(1)
c        write(*,*) ' hpa(1)=',hpa(1)
c ------------------ Get starting points for Gilks -----------------------
         IF (hpa(2) .gt. 0.0) THEN
c        Both points to the left of the mode; push a(2) to the right
30          a(2)=a(2)+step
            CALL EVALA1(a(2), ha(2), hpa(2))
c        If a(2) still to the left of the mode repeat, else continue
            IF (hpa(2).gt.0.0) THEN
               GOTO 30
            ELSE
               a(1)=a(2)-step
            ENDIF
         ELSE
            IF (hpa(1) .le. 0.0) THEN
c        Both points to the right of the mode
c        Insert new point to the left of a(1)
               a(2)=a(1)
               ha(2)=ha(1)
               hpa(2)=hpa(1)
40             a(1)=a(1)-step
               CALL EVALA1(a(1), ha(1), hpa(1))
               IF (hpa(1) .le. 0.0) THEN
                  GOTO 40
c              ELSE
c                 a(2)=a(1)+step
               ENDIF
             ENDIF
         ENDIF
c -----------------------------------------------------------------------
         ifcount=0
50       CALL INITIAL(ns, im, emax, a, ha, hpa, lb, xlb, ub, xub,
     +                ifault, iwv, rwv, evala1)
         IF((ifault.eq.3) .or. (ifault.eq.4)) GOTO 20
C
         CALL SAMPLE(iwv,rwv,evala1,val,ifault,iseed)
         IF(ifault.eq.5) THEN
c          PRINT*, 'ifault =', ifault, ' trying again'
           a(1)=a(1)-step
           a(2)=a(2)+step
           ifcount=ifcount+1
           IF(ifcount.le.maxtries) THEN
             GOTO 50
           ELSE
c             PRINT*, 'Too many tries resetting starting values'
              a(1)=0.0
              a(2)=0.0
              step=1.2d0*step
              GOTO 20
           ENDIF
         ENDIF
         IF(ifault.eq.7) THEN
           a(1)=a(1)-step
           a(2)=a(2)+step
           ifcount=ifcount+1
           IF(ifcount.le.maxtries) THEN
             GOTO 50
           ELSE
c            PRINT*, 'Too many tries resetting starting values'
             a(1)=0.0
             a(2)=0.0
             step=1.2d0*step
             GOTO 20
           ENDIF
         ENDIF
         IF (IFAULT .NE. 0) THEN
c           PRINT*, ' The sampling status is ', IFAULT
         ENDIF
C
         alam=val
C
       return
       end

       subroutine evala1(z1,ahz,ahpz)
c      ming-hui chen, wpi
c      march 6, 1998
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       real*8 z1,ahz,ahpz
       integer 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
       alam=z1
c......calculate log likelihood and its derivative
       sum1=0.0d0
       sum2=0.0d0
       do 300 i=1,n
         ein=real(nbig(i))-d(i)
         sum1=sum1-y(i)**alpha*dexp(alam)*ein
         sum2=sum2-y(i)**alpha*dexp(alam)*ein
         sum1=sum1+(alam-y(i)**alpha*dexp(alam))*d(i)
         sum2=sum2+(1.0d0-y(i)**alpha*dexp(alam))*d(i)
 300   continue
c......from historical data
       sum3=0.0d0
       sum4=0.0d0
       if (ao .eq. 0.0d0) goto 322
       do 330 i=1,no
        ein=real(nobig(i))
        sum3=sum3-yo(i)**alpha*dexp(alam)*ein
        sum4=sum4-yo(i)**alpha*dexp(alam)*ein
        sum3=sum3+(alam-yo(i)**alpha*dexp(alam))*do(i)*ao
        sum4=sum4+(1.0d0-yo(i)**alpha*dexp(alam))*do(i)*ao
 330   continue
c322   ein1=-0.0001*alam**2/2.0d0
c      ein2=-0.0001*alam
 322   ein1=-0.1*alam**2/2.0d0
       ein2=-0.1*alam
c322   ein1=-1.0*alam**2/2.0d0
c      ein2=-1.0*alam
       ahz=sum1+sum3+ein1
       ahpz=sum2+sum4+ein2
       return
       end


       subroutine Galpha( iseed )
c      ming-hui chen, wpi
c      march 6, 1998
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       real*8 ur,ur1,x1,p
       real*8 gam
       real*8 s(3),u(0:2)
       real*8 vh(3),devh(3)
       real*8 v(6),pv(5),temp(11),tmax
       real*8 xopt1,xopt2,amax,xopt
       real*8 start(1),xmin(1),xsec(1),step(1)
       real*8 ynewlo,ysec,reqmin
       integer iseed,nopt
       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 fnalpha,den,dev,ahu,ahl,DRNUNF
c......to find xopt
       xiold=dlog(alpha)
       nopt=1
       konvge=5
       kcount=1000
       reqmin=1.0d-10
       start(1)=xiold
       do 120 k1=1,1
        step(k1)=0.2d0
        xmin(k1)=start(k1)
        xsec(k1)=start(k1)
 120   continue
c......call optim.f
       call nelmin(fnalpha, nopt, start, xmin, ynewlo, reqmin, step,
     1     konvge, kcount, icount, numres, ifault)
       xopt=xmin(1)
       u(0)=0.0d0
       s(2)=xopt
       s(1)=s(2)/2.0d0
       s(3)=3.0d0*s(2)
       amax=den(xopt)
 433   vh(1)=den(s(1))-amax
       devh(1)=dev(s(1))
       vh(3)=den(s(3))-amax
       devh(3)=dev(s(3))
       vh(2)=den(s(2))-amax
       devh(2)=dev(s(2))
       amax1=dev(xopt)
       do 420 k=1,2
        u(k)=( vh(k+1)-vh(k)-s(k+1)*devh(k+1)+s(k)*devh(k) )
     1       / ( devh(k)-devh(k+1) )
 420   continue
c      write(*,*) 's=',s(1),s(2),s(3)
c      write(*,*) 'u=',u(0),u(1),u(2)
c      write(*,*) 'vh=',vh(1),vh(2),vh(3)
c      write(*,*) 'devh=',devh(1),devh(2),devh(3)
c......calculate five probabilities over subintervals:
c...... (0,s_1), (s_1,u_1), (u_1,s_2), (s_2,u_2), (u_2,s_3),(s_3,infinity)
       temp(1)=vh(1)
       temp(2)=vh(1)+devh(1)*(-s(1))
       temp(3)=vh(1)+devh(1)*(u(1)-s(1))
       temp(4)=vh(1)
c
       temp(5)=vh(2)
       temp(6)=vh(2)+devh(2)*(u(1)-s(2))
       temp(7)=vh(2)+devh(2)*(u(2)-s(2))
       temp(8)=vh(2)
c
       temp(9)=vh(3)
       temp(10)=vh(3)+devh(3)*(u(2)-s(3))
       temp(11)=vh(3)
       tmax=0.0d0
       do 432 i=1,11
        if (temp(i) .gt. tmax) tmax=temp(i)
 432   continue
       v(1)=(dexp(temp(1)-tmax)-dexp(temp(2)-tmax))/devh(1)
       v(2)=(dexp(temp(3)-tmax)-dexp(temp(4)-tmax))/devh(1)
       if (devh(2) .eq. 0.0d0) then
        v(3)=dexp(temp(5)-tmax)*(s(2)-u(1))
        v(4)=dexp(temp(5)-tmax)*(u(2)-s(2))
       else
        v(3)=(dexp(temp(5)-tmax)-dexp(temp(6)-tmax))/devh(2)
        v(4)=(dexp(temp(7)-tmax)-dexp(temp(8)-tmax))/devh(2)
       endif
       v(5)=(dexp(temp(9)-tmax)-dexp(temp(10)-tmax))/devh(3)
       v(6)=-dexp(temp(11)-tmax)/devh(3)
       pv(1)=v(1)/(v(1)+v(2)+v(3)+v(4)+v(5)+v(6))
       pv(2)=(v(1)+v(2))/(v(1)+v(2)+v(3)+v(4)+v(5)+v(6))
       pv(3)=(v(1)+v(2)+v(3))/(v(1)+v(2)+v(3)+v(4)+v(5)+v(6))
       pv(4)=(v(1)+v(2)+v(3)+v(4))/(v(1)+v(2)+v(3)+v(4)+v(5)+v(6))
       pv(5)=1.0d0
c      write(*,*) 'v=',(v(i2),i2=1,6)
c      write(*,*) 'pv=',(pv(i2),i2=1,5)
       if ((xopt .eq. 0.0d0) .and. (pv(2) .gt. 0.7)) then
         s(1)=s(1)/2.0d0
         s(2)=s(2)/2.0d0
         s(3)=s(3)/2.0d0
         if ((pv(2) .gt. 0.99d0) .and. (s(2) .le. 0.1d-4)) then
          gam=s(2)
          goto 450
         endif
         goto 433
       endif
       if ((xopt .eq. 0.0d0) .and. (pv(2) .le. 0.1d-5)) then
         s(1)=s(1)*1.5d0
         s(2)=s(2)*1.5d0
         s(3)=s(3)*1.5d0
         goto 433
       endif
c      if (pv(5)-pv(4) .lt. 0.1d-6) then
c        s(3)=(s(3)+u(2))/2.0d0
c        goto 433
c      endif
  430  ur1=DRNUNF()
       ur=DRNUNF()
       if (ur1 .le. pv(2)) then
        if (devh(1) .gt. 0.0d0) then
         ein6=-devh(1)*u(1)
         if (ein6 .gt. -70.0d0) then
          ein7=dexp(ein6)
         else
          ein7=0.0d0
         endif
         x1=u(1)+dlog(ur+(1.0d0-ur)*ein7)/devh(1)
        else
         ein6=devh(1)*u(1)
         if (ein6 .gt. -70.0d0) then
          ein7=dexp(ein6)
         else
          ein7=0.0d0
         endif
         x1=dlog(ur*ein7+(1.0d0-ur))/devh(1)
        endif
       endif
       if ((ur1 .le. pv(4)) .and. (ur1 .gt. pv(2))) then
        if ( devh(2) .gt. 0.0d0 ) then
         x1= u(2)+(dlog(ur+(1.0d0-ur)*dexp((u(1)-u(2))*devh(2)))
     1               / devh(2) )
        else
         x1= u(1)+(dlog(ur*dexp((u(2)-u(1))*devh(2))+1.0d0-ur)
     1               / devh(2) )
        endif
       endif
       if (ur1 .gt. pv(4)) then
        x1=u(2)+dlog(1.0d0-ur)/devh(3)
       endif
       ur=DRNUNF()
       if ( (s(1) .le. x1) .and. (x1 .le. s(3)) ) then
         p= ahl(vh,devh,s,u,x1)-ahu(vh,devh,s,u,x1)
         if ( dlog(ur) .le. p) then
           gam=x1
           goto 450
         endif
       endif
       p= den(x1)-ahu(vh,devh,s,u,x1)-amax
       if (dlog(ur) .le. p) then
           gam=x1
           goto 450
       endif
       goto 430
 450   alpha=gam
       return
       end

       real*8 function fnalpha(gam)
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       real*8 gam 
       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 /vnbig/nbig
       common /vnobig/nobig
c      sum=-0.01d0*gam
c      sum=-0.1d0*gam
       sum=-1.0d0*gam
       do 100 i=1,n
        sum=sum-real(nbig(i))*y(i)**gam*dexp(alam)
     1     +d(i)*( dlog(gam)+(gam-1.0d0)*dlog(y(i)) )
 100   continue
       if (ao .gt. 0.0d0) then
        do 200 i=1,no
         sum=sum-(real(nobig(i))+ao*do(i))*yo(i)**gam*dexp(alam)
     1     +ao*do(i)*( dlog(gam)+(gam-1.0d0)*dlog(yo(i)) )
 200    continue
       endif
       fnalpha=-sum
       return
       end

       real*8 function den(gam)
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       real*8 gam 
       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 /vnbig/nbig
       common /vnobig/nobig
c      sum=-0.01*gam
c      sum=-0.1*gam
       sum=-1.0*gam
       do 100 i=1,n
        sum=sum-real(nbig(i))*y(i)**gam*dexp(alam)
     1     +d(i)*( dlog(gam)+(gam-1.0d0)*dlog(y(i)) )
 100   continue
       if (ao .gt. 0.0d0) then
        do 200 i=1,no
         sum=sum-(real(nobig(i))+ao*do(i))*yo(i)**gam*dexp(alam)
     1     +ao*do(i)*( dlog(gam)+(gam-1.0d0)*dlog(yo(i)) )
 200    continue
       endif
       den=sum
       return
       end

       real*8 function dev(gam)
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       real*8 gam 
       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 /vnbig/nbig
       common /vnobig/nobig
c      sum=-0.01d0
c      sum=-0.1d0
       sum=-1.0d0
       do 100 i=1,n
        sum=sum-real(nbig(i))*y(i)**gam*dexp(alam)*dlog(y(i))
     1     +d(i)*( 1.0d0/gam+dlog(y(i)) )
 100   continue
       if (ao .gt. 0.0d0) then
        do 200 i=1,no
         sum=sum-(real(nobig(i))
     1     +ao*do(i))*yo(i)**gam*dexp(alam)*dlog(yo(i))
     1     +ao*do(i)*( 1.0d0/gam+dlog(yo(i)) )
 200    continue
       endif
       dev=sum
       return
       end
  
       real*8 function ahu(vh,devh,s,u,gam)
       implicit real*8 (a-h,o-z)
       real*8 vh(3),devh(3),gam
       real*8 s(3),u(0:2)
       if (gam .le. u(1)) then
         ahu=vh(1)+(gam-s(1))*devh(1)
       endif
       if ( (u(1) .lt. gam) .and. (gam .le. u(2)) ) then
         ahu=vh(2)+(gam-s(2))*devh(2)
       endif
       if (gam .gt. u(2)) then
         ahu=vh(3)+(gam-s(3))*devh(3)
       endif
       return
       end

       real*8 function ahl(vh,devh,s,u,gam)
       implicit real*8 (a-h,o-z)
       real*8 vh(3),devh(3),gam
       real*8 s(3),u(0:2)
       if ( (s(1) .le. gam) .and. (gam .le. s(2)) ) then
        ahl=( (vh(2)-vh(1))*gam+s(2)*vh(1)-s(1)*vh(2) )
     1        /(s(2)-s(1))
       endif
       if ( (s(2) .lt. gam) .and. (gam .le. s(3)) ) then
        ahl=( (vh(3)-vh(2))*gam+s(3)*vh(2)-s(2)*vh(3) )
     1        /(s(3)-s(2))
       endif
       return
       end
      

       subroutine Gao( iseed )
c      ming-hui chen, wpi
c      march 6, 1998
       implicit real*8 (a-h,o-z)
       parameter( n=284, no=650 )
       real*8 y(n),d(n),x(4,n)
       real*8 yo(no),do(no),xo(4,no)
       real*8 beta(4),alam,alpha
       real*8 ao,alamo,deltao
       real*8 xinew,xiold
       real*8 start(1),xmin(1),xsec(1),step(1)
       real*8 ynewlo,ysec,reqmin
       real*8 ximean,xistd,aLn
       integer konvge,kcount
       integer nopt
       integer 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
       common /vaLn/aLn
       common /vadim/adim
       external DRNUNF,fnao,DRNNOF
       adim=0.0d0
       sum1=0.0d0
       do 301 i=1,no
        adim=adim+real(nobig(i))
        ein=0.0d0
        do 302 j=1,4
         ein=ein+xo(j,i)*beta(j)
 302    continue
        if (ein .gt. -80.0d0) then
         sum1=sum1-dexp(ein)
        endif
         sum1=sum1+do(i)*( ein+dlog(alpha)
     1           +(alpha-1.0d0)*dlog(yo(i))
     1           +alam-yo(i)**alpha*dexp(alam) )
 301   continue
       aLn=sum1
       xiold=dlog(ao/(1.0d0-ao))
       nopt=1
       konvge=5
       kcount=1000
       reqmin=1.0d-10
       start(1)=xiold
       do 420 k1=1,1
        step(k1)=0.2d0
        xmin(k1)=start(k1)
        xsec(k1)=start(k1)
 420   continue
c......call optim.f
       call nelmin(fnao, nopt, start, xmin, ynewlo, reqmin, step,
     1     konvge, kcount, icount, numres, ifault)
       ximean=xmin(1)
c......calculate Hessian matrix
       if (ximean .gt. 0.0d0) then
        d1=dexp(-ximean)/(1.0d0+dexp(-ximean))**2
        d11=-dexp(-2.0d0*ximean)*(1.0d0-dexp(-ximean))
     1        /(1.0d0+dexp(-ximean))**3
        gam=1.0d0/(1.0d0+dexp(-ximean))
       else
        d1=dexp(ximean)/(1.0d0+dexp(ximean))**2
        d11=-dexp(ximean)*(dexp(ximean)-1.0d0)
     1      /(1.0d0+dexp(ximean))**3
        gam=dexp(ximean)/(1.0d0+dexp(ximean))
       endif
       sum1=d11*(aLn+(adim+deltao-1.0d0)/gam)
     1       -(adim+deltao-1.0d0)*d1**2/gam**2
       sum2=-(alamo-1.0d0)*(d11+d1**2/(1.0d0-gam))/(1.0d0-gam)
       sum3=-2.0d0*d1
       sig2inv=-(sum1+sum2+sum3)
       if (sig2inv .lt. 0.0d0) then
         write(*,*) 'error sig2inv is nagtive'
         sig2inv=-sig2inv
       endif
       std=dsqrt(1.0d0/sig2inv)
       sumold=ao*aLn+(adim+deltao-1.0d0)*dlog(ao)
     1    +(alamo-1.0d0)*dlog(1.0d0-ao)+xiold
     1    -2.0d0*dlog(1.0d0+dexp(xiold))
     1  +(xiold-ximean)**2*sig2inv/2.0d0
c.......perform metropolis
        do 440 imetr=1,10
         call rnset( iseed )
          rv=DRNNOF()
         call rnget( iseed )
         xinew=ximean+std*rv
         if (xinew .gt. 0.0d0) then
          gam =1.0d0/(1.0d0+dexp(-xinew))
         else
          gam=dexp(xinew)/(1.0d0+dexp(xinew))
         endif
         sumnew=gam*aLn+(adim+deltao-1.0d0)*dlog(gam)
     1      +(alamo-1.0d0)*dlog(1.0d0-gam)+xinew
     1     -2.0d0*dlog(1.0d0+dexp(xinew))
     1    +(xinew-ximean)**2*sig2inv/2.0d0
         ratio=sumnew-sumold
         if (ratio .ge. 0.0d0) then
          xiold=xinew
          ao=gam
          sumold=sumnew
          ip=ip+1
         else
          call rnset( iseed )
           u=DRNUNF()
          call rnget( iseed )
          if (dlog(u) .le. ratio) then
           xiold=xinew
           sumold=sumnew
           ao=gam
           ip=ip+1
          endif
         endif
 440    continue
       return
       end

       real*8 function fnao(v)
       implicit real*8 (a-h,o-z)
       real*8 v(1)
       real*8 delta,alam,aLn
       COMMON /vdeltao/deltao
       COMMON /valamo/alamo
       common /vaLn/aLn
       common /vadim/adim
       if (v(1) .gt. 0.0d0) then
          gam =1.0d0/(1.0d0+dexp(-v(1)))
       else
          gam=dexp(v(1))/(1.0d0+dexp(v(1)))
       endif
       sum=adim*dlog(gam)+gam*aLn+(deltao-1.0d0)*dlog(gam)
     1       +(alamo-1.0d0)*dlog(1.0d0-gam)+v(1)
     1       -2.0d0*dlog(1.0d0+dexp(v(1)))
       fnao=-sum
       return
       end
