Generate random matrix, subject to equality constraints in R

Question

I would like to generate an m by n random matrix, drawing from a uniform distribution for each element, subject to the constraint that each column must sum to one. I would also like to set the lower and upper bound constraints for each element.

Here is an example of the desired output (although not truly random), mymatrix:

dput(mymatrix)
structure(c(0.5, 1.59365940317881e-17, 1.89403184501105e-18, 
0.5, 8.96677015179615e-19, 1.7326456004408e-18, 0.499999858832886, 
0.499999838867718, 2.22092666609155e-08, 2.69754277003262e-07, 
2.78754819399023e-11, 1.0307976832927e-08, 0.42857144673908, 
0.499999904492316, 3.75555965298514e-08, 4.75046037021098e-08, 
3.3675224155978e-10, 0.0714285633716491, 0.357142848840934, 0.499999941890735, 
7.30974537548529e-08, 2.67022430383198e-08, 2.80475367418637e-10, 
0.142857109188155, 0.28571428969238, 0.499999933476062, 5.02671502622386e-08, 
3.11835918349941e-08, 1.10520352747797e-09, 0.214285694275616
), .Dim = c(6L, 5L), .Dimnames = list(NULL, NULL))

colSums(mymatrix)
#[1] 1 1 1 1 1

max(mymatrix)
#[1] 0.5

In this particular example, note that no element is less than zero or greater than .5; ideally these bounds could be set as variables. I would also like the number of rows and columns generated to be set as a variable. However, the constraint that all columns should sum to 1 should be fixed. Any help is greatly appreciated.

Answer 1

You can use the dirichlet distribution from the MCMCpack library. All alpha parameters should be set to 1.

library(MCMCpack)
mat=rdirichlet(n_columns,rep(1,n_rows))
mat=t(mat) #rotated because you want columns to add to one, not the rows

I'm not sure how you go about setting limits on how big elements can be (like the 0.5 you set above), but if computation time isn't imperative, I would just throw away any columns that have too large of elements. All you'd need to do is run the code

mat=mat[,which(apply(mat,2,max)<0.5)]

To get rid of the unwanted columns.