There is nothing original about this question. It was asked here.
I am just curious about an answer that is beyond my mathematical level.
In one of the simulations appearing in the comments to the OP, the number $e$ is approximated through a formula like this:
$$\Large2 + \mathrm E\left(\frac{1}{\left(\lceil \frac{1}{X\, \sim \,U(0,1)}\rceil-2\right)!} \right)$$
Hopefully this is acceptable mathematical notation for the code (2 + mean(1/factorial(ceiling(1/runif(1e5))-2))).
The question is, what mathematical concept, geometrical approximation, or distribution approximation underpins this formula?
Note that the ceiling changes when its argument passes through an integer. This happens when the uniform variable passes through $1/n$ for some $n$. The ceiling will be $n$ when the uniform variable is between $1/n$ and $1/(n-1)$. The length of this interval is $\frac{1}{n-1}-\frac{1}{n}=\frac{1}{n(n-1)}$, and the probability of a uniform random variable on (0,1) lying in some subinterval of (0,1) is the length of that interval.
So the overall random variable inside the expectation there is equal to $\frac{1}{(n-2)!}$ with probability $\frac{1}{n(n-1)}$, for $n=2,3,\dots$. So the formula for the expectation is
$$\sum_{n=2}^\infty \frac{1}{n(n-1)(n-2)!}=\sum_{n=2}^\infty \frac{1}{n!}.$$
This is the famous Maclaurin series for $e$, except missing the first two terms, which add up to $2$.
Incidentally, this is a pretty terrible way to approximate $e$; a few runs on my machine with 100,000 simulations each gave an accuracy of about $10^{-3}$. Given how Monte Carlo's accuracy scales (the error is essentially proportional to $n^{-1/2}$), this means you should expect to need on the order of $10^{16}$ simulations to get an accuracy of $10^{-9}$. Yet only 13 terms of the sum itself are required for this accuracy.