applying rmultinom in R without iterating over matrix

I want to use rmultinom(), combined with a transition matrix, to generate whole number outputs. However, I can’t figure out how to do it without iterating over the matrix. Here is an example:

a = matrix(runif(16),nrow=4,ncol=4)
a = apply(a,2,FUN = function(x) x/sum(x))

b = c(10,10,10,10)
out = c(0,0,0,0) # initialize
for (i in 1:ncol(a)){
  tmp = rmultinom(1,b,a[,i])
  out = tmp + out
}

a represents a transition matrix, with each column summing to 1. b is a starting vector of integers. The loop iterates along the columns to generate a vector in out that sums to the initial numbers in b. How can I do this without using a loop? The results would be similar to if I multiply a %*% b, but this leaves me with floating point values.

>Solution :

You could do apply and rowSums (this will be stochastic):

library(magrittr)
set.seed(1)
a = matrix(runif(16),nrow=4,ncol=4)
a = apply(a,2,FUN = function(x) x/sum(x))

b = c(10,10,10,10)
out = c(0,0,0,0) # initialize

out <- apply(a, 2, function(.x) rmultinom(1, b, .x)) %>% 
  rowSums()
out
[1] 10 10 12  8
sum(out)
[1] 40

You can also round the answer of the matrix multiplication (which will be deterministic):

out <- round(a %*% b)
out
     [,1]
[1,]   11
[2,]    7
[3,]   11
[4,]   11
sum(out)
[1] 40