program foodgibbs c Bayesian order-restricted inference for a categorized c population: an application to the meal, ready-to-eat c algorithm: gibbs sampling c programmed by ming-hui chen c started in department of mathematical sciences, c worcester polytechnic institute, december 28, 1993 c c this program produces mean shelf-life times c and without constraints on mu(ijk) c use the model: mu(ijk)=theta(k) + error c for whole 12 entrees c July 31, 1995 c implicit real*8 (a-h,o-z) parameter(n=9828) c......parameters 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),sig c......for RIS c real*8 risamu(1:500,1:4,1:9) c real*8 risamu(1:50000,1:4,1:9) real*8 risamu(1:10000,1:4,1:9) c......cumulative varaibles real*8 samu(1:12,1:4,1:9),s2amu(1:12,1:4,1:9) real*8 sta(12,4),s2ta(12,4) real*8 ssig2(12),s2sig2(12) real*8 ssqd(1:12),s2sqd(1:12) real*8 sms(1:12,1:4,1:9),s2ms(1:12,1:4,1:9) real*8 sems(1:4,1:9) real*8 summs real*8 sqsein real*8 sqms(5001,1:4,1:9) c......estimated varaiables real*8 eamu(1:12,1:4,1:9),sdamu(1:12,1:4,1:9) real*8 eta(12,4),sdta(12,4) real*8 esig2(1:12),sdsig2(1:12) real*8 esqd(1:12),sdsqd(1:12) real*8 ems(1:12,1:4,1:9),sdms(1:12,1:4,1:9) real*8 ta(4),delta2 c......vectors for reading data integer ifood(n),itemp(n),itime(n),iscore(n) c......index and constants integer nc(4),indexftptm(1:12,1:4,1:9) integer nrep,iseed,lseed,ind,irep integer nsum, icum c......set constant real*8 time(1:9),temp(1:4) c......store Gibbs-observation variables real*8 avgz(1:12,1:4,1:9) real*8 summu,summuth,sumsig1,sumsig2 integer iout,m integer indexafd(1:4,1:9) c......flag indicator variable for foodid integer indexfd(12) real*8 seqsig2(5000),seqs2del(5001) real*8 seqta(5001,1:4) c real*8 cm(5001),sinvcm,scm,sum2in real*8 anu,sum1,sum2,sumin,do1 c......for marginal densities real*8 d1LOW(5001), d1mean(5001), d1std(5001) real*8 d2UPP(5001),d2LOW(5001),d2mean(5001),d2std(5001) real*8 d3UPP(5001),d3mean(5001),d3std(5001) real*8 aleft, aright real*8 DNORDF real*8 weight(5001),weight1(5001) real*8 ahpdL(1:4,1:9),ahpdU(1:4,1:9) real*8 xac(50001) integer nq data nc/9,9,7,5/ data time/0.0d0,6.0d0,12.0d0,18.0d0,24.0d0,30.0d0, 1 36.0d0,48.0d0,60.0d0/ data temp/4.0d0,21.0d0,30.0d0,38.0d0/ common /para/s common /indic/indexafd external DNORDF c c.......initialize indexfd do 15 i=1,12 indexfd(i)=0 15 continue ntfood=0 m=10000 c......enter the number of iterations for gibbs sampler c print*,'enter the number of iterations for gibbs sampler' c read(*,*) nrep nrep=75000 do 10 i=1,12 indexfd(i)=0 10 continue indexfd(2)=1 c......set initial random number seeds iseed=999999999 c......set constants for gamma gamma(9)=9999999.0d2 gamma(1)=0.0d0 gamma(0)=-9999999.0d2 c......set ta(4) and delta2 ta(1)=0.456979d0 ta(2)=0.428031d0 ta(3)=0.400165d0 ta(4)=0.366237d0 delta2=0.001662d0 anu=40.0d0 c......read data call rdat(n,ifood,itemp,itime,iscore) c......set index function call assind(n,nc,ifood,itemp,itime,indexftptm) do 40 i=1,4 do 50 j=1,nc(i) indexafd(i,j)=indexftptm(1,i,j)-1 50 continue 40 continue c......set the starting values of parameters call starting(sig2,theta,amu) gamma(2)=0.098631d0 gamma(3)=0.19147d0 gamma(4)=0.29905d0 gamma(5)=0.37671d0 gamma(6)=0.52417d0 gamma(7)=0.71926d0 gamma(8)=1.0d0 c call rnset( iseed ) c T_start = cpsec() call RIS(m,nc,indexfd,indexftptm,amu,risamu,iseed) c......burn in the Gibbs markov chain do 90 i=1,500 c......call gibbs sampler call gibbs(n,nc,ifood,itemp,itime,iscore,indexftptm, 1 z,gamma,amu,theta,sqdelta,sig2,indexfd,ntfood,iseed) 90 continue irep=0 nsum=0 icum=0 c......initialize cumulative variables do 260 k=1,12 if (indexfd(k) .eq. 0) goto 260 ssqd(k)=0.0d0 s2sqd(k)=0.0d0 ssig2(k)=0.0d0 s2sig2(k)=0.0d0 do 263 k1=1,4 sta(k,k1)=0.0d0 s2ta(k,k1)=0.0d0 263 continue do 270 i=1,4 do 280 j=1,nc(i) samu(k,i,j)=0.0d0 s2amu(k,i,j)=0.0d0 sms(k,i,j)=0.0d0 s2ms(k,i,j)=0.0d0 280 continue 270 continue 260 continue sinvcm=0.0d0 c do 100 i=1,nrep if ( (i .eq. 1) .or. (i .eq. irep) ) then write(*,*) 'i=',i irep=irep+100 endif c......call gibbs sampler call gibbs(n,nc,ifood,itemp,itime,iscore,indexftptm, 1 z,gamma,amu,theta,sqdelta,sig2,indexfd,ntfood,iseed) c......store information if ( (i .eq. 1) .or. (i .eq. icum) ) then c icum=icum+10 icum=icum+15 if (i .gt. 1) then nsum=nsum+1 do 777 k=1,12 if (indexfd(k) .eq. 0) goto 777 do 778 i1=1,4 do 779 j=1,nc(i1) avgz(k,i1,j)=0.0d0 779 continue 778 continue 777 continue do 401 i1=1,n if (indexfd(ifood(i1)) .eq. 0) goto 401 avgz(ifood(i1),itemp(i1),itime(i1))= 1 avgz(ifood(i1),itemp(i1),itime(i1))+z(i1) 401 continue seqsig2(nsum)=sig2(2) seqs2del(nsum)=sqdelta(2) do 403 j=1,4 seqta(nsum,j)=theta(2,j) 403 continue c seqamu(nsum)=amu(2,3,3) c seqz(nsum)=avgz(2,3,3)/36.0d0 c........store info for density estimation c........density 1 for amu(2,1,1) d1mean(nsum)=(avgz(2,1,1)*sqdelta(2)+ 1 theta(2,1)*sig2(2))/(36.0d0*sqdelta(2)+sig2(2)) d1std(nsum)=dsqrt(sqdelta(2)*sig2(2)/ 1 (36.0d0*sqdelta(2)+sig2(2))) aleft=amu(2,1,3) aleft=max(aleft,amu(2,2,3)) do 415 j=3,4 aleft=max(aleft,amu(2,j,2)) 415 continue d1LOW(nsum)=aleft c........density 2 for amu(2,2,6) d2mean(nsum)=(avgz(2,2,6)*sqdelta(2)+ 1 sig2(2)*theta(2,2))/(36.0d0*sqdelta(2)+sig2(2)) d2std(nsum)=d1std(nsum) aleft=amu(2,2,7) aleft=max(aleft,amu(2,3,6)) d2LOW(nsum)=aleft aright=amu(2,1,6) aright=min(aright,amu(2,2,5)) d2UPP(nsum)=aright c........density 3 for amu(2,4,5) d3mean(nsum)=(avgz(2,4,5)*sqdelta(2)+ 1 sig2(2)*theta(2,4))/(36.0d0*sqdelta(2)+sig2(2)) d3std(nsum)=d1std(nsum) aright=amu(2,3,5) aright=min(aright,amu(2,4,4)) d3UPP(nsum)=aright c c.......calculate cm scm=0.0d0 do 303 j=1,m sum1=0.0d0 sum2=0.0d0 do 305 i1=1,4 if (i1 .eq. 1) then do 307 j1=1,nc(i1) if (indexafd(i1,j1) .eq. 0) goto 307 sum1=sum1+0.5d0*dlog(sqdelta(2)) 1 +(risamu(j,i1,j1)-theta(2,i1))**2 1 /(2.0d0*sqdelta(2)) sum2=sum2+0.5d0*dlog(delta2) 1 +dlog( 1.0d0+(risamu(j,i1,j1)-ta(i1))**2 1 /(anu*delta2) )*(anu+1.0d0)/2.0d0 307 continue else do 309 j1=2,nc(i1) if (indexafd(i1,j1) .eq. 0) goto 309 sum1=sum1+0.5d0*dlog(sqdelta(2)) 1 +(risamu(j,i1,j1)-theta(2,i1))**2 1 /(2.0d0*sqdelta(2)) sum2=sum2+0.5d0*dlog(delta2) 1 +dlog( 1.0d0+(risamu(j,i1,j1)-ta(i1))**2 1 /(anu*delta2) )*(anu+1.0d0)/2.0d0 309 continue endif 305 continue if (sum2-sum1 .gt. -50.0d0) then scm=scm+dexp(sum2-sum1) endif 303 continue cm(nsum)=scm/real(m) sinvcm=sinvcm+1.0d0/cm(nsum) write(*,*) 'nsum=',nsum,' cm(nsum)=',cm(nsum) c c.........for amu and mean score k=2 c.........for sqdelta ssqd(k)=ssqd(k)+sqdelta(k)/cm(nsum) s2sqd(k)=s2sqd(k)+sqdelta(k)**2/cm(nsum) c.........for sig2 ssig2(k)=ssig2(k)+sig2(k)/cm(nsum) s2sig2(k)=s2sig2(k)+sig2(k)**2/cm(nsum) c......for theta do 313 k1=1,4 sta(k,k1)=sta(k,k1)+theta(k,k1)/cm(nsum) s2ta(k,k1)=s2ta(k,k1)+theta(k,k1)**2/cm(nsum) 313 continue c.........for amu and mean score sig=dsqrt(sig2(k)) summu=0.0d0 summuth=0.0d0 c do 310 i1=1,4 if (i1 .eq. 1) then do 320 j=1,nc(i1) if (indexafd(i1,j) .eq. 0) goto 320 samu(k,i1,j)=samu(k,i1,j)+amu(k,i1,j)/cm(nsum) s2amu(k,i1,j)=s2amu(k,i1,j)+ 1 amu(k,i1,j)**2/cm(nsum) c..............cumulate summu and summuth summu=summu+amu(k,i1,j) summuth=summuth+(amu(k,i1,j)-theta(k,i1))**2 c..............for mean score summs=1.0d0 do 315 l=1,8 summs=summs+ 1 (1.0d0-DNORDF( (gamma(l)-amu(k,i1,j))/sig )) 315 continue sqms(nsum,i1,j)=summs sms(k,i1,j)=sms(k,i1,j)+summs/cm(nsum) s2ms(k,i1,j)=s2ms(k,i1,j)+summs**2/cm(nsum) 320 continue else do 325 j=2,nc(i1) if (indexafd(i1,j) .eq. 0) goto 325 samu(k,i1,j)=samu(k,i1,j)+amu(k,i1,j)/cm(nsum) s2amu(k,i1,j)=s2amu(k,i1,j)+ 1 amu(k,i1,j)**2/cm(nsum) c..............cumulate summu and summuth summu=summu+amu(k,i1,j) summuth=summuth+(amu(k,i1,j)-theta(k,i1))**2 c..............for mean score summs=1.0d0 do 316 l=1,8 summs=summs+ 1 (1.0d0-DNORDF( (gamma(l)-amu(k,i1,j))/sig )) 316 continue sqms(nsum,i1,j)=summs sms(k,i1,j)=sms(k,i1,j)+summs/cm(nsum) s2ms(k,i1,j)=s2ms(k,i1,j)+summs**2/cm(nsum) 325 continue endif 310 continue c endif endif 100 continue c T_end = cpsec() call rnget( iseed ) lseed=iseed c T_gibbs = T_end - T_start c c......calculate the weights do 4000 i=1,nsum weight(i)=1.0d0/(cm(i)*sinvcm) 4000 continue c c......obtain estimates c do 340 k=1,12 if (indexfd(k) .eq. 0) goto 340 c.........for sqdelta esqd(k)=ssqd(k)/sinvcm sdsqd(k)=dsqrt( s2sqd(k)/sinvcm - 1 esqd(k)**2 ) c.........for sig2 esig2(k)=ssig2(k)/sinvcm sdsig2(k)=dsqrt( s2sig2(k)/sinvcm 1 - esig2(k)**2 ) c.........for theta do 343 k1=1,4 eta(k,k1)=sta(k,k1)/sinvcm sdta(k,k1)=dsqrt( s2ta(k,k1)/sinvcm 1 -eta(k,k1)**2 ) 343 continue c.........for mu do 350 i=1,4 c.............for mean score q=real(nsum)*0.95d0 nq=nint(q) if (i .eq. 1) then do 360 j=1,nc(i) if (indexafd(i,j) .eq. 0) goto 360 eamu(k,i,j)=samu(k,i,j)/sinvcm sdamu(k,i,j)=dsqrt( s2amu(k,i,j)/sinvcm 1 - eamu(k,i,j)**2 ) ems(k,i,j)=sms(k,i,j)/sinvcm sdms(k,i,j)=dsqrt( s2ms(k,i,j)/sinvcm 1 - ems(k,i,j)**2 ) do 4100 i1=1,nsum xac(i1)=sqms(i1,i,j) weight1(i1)=weight(i1) 4100 continue do 4111 i1=1,nsum-1 do 4112 j1=i1+1,nsum if (xac(i1) .gt. xac(j1)) then temp2=xac(i1) xac(i1)=xac(j1) xac(j1)=temp2 temp1=weight1(i1) weight1(i1)=weight1(j1) weight1(j1)=temp1 endif 4112 continue 4111 continue whpd=1000000.0d0 do 4600 i1=1,nsum-nq sum1=0.0d0 q1=real(i1)/real(nsum) q2=(real(i1)+real(nq))/real(nsum) do 4601 i2=1,nsum if (i2 .gt. 1) then sum1=sum1+weight1(i2-1) endif sum2=sum1+weight1(i2) if ( (sum1 .lt. q1) .and. (sum2 .ge. q1) ) 1 alow1=xac(i2) if ( (sum1 .lt. q2) .and. (sum2 .ge. q2) ) 1 aupp1=xac(i2) 4601 continue wb=aupp1-alow1 if (whpd .gt. wb) then whpd=wb aupp=aupp1 alow=alow1 endif 4600 continue ahpdL(i,j)=alow ahpdU(i,j)=aupp 360 continue else do 365 j=2,nc(i) if (indexafd(i,j) .eq. 0) goto 365 eamu(k,i,j)=samu(k,i,j)/sinvcm sdamu(k,i,j)=dsqrt( s2amu(k,i,j)/sinvcm 1 - eamu(k,i,j)**2 ) c.............for mean score ems(k,i,j)=sms(k,i,j)/sinvcm sdms(k,i,j)=dsqrt( s2ms(k,i,j)/sinvcm 1 - ems(k,i,j)**2 ) q=real(nsum)*0.95d0 nq=nint(q) do 5100 i1=1,nsum xac(i1)=sqms(i1,i,j) weight1(i1)=weight(i1) 5100 continue do 5111 i1=1,nsum-1 do 5112 j1=i1+1,nsum if (xac(i1) .gt. xac(j1)) then temp2=xac(i1) xac(i1)=xac(j1) xac(j1)=temp2 temp1=weight1(i1) weight1(i1)=weight1(j1) weight1(j1)=temp1 endif 5112 continue 5111 continue whpd=1000000.0d0 do 5600 i1=1,nsum-nq sum1=0.0d0 q1=real(i1)/real(nsum) q2=(real(i1)+real(nq))/real(nsum) do 5601 i2=1,nsum if (i2 .gt. 1) then sum1=sum1+weight1(i2-1) endif sum2=sum1+weight1(i2) if ( (sum1 .lt. q1) .and. (q1 .le. sum2) ) 1 alow1=xac(i2) if ( (sum1 .lt. q2) .and. (q2 .le. sum2) ) 1 aupp1=xac(i2) 5601 continue wb=aupp1-alow1 if (whpd .gt. wb) then whpd=wb aupp=aupp1 alow=alow1 endif 5600 continue ahpdL(i,j)=alow ahpdU(i,j)=aupp 365 continue endif 350 continue 340 continue c.......to obtain the standard error do1=real(nsum)/real(m) if (do1 .gt. 0.0d0) goto 1001 do 600 i=1,4 if (i .eq. 1) then do 610 j=1,nc(i) if (indexafd(i,j) .eq. 0) goto 610 write(*,*) 'for se i=',i,' j=',j sumsig1=0.0d0 do 640 i1=1,nsum sumsig1=sumsig1+( real(nsum)* 1 (sqms(i1,i,j)-ems(2,i,j)) 1 /(cm(i1)*sinvcm) )**2 640 continue sumsig1=sumsig1/real(nsum) sumsig2=0.0d0 icount=0 do 645 j1=1,m if ( (j1 .eq. 1) .or. (j1 .eq. icount) ) then icount=icount+100 write(*,*) 'in se j1=',j1 endif sumin=0.0d0 do 650 i1=1,nsum c sum1=0.0d0 sum2=0.0d0 do 605 k1=1,4 if (k1 .eq. 1) then do 607 k2=1,nc(k1) if (indexafd(k1,k2) .eq. 0) goto 607 sum1=sum1+0.5d0*dlog(seqs2del(i1)) 1 +(risamu(j1,k1,k2)-seqta(i1,k1))**2 1 /(2.0d0*seqs2del(i1)) sum2=sum2+0.5d0*dlog(delta2) 1 +dlog( 1.0d0+(risamu(j1,k1,k2)-ta(k1))**2 1 /(anu*delta2) )*(anu+1.0d0)/2.0d0 607 continue else do 609 k2=2,nc(k1) if (indexafd(k1,k2) .eq. 0) goto 609 sum1=sum1+0.5d0*dlog(seqs2del(i1)) 1 +(risamu(j1,k1,k2)-seqta(i1,k1))**2 1 /(2.0d0*seqs2del(i1)) sum2=sum2+0.5d0*dlog(delta2) 1 +dlog( 1.0d0+(risamu(j1,k1,k2)-ta(k1))**2 1 /(anu*delta2) )*(anu+1.0d0)/2.0d0 609 continue endif 605 continue sum2in=sum2-sum1 if (sum2in .gt. -50.0d0) then sqsein=dexp(sum2in) else sqsein=0.0d0 endif c sumin=sumin+(sqsein-cm(i1))* 1 real(nsum)*(sqms(i1,i,j)-ems(2,i,j)) 1 /(cm(i1)**2*sinvcm) 650 continue sumin=sumin/real(nsum) sumsig2=sumsig2+sumin**2 645 continue sumsig2= sumsig2/real(m) write(*,*) 'sumsig1=',sumsig1,' sumsig2=', 1 sumsig2 sems(i,j)=dsqrt( (sumsig1+do1*sumsig2) 1 / real(nsum) ) write(*,*) 'sems(i,j)=',sems(i,j) 610 continue else do 620 j=2,nc(i) if (indexafd(i,j) .eq. 0) goto 620 write(*,*) 'for se i=',i,' j=',j sumsig1=0.0d0 do 740 i1=1,nsum sumsig1=sumsig1+( real(nsum)* 1 (sqms(i1,i,j)-ems(2,i,j)) 1 /(cm(i1)*sinvcm) )**2 740 continue sumsig1=sumsig1/real(nsum) sumsig2=0.0d0 icount=0 do 745 j1=1,m if ( (j1 .eq. 1) .or. (j1 .eq. icount) ) then icount=icount+100 write(*,*) 'in se j1=',j1 endif sumin=0.0d0 do 750 i1=1,nsum c sum1=0.0d0 sum2=0.0d0 do 705 k1=1,4 if (k1 .eq. 1) then do 707 k2=1,nc(k1) if (indexafd(k1,k2) .eq. 0) goto 707 sum1=sum1+0.5d0*dlog(seqs2del(i1)) 1 +(risamu(j1,k1,k2)-seqta(i1,k1))**2 1 /(2.0d0*seqs2del(i1)) sum2=sum2+0.5d0*dlog(delta2) 1 +dlog( 1.0d0+(risamu(j1,k1,k2)-ta(k1))**2 1 /(anu*delta2) )*(anu+1.0d0)/2.0d0 707 continue else do 709 k2=2,nc(k1) if (indexafd(k1,k2) .eq. 0) goto 709 sum1=sum1+0.5d0*dlog(seqs2del(i1)) 1 +(risamu(j1,k1,k2)-seqta(i1,k1))**2 1 /(2.0d0*seqs2del(i1)) sum2=sum2+0.5d0*dlog(delta2) 1 +dlog( 1.0d0+(risamu(j1,k1,k2)-ta(k1))**2 1 /(anu*delta2) )*(anu+1.0d0)/2.0d0 709 continue endif 705 continue sum2in=sum2-sum1 if (sum2in .gt. -50.0d0) then sqsein=dexp(sum2in) else sqsein=0.0d0 endif c sumin=sumin+(sqsein-cm(i1))* 1 real(nsum)*(sqms(i1,i,j)-ems(2,i,j)) 1 /(cm(i1)**2*sinvcm) 750 continue sumin=sumin/real(nsum) sumsig2=sumsig2+sumin**2 745 continue sumsig2= sumsig2/real(m) write(*,*) 'sumsig1=',sumsig1,' sumsig2=', 1 sumsig2 sems(i,j)=dsqrt( (sumsig1+do1*sumsig2) 1 / real(nsum) ) write(*,*) 'sems(i,j)=',sems(i,j) 620 continue endif 600 continue c......write all estimates and other informations 1001 open(unit=15,file='normc_hpd.out',access='sequential', 1 status='unknown') write(15,*) 1 '***************************************************' write(15,*) 1 '* *' write(15,*) 1 '* OUTPUT FOR FOOD SHELF-LIFE MODEL *' write(15,*) 1 '* (INCLUDING MEAN SHELF_LIFE TIME) *' write(15,*) 1 '* (Using Multiple Linear Fitting) *' write(15,*) 1 '* (Model: mu(ijk)=theta_k + error) *' write(15,*) 1 '* ( Whole 12 Entrees ) *' write(15,*) 1 '* ( Fitted Hyperparameters Using data ) *' write(15,*) 1 '* ( Using Arrhenius Model for Shelf Life ) *' write(15,*) 1 '* *' write(15,*) 1 '***************************************************' write(15,*) ' PART I: Mean Scores (Estimates) ' write(15,*) ' (standard Deviations) ' write(15,*) ' ============================== ' 996 format(2x, 9f9.6) 896 format(1x,'(', 9(f9.6),')') 997 format(11x,9f9.6) 897 format(10x,'(', 8(f9.6),')') 998 format(11x,9f9.6) 898 format(10x,'(', 6(f9.6),')') 999 format(11x,9f9.6) 899 format(10x,'(', 4(f9.6),')') do 402 k=1,12 if (indexfd(k) .eq. 0) goto 402 write(15,*) ' FOOD ID=:',k write(15,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~ ' do 400 j=1,4 if (j .eq. 1) then write(15,996) (ems(k,j,l),l=1,nc(j)) write(15,896) (sdms(k,j,l),l=1,nc(j)) c write(15,896) (sems(j,l),l=1,nc(j)) write(15,996) (ahpdL(j,l),l=1,nc(j)) write(15,996) (ahpdU(j,l),l=1,nc(j)) endif if (j .eq. 2) then write(15,997) (ems(k,j,l),l=2,nc(j)) write(15,897) (sdms(k,j,l),l=2,nc(j)) c write(15,897) (sems(j,l),l=2,nc(j)) write(15,997) (ahpdL(j,l),l=2,nc(j)) write(15,997) (ahpdU(j,l),l=2,nc(j)) endif if (j .eq. 3) then write(15,998) (ems(k,j,l),l=2,nc(j)) write(15,898) (sdms(k,j,l),l=2,nc(j)) c write(15,898) (sems(j,l),l=2,nc(j)) write(15,998) (ahpdL(j,l),l=2,nc(j)) write(15,998) (ahpdU(j,l),l=2,nc(j)) endif if (j .eq. 4) then write(15,999) (ems(k,j,l),l=2,nc(j)) write(15,899) (sdms(k,j,l),l=2,nc(j)) c write(15,899) (sems(j,l),l=2,nc(j)) write(15,999) (ahpdL(j,l),l=2,nc(j)) write(15,999) (ahpdU(j,l),l=2,nc(j)) endif 400 continue write(15,*) 1 '---------------------------------------------------------' 402 continue write(15,*) ' PART II: Latent Variable Means (Estimates) ' write(15,*) ' (standard Deviations) ' write(15,*) ' ========================================== ' do 451 k=1,12 if (indexfd(k) .eq. 0) goto 451 write(15,*) 'Food ID=', k write(15,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~ ' do 450 j=1,4 if (j .eq. 1) then write(15,996) (eamu(k,j,l),l=1,nc(j)) write(15,896) (sdamu(k,j,l),l=1,nc(j)) endif if (j .eq. 2) then write(15,997) (eamu(k,j,l),l=2,nc(j)) write(15,897) (sdamu(k,j,l),l=2,nc(j)) endif if (j .eq. 3) then write(15,998) (eamu(k,j,l),l=2,nc(j)) write(15,898) (sdamu(k,j,l),l=2,nc(j)) endif if (j .eq. 4) then write(15,999) (eamu(k,j,l),l=2,nc(j)) write(15,899) (sdamu(k,j,l),l=2,nc(j)) endif 450 continue write(15,*) 1 '---------------------------------------------------------' 451 continue write(15,*) ' PART III: Hyper-Mean of Latent Variable ' write(15,*) ' Means Mu: Theta (Estimates) ' write(15,*) ' (standard Deviations) ' write(15,*) ' ========================================== ' write(15,1003) 1003 format(1x,'Food',1x,' Estimates ',1x,'(Standard Error)') 1004 format(1x,I4,1x,I4,1x,f11.6,1x,'(',f14.6,')') do 454 k=1,12 if (indexfd(k) .eq. 0) goto 454 do 457 k1=1,4 write(15,1004) k,k1,eta(k,k1),sdta(k,k1) 457 continue 454 continue write(15,*) 1 '---------------------------------------------------------' write(15,*) ' PART IV: Variance of Latent Variables ' write(15,*) ' Means: Delta Square (Estimates) ' write(15,*) ' (standard Deviations) ' write(15,*) ' ========================================== ' write(15,1006) 1006 format(1x,'Food ID',1x,' Estimate ',1x,'Standard Error') do 456 k=1,12 if (indexfd(k) .eq. 0) goto 456 write(15,1005) k,esqd(k),sdsqd(k) write(15,*) 'Latent variable variances' write(15,1005) k,esig2(k),sdsig2(k) 1005 format(2x,I6,1x,2f14.6) write(15,*)'=======================' 456 continue write(15,*) 1 '---------------------------------------------------------' write(15,*) ' PART VII: General Information ' write(15,*) ' ============================ ' write(15,*) 'Gibbs sampling runtime =', T_gibbs write(15,*) 'number of Gibbs iterations=',nrep write(15,*) 'Last Random Number seed=',lseed write(15,*) 'nsum=', nsum write(15,*) 'sinvcm=', sinvcm write(15,*) 1 '---------------------------------------------------------' close(15) c......write out for density 2 for amu(2,2,6) open(unit=12,file='den2.dat',access='sequential', 1 status='unknown') do 381 i=1,nsum write(12,1998) i,d2LOW(i),d2UPP(i), 1 d2mean(i),d2std(i),cm(i) 1998 format(1x,I5,1x,5f16.8) 381 continue close(12) c......write out for density 3 for amu(2,4,5) open(unit=14,file='den3.dat',access='sequential', 1 status='unknown') do 382 i=1,nsum write(14,1997) i,d3UPP(i),d3mean(i),d3std(i),cm(i) 1997 format(1x,I5,1x,5f16.8) 382 continue close(14) stop end c c include start.f for setting starting values of parameters include 'springstart.f' c c inclue gibbs subroutine include 'gibbsspring.f' c c include ris generation include 'ris.f' c subroutine assind(n,nc,ifood,itemp,itime,indexftptm) c assign real values to indicated function c indexftptm for food, temp and time c 1: for the possible cell; 2: for nonempty cell integer ifood(n),itemp(n),itime(n),nc(4) integer indexftptm(1:12,1:4,1:9) c......initialize indexftptm do 510 k=1,12 do 515 i=1,4 do 520 j=1,nc(i) indexftptm(k,i,j)=1 520 continue 515 continue 510 continue do 525 i=1,n indexftptm(ifood(i),itemp(i),itime(i))=2 525 continue return end subroutine rdat(n,ifood,itemp,itime,iscore) c data file: entreew.dat c total number of observations: 9828 integer ifood(n),itemp(n),itime(n),iscore(n) open(unit=10, file= 1 'entreew.dat', 1 status='old') do 500 i=1,n read(10,*) ifood(i),itemp(i),itime(i),iscore(i) 500 continue close(10) return end c include 'tnorm.f' c