import numpy as np

# Compute the probability by using the alternating sum;
# Unstable, do not use it
def prob_exact_unstable(n, m):
    P = 0
    for k in range(m):
        if k == 0:
            A = 1
            B = 1
        else:
            A = -A
            B = B * (m-k+1) / k
        C = np.pow(1-k/m, n)
        P += A * B * C
    return P

# Compute the probability by using dynamic programming; stable
def pre_compute_prob_table(N, M):
    P = np.zeros((N+1, M+1))
    P[0][0] = 1
    for n in range(0, N):
        for m in range(1, min(M+1, n+2)):
            P[n+1][m] = P[n][m] + np.pow(1-1/m,n) * P[n][m-1]
    return P

def prob_exact(n, m, P):
    return P[n][m]

# Compute the approximate probability
def prob_approx(n, m):
    return np.pow( 1 - np.exp(-n/m), m )

# Compute the probability with fixed m and a list of n between 0 and N
m = 480
N = 6000
P = pre_compute_prob_table(N, m)
n_list = list(range(N))
prob_exact_list, prob_approx_list = [], []
for n in n_list:
    prob_exact_list.append( prob_exact(n, m, P) )
    prob_approx_list.append( prob_approx(n, m) )

# Plot. The exact formula and approximate formula give very close results
import matplotlib.pyplot as plt
plt.plot(n_list, prob_exact_list, label='prob_exact')
plt.plot(n_list, prob_approx_list, label='prob_approx')
plt.legend()
plt.xlabel('n')
plt.title(f'probability of seeing each of the {m} coupons at least once, as n varies')

# Compare prob_exact_unstable() and prob_exact()
n = 1000
print(f'Unstable: prob_exact_unstable({n}, {m}) = {prob_exact_unstable(n, m)}')
print(f'Stable: prob_exact({n}, {m}) = {prob_exact(n, m, P)}')

Unstable: prob_exact_unstable(1000, 480) = 31960815.770431846
Stable: prob_exact(1000, 480) = 2.8147320313908906e-34