Hi, I am a newbie at R and was thinking of ways to vectorize this loop but can't quite figure out if anything more can be done : Here the dimensions of the arrays are (N=2000, K=2, D=1000): responsibilities : N * K pi: K corp: N*D theta: K*D for (n in 1:N) { for (k in 1:K) { responsibilities[n, k] <- (pi[k]) + (log1p(dmultinom(x=corp[n, ], prob=exp(theta[k, ])))) } denom <- log.sum(responsibilities[n, ]) responsibilities[n,] <- responsibilities[n, ] - denom } Any help is much appreciated! thanks. [[alternative HTML version deleted]]