Nothing
recursions <- function(S,R,yv,Psi,piv,Pi,k,lth,matr,Bm,Cm,bv,mod){
# forward recursion
sS = dim(S)
ns = sS[1]
TT = sS[2]
if(length(sS)==2){
r = 1
S = array(S,c(ns,TT,r))
}else r = sS[3]
b = max(bv)
miss = !is.null(R)
M = array(1,c(k,TT,ns))
if(miss) for(i in 1:ns) for(t in 1:TT) for(j in 1:r) M[,t,i] = M[,t,i]*(Psi[S[i,t,j]+1,,j]*R[i,t,j]+(1-R[i,t,j]))
else for(i in 1:ns) for(t in 1:TT) for(j in 1:r) M[,t,i] = M[,t,i]*Psi[S[i,t,j]+1,,j]
Q = array(0,c(k,TT,ns))
pv = rep(0,ns)
for(i in 1:ns){
q = M[,1,i]*piv
Q[,1,i] = q
for(t in 2:TT){
q = M[,t,i]*(t(Pi[,,t])%*%q)
Q[,t,i] = q
}
pv[i] = sum(q)
}
#log-likelihood
lk = sum(yv*log(pv))
# backward recurion
Qb = array(0,c(k,TT,ns))
for(i in 1:ns){
qb = rep(1,k)
Qb[,TT,i] = qb
for(t in seq(TT-1,1,-1)){
qb = (Pi[,,t+1])%*%(M[,t+1,i]*qb)
Qb[,t,i] = qb
}
}
# posterior probabilities
F1 = Q*Qb
for(i in 1:ns) F1[,,i] = F1[,,i]/pv[i]
F2 = array(0,c(k,k,TT,ns))
QbM = Qb*M
for(i in 1:ns) for(t in 2:TT) F2[,,t,i] = Pi[,,t]*(Q[,t-1,i]%o%QbM[,t,i])/pv[i]
#derivative of the forward recursion
ind = 0
#derivative of Psi
Psid = array(NA,c(b+1,k,r,lth))
for(j in 1:r) Psid[1:(bv[j]+1),,j,] = 0
for(u in 1:k) for(j in 1:r){
indj = 1:(bv[j]+1)
Om = diag(Psi[indj,u,j])-Psi[indj,u,j]%o%Psi[indj,u,j]
D = Om%*%matr[[j]]$Am
for(h in 1:bv[j]){
ind = ind+1
Psid[indj,u,j,ind] = D[,h]
}
}
#derivative of piv
pivd = matrix(0,k,lth)
Om = diag(piv)-piv%o%piv
D = Om%*%Bm
for(j in 1:(k-1)){
ind = ind+1
pivd[,ind] = D[,j]
}
#derivative of Pi
if(mod==0){
Pid = array(0,c(k,k,TT,lth))
for(t in 2:TT) for(u in 1:k){
Om = diag(Pi[u,,t])-Pi[u,,t]%o%Pi[u,,t]
D = Om%*%Cm[,,u]
for(j in 1:(k-1)){
ind = ind+1
Pid[u,,t,ind] = D[,j]
}
}
}
if(mod==1){
Pid = array(0,c(k,k,lth))
for(u in 1:k){
Om = diag(Pi[u,,2])-Pi[u,,2]%o%Pi[u,,2]
D = Om%*%Cm[,,u]
for(j in 1:(k-1)){
ind = ind+1
Pid[u,,ind] = D[,j]
}
}
}
if(mod>1){
Pid = array(0,c(k,k,2,lth))
for(u in 1:k){
Om = diag(Pi[u,,2])-Pi[u,,2]%o%Pi[u,,2]
D = Om%*%Cm[,,u]
for(j in 1:(k-1)){
ind = ind+1
Pid[u,,1,ind] = D[,j]
}
}
for(u in 1:k){
Om = diag(Pi[u,,mod+1])-Pi[u,,mod+1]%o%Pi[u,,mod+1]
D = Om%*%Cm[,,u]
for(j in 1:(k-1)){
ind = ind+1
Pid[u,,2,ind] = D[,j]
}
}
}
#derivative of the forward recursion
Md = array(0,c(k,TT,ns,lth))
for(i in 1:ns) for(t in 1:TT) for(j in 1:r){
if(miss) Md[,t,i,] = Md[,t,i,]+(Psid[S[i,t,j]+1,,j,]*R[i,t,j])*M[,t,i]/(Psi[S[i,t,j]+1,,j]*R[i,t,j]+(1-R[i,t,j]))
else Md[,t,i,] = Md[,t,i,]+Psid[S[i,t,j]+1,,j,]*M[,t,i]/Psi[S[i,t,j]+1,,j]
}
Qd = array(0,c(k,TT,ns,lth))
pvd = matrix(0,ns,lth)
for(i in 1:ns){
Qd[,1,i,] = Md[,1,i,]*piv+M[,1,i]*pivd
for(t in 2:TT){
q = Q[,t-1,i]
if(mod==0) Qd[,t,i,] = Md[,t,i,]*as.vector(t(Pi[,,t])%*%q)+M[,t,i]*(t(Pi[,,t])%*%Qd[,t-1,i,])
for(j in 1:lth){
qd = Qd[,t-1,i,j]
if(mod==0) Qd[,t,i,j] = Qd[,t,i,j]+M[,t,i]*(t(Pid[,,t,j])%*%q)
if(mod==1) Qd[,t,i,j] = Md[,t,i,j]*(t(Pi[,,t])%*%q)+M[,t,i]*(t(Pid[,,j])%*%q)+M[,t,i]*(t(Pi[,,t])%*%qd)
if(mod>1){
if(t<=mod) Qd[,t,i,j] = Md[,t,i,j]*(t(Pi[,,t])%*%q)+M[,t,i]*(t(Pid[,,1,j])%*%q)+M[,t,i]*(t(Pi[,,t])%*%qd)
if(t>mod) Qd[,t,i,j] = Md[,t,i,j]*(t(Pi[,,t])%*%q)+M[,t,i]*(t(Pid[,,2,j])%*%q)+M[,t,i]*(t(Pi[,,t])%*%qd)
}
}
}
for(j in 1:lth) pvd[i,j] = sum(Qd[,TT,i,j])
}
#score
sc = (yv/pv)%*%pvd
# derivative of backward recurion
Qbd = array(0,c(k,TT,ns,lth))
for(j in 1:lth){
for(i in 1:ns){
qbd = rep(0,k)
Qbd[,TT,i,j] = qbd
for(t in seq(TT-1,1,-1)){
qb = Qb[,t+1,i]
if(mod==0) qbd = (Pid[,,t+1,j])%*%(M[,t+1,i]*qb)+(Pi[,,t+1])%*%(Md[,t+1,i,j]*qb)+(Pi[,,t+1])%*%(M[,t+1,i]*qbd)
if(mod==1) qbd = (Pid[,,j])%*%(M[,t+1,i]*qb)+(Pi[,,t+1])%*%(Md[,t+1,i,j]*qb)+(Pi[,,t+1])%*%(M[,t+1,i]*qbd)
if(mod>1){
if(t>=mod) qbd = (Pid[,,2,j])%*%(M[,t+1,i]*qb)+(Pi[,,t+1])%*%(Md[,t+1,i,j]*qb)+(Pi[,,t+1])%*%(M[,t+1,i]*qbd)
if(t<mod) qbd = (Pid[,,1,j])%*%(M[,t+1,i]*qb)+(Pi[,,t+1])%*%(Md[,t+1,i,j]*qb)+(Pi[,,t+1])%*%(M[,t+1,i]*qbd)
}
Qbd[,t,i,j] = qbd
}
}
}
# posterior probabilities
F1d = array(0,c(k,TT,ns,lth))
F1t = Q*Qb
for(j in 1:lth){
F1dj = Qd[,,,j]*Qb+Q*Qbd[,,,j]
for(i in 1:ns) F1d[,,i,j] = F1dj[,,i]/pv[i]-F1t[,,i]*pvd[i,j]/pv[i]^2
}
F2d = array(0,c(k,k,TT,ns,lth))
for(j in 1:lth){
QbMdj = Qbd[,,,j]*M+Qb*Md[,,,j]
for(i in 1:ns) for(t in 2:TT){
if(mod==0) F2d[,,t,i,j] = Pid[,,t,j]*(Q[,t-1,i]%o%QbM[,t,i])/pv[i]+Pi[,,t]*(Qd[,t-1,i,j]%o%QbM[,t,i])/pv[i]+
Pi[,,t]*(Q[,t-1,i]%o%QbMdj[,t,i])/pv[i]-Pi[,,t]*(Q[,t-1,i]%o%QbM[,t,i])*pvd[i,j]/pv[i]^2
if(mod==1) F2d[,,t,i,j] = Pid[,,j]*(Q[,t-1,i]%o%QbM[,t,i])/pv[i]+Pi[,,t]*(Qd[,t-1,i,j]%o%QbM[,t,i])/pv[i]+
Pi[,,t]*(Q[,t-1,i]%o%QbMdj[,t,i])/pv[i]-Pi[,,t]*(Q[,t-1,i]%o%QbM[,t,i])*pvd[i,j]/pv[i]^2
if(mod>1){
if(t <=mod) F2d[,,t,i,j] = Pid[,,1,j]*(Q[,t-1,i]%o%QbM[,t,i])/pv[i]+Pi[,,t]*(Qd[,t-1,i,j]%o%QbM[,t,i])/pv[i]+
Pi[,,t]*(Q[,t-1,i]%o%QbMdj[,t,i])/pv[i]-Pi[,,t]*(Q[,t-1,i]%o%QbM[,t,i])*pvd[i,j]/pv[i]^2
if(t>mod) F2d[,,t,i,j] = Pid[,,2,j]*(Q[,t-1,i]%o%QbM[,t,i])/pv[i]+Pi[,,t]*(Qd[,t-1,i,j]%o%QbM[,t,i])/pv[i]+
Pi[,,t]*(Q[,t-1,i]%o%QbMdj[,t,i])/pv[i]-Pi[,,t]*(Q[,t-1,i]%o%QbM[,t,i])*pvd[i,j]/pv[i]^2
}
}
}
out = list(lk=lk,sc=sc,F1=F1,F2=F2,F1d=F1d,F2d=F2d,M=M,Q=Q,pv=pv)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.