[R] using mapply to avoid loops
David Winsemius
dwinsemius at comcast.net
Wed Oct 14 15:15:41 CEST 2009
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
More information about the R-help
mailing list