Randomized Marking of a Tetrahedron

Nicolás Guarín-Zapata

2020-06-05 22:58

Also available in:

Yesterday (June 4, 2020), Christian Howard posted on Twitter the following question

You are given a tetrahedron τ. For each triangular facet of τ, we uniformly at random mark one of their edges. What is the probability that there exists an edge of τ that is marked twice?

I thought about a little bit but I couldn't find how to count properly. Then, a number popped up in my mind out of the blue: \(2/3\), but I don't know why. So, I decided to run a simulation to check this number.

The right answer is \(51/81\) that is approximately 63%. This calculation is well explained in Christian's blog and with some cool drawings (and memes).

The algorithm

The algorithm is quite simple. I number the edges in each face following a counter-clockwise orientation:

face 0: edge 0, edge 1, edge 2
face 1: edge 0, edge 3, edge 4
face 2: edge 1, edge 5, edge 3
face 3: edge 2, edge 4, edge 5

Then, I take each face and pick a random number from \((0, 1, 2)\) and mark the corresponding edge, and move to the following face. I repeat this process several times and I count the favorable cases and divide them by the number of trials to get an estimate of the probability.

Following is a Python code that follows this idea.

import numpy as np
import matplotlib.pyplot as plt

faces = np.array([
        [0, 1, 2],
        [0, 3, 4],
        [1, 5, 3],
        [2, 4, 5]])


def mark_edges():
    marked_edges = np.zeros((6), dtype=int)
    for face in faces:
        num = np.random.randint(0, 3)
        edge = face[num]
        marked_edges[edge] += 1
    return marked_edges


def comp_probs(N_min=1, N_max=5, ntrials=100):
    prob = []
    N_vals = np.logspace(N_min, N_max, ntrials, dtype=int)
    for N in N_vals:
        cont_marked = 0
        for cont in range(N):
            marked = mark_edges()
            if 2 in marked:
                cont_marked += 1
        prob.append(cont_marked/N)
    return N_vals, prob


#%% Computation
N_min = 1
N_max = 5
ntrials = 100
np.random.seed(seed=2)
N_vals, prob = comp_probs(N_min, N_max, ntrials)

#%% Plotting
plt.figure(figsize=(4, 3))
plt.hlines(0.63, 10**N_min, 10**N_max, color="#3f3f3f")
plt.semilogx(N_vals, np.array(prob), "o", alpha=0.5)
plt.xlabel("Number of trials")
plt.ylabel("Estimated probability")
plt.savefig("prob_tet.svg", dpi=300, bbox_inches="tight")
plt.show()

And we can see the following evolution for different number of trials.

Estimated probabilities for different sample sizes.

Nicolás' blog

The algorithm

Comments