subroutine gibbloop(xx,xxest,ig,smixest,pmixest,
     6              zmixest,xbar,s,n,Mprior,nuprior,rnseed,
     6              maxmix,Kbigest,njmatrix,iiest,
     6              kklow,kkhigh,kkdelta,kkgibbs)
           implicit real*8(a-h,o-z)
           real*8 xx(n),xxest(n,1),xbar(5),s(5)
           real*8 pmixest(maxmix,1),smixest(maxmix,1),zmixest(maxmix,1)
           integer ig(n),njmatrix(maxmix,1),iiest(n,maxmix)

c          /* ******** */
c          /* update A */
c          /* ******** */
          temp = 0.0 
          do 40 j=1,Kbigest
 40          temp = temp + (1.0 / (smixest(j,1) ** 2))
          Amat = rchisq(Kbigest*nuprior,rnseed) / temp 
c          if(kkgibbs.ge.1000) print *, "Amat, kkgibbs = ",kkgibbs
c          /* ******** */
c          /* update p */
c          /* ******** */
          a1 = 0.0 
          do 50 j=1,Kbigest
               a2 = qqqgamma(1.0d0+njmatrix(j,1),rnseed) 
               pmixest(j,1) = a2
 50            a1 = a1 + a2 
          do 60 j=1,Kbigest
 60          pmixest(j,1) = pmixest(j,1) / a1
c          if(kkgibbs.ge.1000) print *, "pmix, kkgibbs = ",kkgibbs
c         /* ***************** */
c         /* update the sigmas */
c         /* ***************** */
          if (Kbigest.eq.1) then
          aa = amat
          do 65 i=1,n
 65          aa = aa + ( (xxest(i,1) - zmixest(1,1)) ** 2)
           smixest(1,1) = dsqrt(aa/rchisq(n + nuprior,rnseed)) 
           endif
           if (Kbigest.ge.2) then
          do 300 j=1,Kbigest
             if (j.eq.1) then
                r = (zmixest(1,1) - zmixest(2,1)) *
     6              (zmixest(1,1) - zmixest(2,1)) /
     6              (2.0 * Mprior*Mprior) 
                endif
             if (j.eq.Kbigest) then
                r = (zmixest(Kbigest,1) - zmixest(Kbigest-1,1)) *
     6              (zmixest(Kbigest,1) - zmixest(Kbigest-1,1)) /
     6              (2.0 * Mprior*Mprior) 
                endif
             if((j.gt.1).and.(j.lt.Kbigest)) then
                r = (((zmixest(j,1) - zmixest(j-1,1)) ** 2) /
     6              (2.0 * Mprior*Mprior) )
     6             +
     6               (((zmixest(j,1) - zmixest(j+1,1)) ** 2) /
     6              (2.0 * Mprior*Mprior)) 
                    endif
               u = ((njmatrix(j,1)-1)*s(j)*s(j))
     6                     +Amat+r+
     6               (njmatrix(j,1)*((xbar(j)-zmixest(j,1))**2))
           cand = dsqrt(u/rchisq(njmatrix(j,1)+ nuprior,rnseed)) 

             if (j.eq.1) then
                    temp1 = (1.0/(cand*cand)) + 
     6                  (1.0/(smixest(2,1)*smixest(2,1))) 
                    temp2 = (1.0/(smixest(1,1)*smixest(1,1))) +
     6                  (1.0/(smixest(2,1)*smixest(2,1))) 
                    ratio = 0.5*dlog(temp1) - 0.5*dlog(temp2) 
               endif
             if (j.eq.Kbigest) then
           temp1 = (1.0/(cand*cand)) + 
     6           (1.0/(smixest(Kbigest-1,1)*smixest(Kbigest-1,1))) 
           temp2 = (1.0/(smixest(Kbigest,1)*smixest(Kbigest,1))) 
     6        +  (1.0/(smixest(Kbigest-1,1)*smixest(Kbigest-1,1))) 
                    ratio = 0.5*dlog(temp1) - 0.5*dlog(temp2) 
                    endif
             if((j.gt.1).and.(j.lt.Kbigest)) then
             temp1 = (1.0/(cand*cand))+
     6                  (1.0/( smixest(j+1,1)*smixest(j+1,1)))
             temp2 = (1.0/(cand*cand))+
     6                  (1.0/( smixest(j-1,1)*smixest(j-1,1))) 
             temp3 = (1.0/( smixest(j,1)*smixest(j,1)) ) +
     6                  (1.0/( smixest(j+1,1)*smixest(j+1,1)))
             temp4 = (1.0/( smixest(j,1)*smixest(j,1)) )+
     6                  (1.0/( smixest(j-1,1)*smixest(j-1,1))) 
                    ratio = 0.5*dlog(temp1) + 0.5*dlog(temp2)
     6                       -0.5*dlog(temp3) - 0.5*dlog(temp4) 
                            endif
          ratio = expp(ratio) 
          u = qqqunif(rnseed) 
          if (u.lt.ratio) smixest(j,1) = cand 
 300      continue
          endif
c          if(kkgibbs.ge.1000) print *, "smix, kkgibbs = ",kkgibbs
c          /* *************** */
c          /* update the mu's */
c          /* *************** */
          if (Kbigest.eq.1) then
          temp = qqqrnorm(rnseed) 
          zmixest(1,1) = xbar(1) + 
     6              (smixest(1,1) * temp / dsqrt(dfloat(n)))
          endif
          if (Kbigest.ge.2) then
          temp1 = smixest(1,1)*smixest(1,1) 
          temp2 = smixest(2,1)*smixest(2,1) 
          v = 2.0/(  (1.0/temp1) + (1.0/temp2) ) 
          a1 = njmatrix(1,1)/temp1 
          a2 = 1.0/(Mprior*Mprior*v) 
          b  = ((xbar(1)*a1) + (zmixest(2,1)*a2))/(a1+a2) 
          c  = 1.0/dsqrt(a1+a2) 
c          print *, " "
c          print *, "temp1,temp2,v,a1,a2,b,c,zmixest(2,1),rnseed=",
c     6              temp1,temp2,v,a1,a2,b,c,zmixest(2,1),rnseed
          zmixest(1,1) = rnorm1(b,c,zmixest(2,1),rnseed) 
c           print *, "zmixest(1,1) = ",zmixest(1,1)
c          print *, " "
          if (Kbigest.ge.3) then 
          do 400 j=2,(Kbigest-1)
               temp1 = smixest(j,1)*smixest(j,1) 
               temp2 = smixest(j-1,1)*smixest(j-1,1) 
               v = 2.0/(  (1.0/temp1) + (1.0/temp2) ) 
               temp3 = smixest(j+1,1)*smixest(j+1,1) 
               w = 2.0/(  (1.0/temp1) + (1.0/temp3) ) 
               a1 = njmatrix(j,1)/temp1 
               a2 = 1.0/(Mprior*Mprior*v) 
               a3 = 1.0/(Mprior*Mprior*w) 
               b  = ((xbar(j)*a1) + (zmixest(j-1,1)*a2)
     6               + (zmixest(j+1,1)*a3))/
     6               (a1+a2+a3) 
               c  = 1.0/dsqrt(a1+a2+a3) 
               zmixest(j,1) = rnorm2(b,c,zmixest(j-1,1),
     6                             zmixest(j+1,1),rnseed) 
c           print *, "j and zmixest(j,1) = ",j,zmixest(j,1)
 400           continue
               endif
          temp1 = smixest(Kbigest,1)*smixest(Kbigest,1) 
          temp2 = smixest(Kbigest-1,1)*smixest(Kbigest-1,1) 
          w = 2.0/(  (1.0/temp1) + (1.0/temp2) ) 
          a1 = njmatrix(Kbigest,1)/temp1 
          a2 = 1.0/(Mprior*Mprior*w) 
          b  = (xbar(Kbigest)*a1 + zmixest(Kbigest-1,1)*a2)
     6                /(a1+a2) 
          c  = 1.0/dsqrt(a1+a2) 
c          print *, "xbar = ",(xbar(j),j=1,Kbigest)
          dddd = zmixest(Kbigest-1,1)
          zmixest(Kbigest,1) = rnorm3(b,c,dddd,
     6                                    rnseed) 
c           print *, "zmixest(Kbigest,1) = ",zmixest(Kbigest,1)
          endif
c          if(kkgibbs.ge.1000) print *, "zmix, kkgibbs = ",kkgibbs
c          /* ************************ */
c          /* update membership vector */
c          /* ************************ */
          if (Kbigest.ge.2) then
          do 500 ii=kklow,kkhigh
               if(ig(ii).eq.1) icand = 2 
               if(ig(ii).eq.Kbigest)icand = Kbigest-1
            if((ig(ii).gt.1).and.(ig(ii).lt.Kbigest)) then
                    u = qqqunif(rnseed) 
                    if(u.lt.0.5d0) icand = ig(ii)-1 
                    if(u.ge.0.5d0) icand = ig(ii)+1 
              endif
       ratio = dlog(pmixest(icand,1))+
     6        zlnorm(xx(ii),zmixest(icand,1),smixest(icand,1))-
     6        dlog(pmixest(ig(ii),1))-
     6        zlnorm(xx(ii),zmixest(ig(ii),1),smixest(ig(ii),1)) 
       if((ig(ii).eq.1).or.(ig(ii).eq.Kbigest))
     6                       ratio = ratio + dlog(0.5d0) 
       if((icand.eq.1).or.(icand.eq.Kbigest))
     6                       ratio = ratio + dlog(2.0d0) 
       ratio = expp(ratio) 
       u = qqqunif(rnseed) 
       if(u.lt.ratio)ig(ii) = icand 
 500   continue
       endif
c          if(kkgibbs.ge.1000) print *, "ig, kkgibbs = ",kkgibbs

c          /* ******************************* */          
c          /* get new values of xbar, s and njmatrix 
c             and iiest based on new groups */
c          /* ******************************* */          
       call summary(xx, n, ig, Kbigest,xbar,s,
     5           njmatrix,iiest,maxmix)
c          if(kkgibbs.ge.1000) print *, "summary, kkgibbs = ",kkgibbs
c      *********************************************************
c      Update kklow and kkhigh
c      *********************************************************
        kklow  = kklow  + kkdelta
        kkhigh = kkhigh + kkdelta
        if (kklow.gt.n) then
              kklow  = 1
              kkhigh = kkdelta
              endif
        if ((kklow.le.n).and.(kkhigh.gt.n)) kkhigh = n
c          if(kkgibbs.ge.1000) print *, "kklow, kkgibbs = ",kkgibbs
        return
        end


c     ***********************************************************
c     The function summary from Larry
c
c     compute number of observations, mean and stan dev of xx
c      grouped by g into at most k groups.
c      xbar, s and njmatrix get modified 
c     ***********************************************************
       subroutine summary(xx, n, ig, Kbigest,xbar,s,
     5           njmatrix,iiest,maxmix)
       implicit real*8(a-h,o-z)
       real*8 xx(n),xbar(5),s(5)
       integer njmatrix(maxmix,1),ig(n),iiest(n,maxmix)
       do 100 j=1,Kbigest
          xbar(j) = 0.
          s(j)    = 0.
          njmatrix(j,1) = 0
          do 80 i=1,n
 80          iiest(i,j) = 0
 100      continue
       do 200 i=1,n
          xbar(ig(i))         = xbar(ig(i)) + xx(i) 
          s(ig(i))            = s(ig(i)) + (xx(i)*xx(i)) 
          njmatrix(ig(i),1)   = njmatrix(ig(i),1) + 1 
          iiest(i,ig(i))      = 1
 200      continue
       do 300 j=1,Kbigest
          if(njmatrix(j,1).gt.0)xbar(j)=xbar(j) / njmatrix(j,1) 
          if(njmatrix(j,1).gt.1)
     6       s(j) = dsqrt( (s(j)-njmatrix(j,1)*xbar(j)*xbar(j))
     6         /(njmatrix(j,1)-1.0) ) 
          if(njmatrix(j,1).le.1) s(j) = 0.
 300      continue
          return
          end


c     ***********************************************************
c     The function expp from Larry
c     ***********************************************************
       double precision function expp(r)
       implicit real*8(a-h,o-z)
       if(r.gt.0.0d0)out = 1.0 
       if(r.lt.-10.0d0)out = 0.0 
       if((r.le.0.0d0).and.(r.ge.-10.0d0))  out = dexp(r) 
       expp = out
       return
       end



c     ***********************************************************
c     Normal Loglikelihood  from Larry
c     ***********************************************************
       double precision function zlnorm(x, zmu, sigma)
       implicit real*8(a-h,o-z)
       out = -0.5*dlog(2*3.141592653d0)
     6         -((x-zmu)*(x-zmu)/(2.0*sigma*sigma))
     6            -dlog(sigma) 
       zlnorm = out
       return
       end


c     ***********************************************************
c     The function rnorm1 from Larry
c     ***********************************************************
       double precision function rnorm1( zmu,sigma,trunc,rnseed)
       implicit real*8(a-h,o-z) 
c       print *, "rnorm1 --> zmu,sigma,trunc,rnseed = ",
c     6           zmu,sigma,trunc,rnseed
       aa = (trunc - zmu) / sigma
       out = pnorm(  (trunc-zmu)/sigma ) 
c       out = (out + 0.00000001d0) / (1.00000002d0)
c       print *, "rnorm1 --> out,aa = ",out,aa
       out = out*qqqunif(rnseed) 
c       print *, "rnorm1 --> out = ",out
       out = zmu+sigma*qnorm(out) 
       if (out.ge.trunc) out = trunc-abs(sigma * qqqrnorm(rnseed))
       rnorm1 = out
       return 
       end


c     ***********************************************************
c     The function rnorm2 from Larry
c     ***********************************************************
       double precision function  rnorm2( zmu,  sigma,  trunc1, 
     6                              trunc2,rnseed)
       implicit real*8(a-h,o-z) 
       u1 =  pnorm(  (trunc1-zmu)/sigma ) 
       u2 =  pnorm(  (trunc2-zmu)/sigma ) 
       out = u1 + (u2-u1)*qqqunif(rnseed) 
c       out = (out + 0.00000001d0) / (1.00000002d0)
       out = zmu+sigma*qnorm(out) 
       rnorm2 = out
       return 
       end


c     ***********************************************************
c     The function rnorm3 from Larry
c     ***********************************************************
       double precision function  rnorm3( zmu,  sigma,  
     6                      trunc,rnseed)
       implicit real*8(a-h,o-z) 
       out = pnorm(  (trunc-zmu)/sigma ) 
       out = out + (1.0-out)*qqqunif(rnseed) 
       sspp = qnorm(out)
       out = zmu+(sigma*sspp)
       if (out.le.trunc) out = trunc+abs(sigma * qqqrnorm(rnseed))
       rnorm3 = out
       return 
       end



c     ***********************************************************
c     Function to generate a uniform[0,1]
c     ***********************************************************
       double precision function  qqqunif(rnseed)
       implicit real*8(a-h,o-z)
       real*8 uu(1)
       call runif(rnseed,1,uu)
       qqqunif = uu(1)
       return
       end
c
c     return(((double) rand())/RAND_MAX) 




c     ***********************************************************
c     Generate a normal random variable, from Larry
c     ***********************************************************
       double precision function  qqqrnorm(rnseed)
       implicit real*8(a-h,o-z)
       real*8 uu(1)
       call rnorm(rnseed,1,uu)
       qqqrnorm = uu(1)
       return
       end
c    static int iset=0 
c    static double gset 
c    double fac, rsq, v1, v2 
c
c    if(iset==0) {
c      do {
c          v1 = 2.0*qqqunif(rnseed)-1.0 
c          v2 = 2.0*qqqunif(rnseed)-1.0 
c          rsq = v1*v1 + v2*v2 
c         } while (rsq >= 1.0 ||rsq == 0.0) 
c      fac = dsqrt (-2.0*dlog(rsq)/rsq) 
c      gset=v1*fac 
c      iset=1 
c      return v2*fac 
c  } else{
c      iset=0 
c      return gset 
c      }
c
c  }


c     ***********************************************************
c     Generate a chi-squared with idf degrees of freedom
c     ***********************************************************
       double precision function  rchisq(idf,rnseed)
       implicit real*8(a-h,o-z)
       out = 0.
       do 100 i=1,idf
          temp = qqqrnorm(rnseed) 
          out  = out + (temp*temp)
 100      continue
       rchisq = out
       return
       end

c     ***********************************************************
c     Normal Density from Larry
c     ***********************************************************
       double precision function dnorm(x)
       implicit real*8(a-h,o-z)
       out = dexp(zlnorm(x,0.0d0,1.0d0))
       dnorm = out 
       return
       end


c     ***********************************************************
c     Normal cdf from Larry
c     ***********************************************************
       double precision function pnorm(x)    
       implicit real*8(a-h,o-z)
       if(x.ge.0.0d0) then
       tt =  1.0/(1.0+ 0.2316419*x) 
       out = 0.319381530 * tt   -
     6      0.356563782 * tt*tt +
     6      1.781477937 * tt*tt*tt -
     6      1.821255978 * tt*tt*tt*tt +
     6      1.330274429 * tt*tt*tt*tt*tt 
       out  = 1.0 - dnorm(x)*out 
       endif
       if(x.lt.0.0d0) then
       z = -x 
       tt =  1.0/(1.0+ 0.2316419*z) 
       out = 0.319381530 * tt   -
     6      0.356563782 * tt*tt +
     6      1.781477937 * tt*tt*tt -
     6      1.821255978 * tt*tt*tt*tt +
     6      1.330274429 * tt*tt*tt*tt*tt 
       out  = 1.0 - dnorm(z)*out 
       out = 1.0 - out 
       endif 
       pnorm = out
       return
       end

c     ***********************************************************
c     Larry's function rgamma
c     ***********************************************************
       double precision function  qqqgamma(a,rnseed)
       implicit real*8(a-h,o-z)
       b = a-1.0 
       c = 3.0*a - 3.0/4.0 
       iaccept =0 
       do 100 kkjj=1,99999
          u = qqqunif(rnseed) 
          v = qqqunif(rnseed) 
          if(u.eq.0.0d0)u=0.00001 
          if(v.eq.0.0d0)v=0.00001 
          if(u.eq.1.0d0)u=0.99999 
          if(v.eq.1.0d0)v=0.99999 
          w = u*(1.0-u) 
          y = dsqrt(c/w)*(u-0.5) 
          x = b+y 
          if(x.ge.0.0d0) then
               z = 64.0*w*w*w*v*v 
               if(z.le.(1.0-2.0*y*y/x)) iaccept=1 
               if(iaccept.eq.0) then
                  if(dlog(z).le.(2.0*(b*dlog(x/b)-y))) iaccept=1 
                  endif
                endif
       if(iaccept.eq.1) go to 200
 100   continue
 200   continue
       qqqgamma = x
       return 
       end



c     ***********************************************************
c       find x such that Pr(Z<x)=q assuming  1-q <=.5 
c     ***********************************************************
        double precision function  qfun( q)
        implicit real*8(a-h,o-z)
        p = 1.0-q 
        y = -dlog(2.0*p) 
        if(p.gt.0.0000001d0) then
          x = dsqrt( ((4.0*y+100.0)*y+205.0)*y*y/
     6               ( ((2.0*y+56.0)*y+192.0)*y+131.0) ) 
        endif
        if(p.le.0.0000001d0) then
        x = dsqrt(   ( ((2.0*y+280.0)*y+572.0)*y)
     6                          /( (y+144.0)*y+603.0)  ) 
        endif
        qfun = x
        return
        end

c     ***********************************************************
c      find x such that Pr(Z<x)=q 
c     ***********************************************************
       double precision function  qnorm(q)
       implicit real*8(a-h,o-z)
       if((q*(1.0-q)).lt.0.00001d0) q = 0.00001*(1.0-q)+0.99999*q 
       if(q.ge.0.5d0) x = qfun(q) 
       if(q.lt.0.5d0) x = -qfun(1.0-q) 
       qnorm = x
       return 
       end




c     ***********************************************************
c     subroutine to fill the first n elements of the real*8 array
c     x with iid U(0,1) random variables
c    
c     The argument dseed is double precison and on input contains
c     a large whole number for a seed. On output it contains a new
c     seed that canbe used as input the next time runif is called
c     ************************************************************
         subroutine runif(dseed,n,x)
      implicit real*8(a-h,o-z)
      real *8 x(n)
      real *8 dseed,svnto5,t31m1
      data svnto5,t31m1/16807.d0,2147483647.d0/
      do 10 i=1,n
        dseed=dmod(svnto5*dseed,t31m1)
 10      x(i)=dseed/t31m1
         return
         end
c

c
c     ***********************************************************
c     subroutine to generate n iid N(0,1) rv's x(1)....x(n)
c     dseed is a double precision seed to be  used in runif
c     ***********************************************************
      subroutine rnorm(dseed,n,x)
      implicit real*8 (a-h,o-z)
      real*8  x(n)
      real*8 dseed
      nn=n-1
      if(2*(n/2).eq.n) nn=n
      if(n.eq.nn) go to 20
c 5    call runif(dseed,2,r)
 5    call runif(dseed,1,r1) 
      call runif(dseed,1,r2) 
      r1=2.*(r1-.5)
      r2=2.*(r2-.5)
      s2=r1*r1+r2*r2
      if(s2.gt.1.) go to 5
      y=dsqrt(-2.*(dlog(s2)/s2))
      x(n)=r1*y
      if(n.eq.1) go to 99
 20   call runif(dseed,nn,x)
      nnb2=nn/2
      do 30 i=1,nnb2  
      i2=2*i
      i1=i2-1
 12   r1=2.*(x(i1)-.5)
      r2=2.*(x(i2)-.5)
      s2=r1*r1+r2*r2
      if(s2.le.1.) go to 15
      call runif(dseed,2,x(i1))
      go to 12
 15   y=dsqrt(-2.*(dlog(s2)/s2))
      x(i1)=y*r1
      x(i2)=y*r2
 30   continue
 99   continue
      return
      end



      subroutine largamma(a,rnseed,x)
      implicit real*8(a-h,o-z)
      real*8 uu(1),vv(1)
      b = a - 1.0
      c = (3.0 * a) - (3.0 / 4.0)
      accept = 0
      do 100 kk=1,9999
       call runif(rnseed,1,uu)
       call runif(rnseed,1,vv)
       u = uu(1)
       v = vv(1)
       if (u.le.0.00001d0) u =0.00001d0
       if (v.le.0.00001d0) v =0.00001d0
       if (u.ge.0.99999d0) u =0.99999d0
       if (v.ge.0.99999d0) v =0.99999d0
          w = u * (1.0-u)
          y = dsqrt(c / w) * (u-0.5)
          x = b+y
          if (x.ge.0) then
               z = 64.0 * w * w * w * v * v
               if (z.le.(1.0 - (2.0 * y * y / x))) accept=1 
               if (accept.lt.0.5) then
          if (dlog(z).le.(2.0 * (b * dlog(x / b)-y))) accept=1  
             endif 
             endif
         if (accept.gt.0.5) goto 200
 100     continue
 200     continue
         return
         end


       subroutine largamm1(a,rnseed,x)
       implicit real*8(a-h,o-z)
       real*8 uu(1)
       if (a.le.0)     print *, "Danger in largamm1, a = ",a
       if (a.gt.1.0d0) print *, "Danger in largamm1, a = ",a
       call runif(rnseed,1,uu)
       u = uu(1)
       if (u.lt.0.0001d0) u = 0.0001d0
       call  largamma(a+1,rnseed,z)
       z = dlog(z) + ( (1.0d0 / a) * dlog(u))
       z = dexp(z)
       return
       end