program testinghpd 
c      a test program for using subroutine hpd.f  
c      input
c         hpd.f
c      ming-hui chen, wpi
c      july 23, 2001
       implicit real*8 (a-h,o-z)
       real*8 x(50000)
       real*8 shape,scale
       real*8 r(1)
       real*8 alpha,alow(2),aupp(2) ! variables used in hpd.f
       integer iseed  
       external DRNGAM  ! an IMSL sobroutine for generating a gamma random variable  

c......read the number of Gibbs iterates
       print*,'Please enter sample size n (<=50000): n='
       read(*,*) n
       print*,'Please enter scale parameter for gamma: scale='
       read(*,*) scale 
       print*,'Please enter shape parameter for gamma: shape='
       read(*,*) shape 
       print*,'Please enter confidence level (between 0 and 1): alpha='
       read(*,*) alpha 

       iseed = 9999999

       do 100 i=1,n
        call rnset( iseed )
         call DRNGAM(1,shape,r)
        call rnget( iseed )
        x(i)=r(1)/scale
 100   continue

       call hpd(n,alpha,x,alow,aupp)

       write(*,900) nint(100.0d0*(1.0d0-alpha)),alow(1),aupp(1)
 900   format( I4, '%hpd interval:',' (',f10.3,',',f10.3,')') 
       write(*,901) nint(100.0d0*(1.0d0-alpha)),alow(2),aupp(2)
 901   format( I4, '%credible interval:',' (',f10.3,',',f10.3,')') 

       stop
       end
      
       include 'hpd.f'