For exactly two more aces and exactly two kings in 24 draws from a
deck of 51 cards (missing an ace), I also get
$$ \frac{{3 \choose 2}{4 \choose 2}{44 \choose 20}}{{51 \choose 24}} = 0.138,$$
(to three places), computed in R as:
18*choose(44, 20)/choose(51,24)
## 0.1380654
Here is a simulation of a million such draws with probabilities
correct to 2 or 3 places.
m = 10^6; nr.ace = nr.kng = numeric(m)
deck = 2:52 # aces are 1,2,3,4; kings 5,6.7,8
for (i in 1:m) {
draw = sample(deck, 24)
nr.ace[i] = sum(match(2:4, draw, nomatch=0)>0)
nr.kng[i] = sum(match(5:8, draw, nomatch=0)>0) }
mean(nr.ace == 2)
## 0.357932
mean(nr.kng == 2)
## 0.387775
mean(nr.ace==2 & nr.kng==2)
## 0.137598 # approx 0.138 as in exact combinatorial result
mean(nr.ace==1 & nr.kng==3)
## 0.092077
AK = as.data.frame(cbind(nr.ace,nr.kng))
table(AK)/m
nr.kng
## nr.ace 0 1 2 3 4
## 0 0.007532 0.035196 0.055101 0.035380 0.007596
## 1 0.026295 0.109556 0.157987 0.092077 0.018387
## 2 0.027684 0.105594 0.137598 0.073657 0.013399
## 3 0.008706 0.030502 0.037089 0.017775 0.002889
The approximate probability of $P(A = 2, K = 2)$ is found in
cell $(2,2)$ of the table. The approximate probability $P(A = 1, K = 3)$ is in cell $(1, 3).$
Related probabilities can also be approximated from the table.
For example, the total probability
$P(A = 2) \approx 0.358$ is found separately in the printout
above and as the total of row 2 of the table.
sum(c(0.027684, 0.105594, 0.137598, 0.073657, 0.013399))
## 0.357932
Its exact
probability (to three places) is
$$ \frac{{3 \choose 2}{48 \choose 23}}{{51 \choose 24}} = 0.358,$$
3*choose(48, 22)/choose(51, 24)
## 0.3578391
However, to get the probability of either 'two aces and two kings' OR
'one ace and three kings', you should add only two entries in
the table $(2,2)$ and $(1,3).$
Addendum: 'Expanded' R code, demonstrating method of counting aces (cards 2 through 4) in 24 draws:
draw = sample(deck, 24); draw # list of 24 cards drawn w/o replacement
## 17 40 46 51 44 30 50 24 4 7 52 31 13 28 25 18 22 42 3 5 43 48 19 8
ace.posn = match(2:4, draw, nomatch=0); ace.posn
## 0 19 9 # card 2 not drawn, 19th was card 3, 9th was card 4
tf.aces = (ace.posn > 0); tf.aces
## FALSE TRUE TRUE # TRUE for each ace found
nr.aces = sum(tf.aces); nr.aces # counts TRUE's
## 2 # count of aces in 'draw' is stored in vector 'nr.ace'
The count of the number of aces in a 'draw' does not depend on which aces or their order. This is done a million times.