On Oct 14, 2009, at 12:58 AM, Stephen Samaha wrote:
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)
(It would be more courteous to offer the error messages you are
currently keeping secret.)
Could it be because you are offering a mix of vectors and scalars to
mapply without properly segregating them? See the help page for mapply
and the moreArgs argument.
--
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
______________________________________________
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.
