Hello, I would like to use mapply to avoid using a loop but for some reason, I
can't seem to get it to work. I've included copies of my code below. The first
set of code uses a loop (and it works fine), and the second set of code
attempts to use mapply but I get a "subscript out of bounds" error. Any
guidance would be greatly appreciated. Xj, Yj, and Wj are also lists, and s2,
TAU, and GAMMA are scalars.
Thank You.
# THIS WORKS USING THE LOOP
for (j in 1:J) {
V.tilde.j <- solve((1/s2)*t(Xj[[j]])%*%Xj[[j]] + solve(TAU))
# Not singular case:
if(round(det(t(Xj[[j]])%*%Xj[[j]]),8)!=0) {
Beta.hat.j <- solve(t(Xj[[j]])%*%Xj[[j]])%*%t(Xj[[j]])%*%Yj[[j]]
V.j <- s2*solve(t(Xj[[j]])%*%Xj[[j]])
Lambda.j <- solve(solve(V.j) + solve(TAU))%*%solve(V.j)
Beta.tilde.j <- Lambda.j%*%Beta.hat.j + (diag(P) - Lambda.j)%*%Wj[[j]]%*%GAMMA
}
# Singular case
else {
Beta.tilde.j <- V.tilde.j%*%((1/s2)*t(Xj[[j]])%*%Yj[[j]] +
solve(TAU)%*%Wj[[j]]%*%GAMMA)
}
BETA.Js[[j]] <- t(rmnorm(1, mean=as.vector(Beta.tilde.j), V.tilde.j))
}
# THIS DOESN'T WORK USING MAPPLY
update.betas <- function(s2,Xj,Yj,TAU,Wj,GAMMA) {
V.tilde.j <- solve((1/s2)*t(Xj[[j]])%*%Xj[[j]] + solve(TAU))
# Not singular case:
if(round(det(t(Xj[[j]])%*%Xj[[j]]),8)!=0) {
Beta.hat.j <- solve(t(Xj[[j]])%*%Xj[[j]])%*%t(Xj[[j]])%*%Yj[[j]]
V.j <- s2*solve(t(Xj[[j]])%*%Xj[[j]])
Lambda.j <- solve(solve(V.j) + solve(TAU))%*%solve(V.j)
Beta.tilde.j <- Lambda.j%*%Beta.hat.j + (diag(P) - Lambda.j)%*%Wj[[j]]%*%GAMMA
}
# Singular case
else {
Beta.tilde.j <- V.tilde.j%*%((1/s2)*t(Xj[[j]])%*%Yj[[j]] +
solve(TAU)%*%Wj[[j]]%*%GAMMA)
}
BETA.Js[[j]] <- t(rmnorm(1, mean=as.vector(Beta.tilde.j), V.tilde.j))
return(Beta.tilde.j)
}
BETA.Js <- mapply(update.betas,s2,Xj,Yj,TAU,Wj,GAMMA,SIMPLIFY=FALSE)
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
