{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 02402 Week 2\n",
    "\n",
    "Welcome to week 2 of 02402 Statistics (PF)\n",
    "\n",
    "Today we will start using Python for statistics and simulation. \n",
    "\n",
    "We will also start using the libraries: Numpy, Matplotlib and Pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Descriptive statistics in Python\n",
    "\n",
    "We start with some descriptive statistics - calculating the mean, average etc from a data sample. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to be able to work with a data type that behaves as a vector. \n",
    "\n",
    "For this we use Numpy arrays from the *Numpy* library. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the Numpy library\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now work with a sample, consisting of 10 measurements of students heights. \n",
    "\n",
    "The 10 observations have the values: 168, 161, 167, 179, 184, 166, 198, 187, 191 and 182"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# store sample data in variable x:\n",
    "x = np.array([168, 161, 167, 179, 184, 166, 198, 187, 191, 182])\n",
    "\n",
    "# print out the values:\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute Statistics (DK: Nøgletal): mean, variance, quantiles and more:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate mean of x (average height of students)\n",
    "np.mean(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# \"mean()\" can also be called as a \"method\" (for the ndarray)\n",
    "x.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Have a look in the online documentation for numpy arrays: https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html\n",
    "\n",
    "The datatype \"ndarray\" (also called a numpy array) has many methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets try some other \"methods\"\n",
    "x.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what about variance? \n",
    "\n",
    "# OBS: we need to remember ddof = 1 in order to calculate the \"sample variance\"\n",
    "\n",
    "x.var(ddof=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why ddof=1? have a look in the documentation for explanation: https://numpy.org/doc/stable/reference/generated/numpy.var.html#numpy.var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# standard deviation (also remember ddof=1 for \"sample standard deviation\")\n",
    "x.std(ddof=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what about the median?\n",
    "x.median()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "theres an error!\n",
    "\n",
    "This is because there is no \"method\" (for numpy arrays) called median. \n",
    "\n",
    "OK - then we call the median() function directly from numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.median(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we could also call .mean(), .min(), .var() etc. directly from numpy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.mean(x))\n",
    "print(np.min(x))\n",
    "print(np.var(x, ddof=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can also get other quantiles (0.50-quantile is the same as the median)\n",
    "np.quantile(x, [0,0.10,0.25,0.50,0.75,0.90,1.00], method='averaged_inverted_cdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " # Python has two equivalent funtions for calculating quantiles: \"percentile\" and \"quantile\"\n",
    "np.percentile(x, [0,10,25,50,75,90,100], method='averaged_inverted_cdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare with sorted data\n",
    "sorted_x = np.sort(x)\n",
    "print(sorted_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the method=\"averaged_inverted_cdf\"  !<br>\n",
    "\n",
    "See the documentaion: https://numpy.org/doc/stable/reference/generated/numpy.percentile.html#numpy.percentile\n",
    "\n",
    "There are many different ways to define percentiles - the definition we use in this book is the same as \"averaged_inverted_cdf\". \n",
    "\n",
    "(we will come back to this later) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize the data with plots:\n",
    "\n",
    "We use the *Matplotlib* library to produce plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the matplotlib.pyplot package \n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Distribution and histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now make a histogram of the sample data\n",
    "plt.hist(x)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Customize your histogram\n",
    "\n",
    "# Make a \"density plot\"\n",
    "plt.hist(x, bins=8, edgecolor='black', color='red', density=True)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('Now the y-axis shows density instead of counts')\n",
    "plt.title('Histogram Example')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# specifying bin-edges:\n",
    "plt.hist(x, bins=[160,165,170,175,180,185,190,195,200], edgecolor='black', color='red', density=True)\n",
    "plt.xlabel('x - now grouped into user-specified bins')\n",
    "plt.ylabel('Density')\n",
    "plt.title('Histogram Example')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histograms are **important** - they show how the data is **distributed** and are often the first choice of visualising a sample <br>\n",
    "\n",
    "Histograms serve as *empirical distributions* (\"empirical pdf\")<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Cumulative distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the \"empirical cumulated density function\" (empirical cdf)\n",
    "plt.ecdf(x)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare with values \n",
    "print(sorted_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can check that the vertical steps in the cumulated distribution plot corresponds to the exact values in the data sample. \n",
    "\n",
    "In the cumulated distribution all detailed information is kept - it is another way to visualise the distribution of data (without having to choose a specific bin-size). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets visualize the \"averaged_inverted_cdf\" method (how we find quantiles):\n",
    "\n",
    "q1, q2, q3 = np.quantile(x, [.25, .50, .75], method='averaged_inverted_cdf')\n",
    "plt.ecdf(x)\n",
    "\n",
    "# (Don't worry to much about how this plot works)\n",
    "for q, p, color in zip(\n",
    "    [q1, q2, q3],\n",
    "    [0.25, 0.50, 0.75],\n",
    "    ['tab:green', 'tab:orange', 'tab:red']\n",
    "):\n",
    "    # horizontal line from y-axis to ECDF\n",
    "    plt.plot([x.min(), q], [p, p], '--', color=color)\n",
    "\n",
    "    # vertical line down to x-axis\n",
    "    plt.plot([q, q], [0, p], '--', color=color)\n",
    "\n",
    "    # mark intersection\n",
    "    plt.plot(q, p, 'o', color=color)\n",
    "\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# printing the 0.90-quantile:\n",
    "print(np.quantile(x, [.90], method='averaged_inverted_cdf'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Boxplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a boxplot\n",
    "plt.boxplot(x)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "see documentation for definition of box and whiskers: \n",
    "\n",
    "https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.boxplot.html#matplotlib.axes.Axes.boxplot\n",
    "\n",
    "Also study example 1.29 in the book, to understand how whiskers can be plotted in different ways. <br>\n",
    "OBS: the \"modified boxplot\" is the default in Python. The \"basic boxplot\" needs the extra input \"whis = (0,100)\". But since the example shown here (in the plot above) has no \"outliers\" it does not make a difference here. In example 1.29 an outlier is added to the data in order to visualise the difference. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Scatter plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We store 10 measurements of **student weights** in variable y:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.array([65.5, 58.3, 68.1, 85.7, 80.5, 63.4, 102.6, 91.4, 86.7, 78.9])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a scatter plot of x (student heights) and y (student weights):\n",
    "plt.scatter(x,y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does the plot tell you?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate covariance:\n",
    "np.cov(x,y, ddof=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the four values?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print only the off diagonal (top right) element:\n",
    "np.cov(x,y, ddof=1)[0,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate correlation\n",
    "np.corrcoef(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are these four values?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do you interpret a correlation of 0.9656 ? Does the number match your expectation from the scatter plot?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DataFrames - Working with real data\n",
    "\n",
    "**For real datasets it is very often relevant to read the data from a file** - e.g. a csv file. \n",
    "\n",
    "We will try \"importing\" a dataset here. You will also do this when you start working with the *project exercises*. \n",
    "\n",
    "In Python we store tabular data in a **\"DataFrame\"** for which we need the *Pandas* library. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the Pandas library\n",
    "import pandas as pd "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read data from a csv file:\n",
    "csv_data = pd.read_csv(\"studentheights.csv\", sep=';')\n",
    "\n",
    "# OBS: for this code to work you need to have the file \"studentheights.csv\" in the SAME FOLDER where \n",
    "# you keep this jupyter notebook !\n",
    "# If the datafile is stored somewhere else, you need to specify the file path, e.g.:\n",
    "# csv_data = pd.read_csv(\"C:/folder/folder/studentheights.csv\", sep=';')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view the first few rows:\n",
    "csv_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print the number of rows in the dataset:\n",
    "print(len(csv_data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try opening the file \"studentheights.csv\" in another program and compare the ourput from Python with what you see in the file. Does the number of rows match?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# try calculating the mean of the heights. You will need to specify which column you want to use:\n",
    "csv_data[\"Height\"].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try also computing other statistics than the mean. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To get an overview of the data in a DataFrame use the \".describe\" method:\n",
    "csv_data.describe(include='all')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check what happens if you do not include the \"include=all\" input. (also consider checking the documentation or ask an AI what the difference is)\n",
    "\n",
    "Add the \"include=all\" and compute the output again:\n",
    "Do you understand *why* the information for the two columns is not the same? What is shown for numerical variables and what is shown for catagorical variables?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The DataFrame has a direct method for making histograms:\n",
    "csv_data.hist(column=\"Height\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The DataFrame has a direct method for making a boxplot:\n",
    "csv_data.boxplot(column=\"Height\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to do a boxplot by gender, we need to include the \"by=..\" argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "csv_data.boxplot(column=\"Height\", by='Gender')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we could also use the \"by=..\" argument for making separate histograms:\n",
    "csv_data.hist(column=\"Height\", by=\"Gender\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**DataFrames can also be defined directly**\n",
    "\n",
    "It *is* possible to type values directly into a *DataFrame*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# typing data directly into a DataFrame (here we make a new DataFrame with three columns):\n",
    "student_data = pd.DataFrame({\n",
    "    'height':  [168, 161, 167, 179, 184, 166, 198, 187, 191, 179],\n",
    "    'weight':  [65.5, 58.3, 68.1, 85.7, 80.5, 63.4, 102.6, 91.4, 86.7, 78.9],\n",
    "    'gender':  ['female', 'female', 'female', 'male', 'male', 'female', 'male', 'male', 'male', 'male']\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the code above and make sure you understand what is going on! This will be important for understanding many examples in the book and also in the exam exercises. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# virew the top 5 rows:\n",
    "student_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is good practice to always have one *observational unit* (e.g. one person) in each row and different *observational variables* (e.g. different information about that person) in the different columns (and not the other way around). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Note on *Libraries*\n",
    "\n",
    "Confused about when to use which libraries?\n",
    "\n",
    "We recommend loading the full list as used in the course and in the book:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Put this in the top of every notebook to make sure you have everything you need:\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.stats.power as smp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OBS: These aliases are also used in the exam - i.e. when we write *np.something*, we asume you know that np is an alias for *numpy*. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: First simulation in Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(8366)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Simmulating a discrete random variable (stochastic variable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the stochastic variabel X with following f(x) (pdf):\n",
    "\n",
    "| X: |0|1|2|3|\n",
    "|:--:|:--:|:--:|:--:|:--:|\n",
    "| f(x): | 0.1 | 0.3 | 0.4 | 0.2 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulate ONE realization of the random variable:\n",
    "result = np.random.choice(a=[0,1,2,3], size=1, replace=True, p=[0.1,0.3,0.4,0.2])  \n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat the code in the cell above. \n",
    "\n",
    "Does the value change? \n",
    "\n",
    "Is the simulation behaving as you expect?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can also \"draw\" many obervation in one go:\n",
    "sample = np.random.choice(a=[0,1,2,3], size=5000, replace=True, p=[0.1,0.3,0.4,0.2])\n",
    "print(sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do these values look as expected?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets make a histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(sample, bins=[-0.5,0.5,1.5,2.5,3.5], edgecolor='black')\n",
    "plt.xticks([0,1,2,3])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try re-making the sample a few times (and re-plotting the histogram).\n",
    "\n",
    "Does the histogram fit with the simulated probabilities?\n",
    "\n",
    "What happens if you make a larger sample (e.g., 1000 observations)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We can compare simulations to theoretical expectations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# E[X]:\n",
    "print(0*0.1 + 1*0.3 + 2*0.4 + 3*0.2)\n",
    "print(sample.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# V[X]:\n",
    "print( (0-1.7)**2*0.1 + (1-1.7)**2*0.3 + (2-1.7)**2*0.4 + (3-1.7)**2*0.2)\n",
    "print(sample.var(ddof=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate rolling a die"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the stochastic variabel X describing rolling a (fair) die. \n",
    "\n",
    "We *assume* that the die is fair and that there is an even probability of getting each possible outcome. \n",
    "\n",
    "The f(x) (pdf) is therefore:\n",
    "\n",
    "| X: |1|2|3|4|5|6|\n",
    "|:--:|:--:|:--:|:--:|:--:|:--:|:--:|\n",
    "| f(x): | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a single roll (try repeating a few times):\n",
    "dice_1 = np.random.choice([1,2,3,4,5,6], size=1) # we do not need to specify p if all probabilities are equal\n",
    "print(dice_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 10 rolls (try repeating a few times):\n",
    "dice_10 = np.random.choice([1,2,3,4,5,6], size=10) \n",
    "print(dice_10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# count the number of each outcome:\n",
    "print(np.bincount(dice_10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What do these counts represent? Can you make sense of them (compare with the values in dice_10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare the counts with the histogram:\n",
    "plt.hist(dice_10, bins=[-0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5], edgecolor='black')\n",
    "plt.xticks([0,1,2,3,4,5,6,7])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Look at Example 2.15 in the book"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of realizations\n",
    "n = 30\n",
    "# Simulate rolls with a fair dice\n",
    "xFair = np.random.choice(range(1, 7), size=n, replace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "what is going on?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets print a few things to understand what has happened:\n",
    "print(range(1, 7))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The range(1,7) just outputs \"range(1,7)\" (???) <br>\n",
    "We need to do something else to print the actual values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = range(1,7)\n",
    "for n in x:\n",
    "  print(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output are values 1,2,3,4,5,6 - like the values on a dice, so this is good! <br>\n",
    "(you can read about the range function here: https://www.w3schools.com/python/ref_func_range.asp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets print the simulated values:\n",
    "print(xFair)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do these values look fine?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also notice how certain details are left out of \"np.random.choice(range(1, 7), size=n, replace=True)\" <br>\n",
    "We do not need to write \"a = \" <br>\n",
    "We do not need to state \"p = \" <br>\n",
    "\n",
    "Read more about this in the documentation: https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html <br>\n",
    "(consider asking an AI to explain it to you, if you do not understand)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now work through the rest of example 2.15:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Probability distributions and scipy.stats "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import scipy.stats library\n",
    "import scipy.stats as stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example with the binomial distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.binom.pmf(k=6, n=6, p=0.70)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare with formula (for x=n=6 the binomial coefficient is just 1):\n",
    "1 * 0.70**6 * (1-0.70)**0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compute the probability of every possible outcome (P(X = 0), P(X = 1), P(X = 2), P(X = 3), P(X = 4), P(X = 5) and P(X = 6)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(stats.binom.pmf(k=[0,1,2,3,4,5,6], n=6, p=0.70))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can also visualise these probabilities (plot the pdf):\n",
    "plt.bar([0,1,2,3,4,5,6], stats.binom.pmf(k=[0,1,2,3,4,5,6], n=6, p=0.70), width=0.1, color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute cdf, inverse cdf, mean, variance, etc\n",
    "\n",
    "Python has many other methods for every distribution. \n",
    "\n",
    "We can also compute the cdf (.cdf), the inverse cdf (.ppf), the Expectation value/the mean (.mean), the variance (.var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the cdf; P(X <= x):\n",
    "print(stats.binom.cdf(k=[0,1,2,3,4,5,6], n=6, p=0.70))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualise the cdf:\n",
    "plt.bar([0,1,2,3,4,5,6], stats.binom.cdf(k=[0,1,2,3,4,5,6], n=6, p=0.70), width=0.1, color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the inverse cdf, in python called \".ppf\" for percent point function. \n",
    "# For instance we can compute the quartiles:\n",
    "stats.binom.ppf(q=[0.25, 0.50, 0.75], n=6, p=0.70)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you visually verify these values from the cdf plot above?\n",
    "\n",
    "hints:<br>\n",
    "find 0.25 on the y-axis and then go to corresponding x-value - this should be Q1<br>\n",
    "find 0.50 on the y-axis and then go to corresponding x-value - this should be Q2 = the median<br>\n",
    "find 0.75 on the y-axis and then go to corresponding x-value - this should be Q3 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the expectation value / the mean:\n",
    "stats.binom.mean(n=6, p=0.70)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the variance:\n",
    "stats.binom.var(n=6, p=0.70)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the standard deviation:\n",
    "stats.binom.std(n=6, p=0.70)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the median:\n",
    "stats.binom.median(n=6, p=0.70)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating random variates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can simulate a random variate - that is a single observation of the random variable - using .rvs:\n",
    "print(stats.binom.rvs(size=1, n=6, p=0.70))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try repreating the code above a few times. \n",
    "\n",
    "What are we simulating?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can also simulate many obersavtions in one go:\n",
    "print(stats.binom.rvs(size=100, n=6, p=0.70))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pernille",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}