# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Homework 1 - Probability Estimation

# ## Intro
#
# ### The notebooks
#
# The hands-on part (תרגילים רטובים) of the homework assignment in this course will be given in the form of Colaboratory notebooks such as this one.
#
# These notebooks combine cells of formatted text, such as this cell, along with Python  code cells, such as this one:

x = 3
print('x is equal to ' + str(x))

# ### Making your own copy of the notebook
#
# When opening this notebook using the [public link](https://colab.research.google.com/github/technion046195/semester_2019_spring/blob/master/assignments/homework_1_probability_estimation.ipynb), you are opening the public shared copy of this notebook in a **viewing only mode**. In this mode you will not be able to make any changes. You should see the status **VIEWING** on the upper right corner of the page.
#
# To make changes to the notebook and be able to save them you will need to create your own copy of this notebook. For that, you will need a Google account.
#
# 1. In the rare case that you don't already have an account, you can create one [here](https://accounts.google.com/SignUp).
# 2. Make sure that you are logged in into your account. If you are not logged in you will have a **Sign in** button on the upper right corner of this page, Click on it and log into your account.
#
# Now, to make a personal copy of this notebook open the **file** menu from the menu bar and select **Save a copy in drive...** . The notebook's name at the top of this page should change to "**Copy of homework_1_probability_estimation.ipynb**" and the status on the upper right corner should change into **EDITING** 
#
# #### Notebooks name and location
#
# By default, the notebook you have just copied is saved to your [Google drive](https://drive.google.com/drive) under a folder called **Colab Notebooks**. You can view this folder by selecting **Locate in Drive** from the **file** menu.
#
# **Optional**: You can rename the notebook by clicking on the notebook's name in the tile of this page.
#
# Now we are ready to start editing 👍

# ### ✍️ Saving and loading notebooks
#
# - To save a notebook press **Ctrl+S** (or select **save** from the **file** menu).
# - To load a saved notebook, open the **Colab Notebooks** folder in your [Google drive](https://drive.google.com/drive), and double click on the saved notebook. (If you did not rename it should appear as "**Copy of homework_1_probability_estimation.ipynb**") Note that in Google Drive you can have more than 1 file with the same name in the same folder.
#
# Before you start, make sure that you can save and load this notebook. Start be replacing the ID numbers in the code cell below  with our own ID numbers. Then save and load this notebook, and make sure that the IDs were indeed saved.

# +
## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
# Replace the IDs bellow with our own
student1_id = '012345678'
student2_id = '012345678'
## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%

print('Hello ' + student1_id + ' & ' + student2_id)
# -

# ### Expending all cells
#
# When you open a new notebook some sections might be hidden (collapsed), showing only the title of the section with a "↳ # cells hidden" note beneath it. Before you start, make show you expand all the sections. You can do so by pressing **Ctrl+[** or by clicking on all the "↳ # cells hidden" lines.
#

# ### Running a cell
#
# To run a specific code cell, such as the one above, select the cell using the mouse or the up/down keys, and press **Shift+Enter**. 
#
# Run the code cell above. You should see the output of this cell beneath it.

# ### ✍️ Your assignments
#
# In these assignments, you will have to follow the instruction in the notebooks and fill in the missing code blocks according to the instructions. The cells containing the instructions will be marked with ✍️. The line of code which you will be expected to fill in will be surrounded by the two following comments:
#
#
# ``` python
# ## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
#
# ## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%
# ```
#
# For example: fill in the following code block to assign the value $x+10$ to the variable $y$ (replace the ... with your code):

# +
## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
# assign x + 10 to y:
y = ...
## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%

print('y is equal to ' + str(y))
# -

#
# ### Running order
#
# The notebook engine lets you run the code cells in any order you prefer, but to make everything work as it was designed to, you need to make sure you run the code cell in the correct order (in which they appear in the notebook). If you want to start over and reset the Python's state (all code cells which ran previously) press **Ctrl+M I** (**Ctrl+M** followed by **I**), or select **Restart runtime...** from the **Runtime** menu.
#
# Restart the runtime and run all 3 code cells above this line.

# ### The last line in a code cell
#
# For convenience, if the last line of a code cell returns some value, it will be appended to the cell's output as well. For example, run the following cell:

## The output of this code cell would be y to the power of 2
print('This is a ordinary print line')
y ** 2  # The result of this line will appear in the output as well

# ## Python
#
# ### Learning Python
#
# The goal of these homework assignments is not to teach you Python, and it is expected that if you have no previous experience in Python it will take you a while to get the hang of it.
#
# We did try however to add some explanation for some of the concepts and the functions which we will be using. This is in order to make the first steps easier, especially for those will some background in Matlab, and to help you get through this assignment without the need to first thoroughly learning Python. 
#
# However, you will probably need every now and then to relay on the web for help with solving bugs and understanding parts of the code. There is a large number of great resources available online for learning Python, and these resources are be much better than any introduction which we can write on the subject. We have listed a few in the "online resources section" below. You are encouraged to check them out whenever you get stuck.
#
# And of course: <http://bfy.tw/N2JF>
#
# ### A Comment About Python Versions
#
# Few versions of Python have been released so far. The latest (and most popular) version today is Python 3 although there a lot of places still using Python 2. When looking at online resources keep in mind that there are some minor syntax changes between the two versions. In this course, we will be using version 3 and these notebooks by using it by default, so in general, you do not have to worry about it.

# ### Python Packages
#
# The Python language by itself offers only the fundamental functionality of a programing language such as basic variables (integers, floats, strings, etc.), if statements, loops, functions, etc. . The strength of the language comes, among others, from the vast amount of packages which can be imported and used for performing advanced tasks.
#
# In this course we will be using the following very common packages:
#
# <dl>
#   <dt>[NumPy](http://www.numpy.org/)</dt>
#   <dd>A package for defining and working with arrays (of various dimension). It contains stuff such as linear algebra operators, complex array indexing, vectorized operations, random numbers generation, Fourier transformations and more. It basically gives Python the basic Matlab like capabilities.</dd>
#
#   <dt>[SciPy](https://docs.scipy.org/doc/scipy-1.2.0/reference/)</dt>
#   <dd>This package extends the NumPy package and adds stuff such as interpolations, some signal processing tools, sparse matrices and more.</dd>
#   
#   <dt>[Matplotlib](https://matplotlib.org/)</dt>
#   <dd>A plotting library. When it comes to plotting graphs python is still behind the curve in comparison to some other programming languages such as Matlab, R and JavaScript especially when it comes to the interactivity of the figure and the ability to edit it. In recent years, [a lot of new promising alternatives](https://mode.com/blog/python-data-visualization-libraries) started to appear, but matplotlib is still the most versatile and the easiest tool to start with. Matplotlib is build to have a plotting interface similar to Matlab's</dd>
#
#   <dt>[Pandas](https://pandas.pydata.org/)</dt>
#   <dd>Pandas is a great tool for significantly simplifying the work with datasets. It does so through the uses of DataFrame objects. Working with DataFrames usually takes quite a while to master, but it is worth the effort, especially when working with data a lot</dd>
#   
#   <dt>[Scikit-Learn](https://scikit-learn.org/stable/)</dt>
#   <dd>A very popular machine learning package which contains implementation for all most of the material covered in this course and much more.</dd>
# </dl>
#
# *Later in this course, we will also look at the [PyTorch](https://pytorch.org/) and [XGBoost](https://xgboost.readthedocs.io/en/latest/) packages*

# ### Importing Packages
#
# The following code cell imports the NumPy, Pandas and Matplotlib packages. The exact syntax which is used to import the packages is not important at this point.

# +
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

## This line makes matplotlib plot the figures inside the notebook
# %matplotlib inline
# -

# ### Usage
#
# For example, the following code cell uses the NumPy package to print the powers of 2 from $2^0$ to $2^{10}$:

# +
## Define an array running for 0 to 10:
powers = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

## Calculate the powers of 2:
print(2 ** powers)
# -

# ### ✍️ Warming up
#
# Edit the code below so that it will print the multiplication table (לוח הכפל). Use the following NumPy functionalities:
#
# - A 2D array is defined in NumPy by using two sets of square brackets (see the definition of *a* below).
# - The transpose of a 2D array *x* is created by *x.T*.
# - Matrix multiplication between to 2D arrays *A* and *B* in NumPy is done using the *@* operator: *A @ B*

# +
a = np.array([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]])

## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
## Fill in the following line so that multiplication_table will store the 11 x 11 multiplication table (of the numbers between 0 and 10).
multiplication_table = ...
## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%

print(multiplication_table)
# -

# Your output should be:
#
# ``` python
# [[  0   0   0   0   0   0   0   0   0   0   0]
#  [  0   1   2   3   4   5   6   7   8   9  10]
#  [  0   2   4   6   8  10  12  14  16  18  20]
#  [  0   3   6   9  12  15  18  21  24  27  30]
#  [  0   4   8  12  16  20  24  28  32  36  40]
#  [  0   5  10  15  20  25  30  35  40  45  50]
#  [  0   6  12  18  24  30  36  42  48  54  60]
#  [  0   7  14  21  28  35  42  49  56  63  70]
#  [  0   8  16  24  32  40  48  56  64  72  80]
#  [  0   9  18  27  36  45  54  63  72  81  90]
#  [  0  10  20  30  40  50  60  70  80  90 100]]
# ```

# ### Some online resourses
#
# If you are looking for a quick rough overview take a llok at this:
# - [A Whirlwind Tour of Python](https://jakevdp.github.io/WhirlwindTourOfPython/): A short overview of the language which is oriented towards data science applications. Another resource from the same author is the [Python data science handbook](https://jakevdp.github.io/PythonDataScienceHandbook/)
#
# If you are looking for a more thorough online course try ine of these:
# - [Python @ Tutorials Point](https://www.tutorialspoint.com/python/): One of the many free online Python courses.
# - [Python for data science](https://www.datacamp.com/courses/intro-to-python-for-data-science): Another online Python course. Free but requires registration.
#
# If you simply want to get started and have some Matlab experience, you can make use of these two web pages:
# - [NumPy for Matlab users 1](https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html): A good place to start with to get a grip of NumPy when coming from matlab.
# - [NumPy for Matlab users 2](http://mathesaurus.sourceforge.net/matlab-numpy.html): Another web page similar to the one above with  a slightly different set of commands.
#
# Once you feel comfortable with the basic Python syntax you can take a look at these short three guides to improve your NumPy, Pandas and Colab  capabilities:
# - [NumPy's quickstart tutorial](https://www.numpy.org/devdocs/user/quickstart.html): A quickstart tutorial on NumPy from NumPy's official website.
# - [10 Minutes to pandas](https://pandas.pydata.org/pandas-docs/stable/getting_started/10min.html): A quickstart tutorial on Pandas from Pandas' official website.
# - [Welcome to Google Colabratory](https://colab.research.google.com/notebooks/welcome.ipynb): The official introduction to Google Colabratory.

# ## The Ultra-Trail du Mont-Blanc
#
# In this assignment, we will examine the distribution of results in The Ultra-Trail du Mont-Blanc race.
#
# Taken from Wikipedia:
#
# >  It takes place once a year on either the last weekend in August or the first weekend of September in the Alps, and follows the route of the Tour du Mont Blanc through France, Italy and Switzerland. It has a distance of approximately 171 kilometres (106 mi), and a total elevation gain of around 10,040 metres (32,940 ft). It is widely regarded as one of the most difficult foot races in the world, and one of the largest with more than 2,500 starters. It is one race during a week-long festival based around Chamonix in France. The races have strict entry and qualification requirements attained by accumulating enough race points through qualifying trail races over the previous two-year period. In 2016 and 2017, 42% and 34% of runners did not finish the UTMB race.
#
# > While the best runners complete the loop in slightly more than 20 hours, most runners take 32 to 46 hours to reach the finish line. Most runners will have to run through two nights in order to complete the race.
#
#
#

# We will start by loading the data using the Pandas package. The data can be as a [CSV](https://en.wikipedia.org/wiki/Comma-separated_values) file [here](https://technion046195.github.io/semester_2019_spring)

# +
data_file = 'https://technion046195.github.io/semester_2019_spring/datasets/utmb.csv'

## Loading the data
dataset = pd.read_csv(data_file)
# -

# *dataset* is now **DataFrame** which contains our data.
#
# The following code cell prints the number of rows in the dataset and the 10 first rows in it:

# +
## Print the number of rows in the data set
number_of_rows = len(dataset)
print('Number of rows in the dataset: N=' + str(number_of_rows))

## Show the first 10 rows
dataset.head(10)
# -

# ### The Data Fields and Types
#
# This dataset contains the results for all the runners which finished the race over all the years in which the race took place. The data fields (the columns) for each runner are as following:
#
# - **Name**: The runner's name.
# - **Year**: The year of the race.
# - **Result**: The runner's result (the time it took him to finish the race) in hours. The race ends after about 47 hours. Therefore, these results are between 0 and 47.
# - **Age group**: The age group to which the runner belongs to: 0 - the youngest (10-22) and 6 - The oldest (80+)
# - **Gender**: The runner's gender: 0-male, 1-female.
# - **Rank**: The overall score rank of the runner in that specific year.

# To extract a column from the DataFrame as a NumPy array we can use the following code:

dataset['Year'].values 

# ### ✍️ The == operator
#
# Correct the code below so that it will count and print the number of runners which have finished the race in 2014.
# - Use the equality operator "a==b" to create a boolean 1D array.
# - Use the summation function *np.sum(x)* to sum over the values of the array.

# +
## Define the array of the years:
years = dataset['Year'].values
           
## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
## Correct the following line:
## Example: the number of runners which have finished the race since 2010 would be: np.sum(years > 2010)
number_of_runners_in_2014 = ...
## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%

print('The number of runners which finished the race in 2014: ' + str(number_of_runners_in_2014))
# -

# ### ✍️ The & operator
#
# Correct the code below so that it will count and print the number of **male** runners which have finished the race in 2014.
# - Use the logical AND operator "a & b" two combine two boolean arrays

# +
## Define the array of the genders:
genders = dataset['Gender'].values

## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
## Fill in the following line:
number_of_male_runners_in_2014 = ...
## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%

print('The number of male runners which finished the race in 2014: ' + str(number_of_male_runners_in_2014))


# -

# ### ✍️ Calculating the histogram
#
# In the cell below we define a function called *calculate_hist*. We would like this function to calculates an histogram based in the data array *x* by dividing the range of *x*s into *n_bins* equal bins.
#
# Complete the code so that it will correctly calculate the histogram:

def calculate_hist(x, n_bins):
    n_samples = len(x)  # Get the number of values in x
    
    ## Define bins
    bins_edges = np.linspace(x.min(), x.max(), n_bins + 1)
    bins_width = np.diff(bins_edges) ## This is equal to running: bins[1:] - bins[:-1]

    ## Initialize the array of histogram values
    hist_values = np.zeros(n_bins)
    
    ## Loop over the bins
    for k in range(n_bins):
        bin_left_edge = bins_edges[k]
        bin_right_edge = bins_edges[k + 1]
      
        ## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
        ## Fill in the following line to count the number of samples in each bin:
        number_of_samples_in_bin = ...
        ## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%
        
        hist_values[k] = number_of_samples_in_bin / n_samples / bins_width[k]
    
    return hist_values, bins_edges, bins_width


# The following code uses the function above to calculate the histogram on the list of runners' running times.

# +
## Define the array of the results:
results = dataset['Result'].values

## Calculate the histogram
n_bins = 150
hist_values, bins_edges, bins_width = calculate_hist(results, n_bins)
# -

# The following code uses the Matplotlib package to plot the histogram.

# +
## Prepare the figure and axes
fig, ax = plt.subplots()

## Plot the histogram's bars
ax.bar(bins_edges[:-1], hist_values, width=bins_width, align='edge')

## Set title and axis labels
ax.set_title('Histogram of Results')
ax.set_ylabel('PDF')
ax.set_xlabel('Result [hours]');
# -

# ### The hist functions in NumPy and Matplotlib
# From here on we will use the NumPy function [numpy.histogram](https://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram.html) to calculate histograms and the Matplotlib function [plt.hist](https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.hist.html?highlight=hist#matplotlib.axes.Axes.hist) to calculate and plot the histograms:

## Calculate and plot the histogram
fig, ax = plt.subplots()
ax.hist(results, bins=n_bins ,density=True)
ax.set_title('Historgram of Results')
ax.set_ylabel('PDF')
ax.set_xlabel('Result [hours]');

# ### ✍️  Indexing
#
# Like in Matlab, we can use an array of booleans *x* to index an array *y* by *y[x]*.
#
# Complete the following code to calculate the mean result of male runners from 2014:

# +
## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
## Fill in the following line so that indices will be a boolean array corresponding to male runners which competed in 2014
## For exmaple,  woman runners since 2012 would be: (years >= 2012) & (genders == 1)
indices = ...
## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%

result_of_male_in_2014 = results[indices]
mean_result = np.mean(result_of_male_in_2014)
print('The mean result of male runners in 2014 is {:.2f} hours'.format(mean_result))
# -

# *In the last line we have used the format function to format the output string. From here on we will be using format for that. You can read about how to use the format function [here](https://pyformat.info/)*

# ### ✍️  Comparing Two Age Groups
#
# Let us now compare the results of two groups of runers:
#
# - **Group 1**: Male runners from age group 1 in 2018
# - **Group 2**: Male runners from age group 3 in 2018
#
# Complete the code bellow to extract the results of these two groups and plots their histograms:

# +
## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%\
## Define the array of age groups (similar to how we have defined years and genders):
age_group = ...

## Define the indices for each group
indices1 = ...
indices2 = ...
## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%

## Extract the results for each group
results1 = results[indices1]
results2 = results[indices2]

## Calculate and plot the histograms
fig, ax = plt.subplots()
ax.hist(results1, bins=40 ,density=True, alpha=0.3, label='Group 1')
ax.hist(results2, bins=40 ,density=True, alpha=0.3, label='Group 2')
ax.set_title('Historgrams for Groups 1 & 2')
ax.set_ylabel('PDF')
ax.set_xlabel('Result [hours]');
ax.legend();
# -

# ### Fitting a parametric distribution model
#
# We will now try to fit a parametric distribution model to the results of the two age groups. We will do so using use the [Beta distribution](https://en.wikipedia.org/wiki/Beta_distribution).
#
# The Beta distribution describes a continues distribution over the finite interval of $\left[0,1\right]$ and has the following PDF:
#
# $$
# p\left(x;a,b\right)=\frac{\Gamma\left(a+b\right)}{\Gamma\left(a\right)\Gamma\left(b\right)}x^{a-1}\left(1-x\right)^{b-1}
# $$
#
# $\Gamma$ here is the [Gamma function](https://en.wikipedia.org/wiki/Gamma_function) and is usually referred to as the extension of the factorial fuction. (since for integers inputs we get that $\Gamma\left(n\right)=\left(n-1\right)!$). The Gamma function appears in the PDFs of various known distributions and is defined as an integral over some integrand. For our needs we will only have to be able to evaluate it for different inputs, and we will do so using SciPy function [scipy.special.gamma](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gamma.html).
#
# $a$ and $b$ are the two parameters of the distribution, which we would like to find.
#
# #### Scaling the distribution
#
# In our case, the results are distributed over the interval $\left[0,47\right]$, so we will have to scale the Beta distribution to fit this range. We can do so by using the following scaled PDF:
#
# $$
# p\left(x;a,b,s\right)=\frac{\Gamma\left(a+b\right)}{s\cdot\Gamma\left(a\right)\Gamma\left(b\right)}\left(x/s\right)^{a-1}\left(1-\left(x/s\right)\right)^{b-1}
# $$
#
# In our case $s=47$.

# ### ✍️  Plotting the Beta Distribution
#
# Fill in the code below to define the function *calc_beta_pdf*, which calculates the Beta distribution's PDF given some $a$ and $b$.
#
# - Use the gamma(x) to calculate $\Gamma\left(x\right)$

# +
## import the gamma function from the SciPy package
from scipy.special import gamma

## Define the function calculating the PDF of the Beta distribution
def calc_beta_pdf(x, a, b, s):
  ## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
  pdf_value = gamma(a + b) / ....
  ## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%
  
  return pdf_value


# -

# The following code plots the PDF for the following values of values of $\left(a,b\right)$:
#
# $$
# \left(0.5,0.5\right),\left(5,1\right),\left(1,3\right),\left(2,2\right),\left(2,5\right),
# $$

# +
## Define the scale parameter s
s = 47

## Define a grid to plot the distribution over.
results_grid = np.arange(0, 47.1, 0.1)  # A grid from 0 to 47 with steps of 0.1

## Prepare the figure and axes
fig, ax = plt.subplots()

## Calculate and plot the PDF for eac hset of parameters
for a, b in [[0.5, 0.5], [5, 1], [1, 3], [2, 2], [2, 5]]:
  beta_pdf = calc_beta_pdf(results_grid, a, b, s)
  ax.plot(results_grid, beta_pdf, label='$a={},b={}$'.format(a, b))

ax.set_title('Beta Distributions')
ax.set_ylabel('PDF')
ax.set_xlabel('Result [hours]');
ax.set_ylim(0, 2.5 / s)
ax.legend();


# -

# ### ✍️  Maximum Likelihood Estimation (MLE)
#
# We will Use the MLE approach to find the optimal parameters $a^*$ and $b^*$ for fitting the beta distribution to the results of group 1 and group2.
#
# To find the optimal parameters which minimizes the *Argmin* use a grid search, i.e., look for the optimal parameters over a selected grid by checking each and every point on the grid. Fill in the code below the define the functions performing the MLE.
#
# - Use np.log(x) to calculate the log of x
# - Use the calc_beta_pdf we have defined above to calculate the log-likelihood
# - Reminder: Use np.sum(x) to calculate the sum over the array x

# +
# Define the log-likelihood function
def beta_log_likelyhood(a, b, s, x):
  ## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
  ## Fill in this line to calculate the log-likelihood. x here is the array contating the data.
  log_likelyhood = 
  ## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%
  
  return log_likelyhood

  
# Define the function which searches for the optimal parameters a & b on a given grid
def beta_mle(a_grid, b_grid, s, x):
  optimal_a = None
  optimal_b = None
  min_value_so_far = np.inf
  for a in a_grid:
    for b in b_grid:
      ## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
      ## Note that this loop searches for the minimal value assigned to current_value.
      ## Write your optimization problem as an argmin expression (and not as an argmax).
      current_value = ...
      ## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%
      if current_value < min_value_so_far:
        optimal_a = a
        optimal_b = b
        min_value_so_far = current_value
  
  return optimal_a, optimal_b


# -

# The following code uses the above functions to find the optimal parameters and plots the estimated PDF over the histogram

# +
## Define the grid
a_grid = np.arange(0.1, 10, 0.1)
b_grid = np.arange(0.1, 10, 0.1)

## Group 1
## =======
## Find optimal MLE parameters
optimal_a1, optimal_b1 = beta_mle(a_grid, b_grid, s, results1)
## Calcualte the PDF
beta_pdf1 = calc_beta_pdf(results_grid, optimal_a1, optimal_b1, s)
## Plot the results
fig, ax = plt.subplots()
ax.hist(results1, bins=40 ,density=True, label='Histogram')
ax.plot(results_grid, beta_pdf1, label='Beta')
ax.set_title('Group 1 - MLE')
ax.set_ylabel('PDF')
ax.set_xlabel('Time [min]')
ax.legend();

## Group 2
## =======
## Find optimal MLE parameters
optimal_a2, optimal_b2 = beta_mle(a_grid, b_grid, s, results2)
## Calcualte the PDF
beta_pdf2 = calc_beta_pdf(results_grid, optimal_a2, optimal_b2, s)
## Plot the results
fig, ax = plt.subplots()
ax.hist(results2, bins=40 ,density=True, label='Histogram')
ax.plot(results_grid, beta_pdf2, label='Beta')
ax.set_title('Group 2 - MLE')
ax.set_ylabel('PDF')
ax.set_xlabel('Time [min]')
ax.legend();
# -

# ### The [scipy.stats.beta](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.beta.html) model
#
# SciPy has a set of object for working with a large range of distributions, one of them is *scipy.stats.beta*.
# A distribution with a given set of parameter can by define by:
# ``` python
# distrib = beta(a, b)
# ```
#
# or with a specific scale and offset:
# ``` python
# distrib = beta(a, b, offset, scale)
# ```
#
# The *distrib* objects can now be used to calculate the PDF using *distrib.pdf*, generate random samples using *distrib.rsv* and more.
# The *scipy.stats.beta* can also be used to estimate the MLE parameters given a dataset by:
# ``` python
# a,  b, offset, scale = beta.fit(data, floc=offset, fscale=scale)
# ```
#
# The following code compares our grid search calculation of the MLE parameters to SciPy's calculation.

# +
## Import the beta distribution model from SciPy
from scipy.stats import beta

## Calcualte the MLE parameters for group 1 using SciPy
optimal_a1_scipy, optimal_b1_scipy, _, _ = beta.fit(results1, floc=0, fscale=s)
## The _ in the line above tells Python to ignore the 3rd and 4th outputs of beta.fit

print('Results using our function: a={:.2f}, b={:.2f}'.format(optimal_a1, optimal_b1))
print('Results using SciPy: a={:.2f}, b={:.2f}'.format(optimal_a1_scipy, optimal_b1_scipy))
# -

# ## The Probability of Group 2 Being Better Then Group 1
#
# We would now like to calculate the following probability:
#
# > The probability that the result of a random runner from group 2 will be better (lower) then the result of a random runner from group 1
#
# We would like to do so based on the two estimated probabilities we have found.
#
# Before you start, write this probability (for yourself, not for submission) as an integral over the PDF of the joint distribution of the two results. We will calculate this integral numerically.

# ### The 2D grid
#
# The following code creates a 2D grid using np.meshgrid (which is similar to Matlab's meshgrid function), and print the upper left corner of the two 2D arrays defining values on the grid.

# +
## Define the 2D grid of results using np.meshgrid (which is similar to Matlab's meshgrid function)
results_grid_2d_1, results_grid_2d_2 = np.meshgrid(np.arange(0, 47.1, 0.1), np.arange(0, 47.1, 0.1))

print('results_grid_2d_1:')
print(results_grid_2d_1[:5, :5])
print()
print('results_grid_2d_2:')
print(results_grid_2d_2[:5, :5])
# -

# ### ✍️  The joint probability
#
# Calculate the PDF of the joint probability of the results of the two runners over the 2D grid of results.

# +
## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
## Calcualte the MLE parameters for group 2 using SciPy
optimal_a2_scipy, optimal_b2_scipy, _, _ = ...

## Define two SciPy beta distribution objects based on the MLE parameters calculated earlier: 
beta_dist1 = beta(optimal_a1_scipy, optimal_b1_scipy, loc=0, scale=s)
beta_dist2 = ...

## Calculate the marginal PDF of the two results over the grid
marginal_pdf1 = beta_dist1.pdf(results_grid_2d_1)
marginal_pdf2 = ...

## Calculate the joint PDF over the grid
joint_pdf = ...
## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%
# -

# ### Plotting the joint distribution

## Plot the joint PDF
from matplotlib import ticker  # ticker is used to select the levels of the conturs map according to a logaritmic scale
fig, ax = plt.subplots()
ax.contour(results_grid_2d_1, results_grid_2d_2, joint_pdf, locator=ticker.LogLocator(numticks=40), colors='black', linewidths=0.2)
img_obj = ax.contourf(results_grid_2d_1, results_grid_2d_2, joint_pdf, locator=ticker.LogLocator(numticks=40), cmap='jet')
plt.colorbar(img_obj);
ax.set_title('The Joint PDF')
ax.set_ylabel('Result of Runner 2')
ax.set_xlabel('Result of Runner 1')
ax.legend();

# ### Calculation the integral
#
# Given a 2D grid $\left\{x_i\right\}$, $\left\{y_j\right\}$ with step sizes of $\Delta x$ and $\Delta y$ and area $S$, we can approximate the integral:
#
# $$
# \underset{S}{\int\int}f\left(x,y\right)dx_1dx_2
# $$
#
# As:
#
# $$
# \sum_{x_i,y_j\in S}f\left(x_i,y_j\right) \Delta x \Delta y
# $$
#
# ✍️ Approximate the desired probability using a numerical approximation of the relevant integral over the joint probability:

# +
dx1 = 0.1
dx2 = 0.1

## %%%%%%%%%%%%%%% Your code here - Begin %%%%%%%%%%%%%%%
## Fill in the following line to define the points on the grid which we would like to integrate over.
## Use results_grid_2d_1 & results_grid_2d_2 to define relevat_area
relevat_area = ...
## %%%%%%%%%%%%%%% Your code here - End %%%%%%%%%%%%%%%%%

prob = np.sum(joint_pdf * relevat_area) * dx1 * dx2

print('The probability that the result of a random runner from group 2 will be better (lower) then the result of a random runner from group 1 is: {:.2f}'.format(prob))
# -

# ## Submission
#
# To submit your code download it as a **ipynb** file from Colab, and upload it to the course's website. You can download this code by selecting **Download .ipynb** from the **file** menu.

StarcoderdataJupyterScriptsDedupFiltered