        program fruit
c       minghui chen, september 2, 1995, purdue university
c       evaluate the bayesian estimations for linear regression
c       model y=b1x1+b2x2+...+bkxk for fruit data
c       input
c          iseed: random number seed
c          n    : iteration size
c          x(10,207): numbers of trees at various ages
c          y(207): fruit produced by total trees x(10,207) 
c          k=10 : dimension of the regression parameters
c       output
c          estsig : estimator of sigma^2
c          sesig  : standard error of estimator of sigma^2
c          estb   : estimators of beta(10)
c          seb    : standard error of estb
c       purpose: to generate (sigma^2, b(10)) by using
c                gibbs sampler
c        .......

       implicit real*8 (a-h,o-z)
       parameter (m=10)
       real*8 x(10,250), y(250) 
       real*8 b(10) 
       real*8 bmac(m,10), smac(m)
       real*8 estb(m,10), s2b(m,10)
       real*8 bbb(10,10000)
       real*8 bbbs(10000),xbar,sdx,rootr,u97_5r
       real*8 xmean1(m),xvar1(m)
       real*8 best(10),bsd(10),bpsr(10),b97_5(10)
       real*8 sigest,sigsd,spsr,s97_5
       real*8 sumsig, sumsig2, sumb(10), sumb2(10)
       real*8 estsig(m), s2sig(m), sig2
       real*8 T_start, T_end, T_hit
       parameter (iprint=3,maxlag=10,nobs=10000)
       integer imean,iseopt,icount,imac,nrep
       real*8  ac(0:maxlag),acv(0:maxlag),seac(maxlag),
     1         xx(nobs), xmean
       real*8  cpo(250),amu(250),ad(250),av(250),ar(250)
       real*8  sum1,sum2,sum3,api,sum4,sumlcpo,af
       external CONST

       iseed=999999999
       api=CONST('PI')
       write(*,*) 'pi=',api
       print *, 'enter iteration size nrep='
       read(*,*) nrep
       open(unit=12,file='start1.data',status='old')
       do  31 i=1,10
         read(12,*) (bmac(i,j),j=1,10),smac(i)
         write(*,996) i,(bmac(i,j),j=1,10),smac(i)
 996     format(1x,I3,10f8.4,f15.3)
 31    continue
       close(12)

c........read in data set betaprgm
       call datain(x,y)

c......mac loop
       T_start = cpsec()
       do 600 imac=1,1
       write(*,*) 'imac=',imac
        do 610 i=1,10
         b(i)=bmac(imac,i)
 610    continue
        sig2=smac(imac)

        call rnset( iseed )
c........warm up markov chain
       icount=0
       do 308 i=1,100
          call ghr(x,y,sig2,b,iseed)
          if ( (i .eq. 1) .or. (i .eq. icount) ) then
          write(*,*) 'warm up i=',i
          icount=icount+100
          endif
 308   continue
c
c.......to generate beta(10), sigma^2
       sumsig = 0.0d0
       sumsig2 = 0.0d0
       do 309 i=1,10
          sumb(i)=0.0d0
          sumb2(i)=0.0d0
 309   continue
c
       icount=0
       do 300 i=1,nrep
          do 305 i1=1,5
           call ghr(x,y,sig2,b,iseed)
 305      continue
          if ( (i .eq. 1) .or. (i .eq. icount) ) then
          write(*,*) 'gibbs i=',i
          icount=icount+100
          endif
          sumsig=sumsig + sig2
          sumsig2=sumsig2+sig2**2
          do 310 j=1,10
             bbb(j,i)=b(j)
             sumb(j)  = sumb(j)  + b(j)
             sumb2(j) = sumb2(j) + b(j)*b(j)
 310      continue
          bbbs(i)=sig2
 300   continue
 
c      to get estimators
       estsig(imac) = sumsig/dfloat(nrep)
       sigest=estsig(imac)
       s2sig(imac) =(sumsig2-dfloat(nrep)*(estsig(imac)**2))
     1               /dfloat(nrep-1)
       sigsd=dsqrt( s2sig(imac) )
       do 320 j = 1, 10
        estb(imac,j) = sumb(j)/dfloat(nrep)
        s2b(imac,j)=(sumb2(j)-dfloat(nrep)*(estb(imac,j)**2))
     1                /dfloat(nrep-1)
        best(j)=estb(imac,j)
        bsd(j)=dsqrt( s2b(imac,j) )
 320   continue

      
       imean=1
       iseopt=1
       do 301 j=1,10
         do 302 i=1,nrep
             xx(i)=bbb(j,i)
 302     continue
         write(*,*) '====================='
         write(*,*) 'imac=',imac
         write(*,*) '====================='
         write(*,*) 'autocorrelation for b(',j,')'
         call dacf(nrep,xx,iprint,iseopt,imean,xmean,maxlag,
     1             acv,ac,seac)
 301    continue
         write(*,*) 'autocorrelation for sig2'
         call dacf(nrep,bbbs,iprint,iseopt,imean,xmean,maxlag,
     1             acv,ac,seac)

 600    continue
       call rnget( iseed )
       T_end = cpsec()
       T_hit=T_end-T_start
       write(*,*) 'calculate cpo'
       sumlcpo=0.0d0
       icount=0
       do 400 k=1,207
        if ( (k .eq. 1) .or. (k .eq. icount) ) then
         write(*,*) 'k=',k
         icount=icount+10
        endif
        sum1=0.0d0 ! for summing cpo
        sum3=0.0d0 ! for mu
        sum4=0.0d0 ! for sig2
        do 410 i=1,nrep
         sum2=0.0d0
         do 420 j=1,10
          sum2=sum2+ bbb(j,i)*x(j,k)
 420     continue
         af=dexp( -(y(k)-sum2)**2/(2.0d0*bbbs(i))
     1            -0.5d0*dlog(2.0d0*api*bbbs(i)) )
         sum1=sum1+1.0d0/af
         sum3=sum3+sum2/af
         sum4=sum4+bbbs(i)/af
 410    continue
        cpo(k)=1.0d0/(sum1/dfloat(nrep))
        amu(k)=cpo(k)*sum3/dfloat(nrep)
        av(k)=cpo(k)*sum4/dfloat(k)
        ad(k)=(y(k)-amu(k))/dsqrt(av(k))
        ar(k)=y(k)-amu(k)
        sumlcpo=sumlcpo+dlog( cpo(k) )
 400   continue
c
c
c......print estimation results
       open( unit=11, file='bres.out',
     1       access='sequential',status='unknown')
       write(11,*) ' Using gibbs sampler'
       write(11,*) ' sample size n          = ', nrep
       write(11,*) ' the initial random number seed =', iseed1
       write(11,*) ' the final   random number seed =', iseed
       write(11,*) ' the initial state  and last state'
       write(11,*) 'runtime = ', T_hit
       write(11,*) 'sum of log cpo = ', sumlcpo
       write(11,*) ' the estimator of sigma^2 with s.e.'
       write(11,997)    sigest,sigsd
 997   format(1x,f16.6,'(',f16.6,')')
       write(11,*) ' the estimators of beta(10) and its s.e.'
       write(11,*)  ' est             s.e. '
       do 350 j =1, 10
       write(11,311) j, best(j), bsd(j) 
 311   format(1x, I4,4f16.7 )
 350   continue
       close(11)

       open( unit=16, file='cpo.dat',
     1       access='sequential',status='unknown')
       do 360 i=1,207
       write(16,1999) i,cpo(i),amu(i),av(i),ad(i),ar(i)
 1999  format(I6,1x,f16.6,1x,f16.6,1x,f16.6,1x,f16.6,1x,f16.6)
 360   continue
       close(16)     
       stop
       end

      include 'gibbs1.f'
      include 'tnorm.f'


        subroutine datain(x,y)
c       read in sample data
c       input
c         x(10,207): total trees
c         y(207)   : total fruit in cartons produced by
c                    trees x(10,207)
c         betaprgm : data file
c       output
c         x, y     : for calculating the bayesian estimation
c       auxiliary file
c         betaprgm : data file
        implicit real*8 (a-h,o-z)
        real*8 x(10,250), y(250)
        open(unit=14, file='betaprgm', status='old')
        do 100 i=1,207
        read(14,*) y(i), (x(j,i),j=1,10)
 100    continue
        close(14)
        return
        end
