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