Random triple with uniform marginals

Question

Output a random triple \$(x,y,z)\$ of real numbers such that

• \$x+y+z=\frac{3}{2}\$
• Each of \$x,y,z\$, taken alone, is uniformly distributed from \$0\$ to \$1\$.

Please explain in your answer what distribution your code produces and why it meets these conditions. In particular, please include a demonstration that with many samples produced by running your code many times, the distribution of each variable alone is approximately uniform.

Allowances for precision

Since computers can't produce actual real numbers, you can approximate them by outputting any of:

• Floats
• Decimals
• Rationals

I've tried to make the rules here as friendly as possible.

Your outputs need to be accurate within \$\pm 10^{-3}\$, that is 3 decimal digits. That means you can imagine that each real-number output from your true distribution has been fudged by at most \$\pm 10^{-3}\$ in any way. It's fine if this makes \$x+y+z\$ be a bit off from \$\frac{3}{2}\$.

You may output rationals as (numerator, denominator) pairs of integers or similar. You may also choose to output just one denominator that applies to all three values, that is, output \$a,b,c\$ and \$d\$ to represent \$\frac{a}{d},\frac{b}{d},\frac{c}{d}\$. Rationals don't have to be in reduced form.

Randomness

You can assume any standard PRNG is in fact random. For instance, if it outputs random floats from 0 to 1, you can assume these are uniformly random reals that are independent from each other, and not worry about any endpoint issues. It's also fine if your code fails in a probability-zero event.

To decide whether a pseudo-random source such as time or memory contents is sufficiently random, consider an empirical test where your code is run many times in sequence to check if the randomness varies sufficiently between runs. You may have these runs be done within your language such as in a loop, or by repeatedly calling the code from outside like with a shell command. You may not assume any intentional time delay or other actions taken between runs.

Answer 1

APL (Dyalog Unicode), 15 bytes^SBCS

1|(-,.5∘-,+⍨)?0

Try it online!

A full program that generates three floating-point numbers \$ X,Y,Z \$.

The distribution

Uses the same method as Luis Mendo's MATL answer, which in turn is from this Math.SE answer. The Math.SE answer presents the following random variables:

$$ X = \text{uniform random value in } (0,1) \\ Y = \begin{cases} X + \frac12, \quad X < \frac12 \\ X - \frac12, \quad X \geq \frac12 \end{cases} \\ Z = \frac32 - (X + Y) $$

I found that this is equivalent to

$$ X = \text{uniform random value in } (0,1) \\ Y = \{ X + 1/2 \}, Z = \{-2X\} $$

which, in turn, has the same distribution with

$$ r = \text{uniform random value in } (0,1) \\ X = \{ -r \}, Y = \{ 1/2 - r \}, Z = \{ 2r \} $$

where the notation \$ \{ x \} \$ means the positive fractional part of \$ x \$.

How the code works

1|(-,.5∘-,+⍨)?0
             ?0  ⍝ r←Generate a random number from (0,1)
  ( ,    ,  )    ⍝ Create a 3-element array from r:
   -             ⍝ -r
     .5∘-        ⍝ 0.5 - r
          +⍨     ⍝ r + r (= 2r)
1|               ⍝ Extract positive fractional part from each number

Test script

⍝ Take 10000 samples of X, Y, Z's
data←f¨⍳10000

⍝ Tests that X+Y+Z = 1.5 in all cases
∧/1.5=+/¨data

⍝ Shows the distribution of X,Y,Z in 10 bins
{(⊂∘⍋⌷⊢)↓,∘≢⌸⍵}¨⌊10×↓⍉↑data
                    ↓⍉↑data  ⍝ Reform the data into X's, Y's, and Z's
                ⌊10×         ⍝ Classify into 10 bins
{             }¨     ⍝ For each of X's, Y's, or Z's,
        ↓,∘≢⌸⍵       ⍝ Gather the list of (bin, count)
 (⊂∘⍋⌷⊢)             ⍝ Sort the bins

⍝ Example result
X:  0 1012  1 1016  2 993   3 1025  4 1002  5 969   6 1002  7 977  8 1049  9 955  
Y:  0 969   1 1002  2 977   3 1049  4 955   5 1012  6 1016  7 993  8 1025  9 1002 
Z:  0 1013  1 944   2 1022  3 1052  4 990   5 980   6 1039  7 979  8 959   9 1022 

---

Alternative method, 21 bytes^SBCS

(3*9),3⊥¨3|(⍳3)+⊂?9/3

Try it online!

This one uses the OP's formulation from the same Math.SE question:

1. Generate a uniform random number \$ X \$ from \$ (0, 1) \$.
2. Take its ternary digits \$ X = 0.X_1 X_2 X_3 \cdots \$ and generate ternary digits of \$ Y \$ and \$ Z \$:

$$ Y_i = (X_i + 1) \% 3 \\ Z_i = (X_i + 2) \% 3 $$

3. Then the resulting numbers \$ X,Y,Z \$ are uniform in range \$ (0, 1) \$, and their sum is

$$ \frac3{3^1} + \frac3{3^2} + \frac3{3^3} + \cdots = \frac32 $$

Theoretically it requires infinite number of digits; the code approximates the scheme with just 9 digits, which makes \$ X+Y+Z=1.499923792 \$. I'm not exactly sure if it meets the OP's spec though.

How the code works

(3*9),3⊥¨3|(⍳3)+⊂?9/3
                 ?9/3  ⍝ Generate 9 ternary digits
         3|(⍳3)+⊂      ⍝ Ternary digits for X,Y,Z
      3⊥¨  ⍝ Convert each list of digits to numerator for the fraction
(3*9),     ⍝ Prepend the denominator

Answer 2

MATL, 18 17 11 bytes

rt.5+yE_v1\

The code is based on this method, and incorporates Bubbler's improvement, which saved 6 bytes.

Try it online! Or see the three marginal histograms for 1000 realizations, using 10 bins on the interval (0,1).

Commented code

r       % Push random number, r, uniformly distributed on (0,1)
t.5+    % Duplicate, add 0.5. Gives r+0.5
yE_     % Duplicate from below, multiply by 2, negate. Gives -2*r
v       % Concatenate the three values into a column vector
1\      % Modulo 1. Implicit display

Answer 3

Wolfram Language (Mathematica), 51 48 bytes

(.5+Sqrt[3/8]*(#-Mean@#))&@RandomPoint@Sphere[]

Try it online!

Theory

According to this post on MathOverflow about a cool unsolved problem, referencing this other post also on MathOverflow about distributions with nice projections, the following procedure should provide a distribution satisfying the demands of this question:

1. Pick a random point on a sphere.
2. Project (orthogonally) to a fixed plane. The support of the distribution is now a disc, but with a non-uniform distribution.
3. Inscribe this disc in an equilateral triangle. Take the barycentric coordinates of our projected random point within this equilateral triangle.
4. These three coordinates now sum to a constant and have a uniform distribution on their support. Rescale if necessary.

The proof that this works fundamentally comes down to Archimedes's "Hat-Box" Theorem, which is the assertion that a single coordinate of a random point on a sphere has a uniform distribution on its support.

In the solution above, I use the plane x+y+z=0 because it's easy to project to it by subtracting off the mean (one-third of the total x+y+z) from each coordinate. Then I simply need to shift and scale appropriately (I can even skip the barycentric coordinates, which are themselves just a shift and scale here), which a bit annoyingly involves Sqrt[3/8]. Perhaps I can save a few bytes by taking an approximate value to the tolerances of the original question.

However, I hope that this answer -- or variants of it -- might be efficient to express in some language.

Practice

You can run a Kolmogorov-Smirnov test to see if the distribution is correct. Try it online to compare the distribution of the first coordinate of this answer (which has the same distribution as the other two coordinates) with the uniform distribution. Put briefly, we want the result of this test to not be too "close to 0" ("statistically significant").

Remarks

Fun fact: if you scale the coordinates so that their sum is 1, this can be used to solve a nice problem from the xkcd forums.

Two generals, who each control the same number of units, each secretly decide how to split their units among three battle fronts. Whoever sends more units to a front wins that front, and whoever wins 2 out of 3 fronts wins overall. Margins of victory are irrelevant. What is the optimal strategy , pure or mixed?

Answer 4

JavaScript (ES6), 38 bytes

Also derived from this method from math.SE used by Luis.

_=>[x=Math.random(),(x+.5)%1,--x*-2%1]

Try it online!

How?

• We pick a uniform random value \$x \in [0,1[\$
• We define \$Y=x+\dfrac{1}{2}\$
• We define \$Z=2-2x=-2(x-1)\$

Then we define:

$$y=\cases{ Y&\text{if $Y<1$}\\ Y-1&\text{if $Y\ge 1$} }$$

and:

$$z=\cases{ Z&\text{if $Z<1$}\\ Z-1&\text{if $Z\ge 1$} }$$

Which can also be written as \$y=Y-d\$ and \$z=Z-(1-d)\$ with:

$$d=\cases{ 0&\text{if $x<1/2$}\\ 1&\text{if $x\ge1/2$} }$$

Leading to:

$$x+y+z=x+\left(x+\dfrac{1}{2}-d\right)+(2-2x-(1-d))=\dfrac{3}{2}$$

Answer 5

R, 38 27 bytes

-10 bytes thanks to Nick Kennedy's suggestion of following Bubbler's trick; -1 extra byte by switching \$X\$ to \$1-X\$.

c(x<-runif(1),.5+x,-2*x)%%1

Try it online!

The footer on TIO includes a Kolmogorov-Smirnov test of uniformity for the three marginals.

Equivalent to the easier to read:

R, 38 bytes

c(x<-runif(1),y<-1-2*x+(x>.5),1.5-x-y)

Try it online!

This uses the same method as other answers, rediscovered independently.

Defines \$X\sim U(0,1)\quad Y = \begin{cases}1-2X\text{ if }X<\frac12\\2-2X\text{ if } X\geq\frac12\end{cases}\quad Z=\frac32 - X - Y\$

This way:

• \$X\sim U(0,1)\$ by definition
• \$Y|X<\frac12\$ and \$Y|X\geq\frac12\$ are both \$U(0,1)\$, hence \$Y\sim U(0,1)\$
• \$X+Y|X<\frac12 = 1-X|X<\frac12\sim U(\frac12, 1)\$ and \$X+Y|X\geq\frac12 = 2-X|X\geq\frac12\sim U(1, \frac32)\$ hence \$X+Y\sim U(\frac12, \frac32)\$ and therefore \$Z\sim U(0,1)\$

---

Alternative method: R, 84 bytes

This is much longer, but uses a very different method so I am also posting it (with thanks to Amic Frouvelle).

`?`=runif
x=?1
y=?1
if(abs(x+y-1)>.5|(y>.5)!=(x>.5)&abs(x-y)<.5)y=1-y
c(x,y,1.5-x-y)

Try it online!

Samples \$(X,Y)\$ from the uniform distribution over this shape:

then lets \$Z=\frac32-X-Y\$. These histograms show that the marginals are uniform:

Answer 6

C++ (clang), ^{206 ... 163} 160 bytes

#import<random>
std::vector<double>f(){std::random_device r;std::mt19937 g(r());auto x=std::generate_canonical<double,53>(g),y=x+(x<.5)-.5;return{x,y,1.5-x-y};}

Try it online!

Uses the C++ library's Mersenne twister engine to generate uniform random numbers with \$53\$ bits of randomness in \$[0,1)\$.
I realise we only need to be accurate within \$3\$ decimal digits, but that's going to be at least \$10\$ bits so might as well go the whole hog with \$53\$ (the number of bits in the mantissa part of a double).
Also uses @Bubbler's method to calculate \$(x, y, z)\$.

Ungolfed:

#include <random>
#include <vector>

std::vector<double> f() {
    std::random_device rand_dev;
    std::mt19937 generator(rand_dev());
    double x = std::generate_canonical<double, numeric_limits<double>::digits>(generator);
    double y = x;
    if (y < 1./2) {
        y += 1./2;
    } else {
        y -= 1./2;
    }
    double z = 3./2 - x - y;
    return {x, y, z};
}

Answer 7

Wolfram Language (Mathematica), 79 bytes

Sqrt[1-r[]^2]{Sin@#,Cos@#}.{{1,1,-2},Sqrt@3{1,-1,0}}/4&[2Pi(r=RandomReal)[]]+.5

Try it online!

A bit longer but using the direct Sqrt[1-r^2] distribution on the 2D disc, then rotating and translating the disc into 3D.

Answer 8

05AB1E, 14 17 bytes

₄L¨₄/Ωx¦1αsD.5+¦)

Port of @Bubbler's approach, so make sure to upvote him!

Try it online.
This approach uses random integers in the range \$(0,1)\$.

The footer of the TIOs (=O,) is to verify the result equals 1.5, which will = (print without popping); O (sum); , (pop and print).

Explanation:

₄L            # Push a list in the range [1,1000]
  ¨           # Remove the last value to make the range [1,1000)
   ₄/         # Divide each value by 1000 to make the range [0.001,1) in increments of 0.001
     Ω        # Pop and push a random value of this list
              # (NOTE: 05AB1E doesn't have a random builtin within the range [0,1))
      x       # Double this value without popping
       ¦      # Remove the first digit (the integer portion)
        1α    # And take the absolute difference with 1
      s       # Swap to take the original random value again
       D      # Duplicate it
        .5+   # Add 0.5 to it
           ¦  # And remove the first digit (the integer portion)
            ) # Then wrap all (three) values on the stack into a list
              # (after which this list is output implicitly as result)

---

₄L₄/3ãʒO3;α_}Ω

Try it online (way too slow for TIO, so will use \$10^{-1}\$ accuracy instead of \$10^{-3}\$ by replacing the ₄s (1000) with Ts (10)).
This approach uses random integers in the range \$(0,1]\$.

This previous 14 bytes answer above unfortunately didn't had marginal distribution. The \$x,y\$ values had the same distribution, but the \$z\$ values were always slightly larger in terms of marginal distribution. Here for example the distribution for the \$x,y,z\$ values with steps of \$0.1\$. These total values will be larger for the actual steps of \$0.001\$ of course, but the \$x,y\$ values will always be the same, and \$z\$ unfortunately always slightly larger.

₄L            # Push a list in the range [1,1000]
  ₄/          # Divide each value by 1000 to make the range [0.001,1] in increments of 0.001
    3ã        # Cartesian product of 3, to create all possible triplets of these values
      ʒ       # Filter these triplets by:
       O      #  Sum the three values
        3;α   #  Take the absolute difference with 3/2
           _  #  And check if it's equal to 0
      }Ω      # After the filter: pop and push a random triplet from this filtered list
              # (which is output implicitly as result)

Answer 9

Ruby, 25 bytes

p x=rand,(x+0.5)%1,-2*x%1

Try it online!

Based on the same method as most other answers.

Answer 10

Jelly, 12 bytes

ȷX÷¤._;Ḥ;N%1

Try it online!

A niladic link that returns a list of 3 floats, accurate to within \$\pm 10^{-3}\$, that is 3 decimal digits, and summing to 1.5. Footer demonstrates the sum.

Uses the same method developed by Bubbler.

Answer 11

Perl 5, 50 bytes

$x=rand;$f=$x<.5;say for$x,$x-.5+$f,3/2-2*$x+.5-$f

Try it online!

This answer takes seriously the first word of the question: output (as a verb). Many other answers here moves the output part outside of the byte count :-)