text: 
  # coding: utf-8
  
  # # Numerical simulations of the Monty-Hall "paradox"
  # This short notebook aims at simulating trials of the so-called [Monty-Hall problem](https://en.wikipedia.org/wiki/Monty_Hall_problem), and thus helping to convince about the result thanks to a numerical results rather than a possibly-unclear proof.
  
  # ----
  # ## Definition of the problem
  # 
  
  # There is $M \geq 3$ doors, and behind only one of them there is a treasure.
  # The goal of the player is to find the treasure, following this game:
  # 
  # 1. The player first chose a door, but does not open it yet,
  # 2. All remaining doors but one are opened, and none of them contain the treasure. The player sees $M - 2$ bad doors,
  # 3. So there is just 2 doors, the one first chosen, and the last one. She knows the treasure is behind one of them,
  # 4. And she has to decide if she wants to **stay on her initial choice**, or **switches to the last door**.
  # 
  # Finally, the chosen door is opened, and the player wins this round if she found the treasure.
  
  # ---
  # The goal of this notebook is to numerically prove that the choice of always switching to the last door is the best one.
  
  # ## Starting our numerical simulations
  # We start by importing some modules, then define a function, `randomAllocation()`, to generate a random allocation of treasures behind the `M` doors.
  # 
  # > *Note:* all the code below is generic for any $M \geq 3$, but `M = 3` is used to keep small and clear visualizations.
  
  # In[1]:
  
  
  import random
  
  
  # In[2]:
  
  
  M = 3
  allocation = [False] * (M - 1) + [True]  # Only 1 treasure!
  assert set(allocation) == {True, False}  # Check: only True and False
  assert sum(allocation) == 1              # Check: only 1 treasure!
  
  
  # Just to check:
  
  # In[3]:
  
  
  allocation
  
  
  # We can generate a random spot for the treasure by simply shuffling (with [`random.shuffle()`](https://docs.python.org/3/library/random.html#random.shuffle)):
  
  # In[4]:
  
  
  def randomAllocation():
      r = allocation[:]
      random.shuffle(r)
      return r
  
  
  # Let us quickly check this function `randomAllocation()`:
  
  # In[5]:
  
  
  for _ in range(10):
      print(randomAllocation())
  
  
  # We need to write a small function to simulate the choice of the door to show to the player, `show()`:
  
  # In[6]:
  
  
  def last(r, i):
      # Select a random index corresponding of the door we keep
      if r[i]:  # She found the treasure, returning a random last door
          return random.choice([j for (j, v) in enumerate(r) if j != i])
      else:     # She didn't find the treasure, returning the treasure door
          # Indeed, the game only removes door that don't contain the treasure
          return random.choice([j for (j, v) in enumerate(r) if j != i and v])
  
  
  # In[7]:
  
  
  for _ in range(10):
      r = randomAllocation()
      i = random.randint(0, M - 1)
      j = last(r, i)
      print("- r =", r, "i =", i, "and last(r, i) =", j)
      print("  Stay on", r[i], "or go to", r[j], "?")
  
  
  # We need a function to simulate the first choice of the player, and a simple choice is to select a uniform choice:
  
  # In[8]:
  
  
  def firstChoice():
      global M
      # Uniform first choice
      return random.randint(0, M - 1)
  
  
  # Now we can simulate a game, for a certain left-to-be-written function `strategy()` that decides to keep or to change the initial choice.
  
  # In[9]:
  
  
  def simulate(stayOrNot):
      # Random spot for the treasure
      r = randomAllocation()
      # Initial choice
      i = firstChoice()
      # Which door are remove, or equivalently which is the last one to be there?
      j = last(r, i)
      assert {r[i], r[j]} == {False, True}  # There is still the treasure and only one
      stay = stayOrNot()
      if stay:
          return r[i]
      else:
          return r[j]
  
  
  # We can simulate many outcome of the game for one strategy, and return the number of time it won (i.e. average number of time it found the good door, by finding `r[i] = True` or `r[j] = True`):
  
  # In[10]:
  
  
  N = 10000
  
  def simulateManyGames(stayOrNot):
      global N
      results = [simulate(stayOrNot) for _ in range(N)]
      return sum(results)
  
  
  # ----
  # ## Comparing two strategies, on many randomized trials
  # We will simulate the two strategies, `keep()` vs `change()`, on $N = 10000$ randomized games.
  
  # ### Keeping our first choice, `keep()`
  
  # In[11]:
  
  
  def keep():
      return True  # True == also stay on our first choice
  
  
  # In[12]:
  
  
  rate = simulateManyGames(keep)
  print("- For", N, "simulations, the strategy 'keep' has won", rate, "of the trials...")
  print("  ==> proportion = {:.2%}.".format(rate / float(N)))
  
  
  # $\implies$ We find a chance of winning of about $\frac{1}{M} = \frac{1}{3}$ for this strategy, which is very logical as only the initial choice matters, and due to the uniform location of the treasure behind the $M = 3$ doors, and the uniform first choice with `firstChoice()`.
  
  # ### Changing our first choice, `change()`
  
  # In[13]:
  
  
  def change():
      return False  # False == never stay, ie always chose the last door
  
  
  # In[14]:
  
  
  rate = simulateManyGames(change)
  print("- For", N, "simulations, the strategy 'change' has won", rate, " of the trials...")
  print("  ==> proportion = {:.2%}.".format(rate / float(N)))
  
  
  # $\implies$ We find a chance of winning of about $\frac{M - 1}{M} = \frac{2}{3}$ for this strategy, which is less logical.
  # 
  # Due to the uniform location of the treasure behind the $M = 3$ doors, and the uniform first choice with `firstChoice()`, we have a $\frac{1}{M}$ chance of finding the treasure from the first time.
  # 
  # 1. If we found it, the last door does not contain the treasure, hence **we loose** as we switch to it,
  # 2. However, if we did not find it, the last door has to contain the treasure, hence **we win** as we switch to it deterministically (i.e. always).
  # 
  # The first case has probability $\frac{1}{M}$, probability of loosing, and the second case has probability $\frac{M - 1}{M}$, probability of loosing.
  # 
  # $\implies$ Conclusion : this strategy `change()` has a chance of winning of $\frac{2}{3}$, far better than the chance of $\frac{1}{3}$ for `stay()`.
  # 
  # > We proved numerically the results given and explained [here on the Wikipedia page](https://en.wikipedia.org/wiki/Monty_Hall_problem#Direct_calculation). Great!
  
  # ### Bernoulli choice
  # We can try a randomized strategy, a simple one can be to follow the decision of a (biased) coin:
  # 
  # 1. Toss a coin (with fixed head probability equals to $p \in [0, 1]$)
  # 2. If head, switch to the last door.
  
  # In[15]:
  
  
  def bernoulli(p=0.5):
      return random.random() < p
  
  
  # In[16]:
  
  
  rate = simulateManyGames(bernoulli)
  print("- For", N, "simulations, the strategy 'bernoulli' has won", rate, " of the trials...")
  print("  ==> proportion = {:.2%}.".format(rate / float(N)))
  
  
  # Now we can try different values for $p$, and plot the resulting chance of winning the game as a function of $p$.
  # Hopefully, it should be monotonic, confirming the result explained above.
  
  # In[17]:
  
  
  import numpy as np
  import matplotlib.pyplot as plt
  
  
  # We generate lots of values for $p$, then a function `stratBernoulli()` to create the strategy described above, for some $p \in [0, 1]$.
  
  # In[18]:
  
  
  values_p = np.linspace(0, 1, 500)
  
  def stratBernoulli(p):
      def stayOrNot():
          return bernoulli(p=p)
      return stayOrNot
  
  
  # Let finally do all the simulations, and store the empirical probability of winning the game when following a Bernoulli strategy.
  # 
  # > This line takes about $4$ minutes on my laptop, it's not that quick.
  
  # In[19]:
  
  
  chance_of_winning = [simulateManyGames(stratBernoulli(p)) / float(N) for p in values_p]
  
  
  # In[20]:
  
  
  plt.figure()
  plt.plot(values_p, chance_of_winning, 'r')
  plt.title("Monty-Hall paradox with {} doors ({} random simulation)".format(M, N))
  plt.xlabel("Probability $p$ of staying on our first choice (Bernoulli strategy)")
  plt.ylabel("Probability of winning")
  plt.ylim(0, 1)
  plt.yticks(np.linspace(0, 1, 11))
  plt.show()
  
  
  # ---
  # ## Examples with $M = 100$ doors
  
  # In[27]:
  
  
  def completeSimu():
      global M
      global N
      allocation = [False] * (M - 1) + [True]  # Only 1 treasure!
      
      def randomAllocation():
          r = allocation[:]
          random.shuffle(r)
          return r
      
      def last(r, i):
          # Select a random index corresponding of the door we keep
          if r[i]:  # She found the treasure, returning a random last door
              return random.choice([j for (j, v) in enumerate(r) if j != i])
          else:     # She didn't find the treasure, returning the treasure door
              # Indeed, the game only removes door that don't contain the treasure
              return random.choice([j for (j, v) in enumerate(r) if j != i and v])
      
      def simulate(stayOrNot):
          # Random spot for the treasure
          r = randomAllocation()
          # Initial choice
          i = firstChoice()
          # Which door are remove, or equivalently which is the last one to be there?
          j = last(r, i)
          stay = stayOrNot()
          if stay:
              return r[i]
          else:
              return r[j]
  
      def simulateManyGames(stayOrNot):
          global N
          results = [simulate(stayOrNot) for _ in range(N)]
          return sum(results)
  
      values_p = np.linspace(0, 1, 300)
      chance_of_winning = [simulateManyGames(stratBernoulli(p)) / float(N) for p in values_p]
      plt.figure()
      plt.plot(values_p, chance_of_winning, 'r')
      plt.title("Monty-Hall paradox with {} doors ({} random simulation)".format(M, N))
      plt.xlabel("Probability $p$ of staying on our first choice (Bernoulli strategy)")
      plt.ylabel("Probability of winning")
      plt.ylim(0, 1)
      plt.yticks(np.linspace(0, 1, 11))
      plt.show()
  
  
  # In[28]:
  
  
  M = 4
  completeSimu()
  
  
  # ### Last plot, for $M = 100$
  # It clearly shows the linear behavior we expected.
  
  # In[29]:
  
  
  M = 100
  completeSimu()
  
  
  # $\implies$ it is now clear that the best strategy is the **always change and select the other door**.
  
  # ----
  # ## Conclusion
  # As we saw it, it is clear from the numerical experiments presented above that the strategy of **always changing our first choice** (`change()`), is statistically better than the strategy of **keeping our first choice** (`keep()`).
  # 
  # > [See this repository for other Python notebook doing numerical simulations](https://github.com/Naereen/notebooks/tree/master/simus/).
  