c/* ****************************************************************** */
c/*                                                                    */
c/*                                                                    */
c/**** Generate sigma_e and sigma_u                                 ****/
c/*                                                                    */
c/*                                                                    */
c/* ****************************************************************** */
       subroutine findsesu(sige,sigu,n,nueps,msamp,nuuu,aeps,
     6              yy,beta0,beta1,xx,ww,auu,dummyx,dummyw,numw,
     6              rnseed)
       implicit real*8(a-h,o-z)
       real*8 xx(n,1),ww(msamp,numw),yy(n,1)
       real*8 dummyx(n),dummyw(msamp,numw)
       sedf       = dfloat(n) + dfloat(nueps)
       sudf       = dfloat(msamp) + msamp + nuuu
       sescale = aeps
       do 100 i=1,n
 100   sescale = sescale + 
     6            ((yy(i,1) - beta0 - (beta1 * xx(i,1))) ** 2)
       suscale = auu
       do 200 i=1,msamp
          do 150 j=1,numw
 150         suscale = suscale + ((ww(i,j) - xx(i,1)) ** 2)
 200         continue
       sige       = sichi(sedf,sescale,rnseed)
       sigu       = sichi(sudf,suscale,rnseed)
       return
       end





c/* ****************************************************************** */
c     Generate an inverse chi random variable
c/* ****************************************************************** */
      double precision function sichi(df,scale,rnseed)
      implicit real*8(a-h,o-z)
      real*8 aa(1)
      dfd2 = df / 2.0
      call largamma(dfd2,rnseed,bb)
      sichi = dsqrt(scale / (2 * bb))
c      call rgamma(rnseed,1,aa,1.0d0,dfd2)
c      aa(1)   = 2 * aa(1)
c      sichi   = dsqrt(scale / aa(1))
      return
      end




c/* ****************************************************************** */
c/*                                                                    */
c/*                                                                    */
c/**** Generate beta_0 and beta_1                                   ****/
c      des       = ones(n,1)~xx;
c      sigma     = (sige .^ 2) .* inv(des' * des);
c      vv        = (chol(sigma))';
c      beta      = (inv(des' * des) * des' * yy) + (vv * rndns(2,1,seed));
c/*                                                                    */
c/*                                                                    */
c/* ****************************************************************** */
      subroutine findbeta(rnseed,n,x,y,sige,designx,
     6                      beta0,beta1,ssbprior)
      implicit real*8(a-h,o-z)
      real*8 x(n,1),y(n,1),designx(n,2)
      real*8 xpy(2,1),xpxinv(2,2),xpx(2,2),betahat(2,1)
      real*8 vv(2,2),vvnorm(2,1),eeer(2,1)
        zkappa = (sige ** 2) / ssbprior
        call zeros(xpx,2,2)
        call zeros(xpy,2,1)
        do 200 i=1,n
        designx(i,1) = 1
        designx(i,2) = x(i,1)
        do 150 j=1,2
        xpy(j,1) = xpy(j,1) + (designx(i,j) * y(i,1))
        do 125 k=1,2
 125       xpx(j,k) = xpx(j,k) + (designx(i,j) * designx(i,k))
 150    continue
 200    continue
c
c       Make the adjustment for the prior in beta_1
        xpy(2,1) = zkappa + xpy(2,1)
        xpx(2,2) = zkappa + xpx(2,2)
c       Get inverse of xpx
        aa = (xpx(1,1) * xpx(2,2)) - (xpx(1,2) ** 2)
        xpxinv(1,1) = xpx(2,2) / aa
        xpxinv(2,2) = xpx(1,1) / aa
        xpxinv(1,2) = -xpx(1,2) / aa
        xpxinv(2,1) = -xpx(2,1) / aa
c        call xinverse(xpx,xpxinv,2,bigdet) 
c       Get betahat
        call zmmultm(xpxinv,xpy,betahat,2,2,1)
c       Do the cholesky decomposition of xpxinv and multiply
c            it by sigma_e to get the matrix square root
c            of sigma_e^2 inv(x'*x)
c       This is a two step process. You first call mchol,
c            then set the top right element of xcholmat  equal to 0,
c            to get the lower trianglar matrix.
c            You then multiply xcholmat by the square root of
c            dcholmat times sige to get the matrix square
c            root vv
c        call mchol(xpxinv,xcholmat,dcholmat,ier)
c        xcholmat(1,2) = 0
c        call zeros(ddchol,2,2)
c        ddchol(1,1) = sige * dsqrt(dcholmat(1))
c        ddchol(2,2) = sige * dsqrt(dcholmat(2))
c        call zmmultm(xcholmat,ddchol,vv,2,2,2)
        vv(1,1) = dsqrt(sige * sige * xpxinv(1,1))
        vv(1,2) = 0.
        vv(2,1) = sige * sige * xpxinv(2,1) / vv(1,1)
        vv(2,2) = dsqrt( (sige * sige * xpxinv(2,2)) - (vv(2,1) ** 2))
c       You are now in a position to get the new estimates
c           of beta, according to the gauss formula
c      beta      = (inv(des' * des) * des' * yy) + (vv * rndns(2,1,seed));
        call rnorm(rnseed,2,vvnorm)
        call zmmultm(vv,vvnorm,eeer,2,2,1)
        beta0 = betahat(1,1) + eeer(1,1)
        beta1 = betahat(2,1) + eeer(2,1)
        return
        end





c/* ****************************************************************** */
c/*                                                                    */
c/*                                                                    */
c/**** Compute the current mean and variance of marginal X          ****/
c/*                                                                    */
c/*                                                                    */
c/* ****************************************************************** */
         subroutine findvxex(kdim,pmix,zmumix,sigmix,k,expx,varx)
         implicit real*8(a-h,o-z)
         real*8 pmix(kdim,1),zmumix(kdim,1),sigmix(kdim,1)
         varx       = 0.
         expx       = 0.
         do 100 j=1,k
            varx       = varx + (pmix(j,1) * (sigmix(j,1) ** 2))
     6             + (pmix(j,1) * ( zmumix(j,1) ** 2))
 100        expx       = expx + (pmix(j,1) * zmumix(j,1)) 
            varx       = varx - (expx ** 2)
            return
            end