First convert P to cumulative probabilities within each row.
An efficient way (without external loops) is:
Q <- matrix(cumsum(t(P)) - rep(0:n, rep(k,n)), byrow=T, ncol=k)
Next, generate runif(n) and compare it with each column of Q, summing
the rows of the comparison matrix (and adding 1):
R <- 1 + (Q < runif(n)) %*% rep(1, k)
R is the result you want. This can easily be adapted to give samples
of specified sizes greater than 1 for each multinomial: if Z is the
vector of n sample sizes, replace "Q < runif(n)" by
Q[rep(1:n, Z), ] < runif(sum(Z))
and then make R into a list by using split() (if Z is nonconstant) or
into an n x Z[1] matrix (if Z is constant).
Dave Krantz
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news