I have come with the following analytical result:
$$ \frac{(52+51+50+49)\times 4!}{C_5^{52}}\approx 0.0023 $$
However a monte carlo simulation in python gives me another answer (that I find very low):
import numpy as np
deck = [i for i in range(52)]
iter = 1000000
victory = 0
for _ in range(iter):
hand = np.random.choice(deck, 5, replace=False)
if np.isin([1, 2, 3, 4], hand).sum() == 4:
victory += 1
print(victory / iter)
1.3e-05
Does someone know where is my mistake?
It's actually $$p=\frac{48}{C_{52}^5}\approx 1.8469e-05$$ because you choose all the aces, and you're left with $48$ choices for the last card. In other words, you can choose 4 aces and a card with 48 different ways.
Note: When I run your simulation, I got $1.9e-5$.
Here, we can have 4 aces and one of any other cards for success. One case is when 4 cards are ace and the last one is not i.e AAAAX. This can happen in 432148 ways; the first card can be any of the 4 aces, the second can be any of the remaining 3 aces, the third card can be either of the remaining 2 aces, the fourth card is the last remaining ace and the fifth card can be any of the 48 non-ace or remaining cards. But the non-ace card can be drawn in any of the five trials (XAAAA, AXAAA, AAXAA, AAAXA), so the total number of ways for success is 5*(432148) or 5!48. Similarly, we can draw 5 cards from 52 cards in 52515049*48 or (52!)/(47!) ways. So the required probability will be 5!\cdot\frac{48!}{52!}=1.8468926×10−5.
