       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
