I am moving this to r-devel.

The problem and solution below posted on r-help could have been
a bit slicker if %*% worked with multidimensional arrays multiplying
them so that if the first arg is a multidimensional array it is mulitplied
along the last dimension (and first dimension for the second arg).
Then one could have written:

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

which is a bit nicer than what had to be done, see below, given that %*% only
works with matrices.

I suggest that %*% be so extended to multidimensional arrays.  Note
that this is upwardly compatible and all existing cases would continue
to work unchanged.

On 7/16/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> The double loop is the same as:
>
>  Tbar[] <- matrix(tarray, 9) %*% t(wt) / rep(wti, each = 9)
>
>
> On 7/16/06, RAVI VARADHAN <[EMAIL PROTECTED]> 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]
> >  }
> >
>

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Reply via email to