e in Python function

Here is a function written in Python. It is used to randomly choose at least one from a list.

import random
import math
def choose_at_least_one(list_items):
    chosen_items = [i for i in list_items if random.random() <= (1/len(list_items))]
    if len(chosen_items) == 0:
        return choose_at_least_one(list_items)
    else:
        return chosen_items

The process is very simple. When the function with a list of length $n$ is called, it goes through all the items in the list and the item is chosen if the random number generated by program is less than $1 \over n $. After that if nothing is chosen, then call the function again. Or else return the list with chosen items of length $r$.

Here is a question: When $n \to \infty$, which value will length $r$ approach?

Let the input list has length $n$. Every time this function is called, the returned list has a length of $l$ ($l \leq n$). Suppose the function is called for $c$ times and finally the returned list has a length of $r$ ($r \leq n$).

Then there are many possibilities:

$c=1$, that is to say, the first time the function is called, it returns a list with length $r$.

$c=2$, that is to say, the first time the function is called, it returns a list with length $0$. The second time the function returns a list with length $r$.

$c=3$, that is to say, the first two time the function is called, it returns a list with length $0$. The second time the function returns a list with length $r$.

...

Suppose for the first $c-1$ times, the function returns lists of length 0, and at the $c$th time it returns a list of length $X_c$. The final list has a length of $Y$. Let the length of returned list of the function be $R$.

$P(X_{c}=r)=(P(R=0))^{c-1} \cdot P(R=r)$

$P(R=r) = ({1\over n})^r(1-{1\over n})^{n-r} {n \choose r}$

$P(R=0) = (1-{1\over n})^n$

$P(X_{c}=r)=(1-{1\over n})^{n(c-1)} \cdot ({1\over n})^r(1-{1\over n})^{n-r} {n \choose r}$

$P(Y=r)=\sum_{c=1}^{\infty}{P(X_c=r)}$

$P(Y=r)=({1\over n})^r (1-{1\over n})^{n-r} {n\choose r} {1\over{1-(1- {1\over n})^n}}$

$E(Y) = \sum_{y=1}^{n}{y \cdot P(Y=y)} \\ ={1\over{1-({1-{1\over n})^n}}}\sum_{y=1}^{n}{y \cdot ({1\over n})^y (1-{1\over n})^{n-y}}{n\choose y}\\ ={1\over{({{n}\over{n-1}})^n}-1}\sum_{y=1}^{n}{(n-1)^{-y}y {n \choose y}}$

$\lim_{n \to \infty} {{1}\over {({n\over {n-1}})^n}-1} \\ = \lim_{n \to \infty} {{1}\over {{1\over ({1-{1\over n}})^n}}-1}\\ ={{1}\over {{1\over \lim_{n \to \infty} ({1-{1\over n}})^n}}-1}\\ ={{1}\over {{1\over {e^{-1}}}-1}} = {1\over{e-1}}$

According to binomial theorem, $(x+y)^n=\sum_{k=0}^{n}{{n \choose k}{x^k} y^{n-k}}$

$\frac{d ((x+y)^n)}{dx}=n(x+y)^{n-1}$

$\frac{d (\sum_{k=0}^{n} {n\choose k}x^k y^{n-k})}{dx}=\sum_{k=1}^{n} k {n \choose k} x^{k-1} y^{n-k}$

$n(x+y)^{n-1}=\sum_{k=1}^{n} k {n \choose k} x^{k-1} y^{n-k}$

Let $x={1\over{n-1}}, y=1$

${n\over{n-1}}({1\over{n-1}}+1)^{n-1}=\sum_{k=1}^{n} {k{n \choose k}({1\over{n-1}})^{k}}$

$\lim_{n \to \infty } ({n\over{n-1}}({1\over{n-1}}+1)^{n-1})=\lim_{n \to \infty } (\sum_{k=1}^{n} {k{n \choose k}({1\over{n-1}})^{k}})$

$\lim_{n \to \infty } ({n\over{n-1}}\cdot (1+{1\over{n-1}})^{n-1})         \\=\lim_{n \to \infty } ({n\over{n-1}})\cdot \lim_{n \to \infty }(1+{1\over{n-1}})^{n-1}\\ =1\cdot\lim_{n \to \infty }  (1+{1\over{n-1}})^{n-1}\\ =\lim_{n \to \infty } (1+{1\over{n-1}})^{n-1}$

Let $t=n-1$. When $n \to \infty$, $t \to \infty$

$\lim_{n \to \infty } (1+{1\over{n-1}})^{n-1}=\lim_{t \to \infty } (1+{1\over{t}})^{t}=e$

$\lim_{n \to \infty } ({n\over{n-1}}({1\over{n-1}}+1)^{n-1})=e$

Therefore,

$\lim_{n \to \infty } ({\sum_{k=0}^{n} {{k {n\choose k}}\over{(n-1)^k}}})=e$

$\lim_{n \to \infty} (E(Y))={{e}\over{e-1}}$

We can test the result with Python:

import random
import math
from math import factorial
import sys
def choose_at_least_one(list_items):
    chosen_items = [i for i in list_items if random.random() <= (1/len(list_items))]
    if len(chosen_items) == 0:
        return choose_at_least_one(list_items)
    else:
        return chosen_items
def test(a,b):
    r = 0
    for i in range(a):
        r+=len(choose_at_least_one(list(range(b))))
    return (math.e / (math.e - 1)) - (r/a)

In [4]: math.e / (math.e - 1)
Out[4]: 1.5819767068693265

In [6]: test(999, 100000)
Out[6]: -0.014619889727270241

In [7]: test(9999, 10000)
Out[7]: -0.005482039005261008