Probability of drawing perfect starting hand

Question

I have a Magic the Gathering deck where an ideal starting hand can potentially let me win on turn 3. What is the probability of drawing an ideal hand? How is this computed, and what is a general algorithm that can be used to compute ideal hand probabilities for an arbitrary deck and ideal hand?

This is really a probability question, no knowledge of Magic required. All you need to know is the deck is 60 cards with some duplicates, and the starting hand size is 7 cards.

The deck:

• 10 Forest
• 6 Swamp
• 4 Terramorphic Expanse
• 4 Llanowar Elves
• 4 Elves of Deep Shadow
• 4 Rites of Spring
• 4 Cultivate
• 4 Blood Speaker
• 4 Howling Wolf
• 4 Skyshroud Sentinel
• 4 Wall of Blossoms
• 2 Infernal Spawn of Evil
• 2 Infernal Spawn of Infernal Spawn of Evil
• 4 Infernius Spawnington III, Esq. (for reference, this is a silver-bordered card so not legal in any competitive format)

The hand:

• 1 Forest
• 1 Llanowar Elves -OR- Elves of Deep Shadow
• 1 Rites of Spring
• 1 Infernius Spawnington III, Esq.
• 3 (anything)

---

If you are a Magic player, this is how the ideal game goes:

• Turn 1: Play a Forest and Llanowar Elves.
• Turn 2: Play Rites of Spring, discard 4 cards that are not Infernius Spawnington, search for four lands including a Swamp, play the Swamp, play Infernius, say "I'm here!" to do 3 damage, attack for 9 damage.
• Turn 3: Attack again.

---

Bonus points Although my original question was just asking about the probability of drawing the ideal starting hand, as @Phillip Kendall pointed out, there are some additional details that could increase the probability of accomplishing the desired effect of a 3rd turn win.

Relevant details:

• Because the first player does not draw on their first turn, if you go second, you effectively have a starting hand of 8 cards. (My first thought to account for this, is to run the same computation with a starting hand of 7 cards or 8 cards, and average the result.)
• The Forest and Elf must be in your starting hand (or drawn on your first turn), but the Rites or Infernius could be drawn on your second turn. (I am not sure where to start to account for this.)
• If you are going second, and your starting hand is not ideal, you can mulligan in hopes of getting an ideal hand, and still have enough cards to discard 4 with Rites of Spring, and achieve the 3rd turn win. (I am not sure how to account for this.)

Answer 1

The relevant probability distribution is the multivariate hypergeometric distribution, which governs random choices without replacement from a mix of items of different types. In this case, you have five relevant types of things:

• 10 things of type A (Forests)
• 8 things of type B (mana Elves)
• 4 things of type C (Rite of Spring)
• 4 things of type D (Infernius, Esq.)
• 60-10-8-4-4 = 34 things of type E (anything else)

You need to calculate the probability of getting at least one of each of types A-D in 7 draws.

What the distribution can tell you is the probability of drawing exactly certain numbers, e.g. exactly one Forest, one elf, one Rite, one Infernius, and three cards that are none of those. That would be

C(10, 1) C(8, 1) C(4, 1) C(4, 1) C(34, 3) / C(60, 7) = 17408 / 877743

where

C(N, k) = N!/k!(N-k)!

You then have to add this up for all possible numbers of the various types of cards you could draw:

• 1 Forest, 1 elf, 1 Rite, 1 Infernius, 3 other
• 2 Forest, 1 elf, 1 Rite, 1 Infernius, 2 other
• 3 Forest, 1 elf, 1 Rite, 1 Infernius, 1 other
• 4 Forest, 1 elf, 1 Rite, 1 Infernius, 0 other
• 1 Forest, 2 elf, 1 Rite, 1 Infernius, 2 other
• ...

I wrote a bit of code to go through the possible configurations - it finds 35 of them - and add up the probability for each. The total is 42304/877743, or 0.04819634.

---

I definitely may have missed something, so for anyone who knows Python and wants to check my work, here's the code snippet:

import itertools as it
from fractions import Fraction
from math import factorial

def C(N, k):
    return Fraction(factorial(N), factorial(k) * factorial(N-k))

# Forest, Elves, Rites, Inferniuses, Other
def P(f, e, r, i, o):
    return C(10, f) * C(8, e) * C(4, r) * C(4, i) * C(34, o) / C(60, 7)

result = sum(P(*t, o=7-sum(t)) for t in it.product(range(1, 5), repeat=4) if sum(t) <= 7)

print(result)
print(float(result))

Answer 2

I've written an online program that will work out the probability by simulation. (Only it doesn't deal with mulligans.)

Your ideal game happens with a probability of around 8.7%. (6.8% as starting, 10.6% as drawing)

Simply change the header to modify input. For your example:

const config = `
DECK
10 forest
8 elf
4 rite
4 infernius
34 other

TURN 1
1 forest
1 elf

TURN 2
1 rite
1 infernius

TRIALS: 1000000
`;

DECK is a list of [qty] [card name].
TURN [N] is a list of cards needed by turn N.
TRIALS: [N] sets how many trials to do.

Then press 'Execute' to run the simulations.

Answer 3

Here is a generalized version of @David Z's approach, in F#:

open System.Collections.Generic

// --- Utility functions ---

let memoize f =
    let cache = Dictionary<_,_>()
    fun c ->
        let exist, value = cache.TryGetValue (c)
        match exist with
        | true -> value
        | _ -> 
            let value = f c
            cache.Add (c, value)
            value

let rec factorial n =
    if n < 0 then failwith ""
    match n with
    | 0 | 1 -> bigint 1
    | _ -> (bigint n) * factorial(n-1)
let factorialFast = memoize factorial
    
let rec product lists =
    match lists with
    | h::[] ->
        List.fold (fun acc elem -> [elem]::acc) [] h
    | h::t ->
        List.fold (fun cacc celem ->
            (List.fold (fun acc elem -> (elem::celem)::acc) [] h) @ cacc
            ) [] (product t)
    | _ -> []

// --- Hypergeometric distribution --- 

// C(n, k) in David Z's answer
let possibleSamples (population, sample) = 
    /* I added this special 0 case, because if you have some weird
       degenerate deck where every card is part of the ideal hand, 
       you can end up with 0 as input, and it breaks computations down stream. */
    if population = 0 || sample = 0 then 1.0 
    else
        let a = factorialFast population |> float
        let b = factorialFast sample * factorialFast (population-sample) |> float
        a / b
let possibleSamplesFast = memoize possibleSamples

// P(f,e,r,i,o) in David Z's answer, except o is computed from [f,e,r,i]
let probabilityOfHandItemCounts (itemFreqs: int list) (deckSize: int) (handSize: int) (handItemCounts: int list) =
    let itemPossibleSamples = 
        handItemCounts 
        |> List.zip itemFreqs 
        |> List.map possibleSamplesFast

    let otherPossibleSamples = possibleSamplesFast(
        deckSize - List.sum itemFreqs,
        handSize - List.sum handItemCounts)

    let possibleHands = possibleSamplesFast(deckSize, handSize)

    (otherPossibleSamples :: itemPossibleSamples) 
    |> List.fold (fun a b -> a * b) 1.0
    |> fun x -> x / possibleHands

// for..in expression in David Z's answer
let handsWithRequiredItems (itemFreqs: int list) (handSize: int) =
    itemFreqs
    |> Seq.map (fun f -> [1..f])
    |> Seq.toList
    |> product
    |> List.filter (fun xs -> List.sum xs <= handSize)

// Final summing operation
let probabilityOfHandWithRequiredItems (itemFreqs: int list) (deckSize: int) (handSize: int) =
    handsWithRequiredItems itemFreqs handSize
    |> List.map (probabilityOfHandItemCounts itemFreqs deckSize handSize)
    |> List.sum

let result = probabilityOfHandWithRequiredItems [10; 8; 4; 4] 60 7
printfn "Probability of ideal Infernius starting hand: %f" result
// Prints 0.048196, same as David Z's result

To test the generalization, here is another case: An old-school degenerate Channel/Fireball deck. One copy of each card is required for a first turn win.

• 10 Channel
• 10 Fireball
• 10 Black Lotus
• 10 Mountain

let result2 = probabilityOfHandWithRequiredItems [10; 10; 10; 10] 40 7
printfn "Probability of ideal Channel/Fireball starting hand: %f" result2
// Prints 0.689246, which feels about right intuitively