subroutine mode2DfB0u(Astart,Y,V,Vo,Vs,W,Wt,dsi,dsr,J,r,tol,maxit,
C   outputs:
     $		newA,lf,lA,iter, d1, 
C   workspace:
     $		deltinv,B,Dgam,WDY,gam,mu,H,us,Dus,
     $   	ttol,ldDg,temprr,tempr,tempJr,tempJ,rt)
C
C This subroutine finds the posterior mode of B0=Vo/(Vo+A) for
C UNIVARIATE (p=1) normal outcomes, assuming a prior distribution
C on B* proportional to B*^{ds/2}. Calls subroutine EMfB0stepu.f
C to update H and A at each step.
C
C
	integer J,r,i,k,maxit,iter, miniter, dsi
	double precision Astart, Y(J), V(J), Vo,Vs, W(J,r),Wt(r,J),dsr, 
     $ 		B(J),WDY(r,1), tol, newA, oldA, lf0,lf, deltinv(J), 
     $  	Dgam(r,r),gam(r,1),mu(J),lA,H,ttol,ldDg,us(J),Dus(J),
     $		temprr(r,r),tempr(r,1),tempJr(J,r),tempJ(J),d, d1,
     $		oldB,newB,rt(3)
C
C
C initialize
	miniter=3
	ttol=1.0
	iter=0
	lf0 = -10000000.0
	lf  = 0.0
	d1 = 0.0
	do while(ttol.gt.tol.and.iter.lt.maxit
     $		.or. iter.lt.miniter)
	do 20 i=1,r
	 WDY(i,1)=0.0
	 do 10 k=1,r
	  Dgam(i,k)=0.0
 10	 continue
 20  	continue
C
C update lf, A
C
	call EMfB0stepu(Astart,Y,V,Vo,Vs,W,Wt,dsi,dsr,J,r,newA,lf,lA,
     $		deltinv,B,Dgam,WDY,gam,mu,H,ldDg,us,Dus,
     $		temprr,tempr,tempJr,tempJ,rt)
C
	ttol=lf-lf0
	newB=Vo/(Vo+newA)
	oldB=Vo/(Vo+Astart)
	d1=(lf-lf0)/(newB-oldB)
	lf0=lf
	Astart=newA
	iter=iter+1
	end do
	return
	end


	subroutine EMfB0stepu(Astart,Y,V,Vo,Vs,W,Wt,dsi,dsr,J,r,newA,lf,
     $		lA, deltinv,B,Dgam,WDY,gam,mu,H,ldDg,us,Dus,
     $		temprr,tempr,tempJr,tempJ,rt)
C
C performs E and M steps to update A (newA) and computes log
C posterior density for B0 value implied by Astart, assuming
C p(B*) \propto B*^{ds/2}
C
	integer J,r,i,k,dsi
	double precision Astart, Y(J), V(J), Vo,Vs,W(J,r),Wt(r,J), 
     $		newA, lf, B(J), deltinv(J), us(J),Dus(J), dsr,
     $		Dgam(r,r), WDY(r,1), gam(r,1), mu(J), lA, H,ldDg,
     $		temprr(r,r),tempr(r,1),tempJr(J,r),tempJ(J),x,
     $		c,d,e,f,rt(3)
	do 200 i=1,J
	 deltinv(i)=1/(V(i)+Astart)
	 B(i)=V(i)*deltinv(i)
	 tempJ(i)=Y(i)*deltinv(i)
C tempJ = Deltinv%*%Y
	 do 100 k=1,r
  	  tempJr(i,k)=deltinv(i)*W(i,k)
C tempJr = Deltinv%*%W
 100     continue
 200    continue
	call mmult(Wt,tempJ,r,J,1,WDY)
	call mmult(Wt,tempJr,r,J,r,Dgam)
	call inverse(Dgam,r,tempr,ldDg)
C ldDg = -log(|Dgam|)  (log det of Dgam inverse)
	call mmult(Dgam,WDY,r,r,1,gam)
	call mmult(W,gam,J,r,1,mu)
 	lA=-ldDg/2
	do 300 i=1,J
 	 x = Y(i) - mu(i)
	 lA = lA + (dlog(deltinv(i)) - deltinv(i)*(x**2))/2
 300	continue
	lf = lA - dsr*dlog(Vs+Astart)/2 + 2*dlog((Vo+Astart)/(Vs+Astart))
C
C  E-step
C
	call mmult(W,Dgam,J,r,r,tempJr)
	do 500 i=1,J
	 x = 0
	 us(i) = (1-B(i))*(Y(i)-mu(i))
  	 do 400 k=1,r
	  x = x + tempJr(i,k)*Wt(k,i)
C x = Wi'%*%Dgam%*%Wi
 400	 continue
	 Dus(i) = (1-B(i))*V(i) + x*(1-B(i))**2
 500   	continue
	H=0.0
	do 600 i=1,J
	 H = H + us(i)**2 + Dus(i)
 600	continue
C H = E(sum(uj^2)|Y,Astart)
C
C  M-step
C
	if(Vo.ne.Vs) then
	 c = -(J+dsi)
	 d = H-J*(Vs+Vo)+4*(Vs-Vo)-dsr*Vo
	 e = H*(Vs+Vo)-J*Vs*Vo
	 f = H*Vs*Vo
	 call cubicrt(c,d,e,f,rt,k)
	 if(k.gt.1) then
	  x = 0.0
	  c = Astart*1000000
   	  do 700 i=1,k
	   if(rt(i).gt.0.0) then
	    d=rt(i)-Astart
	    if(d.lt.0.0) then
	     d=-d
	    end if              
C end if(d.lt.0.0)
 	    if(d.lt.c) then
	     newA = rt(i)
	     c = d
	    end if              
C end if(d.lt.c)
	   end if               
C end if(rt(i).gt.0.0)
 700	  continue              
C end do 700 i=1,k
	  else                  
C if not (k.gt.1)
	  newA=rt(1)
	end if			
C end if(k.gt.1)
	else			
C if not (Vo.ne.Vs) (i.e., Vo=Vs)
	 newA = (H-J*Vo + sqrt((H-J*Vo)**2 + 4*(J+dsr)*H*Vo))/(2*(J+dsr))
	end if
C end if(Vo.ne.Vs)
	return
	end

	subroutine postu(Adraws,N,dsi,dsr,f0,adj,Y,V,Vo,Vs,W,Wt,J,r,
     $		quant,zstar,
C outputs
     $		thetap,thetup,thetlp,Vp,gammap,gamlp,gamup,Dgp,Ap,Aqnt,f1,rat,
     $		thetam,Vm,gammam,Dga,maxrat,sumrat,
C workspace
     $		B,deltinv,WDY,gam,Dgam,mu,tempr,temprr,tempJr,tempJ)

	integer N, J, r, i, k, l, draw, dsi
	double precision Adraws(N), ds, f0(N), adj, Y(J), V(J),  Vo, 
     $		Vs, W(J,r), Wt(r,J), quant(7), thetap(J),thetup(J),
     $		thetlp(J), Vp(J), zstar, dsr,
     $   	gammap(r), gamlp(r),gamup(r),Dgp(r,r), Ap, Aqnt(7), 
     $		f1(N), rat(N), 
     $		thetam(J,N), Vm(J,N), gammam(r,N), Dga(r,r,N), 
     $		maxrat, sumrat, B(J), deltinv(J), WDY(r,1), 
     $		gam(r), Dgam(r,r), mu(J), tempr(r), temprr(r,r), 
     $		tempJr(J,r), tempJ(J), A, lf, lA, ldDg, x, e, q

C
C This subroutine processes UNIVARIATE B0 draws from an approximating 
C density f0, calculating the true posterior density f(B0|Y) for a
C prior distribution p(B*) = B*^(d*/2); B* = V*/(V*+A), the 
C normalized importance weights f(B0|Y)/f0(B0), the conditional 
C means and variances of the thetas and of gamma, given each draw
C B0, and the importance weighted unconditional posterior estimates
C 
C initialize:
	e=2.718282
	do 100 i=1,J
	 thetap(i)=0.0
	 thetlp(i)=0.0
	 thetup(i)=0.0
	 Vp(i)=0.0
 100	continue
	do 300 i=1,r
	 gammap(i)=0.0
	 gamlp(i)=0.0
	 gamup(i)=0.0
	 do 200 k=1,r
	  Dgp(i,k)=0.0
 200	 continue
 300 	continue
	Ap=0.0
	maxrat=1.0
	sumrat=0.0
C
	do 1300 draw = 1,N
	 A = Adraws(draw)
	 do 500 i=1,r
	  gam(i)=0.0
	  WDY(i,1)=0.0
	  do 400 k = 1,r
	   Dgam(i,k) = 0.0
 400 	  continue
 500	 continue
	 do 700 i=1,J
	  deltinv(i)=1/(V(i)+A)
	  B(i)=V(i)*deltinv(i)
	  tempJ(i)=Y(i)*deltinv(i)
C tempJ = Deltinv%*%Y
	  do 600 k=1,r
  	   tempJr(i,k)=deltinv(i)*W(i,k)
C tempJr = Deltinv%*%W
 600      continue
 700     continue
	 call mmult(Wt,tempJ,r,J,1,WDY)
	 call mmult(Wt,tempJr,r,J,r,Dgam)
	 call inverse(Dgam,r,tempr,ldDg)
C Dgam = (W'%*%Deltinv%*%W)^{-1}
C ldDg = -log(|Dgam|)  (log det of Dgam inverse)
	 call mmult(Dgam,WDY,r,r,1,gam)
	 call mmult(W,gam,J,r,1,mu)
 	 lA=-ldDg/2
	 do 800 i=1,J
 	  x = Y(i) - mu(i)
	  lA = lA + (dlog(deltinv(i)) - deltinv(i)*(x**2))/2
 800	 continue
	 lf = lA - dsr*dlog(Vs+A)/2+2*dlog((Vo+A)/(Vs+A))+adj
	 f1(draw) = exp(lf)
	 rat(draw) = f1(draw)/f0(draw)
	 Ap = Ap + rat(draw)*A
	 sumrat=sumrat+rat(draw)
	 if(rat(draw).gt.maxrat) then
	  maxrat = rat(draw)
	 end if
C load gammam and Dga
	 do 1000 i=1,r
	  gammam(i,draw) = gam(i)
	  gammap(i) = gammap(i) + rat(draw)*gam(i)
	  do 900 k=1,r
	   Dga(i,k,draw)=Dgam(i,k)
	   Dgp(i,k) = Dgp(i,k) +rat(draw)*Dgam(i,k)
 900	  continue
 1000	 continue
 	 do 1010 i=1,r
	  gamlp(i)=gamlp(i)+rat(draw)*(gammam(i,draw)-zstar*sqrt(Dga(i,i,draw)))
 	  gamup(i)=gamup(i)+rat(draw)*(gammam(i,draw)+zstar*sqrt(Dga(i,i,draw)))
 1010	 continue
C load thetap and Vp
	 call mmult(W,Dgam,J,r,r,tempJr)
	 do 1200 i=1,J
	  thetam(i,draw) = B(i)*mu(i)+(1-B(i))*Y(i)
	  thetap(i) = thetap(i) + rat(draw)*thetam(i,draw)
	  x=0.0
	  do 1100 k=1,r
	   x = x + tempJr(i,k)*Wt(k,i)
 1100	  continue
C x = Wi'%*%Dgam%*%Wi
	  Vm(i,draw)=(1-B(i))*V(i)+x*B(i)**2
	  Vp(i) = Vp(i) + rat(draw)*Vm(i,draw)
	  thetlp(i)=thetlp(i)+rat(draw)*(thetam(i,draw)-zstar*sqrt(Vm(i,draw)))
	  thetup(i)=thetup(i)+rat(draw)*(thetam(i,draw)+zstar*sqrt(Vm(i,draw)))
 1200	 continue
 1300	continue
C normalize:
	Ap = Ap/sumrat
	do 1400 i=1,J
	 thetap(i)=thetap(i)/sumrat
	 thetlp(i)=thetlp(i)/sumrat
	 thetup(i)=thetup(i)/sumrat
	 Vp(i)=Vp(i)/sumrat
 1400	continue
	do 1500 i=1,r
	 gammap(i) = gammap(i)/sumrat
	 gamlp(i)=gamlp(i)/sumrat
	 gamup(i)=gamup(i)/sumrat
	 do 1450 k=1,r
	  Dgp(i,k)=Dgp(i,k)/sumrat
 1450    continue
 1500	continue
C thetap, gammap and Ap are done; Vp and Dgp contain means of
C conditional variances; must add variances of conditonal means.
C Also set Aqnt.
	x = 0.0
	l = 1
	do 2000 draw=1,N
	 q=quant(l)
	 x = x + rat(draw)/sumrat
	 if(x.gt.q.and.l.lt.8) then
	  Aqnt(l)=Adraws(draw)
	  l=l+1
	 end if
	 do 1600 i=1,J
	  Vp(i)=Vp(i)+(rat(draw)*(thetam(i,draw)-thetap(i))**2)/sumrat
 1600	 continue
	 do 1700 i=1,r
	  tempr(i)=gammam(i,draw)-gammap(i)
 1700 	 continue
	 call mmult(tempr,tempr,r,1,r,temprr)
	 do 1900 i=1,r
	  do 1800 k=1,r
	   Dgp(i,k)=Dgp(i,k)+(rat(draw)*temprr(i,k))/sumrat
 1800	  continue
 1900	 continue
 2000	continue
	return
	end
	 

      subroutine inverse(X,p,u,dlogdet)
C X -- a p x p matrix, (returns X^{-1})
C p -- dimension
C u -- work space of length p
C ifault = an integer number
C dlogdet - returns the logrithm of the determinant of X
C
      integer p,k,i,j,ifault
      double precision X(p,p),u(p),flag, dlogdet
      flag=1.0
      dlogdet = 0.0
      do k=1,p
         dlogdet = dlogdet + dlog(X(k,k))
         call swp(X,p,k,flag,u,ifault)
      end do
      do i=1,p
      do j=1,i
         X(i,j)=-X(i,j)
         X(j,i)=X(i,j)
      end do
      end do
      return
      end

      subroutine swp(X,p,k,flag,u,ifault)
c
C This subroutine performs Gauss sweep on the lower triangular
C      part of a symmetric p x p matrix X.
C flag= 1: swp
C flag=-1: reverse swp
c u: work space
c
      integer p,k,i,j,ifault
      double precision X(p,p),u(p),flag,c
      ifault=0
      if(X(k,k).lt.1.0E-10.and.flag.eq.1.0) then
         ifault=1
      end if
      c=1./X(k,k)
      do 1 j=1,k
      u(j)=X(k,j)
    1 X(k,j)=0.0
      do 2 i=k,p
      u(i)=X(i,k)
    2 X(i,k)=0.0
      u(k)=-1.0*flag
      do 3 i=1,p
      do 3 j=1,i
    3 X(i,j)=X(i,j)-c*u(i)*u(j)
      do 4 i=1,p
      do 4 j=i+1,p
    4 X(i,j)=X(j,i)
      return
      end


        subroutine mmult(A,B,p,q,r,C)
C
C this subroutine calculates the matrix product od A and
C B and puts the result in C. p,q,r are dimensions
       integer p,q,r,i,j,k
       double precision t, A(p,q), B(q,r), C(p,r)
        do 100 i=1,p
         do 200 j=1,r
          t=0.0
          do 300 k=1,q
           t=t+A(i,k)*B(k,j)
  300     continue
          C(i,j)=t
  200    continue
  100   continue
        return
        end


	subroutine cubicrt(e,f,g,h,rt,nrt)
C
C returns the (real) solutions to the equation
C ay^3+by^2+cy+d = 0; (ey^3+fy^2+gy+h = 0?)
C nrt gives number of roots
C rt(1:nrt) gives the roots
C
	double precision a, b, c, d,e,f,g,h, p, q, r,rt(3),dd,thet,pi, rt3
C
	integer nrt
	pi = 3.141593
	rt3 = sqrt(3.0)
	rt(1)=0.0
	rt(2)=0.0
	rt(3)=0.0
	p=f/e
	q=g/e
	r=h/e
C y^3+py^2+qy+r = 0
	a=(3*q-p**2)/3 
 	b=(2*(p**3)-9*p*q+27*r)/27
C x^3+ax+b = 0,  x-p/3 = y
	dd=(b**2)/4+(a**3)/27
	if(dd.eq.0.0) then
	 call cbrt(b/2,e)
	 rt(1)= -2*e-p/3
	 rt(2)= e-p/3
	 nrt=2
	end if
	if(dd.gt.0.0) then
   	 e=-b/2+sqrt(dd)
	 f=-b/2-sqrt(dd)
	 call cbrt(e,g)
	 call cbrt(f,h)
	 rt(1)=g+h-p/3
	 nrt=1
	end if
	if(dd.lt.0.0) then
	 c=-b/2
	 d=sqrt(-dd)
	 r=sqrt(c**2+d**2)
	 if(c.gt.0.0) then
	  thet=datan(d/c)
	 else
	  thet = datan(d/c)+pi
	 end if
	 call cbrt(r,h)
	 rt(1) = h*2*dcos(thet/3)-p/3
 	 rt(2) = h*(rt3*dsin(thet/3)-dcos(thet/3))-p/3
	 rt(3) = -h*(rt3*dsin(thet/3)+dcos(thet/3))-p/3
	 nrt=3
	end if
	return
	end

	
	subroutine cbrt(x,cbrtx)
	double precision x, cbrtx
	if(x.lt.0.0) then
	 cbrtx = -exp(dlog(-x)/3)
	else
	 cbrtx = exp(dlog(x)/3)
	end if
	return
	end