10
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ you mean a2<-a1/sum(a1) ? $\endgroup$
    – ghb
    May 15, 2013 at 10:24
  • 3
    $\begingroup$ 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$)? $\endgroup$
    – whuber ♦
    May 15, 2013 at 11:31
  • 2
    $\begingroup$ 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$). $\endgroup$
    – whuber ♦
    May 15, 2013 at 13:41
  • $\begingroup$ @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). $\endgroup$
    – whuber ♦
    May 16, 2013 at 16:48
  • $\begingroup$ 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. $\endgroup$
    – whuber ♦
    May 21, 2013 at 15:28

6 Answers 6

Reset to default
11
$\begingroup$

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
Improve this answer
$\endgroup$
5
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ 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. $\endgroup$
    – Aaron left Stack Overflow
    May 15, 2013 at 16:48
  • 5
    $\begingroup$ @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. $\endgroup$
    – whuber ♦
    May 15, 2013 at 16:57
  • $\begingroup$ Ah, yes, of course. One of those things that is clear only after someone points it out. Thanks. $\endgroup$
    – Aaron left Stack Overflow
    May 15, 2013 at 17:06
  • 2
    $\begingroup$ R code would be library(gtools); rdirichlet(n,c(a,a,a)) $\endgroup$
    – Aaron left Stack Overflow
    May 15, 2013 at 19:10
2
$\begingroup$

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
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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 $\endgroup$
    – Ray Koopman
    May 17, 2013 at 8:22
  • $\begingroup$ 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) $\endgroup$
    – user603
    May 17, 2013 at 8:44
  • 1
    $\begingroup$ 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. $\endgroup$
    – Ray Koopman
    May 17, 2013 at 16:29
1
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ How, then, do you go about choosing $f$ to make the modes all equal $1/3$? Could you provide a simple example? $\endgroup$
    – whuber ♦
    May 16, 2013 at 4:37
  • 1
    $\begingroup$ How can they be independent if they are constrained to sum to 1? $\endgroup$
    – soakley
    May 16, 2013 at 12:58
  • $\begingroup$ @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. $\endgroup$
    – whuber ♦
    May 16, 2013 at 15:58
1
$\begingroup$

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.

Improve this answer
$\endgroup$
13
  • 2
    $\begingroup$ 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. $\endgroup$
    – whuber ♦
    May 15, 2013 at 13:38
  • $\begingroup$ 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. $\endgroup$
    – soakley
    May 15, 2013 at 14:15
  • $\begingroup$ @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 $\endgroup$
    – dfhgfh
    May 15, 2013 at 14:34
  • $\begingroup$ @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. $\endgroup$
    – whuber ♦
    May 15, 2013 at 15:18
  • $\begingroup$ 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. $\endgroup$
    – soakley
    May 15, 2013 at 15:33
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.