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