Skip to content

rmultinom() -- how \\ via own C code?

3 messages · Martin Maechler, Ben Bolker, Kurt Hornik

#
I've had a need for multinomial "random number generation"
occasionally. And other people too.
The following code is currently in the 
(very small ``not very high importance'') CRAN package normix
 --- which I will rename to "nor1mix" very seen because of a 
    ``name registration'' problem

I want to add "this" (well the functionality) to a standard
package -- "mva" probably.

The reason for my post is to ask about importance for C code
(and a an official API to that).  
1) I think we have no clear precedence of a C interface to
   **multivariate** random numbers, and
2) It might be of too marginal importance.

If you (discussants) conclude that a C interface is not at all
desired, I consider using the R code as it is now -- which is
not optimal --- and most importantly:
If this should ever be changed, the change would hardly be
back-compatible. 
In other words, calling the not-yet existing C interface would
produce different random numbers... 

----

- Opinions ?

- Is there already C code around to do this ``in one step'' ?

Martin Maechler <maechler@stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><




## This is based on rmultz2() from S-news by Alan Zaslavsky & Scott Chasalow;
## in R available from library(combinat) -- but it has
## Arg.names like  rbinom();  returns  n x p matrix
rmultinom <- function(n, size, prob) {

    K <- length(prob) # #{classes}
    matrix(tabulate(sample(K, n * size, repl = TRUE, prob) + K * 0:(n - 1),
                    nbins = n * K),
           nrow = n, ncol = K, byrow = TRUE)
}
#
I use multivariate random values frequently -- multinomial (using code
contributed to the R-list by Ian Wilson), Dirichlet (using ratios of
gamma/sum(gamma)), multivariate normal (using mvrnorm() from MASS, which I
wish were called rmvnorm instead!)  I have functions for them in some of
my home-brewed packages.  I would vote for adding some of them to the
recommended packages.

  Many multinomial distributions have matrices or vectors of arguments: if
one wants to generate multiple variates with different parameters, how
should these arguments be stacked?  e.g. for multinomial deviates, should
one be able to specify an m by p matrix for "prob" where m is a divisor of
n?

  The other minor point is that I usually give d- and r-functions for 
these distributions but not p- and q- (since I don't really want to think 
about defining, much less computing, multidimensional cumulative 
distributions ...)

  Ben
On Mon, 27 Jan 2003, Martin Maechler wrote:

            

  
    
#
I have recently added code to r-devel to generate random two-way tables
with given row and column totals.  This works as

	r2dtable(n, r, c)

where r and c are the row and column totals, and returns a *list* of
length n.  If each draw gives a random 'vector' one can stack the
results into a matrix as does e.g. mvrnorm in MASS, or one could return
a list; in the general case, I think we always want a list of random
objects.
[Although e.g. simulating the distribution of random tables with given
marginals and doing computations on that would be rather nice ...]

Best
-k