C subroutine gibbs(n,nc,ifood,itemp,itime,iscore,indexftptm, 1 z,gamma,amu,theta,sqdelta,sig2,indexfd,ntfood,iseed) c purpose: to generate a gibbs iteration implicit real*8 (a-h,o-z) real*8 z(n),gamma(0:9) real*8 amu(1:12,1:4,1:9),theta(1:12,1:4),sqdelta(1:12) real*8 sig2(1:12) integer ifood(n),itemp(n),itime(n),iscore(n) integer nc(4),indexftptm(1:12,1:4,1:9) integer indexfd(12),ntfood integer iseed c......generate z(n) given amu and gamma call Gz(n,ifood,itemp,itime,iscore,amu,gamma, 1 sig2,indexfd,z,iseed) call Gsqdelta(amu,theta,nc,indexfd,indexftptm, 1 ntfood,sqdelta,iseed) c......generate sig2 call Gsig2(n,nc,z,ifood,itemp,itime,indexfd, 1 indexftptm,amu,sig2,iseed) c......generate theta given amu and delta call Gtheta(amu,nc,sqdelta,indexfd,indexftptm,theta,iseed) c......generate amu given z,theta,sqdelta call Gamu(n,nc,ifood,itemp,itime,iscore,z,indexfd, 1 indexftptm,theta,sqdelta,sig2,amu,iseed) return end subroutine Gz(n,ifood,itemp,itime,iscore,amu,gamma, 1 sig2,indexfd,z,iseed) c generate z(n) given amu and gamma implicit real*8 (a-h,o-z) real*8 z(n),gamma(0:9),sig2(1:12) real*8 amu(1:12,1:4,1:9) real*8 aleft,aright,amean,tn,sig integer ifood(n),itemp(n),itime(n),iscore(n) integer iseed,indexfd(12) logical la,lb la=.false. lb=.false. do 300 i=1,n if (indexfd(ifood(i)) .eq. 0) goto 300 amean=amu(ifood(i),itemp(i),itime(i)) sig=dsqrt(sig2(ifood(i))) if (iscore(i) .eq. 1) then la=.true. aright=-amean/sig tn=YTUVN(aleft,aright,la,lb,iseed) z(i)=amean+tn*sig la=.false. endif if ((iscore(i) .ge. 2) .and. (iscore(i) .le. 8)) then aleft =(gamma(iscore(i)-1)-amean)/sig aright=(gamma(iscore(i))-amean)/sig tn=YTUVN(aleft,aright,la,lb,iseed) z(i)=amean+tn*sig endif if (iscore(i) .eq. 9) then lb=.true. aleft=(gamma(8)-amean)/sig tn=YTUVN(aleft,aright,la,lb,iseed) z(i)=amean+tn*sig lb=.false. endif 300 continue return end subroutine Gsqdelta(amu,theta,nc,indexfd,indexftptm, 1 ntfood,sqdelta,iseed) c generate sqdelta given amu and theta implicit real*8 (a-h,o-z) real*8 amu(1:12,1:4,1:9),theta(1:12,1:4),sqdelta(1:12) real*8 sumin real*8 scaler,shape,r(1) integer nc(4),indexfd(12),ntfood integer indexftptm(1:12,1:4,1:9) integer iseed external DRNGAM do 340 i=1,12 if (indexfd(i) .eq. 0) goto 340 sumin=(amu(i,1,1)-theta(i,1))**2 do 350 j=1,4 do 355 l=2,nc(j) if (indexftptm(i,j,l) .eq. 1) goto 355 sumin=sumin+(amu(i,j,l)-theta(i,j))**2 355 continue 350 continue cc shape=23.0d0/2.0d0-1.0d0+3.0000000001d0 cc scaler=sum/2.0d0+0.047699d0 c shape=23.0d0/2.0d0-1.0d0+15.767646d0 c scaler=sumin/2.0d0+0.6567029d0 c shape=23.0d0/2.0d0-1.0d0+3.2571845d0 !using entrees c scaler=sumin/2.0d0+0.4161029d0 !using entrees shape=23.0d0/2.0d0-1.0d0+2.298d0 !fit using pastries data scaler=sumin/2.0d0+0.0049d0 !fit using pastries data (10*variance) c shape=23.0d0/2.0d0-1.0d0+2.0298d0 !fit using pastries data c scaler=sumin/2.0d0+0.0039d0 !fit using pastries data (100*variance) call DRNGAM(1,shape,r) sqdelta(i)=scaler/r(1) 340 continue return end subroutine Gsig2(n,nc,z,ifood,itemp,itime,indexfd, 1 indexftptm,amu,sig2,iseed) implicit real*8 (a-h,o-z) real*8 amu(1:12,1:4,1:9),sig2(1:12) real*8 z(n) real*8 scaler,shape,r(1) real*8 sum(12) integer nc(4),indexfd(12) integer indexftptm(1:12,1:4,1:9) integer iseed c integer icount(12) integer ifood(n),itemp(n),itime(n) external DRNGAM c shape=23.0d0*36.0d0/2.0d0+4.4198d0 !fit from entree data shape=23.0d0*36.0d0/2.0d0+8.44d0 !fit using pastries data x 10 c shape=23.0d0*36.0d0/2.0d0+2.64d0 !fit using pastries data x 100 do 100 i=1,12 sum(i)=0.0d0 c icount(i)=0 100 continue do 200 i=1,n sum(ifood(i))=sum(ifood(i))+(z(i)- 1 amu(ifood(i),itemp(i),itime(i)))**2 c icount(ifood(i))=icount(ifood(i))+1 200 continue do 300 k=1,12 if (indexfd(k) .eq. 0) goto 300 c scaler=sum(k)/(2.0d0)+17.875216 !fit from entrees data scaler=sum(k)/(2.0d0)+0.42505d0 !fit using pastries data x 10 c scaler=sum(k)/(2.0d0)+0.09392d0 !fit using pastries data x 100 call DRNGAM(1,shape,r) sig2(k)=scaler/r(1) c write(*,*) 'k=',k,'scaler=',scaler,'(828)count=',icount(k) 300 continue return end subroutine Gtheta(amu,nc,sqdelta,indexfd,indexftptm,theta,iseed) c......generate theta given amu and delta implicit real*8 (a-h,o-z) real*8 amu(1:12,1:4,1:9),theta(1:12,1:4),sqdelta(1:12) real*8 sum(4),sigma,amean integer nc(4),indexfd(12) integer indexftptm(1:12,1:4,1:9) integer iseed real*8 aleft,aright,tn c real*8 avgmu,alambda integer num(4) logical la,lb c......initialize la=.false. lb=.false. do 360 i=1,12 if (indexfd(i) .eq. 0) goto 360 sum(1)=amu(i,1,1) num(1)=1 do 362 i1=2,4 sum(i1)=0.0d0 num(i1)=0 362 continue do 365 j=1,4 do 370 l=2,nc(j) if (indexfd(i) .eq. 0) goto 370 if (indexftptm(i,j,l) .eq. 1) goto 370 sum(j)=sum(j)+amu(i,j,l) num(j)=num(j)+1 370 continue 365 continue c avgmu=sum(1)/dfloat(num(1)) c alambda=.0072912d0*10.0d0/(.0072912d0*10.0d0+ c 1 sqdelta(i)/dfloat(num(1))) !using pastries c amean=alambda*avgmu+(1.0d0-alambda)*.53138d0 !using pastries c sigma=dsqrt( (1.0d0-alambda)*.0072912d0*10.0d0 ) !using pastries amean=sum(1)/dfloat(num(1)) sigma=dsqrt(sqdelta(i)/dfloat(num(1))) aleft=(theta(i,2)-amean)/sigma lb=.true. tn=YTUVN(aleft,aright,la,lb,iseed) lb=.false. theta(i,1)=amean+sigma*tn do 367 j=2,4 c avgmu=sum(j)/dfloat(num(j)) c alambda=.0072912d0*10.0d0/(.0072912d0*10.0d0+ c 1 sqdelta(i)/dfloat(num(j))) !using pastries c amean=alambda*avgmu+(1.0d0-alambda)*.53138d0 !using pastries c sigma=dsqrt( (1.0d0-alambda)*.0072912d0*10.0d0 ) !using pastries amean=sum(j)/dfloat( num(j) ) sigma=dsqrt(sqdelta(i)/dfloat(num(j))) if (j .lt. 4) then aleft=(theta(i,j+1)-amean)/sigma aright=(theta(i,j-1)-amean)/sigma else aright=(theta(i,j-1)-amean)/sigma la=.true. endif tn=YTUVN(aleft,aright,la,lb,iseed) theta(i,j)=amean+sigma*tn la=.false. 367 continue 360 continue return end subroutine Gamu(n,nc,ifood,itemp,itime,iscore,z,indexfd, 1 indexftptm,theta,sqdelta,sig2,amu,iseed) c......generate amu givenz,theta,sqdelta implicit real*8 (a-h,o-z) real*8 amu(1:12,1:4,1:9),theta(1:12,1:4),sqdelta(1:12) real*8 z(n) real*8 sumz(1:12,1:4,1:9),sig2(1:12) real*8 aleft,aright,tn real*8 amean,sigma integer ifood(n),itemp(n),itime(n),iscore(n) integer nc(4),indexftptm(1:12,1:4,1:9),indexfd(12) integer iseed logical la,lb c......initialize la=.false. lb=.false. do 380 i=1,12 if (indexfd(i) .eq. 0) goto 380 do 385 j=1,4 do 390 l=1,nc(j) sumz(i,j,l)=0.0d0 390 continue 385 continue 380 continue do 400 i=1,n if (indexfd(ifood(i)) .eq. 0) goto 400 sumz(ifood(i),itemp(i),itime(i))= 1 sumz(ifood(i),itemp(i),itime(i))+z(i) 400 continue do 410 i=1,12 if (indexfd(i) .eq. 0) goto 410 lb=.true. amean=(sumz(i,1,1)*sqdelta(i)+theta(i,1)*sig2(i)) 1 / (36.0d0*sqdelta(i)+sig2(i)) sigma=dsqrt(sqdelta(i)*sig2(i)/(36.0d0*sqdelta(i)+sig2(i))) aleft=amu(i,1,3) aleft=max(aleft,amu(i,2,3)) do 415 j=3,4 aleft=max(aleft,amu(i,j,2)) 415 continue aleft=(aleft-amean)/sigma tn=YTUVN(aleft,aright,la,lb,iseed) lb=.false. c........get amu(i,1,1) amu(i,1,1)=amean+sigma*tn do 420 j=1,4 do 430 l=2,nc(j) c.............get two different means and sigmas if (indexftptm(i,j,l) .eq. 1) goto 430 amean=(sumz(i,j,l)*sqdelta(i)+sig2(i)*theta(i,j))/ 1 (36.0d0*sqdelta(i)+sig2(i)) c.............get the bounds of amu's c.............begin first row if ((j .eq. 1) .and. (l .eq. 3)) then aright=(amu(i,1,1)-amean)/sigma aleft=(amu(i,2,3)-amean)/sigma aleft=max((amu(i,1,6)-amean)/sigma,aleft) endif if ((j .eq. 1) .and. (l .eq. 6)) then aright=(amu(i,1,3)-amean)/sigma aleft=(amu(i,2,6)-amean)/sigma aleft=max((amu(i,1,7)-amean)/sigma,aleft) endif if ( (j .eq. 1) .and. (l .lt. nc(1)) 1 .and. (l .gt. 6) ) then aright=(amu(i,j,l-1)-amean)/sigma aleft=(amu(i,j,l+1)-amean)/sigma aleft=max(aleft,(amu(i,j+1,l)-amean)/sigma) endif if ( (j .eq. 1) .and. (l .eq. nc(1)) ) then aright=(amu(i,j,l-1)-amean)/sigma aleft=(amu(i,j+1,l)-amean)/sigma endif c.............end first row c.............begin second row if ((j .eq. 2) .and. (l .eq. 3)) then aright=(amu(i,1,1)-amean)/sigma aright=min(aright,(amu(i,1,3)-amean)/sigma) aleft=(amu(i,2,4)-amean)/sigma aleft=max((amu(i,3,3)-amean)/sigma,aleft) endif if ((j .eq. 2) .and. (l .eq. 4)) then aright=(amu(i,2,3)-amean)/sigma aleft=(amu(i,2,5)-amean)/sigma aleft=max((amu(i,3,4)-amean)/sigma,aleft) endif if ((j .eq. 2) .and. (l .eq. 5)) then aright=(amu(i,2,4)-amean)/sigma aleft=(amu(i,2,6)-amean)/sigma aleft=max((amu(i,3,5)-amean)/sigma,aleft) endif if ( (j .eq. 2) .and. (l .le. nc(3)) 1 .and. (l .gt. 5) ) then aright=(amu(i,j,l-1)-amean)/sigma aright=min(aright, (amu(i,j-1,l)-amean)/sigma) aleft=(amu(i,j,l+1)-amean)/sigma aleft=max(aleft,(amu(i,j+1,l)-amean)/sigma) endif if ( (j .eq. 2) .and. (l .lt. nc(2)) 1 .and. (l .gt. nc(3)) ) then aright=(amu(i,j,l-1)-amean)/sigma aright=min(aright, (amu(i,j-1,l)-amean)/sigma) aleft=(amu(i,j,l+1)-amean)/sigma endif if ( (j .eq. 2) .and. (l .eq. nc(2)) ) then aright=(amu(i,j,l-1)-amean)/sigma aright=min(aright, (amu(i,j-1,l)-amean)/sigma) la=.true. endif c.............end second row c.............begin third row if ((j .eq. 3) .and. (l .eq. 2)) then aright=(amu(i,1,1)-amean)/sigma aleft=(amu(i,3,3)-amean)/sigma aleft=max((amu(i,4,2)-amean)/sigma,aleft) endif if ( (j .eq. 3) .and. (l .le. nc(4)) 1 .and. (l .gt. 2) ) then aright=(amu(i,j,l-1)-amean)/sigma aright=min(aright, (amu(i,j-1,l)-amean)/sigma) aleft=(amu(i,j,l+1)-amean)/sigma aleft=max(aleft,(amu(i,j+1,l)-amean)/sigma) endif if ( (j .eq. 3) .and. (l .lt. nc(3)) 1 .and. (l .gt. nc(4)) ) then aright=(amu(i,j,l-1)-amean)/sigma aright=min(aright, (amu(i,j-1,l)-amean)/sigma) aleft=(amu(i,j,l+1)-amean)/sigma endif if ( (j .eq. 3) .and. (l .eq. nc(3)) ) then aright=(amu(i,j,l-1)-amean)/sigma aright=min(aright, (amu(i,j-1,l)-amean)/sigma) la=.true. endif c.............end third row c.............begin forth row if ((j .eq. 4) .and. (l .eq. 2)) then aright=(amu(i,1,1)-amean)/sigma aright=min(aright,(amu(i,j-1,2)-amean)/sigma) aleft=(amu(i,4,3)-amean)/sigma endif if ( (j .eq. 4) .and. (l .lt. nc(4)) 1 .and. (l .gt. 2) ) then aright=(amu(i,j,l-1)-amean)/sigma aright=min(aright, (amu(i,j-1,l)-amean)/sigma) aleft=(amu(i,j,l+1)-amean)/sigma endif if ( (j .eq. 4) .and. (l .eq. nc(4)) ) then aright=(amu(i,j,l-1)-amean)/sigma aright=min(aright, (amu(i,j-1,l)-amean)/sigma) la=.true. endif c.............end forth row c........get amu(i,j,l) tn=YTUVN(aleft,aright,la,lb,iseed) la=.false. lb=.false. amu(i,j,l)=amean+sigma*tn 430 continue 420 continue 410 continue return end