Generate three random numbers that sum to 1 in R

I would like to generate three random numbers and then standardize them so that they add up to 1.

I would like to repeat this procedure so that in the long run the mode is .33 for each number.

• you mean a2<-a1/sum(a1) ?
– ghb
May 15 '13 at 10:24
• More information is needed because there are many different solutions. How should the individual numbers be distributed? How should they be correlated? Should there be any restrictions on their values (such as lying between $0$ and $1$)?
– whuber
May 15 '13 at 11:31
• We can take a tip from XKCD (and, earlier, from Dilbert) and generate the three random numbers as $(1/3, 1/3, 1/3)$. They satisfy all the conditions of the question in its current form (taking ".33" to be an approximation of $1/3$).
– whuber
May 15 '13 at 13:41
• @user603 How do you demonstrate your three variables have modes of $1/3$? Numerical analysis indicates their modes are near $0.10586$ and this is supported with a large-scale simulation (of $10^9$ replications).
– whuber
May 16 '13 at 16:48
• The votes to close this as "not a real question" are legitimate. However, the many answers--both correct and incorrect!--contain stimulating and useful insights. It would be nice to find a way to keep this thread open.
– whuber
May 21 '13 at 15:28

The mode is a bit of a red herring. Here is a very simple solution to this problem that circumvent the need to define the mode precisely. I'm surprised it has not been proposed earlier. The constraint on the mode can be easily satisfied by drawing samples from a symmetric distribution and scaling them suitably:

$$(x_i,y_i,z_i)\sim\mathrm{i.i.d.}\;\mathcal{L}(\mu,\sigma)$$ $$(x_i^*,y_i^*,z_i^*)=\left(\frac{x_i}{x_i+y_i+z_i},\frac{y_i}{x_i+y_i+z_i},\frac{z_i}{x_i+y_i+z_i}\right)$$

where $\mathcal{L}(\mu,\sigma)$ is a symmetric distribution (so that the mean, the mode and the median are the same) and chosen such that the probability mass below 0 is 0. For example, picking $\mathcal{L}(\mu,\sigma)$ to be $\mathrm{Beta}(2,2)$:

a1 <- matrix(rbeta(100*3,2,2), nc=3)
a1 <- sweep(a1, 1, rowSums(a1), FUN="/")
colMeans(a1)
# [1] 0.3342165 0.3341534 0.3316301


yielding the desired solution

sum(colMeans(a1))
# [1] 1


If X1, X2, and X3 are i.i.d. Gamma(a) then {X1,X2,X3}/(X1+X2+X3) will be Dirichlet(a,a,a).

If a>1 then the mode will be 1/3. The peak will be sharper for larger values of a.

• I looked this up too, and Wikipedia says that's the mode, but is that really correct? I haven't thought about the theory, but in the simulations I did, it appeared to be less than 1/3. May 15 '13 at 16:48
• @Aaron You and Wikipedia are both right. The modes of the two marginals of a Dirichlet$[n,n,n]$ for $n\gt 1$ are at $(n-1)/(3n-2) \ne 1/3$. The (very interesting) error made here is that the components of the mode of a $d$-dimensional distribution are not necessarily the modes of its marginals.
– whuber
May 15 '13 at 16:57
• Ah, yes, of course. One of those things that is clear only after someone points it out. Thanks. May 15 '13 at 17:06
• R code would be library(gtools); rdirichlet(n,c(a,a,a)) May 15 '13 at 19:10

Here is an approximate numerical answer. It can easily be made more precise.
Let $\{U,V,W\} = {X,Y,Z}/(X+Y+Z)$, where $X,Y,Z$ are i.i.d. with a trapezoidal density on $[0,1]$:

$f(x)=1+a-2ax.$ $U,V,W$ will have identical marginals.
Given a numeric 'a', I used Mathematica to get the cdf of $U$:

F[u_] = Assuming[0 < u < 1, Simplify@Integrate[
Boole[x < u(x+y+z)] f[x] f[y] f[z], {x,0,1},{y,0,1},{z,0,1}]


Differentiating $F$ twice, setting the result to zero, and solving the resulting 7th degree polynomial gave the mode. I used a binary search to refine the value of 'a'. I used exact arithmetic throughout, up to the point of solving the polynomial.

  a      mode

1    .318182
7/8   .322065
13/16  .327099
25/32  .330465
49/64  .332373
97/128 .333376 <-- close enough?
3/4   .334221
1/2   .353738
0    .359187

• Doing the integration symbolically and setting the second derivative to zero when u = 1/3 gives us the desired 'a' as a root of ((5a - 48)a + 81)a = 36, or approximately .758206285092054 May 17 '13 at 8:22
• To generate data from a trapezoidal distribution, the R function trapezoid (package of the same names) requires an upper and a lower mode. Would you know what these are in your parametrization? (i'd like to give your solution a try) May 17 '13 at 8:44
• The parallels in my trapezoid are vertical; in the package they're horizontal, so it's a different distribution. To get a random value from my distribution, generate x,y iid U(0,1); if y > 1+a-2ax then use 1-x else use x. May 17 '13 at 16:29

Analytically: Given a joint pdf for $X$,$Y$, and $Z$ $f_{X,Y,Z}(x,y,z)$, if they are iid, then $f_{X,Y,Z}(x,y,z)=f_X(x)f_Y(y)f_Z(z)$, where $f_X(x)=f_Y(y)=f_Z(z)$. You have to find the pdf $$f_{X,Y,Z}\left(\dfrac{x}{x+y+z}\right)$$ After differentiating and equal to zero, you'll find your mode. Obviosly, mode will depend on $f_{X,Y,Z}(x,y,z)$ and consequently on $f_X(x)$, $f_Y(y)$ and $f_Z(z)$.

Numerically: sample three random numbers from your preferred distribution, standardize them and save them in the $i^{th}$-row an Nx3 matrix. Repeat this procedure N times and plot the frequencies of each column.

The analytical solution is preferred instead of trying to demonstrate it from random samples in R, which would be just a numerical approximation.

• How, then, do you go about choosing $f$ to make the modes all equal $1/3$? Could you provide a simple example?
– whuber
May 16 '13 at 4:37
• How can they be independent if they are constrained to sum to 1? May 16 '13 at 12:58
• @Soakley $X$, $Y$, and $Z$ are not constrained to sum to $1$. There is an implicit assumption that the essential support of $f$ is such that the chance of $X/(X+Y+Z)$ etc. being undefined is zero. For instance, assuming that each of $X$, $Y$, $Z$ has zero probability of being non-positive will do. With that assumption, obviously $X/(X+Y+Z)+Y/(X+Y+Z)+Z/(X+Y+Z)=1$ identically, even when $X$, $Y$, and $Z$ are iid--which is an explicit assumption here. The difficulty with this answer lies in actually carrying out the computations it proposes: except in simple cases, they are intractable.
– whuber
May 16 '13 at 15:58

It is still unclear whether the OP wants a solution with a mode of 0.33 or $1 \over 3$ or a mean with one of those two values. Without knowing the exact need, there are multiple possibilities. [1] and [4] below address the problem of getting a mean of ${1 \over 3},$ while [2] and [3] are for a mode of ${1 \over 3}$.

[1] Generate $U_1, U_2, U_3$ as continuous uniform random variates on $[0,{2 \over 3}]$. Let $A = {{U_1} \over {U_1 + U_2 + U_3}},$ $B = {{U_2} \over {U_1 + U_2 + U_3}},$ and $C = {{U_3} \over {U_1 + U_2 + U_3}}.$

[2] Let $X=(1/3)*W_1 + (1/6),$ $Y = (1/3)*W_2 + (1/6),$ and $Z = 1 - X - Y,$ where $W_1$ and $W_2$ are continuous uniform on $[0,1].$ $Z$ has a different distribution than $X$ or $Y,$ but I think this will meet the original mode requirement.

[3] In an effort to produce a more intuitive version, let $R_1$ be right triangular with left endpoint at zero and mode at ${1 \over 3} .$ Let $L_1$ be left triangular with mode at $1 \over 3$ and right endpoint at ${2 \over 3} .$ Then $R_1 + L_1$ has a unique mode at ${2 \over 3},$ and if we define $Q = 1 - (R_1 + L_1)$ then $Q$ has a unique mode at ${1 \over 3}.$ The pdf of $Q$ is given below in the comments.

[4] Trying to be clever, simple, and elegant. Generate 2 independent uniform[0,1] realizations. These 2 points divide the interval from 0 to 1 into 3 pieces. Use the lengths of these 3 pieces as the desired variates. Note how this generalizes intuitively to any sum and any number of random variables. Each variate is identically distributed (another right triangular distribution). Every pairwise correlation is $-{1 \over 2}.$ However, like approach [1], the mean is ${1 \over 3},$ but the mode here is not ${1 \over 3}.$ As whuber noted, it is at zero.

• This is incorrect: the modes are all at the root of $6 x^4-27 x^3+33 x^2-17 x+3$ lying near $0.359187$ and are not equal to $1/3$. This error is interesting in how it reveals one of the subtleties of the question.
– whuber
May 15 '13 at 13:38
• Using a simulation with 10,000 variates, I'm seeing results much closer to $1/3$ than your root for about 200 iterations. But I am calculating the average, not a mode. May 15 '13 at 14:15
• @Whber : Could you explain us how you reach the equation ? And do you have a practical solution example as soakley's ? piece of code maybe ? Thx May 15 '13 at 14:34
• @D.Khireche My result is exact. It is difficult to identify a mode using simulations; even with $10^7$ iterations, the mode can only be pinned down to the $0.36$ range. $200$ or even $10000$ iterations is nowhere near enough to say anything about the mode. Here's R code (2 seconds total): n <- 10^7; x <- matrix(runif(3*n), nrow=n); sim <- x/as.vector(x%*%matrix(1,3)); apply(sim*1000, 2, function(y) (-1/2+which.max(tabulate(ceiling(y))))/1000). I won't offer a "practical solution" until this question has been clarified--until then we're just guessing what is being asked for.
– whuber
May 15 '13 at 15:18
• I have been able to identify the mode as not being 1/3 even with my 200 iterations. That is why I've modified my answer and proposed another one. May 15 '13 at 15:33

Here is a simple approach. Generate $X \sim \mathrm{Unif} [0,{2 \over 3}].$ Let $$Y = \begin{cases} X+{1/3} \ , & \text{if} \ X \le {1/3} \\ X-{1/3} \ , & \text{if} \ {1/3} \lt X \le {2/3} \end{cases}$$

Let $Z = 1 - X - Y.$

Then it's not too hard to show that $X,Y,$ and $Z$ are identically distributed with $\mathrm{Unif} [0,{2 \over 3}]$ distributions. So each has mean and median of ${1 \over 3}$ and has a mode there as well.

Additionally, all pairwise correlations are $-{1 \over 2},$ and only one call to a uniform random generator is needed to get the three variates.