1! Importance sampling of non-gaussian model 2 3subroutine isample(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, dist, & 4p, n, m, r, theta, maxiter,rankp,convtol, nnd,nsim,epsplus,etaplus,& 5aplus1,c,tol,info,antithetics,w,sim,simwhat,simdim, expected, htol) 6 7 implicit none 8 9 integer, intent(in) :: p,m, r, n,nnd,antithetics,nsim,simwhat,simdim& 10 , rankp, expected 11 integer, intent(in), dimension(p) :: dist 12 integer, intent(in), dimension(n,p) :: ymiss 13 integer, intent(in), dimension(5) :: timevar 14 integer, intent(inout) :: maxiter,info 15 integer :: t, j,i,info2 16 double precision, intent(in) :: convtol,tol 17 double precision, intent(inout) :: htol 18 double precision, intent(in), dimension(n,p) :: u 19 double precision, intent(in), dimension(n,p) :: yt 20 double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt 21 double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt 22 double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv 23 double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt 24 double precision, intent(in), dimension(m) :: a1 25 double precision, intent(in), dimension(m,m) :: p1,p1inf 26 double precision, intent(in),dimension(nsim) :: c 27 double precision, intent(inout), dimension(p,n,nsim) :: epsplus 28 double precision, intent(inout), dimension(r,n,nsim) :: etaplus 29 double precision, intent(inout), dimension(m,nsim) :: aplus1 30 double precision, intent(inout), dimension(n,p) :: theta 31 double precision, dimension(p,p,n) :: ht 32 double precision, intent(inout), dimension(simdim,n,3 * nsim * antithetics + nsim) :: sim 33 double precision, dimension(p,(3 * nsim * antithetics + nsim)*(5-simwhat)) :: tsim 34 double precision, dimension(n,p) :: ytilde 35 double precision, dimension(n) :: tmp 36 double precision, dimension(3 * nsim * antithetics + nsim) :: w 37 double precision :: diff 38 double precision, external :: ddot 39 double precision :: lik 40 41 external approx, simgaussian 42 43 ht=0.0d0 44 45 ! approximate 46 call approx(yt, ymiss, timevar, zt, tt, rtv, ht, qt, a1, p1,p1inf, p,n,m,r,& 47 theta, u, ytilde, dist,maxiter,tol,rankp,convtol,diff,lik,info, expected, htol) 48 49 if(info .ne. 0 .and. info .ne. 3) then 50 return 51 end if 52 53 info2 = 0 54 ! simulate signals 55 call simgaussian(ymiss,timevar, ytilde, zt, ht, tt, rtv, qt, a1, p1, & 56 p1inf, nnd,nsim, epsplus, etaplus, aplus1, p, n, m, r, info2,rankp,& 57 tol,sim,c,simwhat,simdim,antithetics) 58 59 if(info2 /= 0) then 60 info = info2 61 return 62 end if 63 64 ! compute importance weights 65 66 w=1.0d0 67 68 if(simwhat.EQ.5) then 69 do j=1,p 70 select case(dist(j)) 71 case(2) !poisson 72 tmp = exp(theta(:,j)) 73 do t=1,n 74 if(ymiss(t,j) .EQ. 0) then 75 76 w = w*exp(yt(t,j)*(sim(j,t,:)-theta(t,j))-& 77 u(t,j)*(exp(sim(j,t,:))-tmp(t)))/& 78 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-sim(j,t,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 79 80 end if 81 end do 82 case(3) !binomial 83 tmp = log(1.0d0+exp(theta(:,j))) 84 do t=1,n 85 if(ymiss(t,j) .EQ. 0) then 86 87 w = w*exp( yt(t,j)*(sim(j,t,:)-theta(t,j))-& 88 u(t,j)*(log(1.0d0+exp(sim(j,t,:)))-tmp(t)))/& 89 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-sim(j,t,:))**2 -(ytilde(t,j)-theta(t,j))**2)) 90 91 end if 92 end do 93 case(4) ! gamma 94 tmp = exp(-theta(:,j)) 95 do t=1,n 96 if(ymiss(t,j) .EQ. 0) then 97 w = w*exp( u(t,j)*(yt(t,j)*(tmp(t)-exp(-sim(j,t,:)))+theta(t,j)-sim(j,t,:)))/& 98 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-sim(j,t,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 99 end if 100 end do 101 case(5) !negbin 102 tmp = exp(theta(:,j)) 103 do t=1,n 104 if(ymiss(t,j) .EQ. 0) then 105 w = w*exp(yt(t,j)*(sim(j,t,:)-theta(t,j)) +& 106 (yt(t,j)+u(t,j))*log((u(t,j)+tmp(t))/(u(t,j)+exp(sim(j,t,:)))))/& 107 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-sim(j,t,:))**2 -(ytilde(t,j)-theta(t,j))**2)) 108 109 end if 110 end do 111 end select 112 end do 113 114 else 115 do j=1,p 116 select case(dist(j)) 117 case(2) !poisson 118 tmp = exp(theta(:,j)) 119 do t=1,n 120 if(ymiss(t,j) .EQ. 0) then 121 do i=1,3 * nsim * antithetics + nsim 122 tsim(j,i) = ddot(m,zt(j,:,(t-1)*timevar(1)+1),1,sim(:,t,i),1) 123 end do 124 w = w*exp(yt(t,j)*(tsim(j,:)-theta(t,j))-& 125 u(t,j)*(exp(tsim(j,:))-tmp(t)))/& 126 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-tsim(j,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 127 128 end if 129 end do 130 case(3) !binomial 131 tmp = log(1.0d0+exp(theta(:,j))) 132 do t=1,n 133 if(ymiss(t,j) .EQ. 0) then 134 do i=1,3 * nsim * antithetics + nsim 135 tsim(j,i) = ddot(m,zt(j,:,(t-1)*timevar(1)+1),1,sim(:,t,i),1) 136 end do 137 w = w*exp( yt(t,j)*(tsim(j,:)-theta(t,j))-& 138 u(t,j)*(log(1.0d0+exp(tsim(j,:)))-tmp(t)))/& 139 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-tsim(j,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 140 141 end if 142 end do 143 case(4) ! gamma 144 tmp = exp(-theta(:,j)) 145 do t=1,n 146 if(ymiss(t,j) .EQ. 0) then 147 do i=1,3 * nsim * antithetics + nsim 148 tsim(j,i) = ddot(m,zt(j,:,(t-1)*timevar(1)+1),1,sim(:,t,i),1) 149 end do 150 w = w*exp( u(t,j)*(yt(t,j)*(tmp(t)-exp(-tsim(j,:)))+theta(t,j)-tsim(j,:)))/& 151 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-tsim(j,:))**2 -(ytilde(t,j)-theta(t,j))**2)) 152 end if 153 end do 154 case(5) !negbin 155 tmp = exp(theta(:,j)) 156 do t=1,n 157 if(ymiss(t,j) .EQ. 0) then 158 do i=1,3 * nsim * antithetics + nsim 159 tsim(j,i) = ddot(m,zt(j,:,(t-1)*timevar(1)+1),1,sim(:,t,i),1) 160 end do 161 w = w*exp(yt(t,j)*(tsim(j,:)-theta(t,j)) +& 162 (yt(t,j)+u(t,j))*log((u(t,j)+tmp(t))/(u(t,j)+exp(tsim(j,:)))))/& 163 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-tsim(j,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 164 end if 165 end do 166 end select 167 end do 168 169 170 end if 171 172 173end subroutine isample 174