R help - Jul 2006 - Manipulation involving arrays

Hi,

I have the following piece of code that is part of a larger function.  This
piece is the most time consuming part of the function, and I would like to make
this a bit more efficient.  Could anyone suggest a way to do this faster?

In particular, I would like to replace the nested "for" loop with a
faster construct.  I tried things like "kronecker" and
"outer" combined with apply, but couldn't get it to work.


Here is a sample code:

 ##########################
 n <- 120
 sigerr <- 5
 covmat <- diag(c(8,6,3.5))
 mu <- c(105,12,10)
 mcsamp <- 10000
 
 Tbar <- array(0, dim=c(3,3,n))
 
 # theta is a mcsamp x 3 matrix
 theta <- mvrnorm(mcsamp, mu = mu, Sigma = covmat)
 
 wt <- matrix(runif(n*mcsamp),n,mcsamp) 
 wti <- apply(wt,1,sum)
 
 tarray <- array(apply(theta,1,function(x)outer(x,x)),dim=c(3,3,mcsamp))
 
 for (i in 1:n) {
 for (k in 1:mcsamp) {
 Tbar[,,i] <- Tbar[,,i] + wt[i,k] * tarray[,,k]
 } 
 Tbar[,,i] <- Tbar[,,i] / wti[i]
 }

###############################################

Thanks very much,
Ravi.

The double loop is the same as:

 Tbar[] <- matrix(tarray, 9) %*% t(wt) / rep(wti, each = 9)


On 7/16/06, RAVI VARADHAN <rvaradhan at jhmi.edu>
wrote:
> Hi,
>
> I have the following piece of code that is part of a larger function.  This
piece is the most time consuming part of the function, and I would like to make
this a bit more efficient.  Could anyone suggest a way to do this faster?
>
> In particular, I would like to replace the nested "for" loop with
a faster construct.  I tried things like "kronecker" and
"outer" combined with apply, but couldn't get it to work.
>
>
> Here is a sample code:
>
>  ##########################
>  n <- 120
>  sigerr <- 5
>  covmat <- diag(c(8,6,3.5))
>  mu <- c(105,12,10)
>  mcsamp <- 10000
>
>  Tbar <- array(0, dim=c(3,3,n))
>
>  # theta is a mcsamp x 3 matrix
>  theta <- mvrnorm(mcsamp, mu = mu, Sigma = covmat)
>
>  wt <- matrix(runif(n*mcsamp),n,mcsamp)
>  wti <- apply(wt,1,sum)
>
>  tarray <- array(apply(theta,1,function(x)outer(x,x)),dim=c(3,3,mcsamp))
>
>  for (i in 1:n) {
>  for (k in 1:mcsamp) {
>  Tbar[,,i] <- Tbar[,,i] + wt[i,k] * tarray[,,k]
>  }
>  Tbar[,,i] <- Tbar[,,i] / wti[i]
>  }
>