<jupyter_start><jupyter_text>Monte Carlo Example[![Latest release](https://badgen.net/github/release/Naereen/Strapdown.js)](https://github.com/eabarnes1010/course_objective_analysis/tree/main/code)[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/eabarnes1010/course_objective_analysis/blob/main/code/monte_carlo_example.ipynb)An illustration of using Monte Carlo techniques.  The problemThis example is based on the "Monte Carlo simulation" example in Chapter 1 of the lecture notes. The idea is as follows...In January (31 days), the maximum daily temperature was 2.2 standard deviations from the climatological mean temperature. If we assume that daily temperatures are normally distributed, how rare is it to have a maximum of 2.2$\sigma$ or greater in 31 daily samples?  The approachWhile we have lots of nice statistical tests for the mean, and even the variance, we've ventured outside of this into the realm of maxima! So what do we do? Since the maximum is a pretty common quantity, there are in fact tests we could use. However, we could also take a Monte Carlo approach (which is what we will do given the title of this example) to solve our problem.The idea is that since we _know_ the underlying population is normal, we can just make some "fake" data and play around with it to find out how rare a 2.2$\sigma$ maximum really is. The first step, of course, is to import our packages.<jupyter_code>#.............................................
# IMPORT STATEMENTS
#.............................................
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

#.............................................
# PLOTTING COMMANDS 
#.............................................
import matplotlib as mpl

# set figure defaults
mpl.rcParams['figure.dpi'] = 150
plt.rcParams['figure.figsize'] = (12.0/2, 8.0/2)<jupyter_output><empty_output><jupyter_text>Setup our variables. The "drawn_max" variable is where we set the maximum we actually got in our 31 samples, and "sample_length" sets the number of values in our sample (i.e. N).<jupyter_code>drawn_max = 2.2
sample_length = 31<jupyter_output><empty_output><jupyter_text>Next, we are going to make a long list of values drawn from a random normal distribution (call this data set "Z"). You can think of this as our _popoulation bucket_. We will also initialize an empty array "M" where we will ultimately store the maximum values from many many different samples of length N=31.<jupyter_code>Z = np.random.randn(1_000_000,1)
M = np.empty([100_000, 1])<jupyter_output><empty_output><jupyter_text>Next we do the heavy lifting. The idea is that since we know our underlying data is normal, we are going to randomly grab N=31 values from this huge list (Z), calculate the maximum in the sample, save it, then rinse and repeat. We will do this 10,000,000 times.<jupyter_code>for iloop in range(np.size(M)):
    ip = np.random.randint(low=0,high=Z.shape[0],size=sample_length) 
    M[iloop] = np.max(Z[ip])<jupyter_output><empty_output><jupyter_text>Ok - so what do we have now? Well, we now have 100,000 maximum values from samples of N=31. So, we can plot the distribution of these maxima (M), and see what they look like! We can then see how rare our actual maximum that we actually got (i.e. drawn_max = 2.2) is. Check it out.<jupyter_code>plt.figure()
bin_width = .1
bin_list = np.arange(-1,10,bin_width)
n, bins = np.histogram(M, bins=bin_list, density=False)
plt.plot([drawn_max, drawn_max],[0, .12],color='red',linewidth=3,linestyle='--', label='actual Jan. maximum')
plt.bar(bins[0:-1],n/float(len(M)),bin_width, facecolor='blue', alpha=0.4)
plt.plot(bins[0:-1]+bin_width/2,n/float(len(M)),color='blue', label='random samples')

plt.xlabel('maximum value')
plt.ylabel('frequency')
titlename = 'Maximum value for sample of length N = ' + str(sample_length) 
plt.title(titlename)
plt.legend(fontsize=8)

plt.show()<jupyter_output><empty_output><jupyter_text>The blue shading and line denotes the frequency of the maxima we got from the 100,000 samples of length N=31 (drawn from a normal don't forget). The red line denotes what we _actually_ got from our real data. How rare is it? Well, just looking at the plot it is pretty clear that a maximum of 2.2 is not rare at all. Of course, we can calculate the number of values in M above 2.2 to see the probability.<jupyter_code>100.*float(sum(M>drawn_max)/np.size(M))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13340732519422863, 'fraction_numerical': 0.023307436182019976, 'mean_word_length': 4.161512027491409, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 26, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '45748', 'license': 'MIT', 'n_tokens_mistral': 1432, 'n_tokens_neox': 1326, 'n_words': 532, 'path': 'code/monte_carlo_example.ipynb', 'repo_name': 'eabarnes1010/course_objective_analysis', 'size': '4505'}