<jupyter_start><jupyter_text>Monte Carlo Modelinghttp://www.analyzemath.com/statistics/probability_questions.html<jupyter_code>import numpy as np
from scipy.stats import binom, itemfreq
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Monte Carlo simulation involves using randomness and iteration to replicate some process.  We've seen many random number generators before.<jupyter_code># Play along!
# In what ways have we already generated random numbers in this course?
One way to get a handle on what Monte Carlo is all about is to simulate some things we already know:<jupyter_output><empty_output><jupyter_text>Simple calculations<jupyter_code># What percentage of the numbers between 1 and 100 are greater than 55?

# This is directly answerable:
len(np.arange(56,101)) / len(np.arange(1,101))
# But we could also simulate this:
k = 100
x = np.random.randint(1,101,k)
np.sum(x > 55) / k
# Play along!  Simulate rolling a six-sided die 100 times.  How many times did you get an even number?
# Play along!  Using Monte Carlo, solve:
# Two coins are tossed, find the probability that two heads are obtained. 
# Play along!
# Two dice are rolled, find the probability that the sum is between 3 and 11 (inclusive)<jupyter_output><empty_output><jupyter_text>With means and standard deviations<jupyter_code># For easy problems, we might just run the direct calculations, but when the situation gets even a little messy, randomness 
# can be very effective.

# I have a sample of 10,000 rats with a mean weight of 230 grams (SD = 25).  What percentage of my rats are above 250 grams?

k = 10000
x = np.random.normal(loc=230,scale=25,size=k)
np.sum(x>250) / k
# I have invented a growth serum that boosts rat weight by M=15.5, SD=22.12
# Give the growth serum to the rats!  (Muhahaha!) and tell me how many are over 250 grams?

k = 10000
rats = np.random.normal(loc=230,scale=25,size=k)
serum = np.random.normal(15.5,12.12,k)
rats = rats + serum
np.sum(rats>250) / k
# These rats are too expensive!  Experiment with k to find out the minimum you think I neeed before the output gets unreliable.
# You'll need to repeat each observations a few times and watch the output to be sure.

k = 50  # Adjust this number at look at the result.  How many rats do I need for a stable, reliable result?
rats = np.random.normal(loc=230,scale=25,size=k)
serum = np.random.normal(15.5,12.12,k)
rats = rats + serum
np.sum(rats>250) / k

# Play along!  Write your answer here, along with a justification.<jupyter_output><empty_output><jupyter_text>random.choice<jupyter_code># One of the nice ways to simulate a random process is with random.choice:

# This will randomly select one of these in the list.  Rerun it to see!

np.random.choice(['Fred','Velma','Shaggy','Scooby','Daphne'])
# You can select multiple people:

selection = np.random.choice(['Fred','Velma','Shaggy','Scooby','Daphne'],size=3)
print(selection)
# You can select multiple people, and you can indicate whether you want to replace (to avoid duplicates)

selection = np.random.choice(['Fred','Velma','Shaggy','Scooby','Daphne'],size=3,replace=False)
print(selection)
# If you aren't replacing, you can't pick a sample larger than the length of your list.
selection = np.random.choice(['Fred','Velman','Shaggy','Scooby','Daphne'],size=6,replace=False)
print(selection)
# Simulate a deck of cards:

suit = ['Clubs','Spades','Diamonds','Hearts']
value = list('23456789')
value.extend(['10','Jack','Queen','King','Ace'])
draw = (np.random.choice(value),np.random.choice(suit))
print('You just drew the',draw[0],'of',draw[1],'.')
# Let's simulate the WHOLE deck using a Pandas Dataframe and np.tile

values = np.tile(value,(4,1))
print(values)
suits = np.tile(suit,(13,1))
print(suits)
print(np.shape(values))
print(np.shape(suits))

# Play along!  Given this information, how can we put our suits and values together?
print(values.T.flatten())
print(suits.flatten())
deck = pd.DataFrame({'suit':suits.flatten(),'value':values.flatten()})
print(deck)
# Play along!  Write out that deck for future use.

# Lets add some attribute data to our deck:

deck['color'] = 'black'
deck.loc[(deck.suit=='Diamonds') | (deck.suit=='Hearts'),'color'] = 'red'
deck.head()
# Let's draw a card!

deck.sample()
# Draw 10,000 cards!  Replacement is set to False by default in a dataframe

k = 10000
my_sample = deck.sample(k,replace=True)
my_sample.head()
# Play along!  How many times did you draw a black card?
# Play along!  How many times did you draw the 3 of diamonds?  How many would you expect?
# Play along! How many queens did you get?<jupyter_output><empty_output><jupyter_text>Replacement<jupyter_code># Suppose we wanted to investigate the odd of getting a flush in poker:

hand = deck.sample(5,replace=False)
print(hand)
# How would I know if I had a flush?  I can use the "unique" method to see if I have just one suit:

hand['suit'].unique()
# Or I could use value_counts to get the number:

hand['suit'].value_counts()
# On this first, cut, imagaine we just want to see if we got a flush or not:

hand = deck.sample(5,replace=False)
len(hand['suit'].unique()) == 1
# We COULD simulate this with a loop:

k = 10000
is_a_flush = np.zeros(k)
for i in range(k):
    hand = deck.sample(5,replace=False)
    is_a_flush[i] = len(hand['suit'].unique()) == 1
np.sum(is_a_flush)
# But you might notice this takes a while!
# Play along!  How long does it take?  
# We could also imagine our deck of cards as numbers from 0 to 51
idx = np.arange(52)

# And then we could shuffle that:
np.random.shuffle(idx)

# And see that the "deck" has been shuffled.
print(idx)

# We can recover our shuffled deck from pandas this way:
deck.loc[idx]
# If we only wanted the first 5 cards as our "hand" we could do this:

# Our population is the 52 cards (0 to 51) and we're asking for 5 of them.
idx = np.random.choice(52,5,replace=False)
hand= deck.loc[idx]
print(hand)
# We could also only pay attention to the suit information:

idx = np.random.choice(52,5,replace=False)
my_suits = deck.loc[idx,'suit']
print(my_suits)

# Solving this with idx-based methods:
k = 10000
is_a_flush = np.zeros(k)
for i in range(k):
    idx = np.random.choice(52,5,replace=False)
    hand= deck.loc[idx]
np.sum(is_a_flush)
# Play along!  Can you figure out how to discard the cards that aren't a flush (up to three), and draw three new cards?
# What's the probability of getting a flush if you can discard up to three cards? (this is a hard one!)<jupyter_output><empty_output><jupyter_text>Simulating a process.  REAL shuffling. Shuffling a deck of cards involves splitting the deck into two piles, and then having cards drop, more or less randomly, from either pile.  Here, I split the piles according to the binomial distribution, and then choose a card from one pile or the other proportionally based on its height (i.e., the number of cards in it).Further Reading:https://github.com/thomaspingel/riffle-shuffle/blob/master/A%20Quick%20Look%20at%20Shuffling.ipynb<jupyter_code>def riffle(deck,n_times=1):
    
    n_cards = deck.size
    
    for j in range(n_times):
        # Split the deck into two piles (left and right) according to the binomial distribution at position "pos"
        pos = binom(n_cards,.5).rvs(1)[0]
        left, right = deck[:pos].copy(), deck[pos:].copy()
        
        # Initialize the deck.  Cards will be placed back here.
        deck[:] = -1

        # The probability of grabbing the card from left or right
        # is proportional to the size of that pile
        for i in range(n_cards):
            p = np.random.rand()
            thresh = left.size / (left.size + right.size)
            if p < thresh:
                deck[i] = left[0]
                left = np.delete(left,0)
            else:
                deck[i] = right[0]
                right = np.delete(right,0)
# Simulate one shuffle

deck = np.arange(52)
riffle(deck)
print(deck)
# And visualize the "rising sequences" in the shuffled deck as an visual indicator of randomness.  An unshuffled deck would be 
# one continuous rising sequence, and a well shuffled deck would appear like random dots (the starting position is unrelated)
# to the ending position.

fix, ax = plt.subplots(figsize=(10,4))
plt.subplot(121)
deck = np.arange(52)
riffle(deck)
plt.plot(deck,'.')
plt.title('One Shuffle')

plt.subplot(122)
deck = np.arange(52)
riffle(deck,100)
plt.plot(deck,'.')
plt.title('One Hundred Shuffles')
# Next, we'll run a simulation look at where the top card ends up.  On a well-shuffled deck, the top card
# should be equally likely to appear anywhere.  We'll use numpy's "where" to find out where it landed.

# Rerun this code block to see where the top card ends up with different shuffles.

deck = np.arange(52)
riffle(deck)
print(deck)
np.where(deck==0)
# Run a simulation of twenty decks, shuffled 20 times each and tracking where the top card ended up.
# In a well-shuffled deck, the top card should be equally likely to appear anywhere.

n_shuffles = 20
n_decks = 20
top_card = 255*np.ones((n_decks,n_shuffles),dtype=np.uint8)
for i in range(n_decks):
    deck = np.arange(52)
    for j in range(n_shuffles):
        riffle(deck)
        top_card[i,j] = np.where(deck==0)[0]
        
# Visualize the results.  Notice most of the time at the beginning, the top card ends up near the top (position 0),
# but with enough shuffles, it's essentially random.
plt.imshow(top_card,interpolation="none")
plt.xlabel('shuffles')
plt.ylabel('deck')
plt.colorbar()<jupyter_output><empty_output><jupyter_text>Finally, we can track the overall distribution of of top cards over many, many shuffles.  Think of this like the image above, but where the number of decks is much larger (say, 1000).<jupyter_code>n_shuffles = 16
n_decks = 1000
top_card = 255*np.ones((n_decks,n_shuffles),dtype=np.uint8)
for i in range(n_decks):
    deck = np.arange(52)
    for j in range(n_shuffles):
        riffle(deck)
        top_card[i,j] = np.where(deck==0)[0]
# Visualize the results; we'll turn it on its side.
fig, ax = plt.subplots(figsize=(10,3))
plt.imshow(top_card.T,interpolation="none")
plt.xlabel('deck')
plt.ylabel('shuffles')
plt.colorbar()
# To see the pattern, we'll construct a histogram for each shuffle, so we can track the likelihood of where
# the top card ends up after each shuffle.

for i in range(15):
    plt.subplot(4,4,i+2)
    plt.hist(top_card[:,i],bins=52)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11032707161247259, 'fraction_numerical': 0.023934394965195005, 'mean_word_length': 4.364705882352941, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 37, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '17536', 'license': 'MIT', 'n_tokens_mistral': 3674, 'n_tokens_neox': 3402, 'n_words': 1263, 'path': '09. Intro to Simulation in Python.ipynb', 'repo_name': 'thomaspingel/intro_to_geospatial_python', 'size': '10487'}