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