Collisions during random insertion unto given domain Here's a satisfying one. Assume a list of maximum length L which starts out empty. For the first element we pick a random number within the domain [0 ... L), we check if the list already contains it (it obviously doesn't when the list is empty), and if it's not already present then we populate the first empty element with that random number. If the number is present in the list (which will happen during later iterations) then we count a collision, we choose another random number, we check it again, and we try to populate the same position in the list. We iterate until the list is complete, counting collisions along the way.
Counting the number of collisions for each element in the list, and graphing that against the list's occupation level, we end up with the following graph:

The x axis represents the list's occupation degree, and the y axis represents the average number of collisions for each attempt to choose a collision-free number. Note that the y axis is logarithmic – an exponential trend line would be represented as a straight line.
I find a lot of satisfying things in this graph: the shape of the curve, the symmetry of the y range (10E-3 ... 10E+3), and the fact that the sweet spot where there's one collision per number selection is at roughly 33%.
Having said that, everything here is based on my empirical attempts to work out how this works – I averaged multiple iterations of computer-generated simulations to get to this graph. I have tested it with both [0...1000) and [1...10000), and the results were basically identical.
I would appreciate a confirmation that my empirical results can be validated using a formal mathematical approach, rather than my naïve computer simulations.
 A: Your problem is a manifestation of the classical occupancy distribution: This question relates to the classical occupancy distribution in statistical theory (for details on the properties and computation of this distribution, see e.g. O'Neill 2020).  Taking $n \geqslant 1$ random draws with replacement from a list with $L$ distinct elements, we let $1 \leqslant K \leqslant \min(n,L)$ be the number of distinct elements on the list (so there have been $n-K$ "collisions").  The random variable $K$ has mass function:
$$\mathbb{P}(K=k | n,L) = \frac{(L)_k \cdot S(n,k)}{L^n} \quad \quad \quad \text{for all } 1 \leqslant k \leqslant \min(n,L),$$
where $(m)_k = \prod_{i=0}^{k-1} (m-i)$ are falling factorials and $S(n,k)$ are the Stirling numbers of the second kind.  The expected value of this random variable is:
$$\mathbb{E}(K) = L \Big( 1 - \frac{1}{L} \Big)^n.$$
The "occupation degree" in your analysis is the value $P \equiv K/L$ and the "number of collisions" is the value $C \equiv n-K$, both of which are random variables that are functions of $K$.  However, for any value of $K$, these quantities should be related deterministically by $C = n-PL$, so it is not very useful to graph them against each other.

Expected "occupation degree" for the list: A simpler way to view the analysis is to ignore the number of collisions and use the fact that the classical occupancy distribution leads to:
$$\log \mathbb{E}(P) = n \log \Big( 1 - \frac{1}{L} \Big).$$
From this equation, we see that if you try plotting the log-occupancy-degree in your simulations against the number of total selections $n$, you should observe a roughly linear relationship between these quantities (subject to some random variation).  Another useful thing to note here is that, when $n$ and $L$ are large, the occupancy distribution is closely approximated by the normal distribution.
