Probability of beginning with 3 lands on initial hand in Magic: the Gathering I am a Magic: the Gathering player and I want to understand the logic and the statistics behind the game.
Let me explain the basics:
A deck contains 60 cards and there is about a third part of those cards need to be land type. (This proportion varies with the type of deck.)
My first question is how to calculate the probability of a starting hand that contains Y (by default 7) cards to have X lands on it (1, 2, 3, etc...).
I assume that after that, with 1 draw each turn the probability of drawing a land type is total_lands - lands in hand divided by the numbers of cards left as it is an independent event, but I might be wrong.
EDIT: Thanks to approbed solution I did the computations and they seem to be well:
from scipy.special import comb
import numpy as np


def compute(desired_lands_in_starting_hand: int, total_land_cards: int, hand_size: int, deck_size: int = 60):
    print(f" Composition of experiment ---> \nYou want {desired_lands_in_starting_hand} lands in {hand_size} cards starting hand")
    print(f"With a total of {total_land_cards} land cards in your {deck_size} cards deck \n")
    k, N, K, n = desired_lands_in_starting_hand, deck_size, total_land_cards, hand_size 
    probability = (comb(K, k, exact=True) * comb(N-K, n-k, exact=True)) / comb(N, n, exact=True)
    print(f" Resulting probability is ---> {np.round(probability, 4) * 100} %")

compute(3, 23, 7, 60)

 A: I remember dabbling with a similar question in the past!
You are looking for the probability of success of WITHOUT replacement: the hyper-geometric distribution Wikipedia link.
You can obtain the probability distribution of getting exactly $k$ successes over $n$ trials (size of your hand) in a population $N$ (your deck size) with $K$ positive outcomes (number of lands). You can do this with suitable programming languages (R, python, Matlab, Stata, Julia...)
If you are not familiar with programming, I think you can find some hypergeometric calculators online as well. I don't want to give you a solution, as I think this is a funny problem to get your hands dirty!
Here is a nice article to dive even deeper: link
A: EDIT: Seeing the other answer and plugging in the numbers, my result is completely wrong. I'll leave this answer here just so that someone may possibly explain where I went wrong!
I'm nowhere near a pro statistician, and there are probably better explanations and smarter ways of doing this, but this is my attempt.
Let's say one third of your cards are lands, so out of 60 cards you have 20 lands. In an initial hand of 7 cards, there are 35 combinations of the positions in which the 3 lands can be found:

For example, you could draw 3 lands in the first three cards, in which case the probability would be:

The first term tells that we have 20 available lands out of 60 cards. Once we draw the first card and it is a land, then for the second card we have 19 lands available and 59 cards in the deck, and so on. This case is the least probable case because of the amount of other cards left in the deck. By analogy we can get the case with the largest probability, which is the case where the last three cards are lands:

And so there are 33 more combinations which give different probabilities. Now, I tried this in excel and saw that the probabilities are log-normally distributed. Here someone with a better understanding can give a theoretical explanation why it is so. Anyhow, this property tells us that we can get the median of all the 35 values by taking the geometric mean:

If the probabilities are normally distributed in the logarithmic space, then we can find the mean by taking the average of the logarithm of the minimum and maximum probabilities, or equivalently the geometric mean of the log-normally distributed values:

This is practically the same value as the calculated median from all the 35 values I got in excel, with the difference starting at the 6th decimal place. However there is a difference in the 5th decimal place between calculating the median using excel median() function and taking the geomean() of the values. Now I'm guessing, but this may be caused by the difference between discrete and continuous distributions.
I hope this makes sense and that I didn't break too much statistical rules with my assumptions! Someone correct me if I'm wrong.
