subroutine gibbs(iseed)
c      gibbs sampling for logistic regression model
c      ming-hui chen, 5/24/99, at wpi
       implicit real*8 (a-h,o-z)
       parameter (n=183)
       real*8  y(n),x(8,n)
       real*8  beta(7)
       integer iseed,idum
       COMMON /vecy/y
       COMMON /vecx/x
       COMMON /vecbeta/beta
c......generate  beta given delta
       call Gbeta(iseed)
       return
       end

       subroutine Gbeta(iseed)
c      generate beta given delta
c      use Gilks log-concave sampling
C *************************** MODEL PARAMETERS *****************************
C -------------------------- Commonly used variables ---------------------
       implicit real*8 (a-h,o-z)
       parameter (n=183)
       PARAMETER(maxtries=15,im=2,NS=10,np=5,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),x(8,n)
       real*8  beta(7)
       integer ifault,iwv(NS+7+in),ifcount
       integer iseed,idum
       logical lb,ub
       COMMON /vecy/y
       COMMON /vecx/x
       COMMON /vecbeta/beta
       COMMON /dummy/idum
       external eval
c     PRINT*, 'Generating betas ...'
      emax=60.0d0
      DO 80 j=1,np
         idum=j
c        write(*,*) 'idum=',idum
c        write(*,*) 'b=',b
c        write(*,*) 'beta(1)=',beta(1)
c        write(*,*) 'beta(2)=',beta(2)
c        Two starting points (step is one std. errors away)
         if (j .eq. 1) then
          step=3.5d0
          a(1)=beta(1)-step
          a(2)=beta(1)+step
         endif
         if (j .eq. 2) then
          step=3.0d0
          a(1)=beta(2)-step
          a(2)=beta(2)+step
         endif
         if (j .eq. 3) then
          step=3.0d0
          a(1)=beta(3)-step
          a(2)=beta(3)+step
         endif
         if (j .eq. 4) then
          step=3.5d0
          a(1)=beta(4)-step
          a(2)=beta(4)+step
         endif
         if (j .eq. 5) then
          step=3.5d0
          a(1)=beta(5)-step
          a(2)=beta(5)+step
         endif
         if (j .eq. 6) then
          step=3.0d0
          a(1)=beta(6)-step
          a(2)=beta(6)+step
         endif
         if (j .eq. 7) then
          step=3.0d0
          a(1)=beta(7)-step
          a(2)=beta(7)+step
         endif
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=183 )
      real*8  y(n),x(8,n)
      real*8  beta(7)
      real*8  sum1,sum2
      real*8  ein
      integer iseed,idum
      COMMON /vecy/y
      COMMON /vecx/x
      COMMON /vecbeta/beta
      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=beta(1)
       do 310 j=2,7
        ein=ein+beta(j)*x(j,i)
 310   continue
       if (ein .gt. 0.0d0) then
         sum1=sum1+ein*(y(i)-1.0d0)
         if (idum .eq. 1) then
          sum2=sum2+y(i)
         else
          sum2=sum2+x(idum,i)*y(i)
         endif
         if (ein .gt. 70.0d0) then
          if (idum .eq. 1) then
            sum2=sum2-1.0d0
          else
            sum2=sum2-x(idum,i)
          endif
         endif
         if (ein .le. 70.0d0) then
          sum1=sum1-dlog(1.0d0+dexp(-ein))
          if (idum .eq. 1) then
           sum2=sum2-1.0d0/(1.0d0+dexp(-ein))
          else
           sum2=sum2-x(idum,i)/(1.0d0+dexp(-ein))
          endif
         endif
       else
         sum1=sum1+ein*y(i)
         if (idum .eq. 1) then
          sum2=sum2+y(i)
         else
          sum2=sum2+x(idum,i)*y(i)
         endif
         if (ein .gt. -70.0d0) then
          sum1=sum1-dlog(1.0d0+dexp(ein))
          if (idum .eq. 1) then
           sum2=sum2-dexp(ein)/
     1          (1.0d0+dexp(ein))
          else
           sum2=sum2-x(idum,i)*dexp(ein)/
     1          (1.0d0+dexp(ein))
          endif
         endif
       endif
 300  continue

c
      ahz=sum1
      ahpz=sum2
      return
      end