rep(1:13, 4)  # Deck of cards

sample(rep(1:13, 4), size = 52) # Shuffled once

rle(sample(rep(1:13, 4), size = 52))$lengths # Patterns

rle(sample(rep(1:13, 4), size = 52))$lengths > 1  # Pairs

pairs <- sapply(1:1e6,
                function(x)
                    sum(rle(sample(rep(1:13, 4),
                                   size = 52))$lengths > 1))

table(pairs)

mean(pairs)

p.hat <- mean(pairs > 0)

sqrt(p.hat * (1 - p.hat) / 1e6)