[R] generate random numbers that sum up to 1

Grant Izmirlian izmirlian at nih.gov
Mon Oct 16 16:57:58 CEST 2006


So,  Alberto, you didn't see my post?  If Y has d independent components that 
are gamma distributed with common rate and shapes A_1, A_2, ..., A_d, then X, 
given by the components of Y divided by their sum has distribution 
Dirichlet(A_1, A_2, ..., A_d).  If you want Uniform on the d-simplex, then use 
A_1 = A_2 = ... = A_d = 1 (just as Duncan said)


original message:
--------------------------------
Duncan Murdoch's definition is _the_ only one that I know.  X is Uniform on A 
means  E phi(X) = \int_A phi(x)  dx / \int_A dx, so that the probability 
density is equal to 1/ \int_A dx everwhere on the set A.  

By the way, another way to simulate X ~ Dirichlet(A1, A2, ..., Ad) 
is to generate d independent gamma variables having equal rate parameter
(doesn't matter, so why not 1) and shape parameters  A1, A2, ..., Ad
Then the vector of components divided by their sum is the desired 
Dirichlet:


n <- 100000
d <- 3  # for three numbers that add to one ( the unit simplex in R^3)
A <- rep(1, 3)  # for uniform
X <- matrix(0, n, d)
for (k in 1:3)  X[,k] <- rgamma(n, shape=A[k], rate=1)
S <- X %*% rep(1, d)
Y <- X/S

Present example will simulate n independant 3 vectors, each having 
non-negative components summing to 1, and having a distribution
assigning equal mass to every possible value.

Changing d and the components of A will provide an arbitrary Dirichlet
on the unit simplex in R^d

Grant Izmirlian

NCI




>> Duncan Murdoch wrote "Another definition of uniform is to have equal
>> density for all possible vectors; the Dirichlet distribution with
>> parameters (1,1,1) would give you that. "



More information about the R-help mailing list