Euler's number, $e$, is a fundamental mathematical constant approximately equal to 2.71828, and it is the base of the natural logarithm. This section introduces a unique probabilistic method for estimating $e$ using a sequence of random numbers. We explore the conceptual underpinnings of this method, provide a step-by-step derivation, and present implementation codes in MATLAB and Python.
Consider a sequence of random numbers, $U_1, U_2, \ldots$, where each number is independently and uniformly drawn from the interval $[0, 1]$. We define $N$ as the position in the sequence where a number is first greater than its immediate predecessor. Formally,
$$ ⁍ $$
The probability $P(N>n)$ represents the probability that the first n numbers are in non-increasing order (i.e., each number is less than or equal to its predecessor). This is because for $N$ to be greater than $n$, the condition $U_n>U_{n−1}$ must not be satisfied for any of the first $n$ comparisons.
Given $n$ numbers, there are $n!$ (factorial of n) ways to order them. Since we assume the sequence of numbers is random, every ordering is equally likely. The specific ordering where $U_1⩾U_1⩾⋯⩾U_n$ is just one of these $n!$ possibilities. Therefore, the probability of this specific ordering occurring is $1/n!$ . Thus,
$$ ⁍ $$
$P(N = n)$, the probability that the first occurrence of a number being greater than its immediate predecessor is exactly at position $n$, is derived as:
$$ ⁍ $$
The expected position $E(N)$ can then be calculated by summing over all $n$, weighted by their probabilities:
$$ ⁍ $$
where the equality is Maclaurin series expansion of $e^x$ evaluated at $x=1$.
Thus, we can create many sequences of U(0,1) r.v.s and in each sequence find the index of the first value that is greater from the previous values, then take the sample mean of all these indices. The following python and MATLAB code below are doing just that.
% Probabilistic method for estimating e in MATLAB
N = 1e6; % Number of trials
counts = zeros(1, N);
for i = 1:N
prev = rand; % Initial random number
n = 2; % Start from the second number
while true
current = rand; % Generate next random number
if current > prev
counts(i) = n;
break;
end
prev = current;
n = n + 1;
end
end
estimated_e = mean(counts);
fprintf('Estimated e: %f\\n', estimated_e);