1! Subroutine for computation of the approximating Gaussian model for non-Gaussian models 2 3subroutine approx(yt, ymiss, timevar, zt, tt, rtv, ht, qt, a1, p1,p1inf, p, n, m, r,& 4theta, u, ytilde, dist, maxiter, tol, rankp, convtol, diff, lik, info, expected, htol) 5 6 implicit none 7 8 integer, intent(in) :: p,m, r, n,rankp, expected 9 integer, intent(in), dimension(n,p) :: ymiss 10 integer, intent(in), dimension(5) :: timevar 11 integer, intent(in), dimension(p) :: dist 12 integer, intent(inout) :: maxiter,info 13 integer :: i, k,tvrqr,kk,jt,dt 14 double precision, intent(in) :: tol,convtol 15 double precision, intent(inout) :: htol 16 double precision, intent(in), dimension(n,p) :: u 17 double precision, intent(in), dimension(n,p) :: yt 18 double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt 19 double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt 20 double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv 21 double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt 22 double precision, intent(in), dimension(m) :: a1 23 double precision, intent(in), dimension(m,m) :: p1,p1inf 24 double precision, intent(inout), dimension(n,p) :: theta 25 double precision, intent(inout), dimension(n,p) :: ytilde 26 double precision, intent(inout), dimension(p,p,n) :: ht 27 double precision, intent(inout) :: diff 28 double precision, dimension(m,r) :: mr 29 double precision, dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr 30 double precision, intent(inout) :: lik 31 double precision, external :: ddot 32 integer, external :: finitex 33 double precision dev, devold 34 double precision, dimension(n,p) :: thetanew, thetaold 35 double precision, dimension(p,n) :: ft,finf 36 double precision, dimension(m,p,n) :: kt,kinf 37 38 external dgemm, pytheta, pthetafirst, approxloop, pthetarest 39 40 !compute rqr 41 tvrqr = max(timevar(4),timevar(5)) 42 do i=1, (n-1)*tvrqr+1 43 call dgemm('n','n',m,r,r,1.0d0,rtv(:,:,(i-1)*timevar(4)+1),m,qt(:,:,(i-1)*timevar(5)+1),r,0.0d0,mr,m) 44 call dgemm('n','t',m,m,r,1.0d0,mr,m,rtv(:,:,(i-1)*timevar(4)+1),m,0.0d0,rqr(:,:,i),m) 45 end do 46 47 ! compute logp(theta) for the first time, no need to compute kt/kinf and ft/finf in successive calls 48 ! in case of totally diffuse initialization term p(theta)=0 for all theta 49 if(rankp .NE. m) then 50 call pthetafirst(theta, timevar, zt, tt, rqr, a1, p1, p1inf, p, m, n, devold, tol,rankp,kt,kinf,ft,finf,dt,jt) 51 end if 52 thetaold = theta 53 devold = -huge(devold) 54 k=0 55 do while(k .LT. maxiter) 56 57 k=k+1 58 ! compute new guess thetanew 59 call approxloop(yt, ymiss, timevar, zt, tt, rtv, ht, qt, rqr, tvrqr, a1, p1,p1inf, p,n,m,r, & 60 theta, thetanew, u, ytilde, dist,tol,rankp,lik, expected) 61 ! and log(p(theta|y)) 62 call pytheta(thetanew, dist, u, yt, ymiss, dev, p, n) 63 if(rankp .NE. m) then 64 call pthetarest(thetanew, timevar, zt, tt, a1, p, m, n, dev, kt,kinf,ft,finf,dt,jt) 65 end if 66 !non-finite value in linear predictor or muhat 67 if(finitex(sum(thetanew, MASK = ymiss.EQ.0)).EQ.0 .OR. & 68 finitex(maxval(exp(thetanew), MASK = ymiss.EQ.0)).EQ.0 ) then 69 if(k .GT. 1) then 70 kk = 0 71 do while(finitex(sum(thetanew, MASK = ymiss.EQ.0)).EQ.0 .OR. & 72 finitex(maxval(exp(thetanew), MASK = ymiss.EQ.0)).EQ.0) 73 kk = kk + 1 74 if(kk .GT. maxiter) then 75 info = 1 76 return 77 end if 78 !backtrack 79 theta = 0.5d0*(thetaold+theta) 80 call approxloop(yt, ymiss, timevar, zt, tt, rtv, ht, qt, rqr, tvrqr, a1, p1,p1inf, p,n,m,r, & 81 theta, thetanew, u, ytilde, dist,tol,rankp,lik, expected) 82 83 call pytheta(thetanew, dist, u, yt, ymiss, dev, p, n) 84 if(rankp .NE. m) then 85 call pthetarest(thetanew, timevar, zt, tt, a1, p, m, n, dev, kt,kinf,ft,finf,dt,jt) 86 end if 87 end do 88 else !cannot correct step size as we have just began 89 info = 1 90 return 91 end if 92 end if 93 94 if(finitex(dev).EQ.0) then !non-finite value of objective function 95 if(k .GT. 1) then 96 kk = 0 97 do while(finitex(dev).EQ.0) 98 kk = kk + 1 99 if(kk .GT. maxiter) then !did not find valid likelihood 100 info = 2 101 return 102 end if 103 104 theta = 0.5d0*(thetaold+theta) 105 call approxloop(yt, ymiss, timevar, zt, tt, rtv, ht, qt, rqr, tvrqr, a1, p1,p1inf, p,n,m,r, & 106 theta, thetanew, u, ytilde, dist,tol,rankp,lik, expected) 107 108 call pytheta(thetanew, dist, u, yt, ymiss, dev, p, n) 109 if(rankp .NE. m) then 110 call pthetarest(thetanew, timevar, zt, tt, a1, p, m, n, dev, kt,kinf,ft,finf,dt,jt) 111 end if 112 113 end do 114 else !cannot correct step size as we have just began 115 info = 2 116 return 117 end if 118 end if 119 120 121 ! decreasing deviance 122 if((dev - devold)/(0.1d0 + abs(dev)) .LT. -convtol .AND. k .GT. 1) then 123 kk = 0 124 do while((dev - devold)/(0.1d0 + abs(dev)) .LT. convtol .AND. kk .LT. maxiter) 125 kk = kk + 1 126 ! previous theta produced too 'big' thetanew 127 ! new guess by halving the last try 128 theta = 0.5d0*(thetaold+theta) 129 call approxloop(yt, ymiss, timevar, zt, tt, rtv, ht, qt, rqr, tvrqr, a1, p1,p1inf, p,n,m,r, & 130 theta, thetanew, u, ytilde, dist,tol,rankp,lik, expected) 131 132 call pytheta(thetanew, dist, u, yt, ymiss, dev, p, n) 133 if(rankp .NE. m) then 134 call pthetarest(thetanew, timevar, zt, tt, a1, p, m, n, dev, kt,kinf,ft,finf,dt,jt) 135 end if 136 137 end do 138 end if 139 140 diff = abs(dev - devold)/(0.1d0 + abs(dev)) 141 if(diff .LT. convtol) then !convergence 142 theta=thetanew 143 info=0 144 exit 145 else 146 thetaold=theta 147 theta=thetanew 148 devold = dev 149 end if 150 end do 151 if(maxiter.EQ.k) then 152 info=3 153 end if 154 maxiter=k 155 156 ! check if the the approximation lead to potentially zero signal-to-noise ratio 157 if(maxval(ht) > htol) then 158 info = -5 159 htol = maxval(ht) 160 return 161 end if 162 163end subroutine approx 164