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:

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.)

3 Answers 3


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


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)

  • 2
    Mathematically, this is correct for "what is the chance of drawing those 7 cards" but it's not the right answer to "what is the chance this deck kills on Turn 3" because it would be possible to draw Spawnington (or other crucial cards) on Turn 2. Jul 8, 2022 at 11:02
  • To account for the 2nd turn draw, couldn't you just swap the 7 for an 8? Also, if you go 2nd, then there is a 9th draw too.
    – JamesFaix
    Jul 8, 2022 at 15:13
  • 2
    @JamesFaix No, because you must have a Forest and an Llanowar Elves in hand on Turn 1. With regards to play/draw, you should probably clarify your question as to the exact problem you're trying to solve - do you actually care about "the perfect 7 card draw", "the chance that your can do 20 damage by the end of your turn 3" or something else? Jul 8, 2022 at 15:35
  • Where is the number 4 coming from in range(1, 5) and repeat=4? Is that because there are 4 classes of items we care about?
    – JamesFaix
    Jul 8, 2022 at 16:26
  • Something like that, yeah. Specifically, it's because there are 4 classes of items that the sample (your opening hand) can contain 1, 2, 3, or 4 copies of in order to be "ideal". I took advantage of a little coincidence there, the fact that no ideal hand can contain more than 4 copies of a single card; without taking that into account, the call would have to be it.product(range(1, 11), range(1,9), range(1,5), range(1,5)).
    – David Z
    Jul 8, 2022 at 18:00

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 = `
10 forest
8 elf
4 rite
4 infernius
34 other

1 forest
1 elf

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.

  • Very cool! I hadn't even thought of an empirical approach. I would prefer a theoretical explanation though.
    – JamesFaix
    Jul 11, 2022 at 19:47
  • Thank you! Calculating the actual probability is definitely more useful, but is often a lot more difficult to do when you have lots of variables. Jul 11, 2022 at 21:44

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)

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 
        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 = 
        |> 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) =
    |> 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.

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

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .