c/* ****************************************************************** */
c   This is a setup routine
c/* ****************************************************************** */
       subroutine findsig8(sig12,sig21,sig11,beta1,sig,sige,sigu)
       implicit real*8(a-h,o-z)
       real*8 sig11(3,3),sig12(1,3),sig21(3,1)
       sig12(1,1) = beta1 * sig * sig
       sig12(1,2) = sig * sig 
       sig12(1,3) = sig * sig
       sig21(1,1) = sig12(1,1)
       sig21(2,1) = sig12(1,2)
       sig21(3,1) = sig12(1,3)
c
       sig11(2,3) = 0
       sig11(3,2) = 0
       sig11(1,1) =  ((beta1 * sig) ** 2) + (sige ** 2)
       sig11(2,2) = (sigu ** 2) 
       sig11(3,3) = (sigu ** 2) 
       sig11(1,2) = beta1 * (sig ** 2) 
       sig11(1,3) = beta1 * (sig ** 2) 
       sig11(2,1) = sig11(1,2)
       sig11(3,1) = sig11(1,3)
       do 200 j=1,2
       do 150 k=1,2
       sig11(j+1,k+1) = sig11(j+1,k+1) +   (sig ** 2)
150    continue
200    continue
       return
       end
c/* ****************************************************************** */
c/*                                                                    */
c/*                                                                    */
c/**** This is the place where you generate the new X's             ****/
c/**** given the Roeder Wasserman indicators                        ****/
c      vvthis is a standard normal random matrix passed to the
c         program.
c
c      Note that this program assumes that numw = 2
c
c      This program can be sped up a lot if you use the full
c         power of the cycling phenomenon, as in Larry's
c         program. For example, if cycper = .10, then this
c         particular subroutine should take only 1/10th the
c         time than if cycper=1. This is a nontrivial difference,
c         and I need to program it in.  There is only one
c         hassle, and that is that there are actually
c         two cases depending on whether W is or is not
c         observed, and this has to be taken into account
c         in the program
c
c/*                                                                    */
c/* ****************************************************************** */
        subroutine findxxxx(yy,ww,sige,sigu,n,msamp,numw,beta0,beta1,
     6               pmix,zmumix,sigmix,iiest,vvthis,Kbigest,
     6               xxest,kklow,kkhigh,maxmix,xxest1)
       implicit real*8(a-h,o-z)
       real*8 yy(n,1),ww(msamp,numw),vvthis(n),xxest(n,1)
       real*8 pmix(maxmix,1),sigmix(maxmix,1),zmumix(maxmix,1)
c      local matrices specific to this program
       real*8 xxest1(n,1),c0(3,1)
       real*8 sig12(1,3),sig21(3,1),sig11(3,3),sig11inv(3,3)
       real*8 aaa(1,3),bbb(1,1)
       integer iiest(n,maxmix)
       do 20 i=kklow,kkhigh
 20       xxest(i,1) = 0.
       do 9900 mmm=1,Kbigest
c      Get the current group mean and standard deviation
       zmu        = zmumix(mmm,1)
       sig        = sigmix(mmm,1)
c      Get the s.d. of X given Y
       sigcond1   = dsqrt( ((sig * sige) ** 2)
     6               / ( ((beta1 * sig) ** 2) + (sige ** 2)) )
c      Get the s.d. of X given (Y,W_1,W_2) and other important
c          variables necessary in this calculation
       c0(1,1) = beta0 + (beta1 * zmu)
       c0(2,1) = zmu
       c0(3,1) = zmu
       call findsig8(sig12,sig21,sig11,beta1,sig,sige,sigu)
       call xinverse(sig11,sig11inv,3,bigdet)  
       call zmmultm(sig12,sig11inv,aaa,1,3,3) 
       call zmmultm(aaa,sig21,bbb,1,3,1)      
       sigcond2   = dsqrt( (sig ** 2)  - bbb(1,1) )
c      if kklow > msamp, then you are on a subset where
c            only Y is observed       
       if (kklow.gt.msamp) then
       do 100 i=kklow,kkhigh
 100      xxest1(i,1) = zmu + 
     6         (beta1 * sig * sig * 
     6           (yy(i,1) - beta0 - (beta1 * zmu))
     6               / ( ((beta1 * sig) ** 2) + (sige ** 2)))
     6                             + (sigcond1 * vvthis(i))
       endif
c      if kkhigh <= msamp, then you are on a subset
c           where W is always observed
       if (kkhigh.le.msamp) then
       do 300 i=kklow,kkhigh
 300   xxest1(i,1)  = zmu 
     6                     + (aaa(1,1) * (yy(i,1)-c0(1,1)))
     6                     + (aaa(1,2) * (ww(i,1)-c0(2,1)))
     6                     + (aaa(1,3) * (ww(i,2)-c0(3,1)))
     6                     + (sigcond2 * vvthis(i))
       endif
c      if kklow <= msamp and kkhigh > msamp, then you have
c          to generate both sets, because your subset
c          contains both kinds of data
       if ((kklow.le.msamp).and.(kkhigh.gt.msamp)) then
       do 325 i=(msamp+1),kkhigh
          xxest1(i,1) = zmu + 
     6         (beta1 * sig * sig * 
     6           (yy(i,1) - beta0 - (beta1 * zmu))
     6               / ( ((beta1 * sig) ** 2) + (sige ** 2)))
     6                             + (sigcond1 * vvthis(i))
 325      continue
       do 350 i=kklow,msamp
 350      xxest1(i,1)  = zmu 
     6                     + (aaa(1,1) * (yy(i,1)-c0(1,1)))
     6                     + (aaa(1,2) * (ww(i,1)-c0(2,1)))
     6                     + (aaa(1,3) * (ww(i,2)-c0(3,1)))
     6                     + (sigcond2 * vvthis(i))
       endif
c      Now generate the X's
       do 400 i=kklow,kkhigh
 400      xxest(i,1) = xxest(i,1) + (xxest1(i,1) * iiest(i,mmm))
9900   continue
       return
       end