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