Suppose there're two players $A$, $B$ playing a dice game, $A$ has a normal dice whose faces are numbered 1 to 6, $B$, on the other hand, has a (regular) icosahedron dice with faces numbered 1 to 20. Each round, $A$ rolls its dice 3 times and get the sum as his score, whereas $B$ rolls only once and get the result as his score, and the one with the greater score wins and gets 1 point. If scores are equal, then they both get 0.5 points. Is this a fair game?
A quick calculation shows that the expected scores for both players are 3.5, so it looks like it's perfectly fair. But to my surprise, when I was doing a Monte Carlo in Python, the programmed told me, out of 10^5 simulations, $A$ has only around 48.x% of the total 10^5 points, and such a 1.x% disadvantage was consistent across several tests.
Now move on to a more complicated version. This time let's introduce a new player $C$, who plays by exactly the same rule as $B$ does. Again, each round the player with the highest score wins; if two player has equal scores that are higher than the third player's, then they both get 1/2 points; and if all three have equal scores then they all get 1/3 points. Is this a fair game?
I thought with a firm "yes" for the same reasons as previously. But the simulation results totally astounded me, much more than the previous 1.x% discrepancy: now in the 3 player game, out of 10^5 simulations, the player $A$ only gets about 27.x% of the points, whereas $B$ and $C$ both get about 36.x% - $A$ has an almost 9% disadvantage!
I know to get everything clear it's best to explicitly compute the probability distribution case by case, and admittedly it's not hard in any technical sense. But, before doing that, is there any easy explanation about the discrepancies between intuition and reality as shown above? Thanks!
Attached below is the code:
import random import numpy as np import pandas as pd def get_A_score(): return sum(random.randint(1, 6) for i in range(3)) def get_B_score(): return random.randint(1, 20) if __name__ == '__main__': n_sim = 100000 A_wins = 0.0 for i in range(n_sim): if get_A_score() > get_B_score(): A_wins += 1 elif get_A_score() == get_B_score(): A_wins += 1/2 # in case of equal numbers print(A_wins / n_sim) print() wins = np.array([0.0, 0.0, 0.0]) for i in range(n_sim): A_score = get_A_score() B_score = get_B_score() C_score = get_B_score() scores = np.array([A_score, B_score, C_score]) winner_index = np.argwhere(scores == np.amax(scores)) n_winners = len(winner_index) wins[winner_index] += 1 / n_winners print(wins / n_sim)