# In R, how to sample from the output of combn(a,b) if the “a choose b” is too large?

I would like to sample 1 million vectors from the set of vectors

$$\mathbf{v} = \{v_1,v_2,...,v_n \}$$

subject to

• $\sum_{i=1}^n v_i = 1$
• $\forall i: v_i \ge 0$
• each $v_i$ has at most 2 decimal places

As an example I have $n = 10$. Effectively, I am choosing to divide $100$ items into $10$ portion such that there is no minimum number of items that can be in any portion (meaning $0$ item in a portion is possible).

This is an interesting combinatorics (albeit elementary) problem. Suppose a rich man wants to divide his \$100 billion wealth amongst his 10 sons. Suppose he wants each son to get at least \$1 billion and the *i*th son should receive $n$ billion where $n$ is an integer greater than 0. Imagine he has stacked his money in \$1 billion piles in a row with some gap between each pile. Then he would have 99 gaps between the piles. To divide his money he has to simply choose 9 cut points from the 99 gaps. So it's 99 choose 9. It general if he has \$c billion and d sons, it would be c-1 choose d ways to divide his wealth. Now suppose he wants to divide his wealth so that there is no minimum and the *i*th son could get $n_i$ billion where $n_i$ is an integer greater than or equal to 0. How many ways are there to do this? He was clever and asked his friend to lend him \$10 billion and he also stacked them into a row together with his other billions. Now there are 110 piles of \$1 billion dollars in a row and there are 109 gaps. He chooses 9 cut points from these 109 gaps. Now he has divided \$110 billion into 10 portions each with at least \$1 billion. Now he retrieves \$1 billion from each portion, and returns the \$10 to his friend. Now the original \$100 billion is divided into 10 piles in such a way that there is no minimum amount! So$109$choose$9$would give me the combinations. In R, combn(a,b) would give me all the cuts points from all possible combinations "a choose b". However running combn(109,9) is not possible on my computer. In R, is there an efficient way to sample from the output of combn? • Your question is unclear. It sounds like you want your v to be a vector of random probabilities, ie they should be nonnegative and sum to 1. Is this correct? – Hong Ooi Aug 14 '13 at 10:46 • why is it taking 9 from 109? The sequence { 0.00, 0.01, 0.02 to 1.00 } has 101 values. I would have thought 9 from 101, but I must be missing something. – TooTone Aug 14 '13 at 13:41 • To answer your final question, try something like replicate(10,sample(109,10)). Though you would need to check for duplicates, they only occur in this example roughly once every 43 trillion replications. I'm not sure how this answers your first part of the question though. – James Aug 14 '13 at 13:41 • To compute combn(a,b) when it is too big, you could try to modify the combn() function to get a function dealing with big integers provided by the gmp package. – Stéphane Laurent Aug 14 '13 at 15:25 • I would try something like N=10; sapply( 1:N, function(i) { diff( c(0, sort(sample(seq(0, 1, by=0.01), 9, replace=TRUE)), 1) ) } ). This is a similar approach to @James. However I think that sampling with replacement is required, because you say zero items in a portion is possible. – TooTone Aug 14 '13 at 15:55 ## 3 Answers There are a couple of good answers here already, but none which directly implement the problem as suggested in the question. The steps to do this are: • Sample 9 out of the 109 numbers from 0.01 to 1.09: sample(seq(0.01, 1.09, by=0.01), 9) • Find the differences between them: • Sort them: sort • Place them in a vector between the minimum and maximum possible values: c(0, ..., 1.1) • Difference: diff • Remove 0.01 from each resulting value, so that 0s are possible: - 0.01 To do this N times, the code can be wrapped up in sapply. The resulting code is: N = 1e6 y <- sapply( 1:N, function(i) { diff( c(0, sort(sample(seq(0.01, 1.09, by=0.01), 9)), 1.1) ) - 0.01 } ) mean(y) hist(y, breaks=seq(-0.005, 1.005, by=0.01), freq=FALSE, main="Histogram of partial sums", xlab="Value") Taking the mean and producing a histogram is suggested by @whuber and yields a mean of 0.1 and the following histogram, consistent with his answer. On my latop (core-i7, 2.3Ghz), wrapping a system.time call around y->... gave an elapsed time of 124s, i.e. just over 2 minutes. • A big +1: that is a clear, simple implementation indeed. For larger partitions it is very efficient. Well done! (I would retract my solution as being overly complicated except that I recall it implements a generalization of the problem, so maybe it would have value for someone.) – whuber Aug 15 '13 at 21:47 • @whuber Thanks! -- although a large part of the credit goes to the OP for suggesting the elegant cutpoint approach. Also I agree your approach should remain -- arguably it is the only one that belongs here not on stackoverflow. Two comments: (1) interesting that your approach isn't much slower than mine; (2) your approach of plotting your results was extremely helpful, as it illustrated that a second approach I tried wasn't the optimization that it appeared to be. – TooTone Aug 15 '13 at 21:58 • Use replicate is better than sapply? – xiaodai Aug 16 '13 at 4:14 Based on the edit it sounds like what is desired is to generate, independently and uniformly at random, a partition of the$N=100$-vector$(1/100, 1/100, \ldots, 1/100)$into$k=10$pieces and to output the sums of elements in each piece. More generally we might ask how to do this for an arbitrary vector$x$. For example, the solution below can be used in the form set.seed(3) z <- rpart(10^(8:0), 6) sapply(z$sums, function(i) sprintf("%09d", i)); z$cuts to generate the$9$-vector$(100000000, 10000000, 1000000, 100000, 10000, 1000, 100, 10, 1)$and partition it into six (possibly empty) parts. Its output--which although random is reproducible by setting the seed to$3$--is [1] "000000000" "100000000" "011100000" "000011100" "000000000" "000000011" [1] 0 1 4 7 7 This partition was generated by cutting the vector at positions$0, 1, 4, 7,$and$7$. Visually this can be represented by placing vertical bars at both ends of the vector and the (random) bars inside to designate the cuts, $$(\vert\ \vert x_1 \vert x_2, x_3, x_4 \vert x_5, x_6, x_7 \vert \vert x_8, x_9\ \vert).$$ The partial sums are$0$(for the initial empty cut between the first two bars),$100000000$for the unary sum$x_1$,$011100000$for the sum$x_2+x_3+x_4$,$000011100$for the sum$x_5+x_6+x_7$,$0$for the empty sum between two successive bars following position$7$, and$000000011$for$x_8+x_9$. The method to generate these is a standard mechanism for generating random permutations of length$k-1$(the for loop below): to permute an$n$-vector$a$,$a_i$is swapped with$a_j$where$j$is uniformly selected in the range$[i, n]$inclusive. This is done for$i=1, 2, \ldots, k-1$. The prefix$(a_1, a_2, \ldots, a_{k-1})$contains the desired random permutation. The rest is just reprocessing the data into the desired form of output. To accommodate the possibility of adjacent cuts (leading to some partial sums of zero), we use a well-known method of considering the$k-1$internal vertical bars as if they, too, were elements of$x$. This extends$x$from an$N$-vector to an$n=N+k-1$-vector. The random permutation selects the$k-1$positions among$1, \ldots, n$at which the internal vertical bars will be placed. The post-processing subtracts$1$from the first index,$2$from the second, ..., and$k-1$from the last so that they point to the positions within$x$itself that occur just before each cut. To represent the two outside vertical bars it adjoins$0$and$N$to these indexes (but does not explicitly output them). A difference of cumulative sums (via cumsum and diff) is an R-centric way to obtain the partial sums within each partition. As a check, let's generate (say)$10,000$such partitions and look at histograms of (a) the$10,000 \times 10=100,000$partial sums in the original problem and (b) the cut indices. We expect the former to average exactly$10 \times 1/100 = 1/10$and the latter to be uniformly distributed among all possible values$0, 1, \ldots, 100$. First the code to do it: x <- rep(1/100, 100); k <- 10 par(mfrow=c(1,2)) set.seed(83) sample <- replicate(1e4, unlist(rpart(x, k)$sums))
hist(sample, breaks=seq(-0.5, 100.5, 1)/100, freq=FALSE,
main="Histogram of partial sums", xlab="Value")
mean(as.vector(sample))

set.seed(83)
sample <- replicate(1e4, unlist(rpart(x, k)$cuts)) hist(sample, breaks=seq(-0.5, 100.5, 1), main="Histogram of cutpoint indices", xlab="Index") mean(as.vector(sample)) Now the results: The figures look appropriate. The left hand histogram in particular provides a reference for any other attempted solution: if it does not produce the same distribution of partial sums (up to sampling error), then that solution is not generating uniformly distributed partitions. By the way, although this code is not well-suited for R implementation, it has a reasonable speed, taking about$N\times k / 50000$seconds on this machine (two seconds for the example). That would be three to four minutes to solve the problem posed in the question with$N=10^6$iterations and$k=10$. Because these partitions can be generated one at a time, there is no RAM limitation. ### R code # # Generate a random partition of vector x into k>0 (possibly empty) pieces # and return the vector of k partial sums of x within those pieces as well as # the k-1 cutpoints between them. (Each partial sum starts *after* the # preceding cut and ends at the current cut.) # # Output is a list with the sums and the corresponding cutpoints. # rpart <- function(x, k) { if (k <= 1) return(list(sums=sum(x), cuts=numeric(0))) n <- length(x)+k-1 a <- 1:n for (i in 1:(k-1)) { j <- sample(i:n, 1) y <- a[i]; a[i] <- a[j]; a[j] <- y } cuts <- c(0, sort(a[1:(k-1)]) - (1:(k-1)), length(x)) psums <- c(0, cumsum(x)) return(list(sums=diff(psums[cuts+1]), cuts=cuts[2:k])) } To hold the set of vectors that combn(109,9) would return if it could, by my estimate you would need upward of 32 petabytes of memory, which exceeds how much R can allocate. Are the stated criteria the only ones that matter to you? Or do the elements of each vector also have to be unique? I went ahead and assumed yes. If so, this should do it (in probably less than half an hour, I only tried it for 1e5 vectors). foo<-replicate(1e6,{oo<-0;while(sum(oo)!=1||length(unique(oo))!=length(oo)) oo<-round((ii<-runif(10))/sum(ii),2);oo}) Obviously, for vectors of 20 elements, substitutes runif(20) and so on. Now, if you don't care about uniqueness, the following will run much faster (a little over a minute on my computer): foo<-replicate(1e6,{oo<-0;while(sum(oo)!=1) oo<-round((ii<-runif(10))/sum(ii),2);oo}) If you don't care about them adding up to exactly 1, the following should also run faster: foo<-replicate(1e6,{oo<-0;while(length(unique(oo))!=length(oo)) oo<-round((ii<-runif(10))/sum(ii),2);oo}) And if you just want real numbers that add up to 1, and put in the stipulation of two decimal places because you thought it would make the task computationally easier, it doesn't. The easiest and fastest way to generate vectors that each add up to 1 is this (21 seconds on my computer): foo<-replicate(1e6,(ii<-runif(10))/sum(ii)) Actually, this is even faster (less than two seconds): foo <-matrix(runif(1e7),nrow=10); foo<-foo/rbind(colSums(foo))[rep(1,10),] ...the uniqueness and summing up to exactly 1 take care of themselves when you relax the two decimal places requirement. By the way, in all the above examples, foo is a matrix, each of whose columns is a vector$v_k$. I'm sorry if the question I answered was not the one you meant to ask. Most people on the StackOverflow sites will generally just flame you for not providing all relevant information (as if it's always obvious what is relevant to answering a question and waht isn't). So, I believe that it's more constructive to instead fill in the blanks with what I would reasonably expect you to mean, and answer that question, thus either guessing right or sending you a positive message to refine the question. Now I have a question, for my own curiousity. What procedure or simulation could possibly require numbers that fit these criteria? Are you trying to empirically get the probability density function of a uniform distribution or something? Because there's probably an easier way. • Could you describe what you think your code does? Apart from the rounding to two decimal places, I cannot see any connection between it and the question (which may largely be the fault of the question, but some clarification would be gladly received). – whuber Aug 14 '13 at 22:00 • It generates a vector random numbers from a uniform distribution over (0,1) and then divides the vector by its sum. Voila, they sum to 1 and are greater than 0. But when I round them to 2 digits (I'm SO curious why someone would want to), they are no longer guaranteed to be unique and sum to 1. So I have to introduce a while statement to check for that and not return until both conditions are satisfied. I think OP is misunderstanding the purpose of combn, like I did when I first encountered that function. – f1r3br4nd Aug 15 '13 at 0:00 • You can use a Dirichlet distribution for the continuous version and a multinomial distribution for the discrete version of that: see cran.r-project.org/web/views/Distributions.html for a list of packages that provide these distributions. But it's not at all clear that either is what the O.P. is looking for. – whuber Aug 15 '13 at 1:42 • In credit risk scoring models practitioners often speak of "weightage", which is nothing but$std(coef_i*raw_i) / (\sum_i{coef_i*raw_i})$where$raw_i$is the empirical data for the ith explanatory variable in the model and$coef_i$is the coefficient you want to apply to that variable. Now imagine if the$raw_i$has been standardised so that they have standard deviation of 1 and a mean of 0. Determining how to assign the$coef_i$is simple - simply choose$coef_i\$ such that they sum to 1. In scoring we are interesting the predictive power and we do not want to overfit, so limit to 2 decimals. – xiaodai Aug 15 '13 at 2:39
• In the question, it was suggested that the problem was equivalent to choosing cutpoints. I thought this was quite an elegant approach (and computationally quite tractable -- several comments, including myself suggested a few lines of R code to do this). Why is your approach preferred? – TooTone Aug 15 '13 at 9:26