Re: [S] How do you quickly generate multinomials with different probabilities?

Dave Krantz (dhk@paradox.psych.columbia.edu)
Sun, 14 Jun 1998 12:10:37 -0400


To restate the problem, suppose P is an n x k matrix with non-negative
entries whose rows sum to 1; how can one generate a sample of size n,
where the ith element has value 1, 2, ..., k with probabilities given
by the ith row of P?

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