subroutine gibbslog(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 Gbetalog(iseed) return end subroutine Gbetalog(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) real*8 start(10),xmin(10),xsec(10),step1(10) real*8 ynewlo,ysec,reqmin 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 COMMON /vecjmac/jmac external evallog,fnbeta c PRINT*, 'Generating betas ...' emax=60.0d0 c if ((jmac .eq. 1) .or. (jmac .eq. 3)) then nopt=np konvge=5 kcount=1000 reqmin=1.0d-10 do 420 k1=1,np start(k1)=beta(k1) step1(k1)=0.2d0 xmin(k1)=start(k1) xsec(k1)=start(k1) 420 continue c.......call optim.f call nelmin(fnbeta, nopt, start, xmin, ynewlo, reqmin, step1, 1 konvge, kcount, icount, numres, ifault) do 421 j=1,np beta(j)=xmin(j) 421 continue c endif 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.0d0 a(1)=beta(1)-step a(2)=beta(1)+step endif if (j .eq. 2) then step=2.5d0 a(1)=beta(2)-step a(2)=beta(2)+step endif if (j .eq. 3) then step=1.5d0 a(1)=beta(3)-step a(2)=beta(3)+step endif if (j .eq. 4) then step=1.5d0 a(1)=beta(4)-step a(2)=beta(4)+step endif if (j .eq. 5) then step=1.5d0 a(1)=beta(5)-step a(2)=beta(5)+step endif if (j .eq. 6) then step=2.0d0 a(1)=beta(6)-step a(2)=beta(6)+step endif if (j .eq. 7) then step=1.5d0 a(1)=beta(7)-step a(2)=beta(7)+step endif 20 CALL evallog(a(1), ha(1), hpa(1)) CALL evallog(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 EVALLOG(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 EVALLOG(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, evallog) IF((ifault.eq.3) .or. (ifault.eq.4)) GOTO 20 C CALL SAMPLE(iwv,rwv,evallog,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 real*8 function fnbeta(zz) 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 zz(5) real*8 ein integer iseed,idum COMMON /vecy/y COMMON /vecx/x do 50 j=1,5 beta(j)=zz(j) 50 continue c.....calculate log likelihood and its derivative sum1=0.0d0 do 300 i=1,n ein=beta(1) do 310 j=2,5 ein=ein+beta(j)*x(j,i) 310 continue ein1=dexp(ein) if (y(i) .eq. 0.0d0) then sum1=sum1-ein1 endif if (y(i) .eq. 1.0d0) then sum1=sum1+dlog(1.0d0-dexp(-ein1)) endif 300 continue c fnbeta=-sum1 return end include 'optim1.f' 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 evallog(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 ein1=dexp(ein) if (y(i) .eq. 0.0d0) then sum1=sum1-ein1 sum2=sum2-x(idum,i)*ein1 endif if (y(i) .eq. 1.0d0) then sum1=sum1+dlog(1.0d0-dexp(-ein1)) ein3=1.0d0-dexp(-ein1) ein4=ein-ein1 ein5=dexp(ein4) sum2=sum2+x(idum,i)*ein5/ein3 endif 300 continue c ahz=sum1 ahpz=sum2 return end include 'gilks1.f'