program threelevel
c       for three level cumulative generalized linear model
c       mgmc algorithm
c       ming-hui chen
c       discussion paper joint with chuanhai liu
c       may 9, 1998
c
        implicit real*8 (a-h,o-z)
        parameter( n=8 )
        parameter (iprint=3,maxlag=20,nobs=50000)
        real*8 z(n),x(n),beta(2),gam(2)
        real*8 y(n)
        real*8 seqgam(50000),seqbeta(2,50000)
        real*8 tau
        real*8 sbeta(2),s2beta(2),ebeta(2),sebeta(2)
        real*8 egam,segam,sgam,s2gam
        real*8  ac(maxlag+1),acv(maxlag+1),seac(maxlag+1),
     1          xac(nobs+1), xacmean
        integer iseed,nrep
        integer imean,iseopt
        common  /vbeta/beta
        common  /vgam/gam
        common  /vx/x
        common  /vy/y
        common  /vz/z
        common  /vtau/tau
        external dacf

        tau=0.001d0
c.......data set
        y(1)=3.0d0
        x(1)=1.0d0
        y(2)=2.0d0
        x(2)=0.0d0
        y(3)=3.0d0
        x(3)=1.0d0
        y(4)=1.0d0
        x(4)=0.0d0
        y(5)=3.0d0
        x(5)=1.0d0
        y(6)=1.0d0
        x(6)=0.0d0
        y(7)=3.0d0
        x(7)=1.0d0
        y(8)=3.0d0
        x(8)=0.0d0
                                                                      
        print*,'enter number of interations nrep='
        read(*,*) nrep

        imean=1
        iseopt=1

c........set initial value
        iseed=9999999
        gam(1)=0.0d0
        gam(2)=1.0d0
        beta(1)=0.0d0
        beta(2)=1.0d0
        call rnset( iseed )

c.......warm up gibbs
        icount=0
        sbeta(1)=0.0d0
        sbeta(2)=0.0d0
        s2beta(1)=0.0d0
        s2beta(2)=0.0d0
        do 100 i=1,5000
         call gibbs( iseed )
         if ( (i .eq. 1) .or. (i .eq. icount)) then
          write(*,*) 'warm up: i=',i
          write(*,800) beta(1),beta(2)
 800      format('beta=',2f16.6)
          write(*,801) gam(2)
 801      format('gam(2)=',f16.6)
          icount=icount+100
         endif
 100    continue
        icount=0
        do 110 i=1,nrep
           call gibbs( iseed )
           if ( (i .eq. 1) .or. (i .eq. icount)) then
          write(*,*) 'Gibbs: i=',i
          write(*,800) beta(1),beta(2)
          write(*,801) gam(2) 
          icount=icount+100
         endif
         do 120 j=1,2 
          seqbeta(j,i)=beta(j)
          sbeta(j)=sbeta(j)+beta(j)
          s2beta(j)=s2beta(j)+beta(j)**2
 120     continue
         seqgam(i)=gam(2)
         sgam=sgam+gam(2)
         s2gam=s2gam+gam(2)**2
 110    continue
c........to obtain posterior estimates
         do 140 j=1,2
          ebeta(j)=sbeta(j)/real(nrep)
          sebeta(j)=dsqrt( (s2beta(j)-real(nrep)*ebeta(j)**2)
     1                     /real(nrep-1) )
 140     continue
         egam=sgam/real(nrep)
         segam=dsqrt( (s2gam-real(nrep)*egam**2)
     1                     /real(nrep-1) )
         open(unit=14,file='ca_mg.dat',access='sequential',
     1       status='unknown')
         do 500 i=1,nrep
           write(14,990) i,seqbeta(1,i),seqbeta(2,i),seqgam(i)
 990       format(1x,I6,3f16.6)
 500     continue
         close(14)
c........compute autocorrelation          
c
c........for beta(1)
         write(*,*) 'for beta(1)'
         do 300 i=1,nrep
           xac(i)=seqbeta(1,i)
 300     continue
         imean=1
         iseopt=1
         call dacf(nrep,xac,iprint,iseopt,imean,xacmean,maxlag,
     1             acv,ac,seac)
c........for beta(2)
         write(*,*) 'for beta(2)'
         do 310 i=1,nrep
           xac(i)=seqbeta(2,i)
 310     continue
         call dacf(nrep,xac,iprint,iseopt,imean,xacmean,maxlag,
     1             acv,ac,seac)
c........for beta(1)+beta(2)
         write(*,*) 'for beta(1)+beta(2)'
         do 320 i=1,nrep
           xac(i)=seqbeta(2,i)+seqbeta(1,i)
 320     continue
         call dacf(nrep,xac,iprint,iseopt,imean,xacmean,maxlag,
     1             acv,ac,seac)
c........for gam
         write(*,*) 'for gam'
         do 330 i=1,nrep
           xac(i)=seqgam(i)
 330     continue
         call dacf(nrep,xac,iprint,iseopt,imean,xacmean,maxlag,
     1             acv,ac,seac)


        open(unit=12,file='ca_mg.out',access='sequential',
     1       status='unknown')
        write(12,*) 'nrep=',nrep
        write(12,*) 'ebeta(1)=',ebeta(1),' se=',sebeta(1)
        write(12,*) 'ebeta(2)=', ebeta(2),' se=',sebeta(2)
        write(12,*) 'egam=',egam,' se=',segam
        close(12)

        stop
        end



        subroutine gibbs( iseed )
        implicit real*8 (a-h,o-z)
        parameter( n=8 )
        real*8 z(n),x(n),beta(2),gam(2)
        real*8 y(n)
        real*8 tau
        integer iseed 
        common  /vbeta/beta
        common  /vgam/gam
        common  /vx/x
        common  /vz/z
        common  /vy/y
        common  /vtau/tau
c.......generate z
        call Gz( iseed )
c.......generate gam 
        call Ggam( iseed )
c.......generate b(i)
        call Gbeta( iseed )
c.......adjusted step
        call Gmg( iseed )
c.......ca-adjusted step
        call Gca_mg( iseed )
        return
        end


        subroutine Gz( iseed )
        implicit real*8 (a-h,o-z)
        parameter( n=8 )
        real*8 z(n),x(n),beta(2),gam(2)
        real*8 y(n)
        real*8 tau
        real*8 aleft,aright,amean,tn,sig
        integer iseed
        logical la,lb
        common  /vbeta/beta
        common  /vgam/gam
        common  /vx/x
        common  /vy/y
        common  /vz/z
        common  /vtau/tau
        external YTUVN
        la=.false.
        lb=.false.

        call rnset( iseed )
        do 300 i=1,n
         amean=beta(1)+beta(2)*x(i)
         sig=1.0d0
         if (y(i) .eq. 1.0d0) then
          la=.true.
          aright=-amean/sig
          tn=YTUVN(aleft,aright,la,lb,iseed)
          z(i)=amean+tn*sig
          la=.false.
         endif
         if (y(i) .eq. 2.0d0) then
          aleft=-amean/sig
          aright=(gam(2)-amean)/sig
          tn=YTUVN(aleft,aright,la,lb,iseed)
          z(i)=amean+tn*sig
         endif
         if (y(i) .eq. 3.0d0) then
          lb=.true.
          aleft=(gam(2)-amean)/sig
          tn=YTUVN(aleft,aright,la,lb,iseed)
          lb=.false.
          z(i)=amean+tn*sig
         endif
 300    continue

        call rnget( iseed )
        return
        end


        subroutine Ggam( iseed )
        implicit real*8 (a-h,o-z)
        parameter( n=8 )
        real*8 z(n),x(n),beta(2),gam(2)
        real*8 y(n)
        real*8 tau
        real*8 scale,shape,r(1)
        real*8 sum
        real*8  alow,aupp
        real*8  amaxz(1:3),aminz(1:3),u
        integer iseed
        common  /vbeta/beta
        common  /vgam/gam
        common  /vx/x
        common  /vy/y
        common  /vz/z
        common  /vtau/tau
        external DRNUNF
        index1=0
        index2=0
        do 400 i=1,n
         if (y(i) .eq. 3.0d0) then
          index1=index1+1
          if (index1 .eq. 1) then
           aupp=z(i)
          else
           if (aupp .gt. z(i)) aupp=z(i)
          endif
         endif
         if (y(i) .eq. 2.0d0) then
          index2=index2+1
          if (index2 .eq. 1) then
           alow=z(i)
          else 
           if (alow .lt. z(i)) alow=z(i)
          endif
         endif
 400    continue
        call rnset( iseed )
          u=DRNUNF()
        call rnget( iseed )
        gam(2)=alow+u*(aupp-alow)
        return
        end


        subroutine Gbeta( iseed )
        implicit real*8 (a-h,o-z)
        parameter( n=8 )
        real*8 z(n),x(n),beta(2),gam(2)
        real*8 y(n)
        real*8 tau
        real*8 scale,shape,r(1)
        real*8 sum
        real*8 xvec(2)
        real*8 xx(1:2,1:2),Bm(1:2,1:2)
        real*8 Binv(1:2,1:2)
        real*8 Binv1(1:2,1:2)
        real*8 Siginv(1:2,1:2)
        real*8 bhat(2),bhat1(2)
        real*8 tol,Rsig(1:2,1:2),RV(1,2)
        integer iseed,irank,lda,ldr,nr
        common  /vbeta/beta
        common  /vgam/gam
        common  /vx/x
        common  /vy/y
        common  /vz/z
        common  /vtau/tau
        external DMACH,DRNMVN
c.......initialize
        do 180 i=1,2
         do 185 j=1,2
          Bm(i,j)=0.0d0
 185     continue
 180    continue
        do 190 i=1,2
         bhat(i)=0.0d0
         Bm(i,i)=tau
 190    continue
        do 200 i=1,n
c........calculate xx(2,2)
         xvec(1)=1.0d0
         xvec(2)=x(i)
         do 209 i1=1,2
          do 210 j1=1,2
           xx(i1,j1)=xvec(i1)*xvec(j1)
 210      continue
 209     continue
         do 250 j=1,2
          do 255 j1=1,2
            Bm(j,j1)=Bm(j,j1)+xx(j,j1)
 255      continue
 250     continue
         do 220 j=1,2
           bhat(j)=bhat(j)+xvec(j)*z(i)
 220     continue
 200    continue
c.......calculate Binv
        call DLINDS(2,Bm,2,Binv,2)
c.......calculate bhat1
        do 270 i=1,2
         bhat1(i)=0.0d0
         do 275 j=1,2
          bhat1(i)=bhat1(i)+Binv(i,j)*bhat(j)
          Binv1(i,j)=Binv(i,j)
 275     continue
 270    continue
c.......do Cholesky factorization
        lda=2
        ldr=2
        tol=100.0d0*DMACH(4)
        call DCHFAC(2,Binv1,lda,tol,irank,Rsig,ldr)
c.......generate multivariate normal
        nr=1
        call rnset( iseed )
         call DRNMVN(nr,2,Rsig,2,RV,1)
        call rnget( iseed )
        do 300 j=1,2
         beta(j)=RV(1,j)+bhat1(j)
 300    continue
        return
        end

        subroutine Gmg( iseed )
        implicit real*8 (a-h,o-z)
        parameter( n=8 )
        real*8 z(n),x(n),beta(2),gam(2)
        real*8 y(n)
        real*8 tau
        real*8 scale,shape,r(1)
        integer iseed
        common  /vbeta/beta
        common  /vgam/gam
        common  /vx/x
        common  /vy/y
        common  /vz/z
        common  /vtau/tau
        external DRNGAM
        shape=real(n+1+2)/2.0d0
        sum=0.0d0
        do 400 i=1,n
         sum=sum+(z(i)-beta(1)-beta(2)*x(i))**2
 400    continue
        scale=sum/2.0d0+tau*(beta(1)**2/2.0d0
     1                 +beta(2)**2/2.0d0)
c
        call rnset( iseed )
c
          call DRNGAM(1,shape,r)
c
        call rnget( iseed )
c
        g=dsqrt(r(1)/scale)
c.......adjusted step
        beta(1)=g*beta(1)
        beta(2)=g*beta(2)
        gam(2)=g*gam(2)
        do 410 i=1,n
         z(i)=g*z(i)
 410    continue
        return
        end

        subroutine Gca_mg( iseed )
        implicit real*8 (a-h,o-z)
        parameter( n=8 )
        real*8 z(n),x(n),beta(2),gam(2)
        real*8 y(n)
        real*8 tau
        real*8 scale,shape,r(1)
        real*8 zstar(4)
        integer iseed
        logical la,lb
        common  /vbeta/beta
        common  /vgam/gam
        common  /vx/x
        common  /vy/y
        common  /vz/z
        common  /vtau/tau
        external DRNNOF,YTUVN
        la=.false.
        lb=.false.
        sum=0.0d0
        num=0
        sig2o=1.0d0/tau
c       alpha=dsqrt(2.0d0)
        alpha=1.0d0
        do 100 i=1,n
         if (x(i) .eq. 1.0d0) then
          sum=sum+z(i)
          num=num+1
         endif
 100    continue
        tstar=sum/real(num)
        ind=0
        do 110 i=1,n
         if (x(i) .eq. 1.0d0) then
          ind=ind+1
          zstar(ind)=z(i)-tstar
         endif
 110    continue
        do 120 j=1,4
         if (j .eq. 1) then
          alow=gam(2)-zstar(j)
         else
          if (alow .lt. gam(2)-zstar(j)) then
            alow=gam(2)-zstar(j)
          endif
         endif
 120    continue
c.......generate T
        amean=alpha*beta(1)
        std=dsqrt(0.25d0+sig2o)
        lb=.true.
        aleft=(alow-amean)/std
        tn=YTUVN(aleft,aright,la,lb,iseed)
        lb=.false.
        T=amean+std*tn
c.......generate amu=beta(1)+beta(2)
        amean=alpha*beta(1)+
     1     sig2o*(T-alpha*beta(1))/(sig2o+0.25d0) 
        std=sig2o-sig2o**2/(sig2o+0.25)
        std=dsqrt(std)
        call rnset( iseed )
         tn=DRNNOF()
        call rnget( iseed )
        amu=amean+tn*std
c.......adjusting stage
        beta(2)=amu-beta(1)
        ind=0
        do 210 i=1,n
         if (x(i) .eq. 1.0d0) then
          ind=ind+1
          z(i)=zstar(ind)+T
         endif
 210    continue
        return
        end

        include 'tnorm.f'