Drawing random lines in a cylinder - How does the distribution of angles look like?

Question

We have a 3D cylinder with a radius of $r = 1$ and a height of $h = 2$.

Now, we pick two random points on the edges of the cylinder - one for the top and one for the bottom circle. By connecting the two dots with a line what we can do next is calculate the angle between the line and the vertical z-axis.

We repeat this procedure $N$-times and collect all the random angles.

What can we say for sure?

• The minimum angle should be $0$°, the case where we pick two points laying on top of each other giving us a parallel line to the z-axis.
• The maximum angle should be $45$° - which is decided for us by $r$ and $h$.

---

My question here is what distribution of angles we should get. Running an experiment will give me the following distribution:

Which is a bit counter intuitive - at least to me - I expected a uniform distribution to be honest. So, I do hope I did not make a mistake in my code, which is why I provide the experiment code below.

In any case I would be interested in whether this distribution makes sense and if it's possible to explain it.

One thing that surely is important: np.random.rand() draws uniform values $x \in [0,1]$ so the points on the edges (v and w) should be uniform too.

%matplotlib inline  # Jupyter Notebooks
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.cm as cm
import seaborn as sns
import numpy as np
import pandas as pd

def angle(a, b):    
    a_u = a / np.linalg.norm(a)
    b_u = b / np.linalg.norm(b)    
    c = np.clip(np.dot(a_u, b_u), -1, 1)    
    return np.degrees(np.arccos(c))

# Create the random points on the cylinder's edges
N = 5000
a = 2 * np.pi * np.random.rand(N) 
b = 2 * np.pi * np.random.rand(N)
v = np.vstack([np.cos(a), np.sin(a), 2*np.ones(N)]).T
w = np.vstack([np.cos(b), np.sin(b), np.zeros(N)]).T
df_v = pd.DataFrame(v, columns=['x', 'y', 'z'])
df_w = pd.DataFrame(w, columns=['x', 'y', 'z'])
# Get the connecting lines
df_u = df_v - df_w

# Calculate the angles of the lines w.r.t. the z-axis
e3 = np.asarray([0, 0, 1])
ds_angles = df_u.apply(lambda x: angle(x, e3), 1)

# Plotting the distribution of angles
ax = plt.figure(figsize=(8, 4.5)).gca()
sns.distplot(df_u.apply(lambda x: angle(x, e3), 1), ax=ax, kde=False)
plt.xlabel('Degree')
plt.show()

Answer 1

Let $X$ denote the point in the circle below, and $Y$ denote the point in the circle above. For any given $X$, with probability $1/2$ we would have that $Y$ lies in the semi-circle whose midpoint is directly above $X$.

This means with probability $1/2$, the angle lies between $0$ and $\arccos(2/\sqrt6)\simeq 35.26^{\circ}$. Hence, the assumption that the distribution of the angles would be uniform is flawed, and about half the time the angles should lie between $35.26^{\circ}$ and $45^{\circ}$ (which is far from half of the possible angles).

---

Here's an attempt at calculating things more explicitly. Let $X, Y$ be uniform in the circle $C=\left\{(x,y,0)\in\mathbb R^3\,|\,x^2+y^2=1\right\}$. Your random line has direction $\big(Y+(0,0,2)\big) - X$. With $W = Y-X$, $h=(0,0,2)$ and $Z$ the cosine of the angle between that line an the $z$-axis, we have

$$Z = \frac{\langle W+h, (0,0,1)\rangle}{\lVert W+h\rVert} = \frac2{\lVert W+h\rVert}.$$

This implies

$$Z^2 = \frac{4}{4 + {\lVert W\rVert}^2} = \frac1{1+{\left\lVert\frac{W}2\right\rVert}^2}.$$

If $\theta_1,\theta_2 \sim \text{Unif}([0,2\pi])$, then $X\sim \big(\cos(\theta_1),\sin(\theta_1)\big)$ and $Y\sim \big(\cos(\theta_2),\sin(\theta_2)\big)$ so that $W\sim\big(\cos(\theta_1)-\cos(\theta_2),\sin(\theta_1)-\sin(\theta_2)\big)$. Hence,

\begin{align} {\left\lVert\frac{W}2\right\rVert}^2 &\sim \frac14\left( \underbrace{(\cos\theta_1)^2+(\cos\theta_2)^2}_1-2\cos\theta_1\cos\theta_2 + \underbrace{(\sin\theta_1)^2+(\sin\theta_2)^2}_1-2\sin\theta_1\sin\theta_2 \right) \\&= \frac14\left( 2 - 2\cos\theta_1\cos\theta_2 - 2\sin\theta_1\sin\theta_2 \right) \\&= \frac12\big( 1 - (\cos\theta_1\cos\theta_2 + \sin\theta_1\sin\theta_2) \big) \\&= \frac12\big( 1 - \cos(\theta_1-\theta_2) \big) = \frac12 - \frac{\cos(\theta_1-\theta_2)}2 . \end{align}

Therefore

$$Z^2 \sim \frac{2}{3-\cos(\theta_1-\theta_2)}.$$

Observing that $Z\geq 0$, it follows that the angle $\alpha$ is distributed as

$$\alpha \sim \arccos\left(\sqrt{\frac{2}{3-\cos(\theta_1-\theta_2)}}\right),$$

where $\theta_1,\theta_2 \sim \text{Unif}([0,2\pi])$. From here on, we can calculate $\mathbb P(\alpha \leq x)$ for $0<x<\pi/4$ as some area in the square $[0,2\pi]^2$ contained in the $\theta_1$-$\theta_2$ plane.

Alternatively, we can observe that $\Theta = \theta_1-\theta_2$ has pdf with support in $[-2\pi,2\pi]$ given by

$$f_{\Theta}(z) = \left\{\begin{array}{cccc} \frac{2\pi+z}{4\pi^2}&&&-2\pi\leq z <0\\ \frac{2\pi-z}{4\pi^2}&&&0\leq z \leq 2\pi \end{array}\right. $$

Then, some algebraic manipulations show that

$$\mathbb P(\alpha \leq x) = \mathbb P\left(\cos(\Theta) \geq 3-\frac{2}{{\cos(x)}^2}\right) .$$

Using the pdf above, we can calculate this explicity in terms of integrals that simplify to

$$\mathbb P(\alpha \leq x) = \frac{1}{\pi}\,\arccos\left(3-\frac{2}{{\cos(x)}^2}\right) .$$

The density function $f_\alpha$ of $\alpha$ has support in $[0,\pi/4]$, and according to Wolfram we have (after simplifications considering its support) that it equals

$$f_{\alpha}(z) = \frac{2}{\pi}\,\frac{1}{\cos(z)\,\sqrt{\cos(2z)}},$$

which is way simpler than I expected, to be honest. Observe that it's indeed a probability density function, in that it is non-negative and its integral is $1$.

Your histogram does somewhat resemble the area under $f_{\alpha}$, as you can verify here, so this looks more like a case of 'wrong intuition' than of 'wrong code'.

---

I must point out that this sheds light to yet another rather curious identity involving $\pi$:

$$\int_{\sqrt{2}/2}^1\,\,\,\frac{1}{u\,\sqrt{2u^2-1}\,\sqrt{1-u^2}}\,\,\,du = \frac\pi2$$ $${}$$ $$\int_{1/2}^1\,\,\,\frac{1}{u\,\sqrt{2u-1}\,\sqrt{1-u}}\,\,\,du = \pi$$

Answer 2

It's possible to find the closed form for the cdf. The angle is $$A = \arccos \frac {((\cos v, \sin v, 1) - (\cos u, \sin u, -1)) \cdot (0, 0, 1)} {| (\cos v, \sin v, 1) - (\cos u, \sin u, -1) |} = \arccos \sqrt {\frac 2 {3 - \cos(v - u)}}.$$ Switching to $\xi = v + u, \eta = v - u$, the cdf is $$\operatorname P(A < x) = \int_0^{2\pi} \int_\eta^{4\pi - \eta} \frac 1 {(2 \pi)^2} [A < x] d\xi d\eta \overset * = \\ \int_0^\pi \frac 2 {(2 \pi)^2} (4\pi - 2\pi) [A < x] d\eta = \\ \frac 1 \pi \int_0^\pi [\eta \lt \arccos(3 - 2 \sec^2 x)] d\eta = \\ \frac 1 \pi \arccos(3 - 2 \sec^2 x),$$ where $[\text {condition}]$ is $1$ if the condition is true and $0$ otherwise. For $*$, $f$ is symmetric around $\pi$, therefore $\int_0^{2\pi} \eta f(\eta) d\eta = 2 \pi \int_0^\pi f(\eta) d\eta$.