notebooks/08_module_numpy_and_scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 8 - NumPy and SciPy: introduction to statistics in python <a id='0'></a>\n",
    "--------------------------------------------------------------------------------------------------\n",
    "\n",
    "## Table of Content <a id='toc'></a>\n",
    "\n",
    "<br>\n",
    "\n",
    "[**NumPy - The fundamental package for scientific computing in Python**](#0)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**Introduction**](#1)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Importing NumPy](#1.1)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**The heart of NumPy: the array**](#2)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**Creating NumPy arrays**](#3)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Creating arrays from lists](#3.1)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Creating homogeneous-filled arrays](#3.2)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Create an identity matrix](#3.3)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Array with random numbers following a uniform distribution between 0.0 and 1.0](#3.4)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Evenly spaced arrays: the np.arange function](#3.5)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Reading and writing arrays from/to files](#3.6)  \n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**Array reshaping**](#4)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**Accessing data in NumPy arrays**](#5)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Indexing](#5.1)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Slicing](#5.2)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Comparison operations](#5.3)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Iterating over arrays](#5.4)  \n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**Modifying array values**](#6)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Assigning new values to an existing index](#6.1)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Adding rows and columns to an existing array](#6.2)  \n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**To copy or not to copy?**](#7)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**NumPy functions**](#8)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**Random numbers in NumPy**](#9)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**Linear algebra built-in capabilities**](#10)\n",
    "\n",
    "<br>\n",
    "\n",
    "[**SciPy.stats and statistics in python**](#11)\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**Manipulation of random distributions**](#12)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Drawing random numbers: rvs](#12.1)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Looking up the quantiles and probability density functions](#12.2)  \n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;[**Statistical tests**](#13)  \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[Statistical modelling and regression](#13.1)\n",
    "\n",
    "<br>\n",
    "\n",
    "[**Exercises**](#14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "# NumPy - The fundamental package for scientific computing in Python <a id='0'></a>\n",
    "-------------\n",
    "\n",
    "\n",
    "## Introduction <a id='1'></a>\n",
    "-------------------\n",
    "\n",
    "**[NumPy](www.numpy.org)** is the fundamental package for scientific computing with Python.  \n",
    "Its **highlights** are:\n",
    "* A powerful N-dimensional array object.\n",
    "* Efficient, broadcasting functions.\n",
    "* Tools for integrating C/C++ and Fortran code.\n",
    "* Useful linear algebra, Fourier transform, and random number generation capabilities.\n",
    "\n",
    "**Why use NumPy ?**\n",
    "* **NumPy arrays are faster** than standard python lists.\n",
    "* NumPy functions allow to **perform many operations with much less code**.\n",
    "* NumPy functions are **faster than native Python implementation**.\n",
    "* Great **collection of mathematical functions** available (numpy, scipy, sympy).\n",
    "\n",
    "### Importing NumPy <a id='1.1'></a>\n",
    "\n",
    "`numpy` is generally imported under the alias **`np`**.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## The heart of NumPy: the array <a id='2'></a>\n",
    "------------------------------------\n",
    "\n",
    "* NumPy's main object is the homogeneous **multidimensional array** (all elements in the array must\n",
    "  of the same type).\n",
    "* Numpy arrays are created with **`np.array()`**.\n",
    "* While a NumPy array object might be reminiscent of a standard python `list`, NumPy arrays are much\n",
    "  more efficient for operations with large numerical data and in general outperform standard Python lists.\n",
    "  \n",
    "<br>\n",
    "\n",
    "**Example:** create a new numpy array, here of dimension 2 (rows) x 3 (columns)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array = np.array([[1.,2.,3.],[4.,5.,6.]])\n",
    "\n",
    "print(\"my_array:\\n\", my_array)\n",
    "print(\"my_array dimensions:\", my_array.shape)          # Number or rows and columns in the array.\n",
    "print(\"my_array number of elements:\", my_array.size)   # Total number of elements in the array.\n",
    "print(\"my_array type of elements\", my_array.dtype)     # Type of data stored in the array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "One great functionality of NumPy arrays is that it is painfully easy to **perform an operation over the entirety of the elements in the array**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_result = my_array * 3   # Multiply by 3 all elements in the array.\n",
    "print(my_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> To perform the same operation with **native python** would take a lot more effort and code:\n",
    ">\n",
    "> ```python\n",
    "> my_list = [[1.,2.,3.],[4.,5.,6.]]\n",
    "> my_result_2 = [[0,0,0],[0,0,0]]\n",
    "> for i in range(len(my_list)):\n",
    ">     for j in range(len(my_list[i])):\n",
    ">         my_result_2[i][j] = my_list[i][j] * 3\n",
    ">\n",
    "> print(my_result_2)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Numpy is also able to make the most out of the constraint of homogeneity in the array data to provide **amazing speed-ups**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from time import time\n",
    "\n",
    "native_data = tuple(range(10**7))\n",
    "numpy_data = np.array(native_data)\n",
    "\n",
    "t0 = time()\n",
    "numpy_data *= 3\n",
    "numpy_time = time() - t0\n",
    "\n",
    "t0 = time()\n",
    "native_data = [ x*3 for x in native_data ]\n",
    "native_time = time() - t0\n",
    "\n",
    "print(\"native timing :\", native_time)\n",
    "print(\"numpy timing :\", numpy_time)\n",
    "print(\"numpy acceleration factor:\" , native_time // numpy_time )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## Creating NumPy arrays <a id='3'></a>\n",
    "---------------------------------\n",
    "\n",
    "### Creating arrays from lists (and more generally from iterables) <a id='3.1'></a>\n",
    "\n",
    "The simplest way to create a new NumPy array is to use the **`np.array()`** method that accepts an iterable as argument, e.g. a `list`, `tuple` or `generator`.\n",
    "\n",
    "* **Create one-dimensional** arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "my_array = np.array([1, 2, 3])          # Create an array from a list...\n",
    "my_array_2 = np.array((7, 8, \"nine\"))   # ...from a tuple\n",
    "my_array_3 = np.array(range(10))        # ...from a generator (here a range object).\n",
    "\n",
    "print(\"object type:\", type(my_array))\n",
    "print(\"array value:\", my_array)\n",
    "print(\"array value:\", my_array_2)\n",
    "print(\"array value:\", my_array_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note how `np.array((7, 8, \"nine\"))` has converted the values of `7` and `8` to strings.\n",
    "> This is because all values in a NumPy array must be of the same type."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "* **Create N-dimensional** arrays: NumPy also allows creating 2-(or more)dimensional arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of a 2-dimensional array:\n",
    "my_array = np.array([[1,2,3],[4,5,6]])\n",
    "\n",
    "print(\"object type:\", type(my_array))\n",
    "print(\"array value:\\n\", my_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of a 3-dimensional array:\n",
    "my_array = np.array([[range(10), range(10)], [range(10), range(10)], [range(10), range(10)]])\n",
    "my_array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Creating homogeneous-filled arrays <a id='3.2'></a>\n",
    "\n",
    "The following functions generate arrays filled with the same value. The **`<shape>`** argument is a tuple/list (iterable) that indicates the dimensions of the array (e.g. rows and columns).\n",
    "* Array **filled with zeroes**: `np.zeros(<shape>)`.\n",
    "* Array **filled with ones**: `np.ones(<shape>)`.\n",
    "* Array **filled with a specific value**: `np.full(<shape>, value)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2-dimensional array filled with 0s.\n",
    "my_array = np.zeros((3, 5))\n",
    "print(my_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3-dimensional array filled with 1s.\n",
    "my_array = np.ones((2, 4, 2))\n",
    "print(my_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2-dimensional array filled with a specific value.\n",
    "my_array = np.full((2, 3), 42)\n",
    "print(my_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Create an identity matrix <a id='3.3'></a>\n",
    "\n",
    "* **`np.eye()`**: `np.eye(<number of rows>, <number of columns>, <position of diagonal>)`.\n",
    "* The 3rd argument, **`<position of diagonal>`** specifies where the diagonal with `1` values should be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array = np.eye(4,4,0)\n",
    "print(\"Diagonal at index position 0:\\n\", my_array)\n",
    "\n",
    "my_array = np.eye(4,4,1)\n",
    "print(\"\\nDiagonal at index position 1:\\n\", my_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Array with random numbers following a uniform distribution between 0.0 and 1.0  <a id='3.4'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array = np.random.rand(2,5)\n",
    "print(my_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Evenly spaced arrays: the `np.arange` function <a id='3.5'></a>\n",
    "\n",
    "**`np.arange()`** allows to generate one-dimensional arrays filled with **evenly spaced numbers**.\n",
    "* Accepts custom `start`, `stop`, and `step` (increment) values: `np.arange(start, stop, step)`.\n",
    "* By default, arrays are filled with integers starting at `0`.\n",
    "* The stop (end) value is *not* included (as usual in python).\n",
    "* Unlike the built-in `range()` function, `np.arange()` supports **float** `start`/`stop` points and float increments (`step`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array_1 = np.arange(6)             # Integers from 0 to 5 - the stop value is not included.\n",
    "my_array_2 = np.arange(6, 21, 3)      # Every 3rd number from 6 to 21, 21 not included.\n",
    "my_array_3 = np.arange(1.1, 7, 1.5)   # Supports float values! Here between 1.1 and 7, incremented by 1.5.\n",
    "\n",
    "print(my_array_1)\n",
    "print(my_array_2)\n",
    "print(my_array_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Reading and writing arrays from/to files <a id='3.6'></a>\n",
    "\n",
    "* **`numpy.loadtxt()`**: read data from a (tabulated) text file.\n",
    "\n",
    "> `numpy.loadtxt(<file name>, dtype=<type of data in file>, comments='#', delimiter=None, converters=None,\n",
    "   skiprows=0, usecols=None, unpack=False, ndmin=0)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Print the content of the file.\n",
    "# Note that we use \"!\" to call a shell command (here \"cat\"). This is a feature of Jupyter Notebooks,\n",
    "# but it will not work on Windows machines.\n",
    "!cat data/test.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: we skip the 1st row of the file, because it contains column names.\n",
    "my_array = np.loadtxt(\"data/test.tab\", dtype=float, delimiter=\" \", skiprows=1)\n",
    "print(my_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Notes:**\n",
    "> * There is a more sophisticated function for handling files with complex formatting : **`numpy.genfromtxt()`**.\n",
    "> * For complex cases, the use of other libraries or custom parsers is required."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "* **`numpy.savetxt()`**: write data to a text file.\n",
    "> `numpy.savetxt(fname, X, fmt='%.18e', delimiter=' ', newline='\\n', header='', footer='', comments='# ')`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array = np.array([[1,2,3], [4,5,6]])\n",
    "np.savetxt(\n",
    "    \"test_out.tab\", \n",
    "    my_array, \n",
    "    fmt=\"%d\", \n",
    "    delimiter=\"\\t\",\n",
    "    header=\"# my sample\\n#test\",\n",
    "    comments=\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the content of the file we just created (note: this command will not work on Windows).\n",
    "!cat test_out.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "### Array reshaping <a id='4'></a>\n",
    "\n",
    "* **`numpy.reshape()`**: changes the shape of an array without changing its data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1-dimensional array.\n",
    "my_array = np.arange(1.1, 6, 0.7)\n",
    "print(my_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reshape to a 2D array.\n",
    "print(np.reshape(my_array, (2, 4)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reshape to a 3D array.\n",
    "print(np.reshape(my_array, (2, 2, 2), order=\"C\"))  # order=\"C\" is the default value.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## Accessing data in NumPy arrays <a id='5'></a>\n",
    "------------------------------------------------\n",
    "\n",
    "### Indexing <a id='5.1'></a>\n",
    "\n",
    "* **One-dimensional** arrays are accessed in the same way as standard lists in python.\n",
    "* **Multi-dimensional** arrays are accessed by indices along every axis. These indices are separated by comas and given in square brackets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Create a new 3x3 Numpy array.\n",
    "my_array = np.arange(1., 10.).reshape((3, 3))\n",
    "print(my_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Accessing a single element in standard way (list-like syntax):\n",
    "print(my_array[1][0])\n",
    "\n",
    "# Accessing a row in the standard way:\n",
    "print(\"row:\", my_array[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Accessing an element using the NumPy syntax**: in NumPy indices are given in square brackets and separated by coma."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Accessing a single element in numpy syntax:\n",
    "print(my_array[1, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Slicing <a id='5.2'></a>\n",
    "\n",
    "* Accessing a subset of an array is done using the **`[rows to select, cols to select]`** syntax.\n",
    "* Passing **`:`** selects all rows or columns.\n",
    "* Remember that indexing is zero-based (i.e. the first row/column has index 0, not 1), and that\n",
    "  the end index of a slice is not included.\n",
    "\n",
    "<br>\n",
    "\n",
    "**Example 1:** accessing a single row or column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(my_array, \"\\n\")\n",
    "\n",
    "# Accessing a single row.\n",
    "print(\"Second row:\", my_array[1,:])\n",
    "\n",
    "# Accessing a single column.\n",
    "print(\"Third column:\", my_array[:, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "**Example 2:** accessing a 2-dimensional subset of an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All rows from the second row to the end, and all columns\n",
    "# from the first one to the third one (not included).\n",
    "print(my_array[1:,:2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# First and 3rd row, and all columns.\n",
    "print(my_array[[0,2],:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "* Similarly to slicing lists/tuples, the **`start:stop:step`** notation can also be used when accessing array elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array = np.arange(1., 21.).reshape((4,5))\n",
    "print(my_array)\n",
    "\n",
    "# Accessing a subset with a step argument.\n",
    "print(\"subset:\\n\", my_array[:,0::2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Comparison operations <a id='5.3'></a>\n",
    "\n",
    "When applying **comparison operators**, such as `>`, `<` or `==`, an array of boolean values of the same shape than the input array is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array = np.arange(1., 10.).reshape((3,3))\n",
    "\n",
    "print(\"my_array:\\n\", my_array, \"\\n\")\n",
    "print(\"my_array > 5:\\n\", my_array > 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The **`.all()`** and **`.any()`** methods of arrays can be used to evaluate whether all values in\n",
    "  a boolean array are `True` or `False` (in non-boolean arrays, these methods can also be used - they test\n",
    "  the *truthyness* of values in the array)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = my_array > 5\n",
    "\n",
    "print(\"All values in the array are greater than 5:\", results.all()) \n",
    "print(\"There is at least one value in the array greater than 5:\", results.any()) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Boolean arrays can also be used to **filter an array** (extract values) that match a given criteria."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A boolean array is used to select values >= 5.\n",
    "print(\"Values which are greater or equal to 5:\", my_array[my_array >= 5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Iterating over arrays <a id='5.4'></a>\n",
    "\n",
    "* When iterating over a multi-dimensional NumPy array with a `for` loop, the iteration is done over rows.\n",
    "* To iterate over **individual elements**, the **`np.nditer()`** can be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The standard for loop iterates over rows.\n",
    "for row in my_array:\n",
    "    print(row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Iterating over all elements in standard way.\n",
    "for x in my_array:\n",
    "    for y in x:\n",
    "        print(y, \" \",end=\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Iterating over all elements in numpy way.\n",
    "for x in np.nditer(my_array):\n",
    "    print(x, \" \", end=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## Modifying array values <a id='6'></a>\n",
    "--------------------------------\n",
    "\n",
    "### Assigning new values to an existing index <a id='6.1'></a>\n",
    "\n",
    "NumPy arrays are mutable, and their values can be modified.\n",
    "* To modify a value, simply assign a new value to an existing index (or range) position.\n",
    "\n",
    "**Example 1:** modify a single value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new array.\n",
    "my_array = np.arange(1,5).reshape((2,2))\n",
    "print(my_array, \"\\n\")\n",
    "\n",
    "# Modify the value in row 0, column 1:\n",
    "my_array[0,1] = 5\n",
    "print(\"After changing one element:\\n\", my_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "* **Example 2:** modify an entire row or column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new array.\n",
    "my_array = np.arange(1,5).reshape((2,2))\n",
    "print(my_array, \"\\n\")\n",
    "\n",
    "# Modify values in the 1st row.\n",
    "my_array[0,:] = np.array([5,6])\n",
    "print(\"After changing one row:\\n\", my_array, \"\\n\")\n",
    "\n",
    "# Modify values in the 2nd column.\n",
    "my_array[:,1] = np.array([7,8])\n",
    "print(\"After changing one column:\\n\", my_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Adding rows and columns to an existing array <a id='6.2'></a>\n",
    "\n",
    "* **`np.append()`**: append a row (`axis=0`) or column (`axis=1`) to an array.\n",
    "* **`np.insert()`**: add a row (`axis=0`) or column (`axis=1`) at a specific position in the array.\n",
    "* **`np.concatenate()`**: bind multiple arrays by rows or columns.\n",
    "\n",
    "<br>\n",
    "\n",
    "**Example 1:** add a row or column with **`np.append()`**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new array\n",
    "my_array = np.arange(0,6).reshape((2,3))\n",
    "print(my_array, \"\\n\")\n",
    "\n",
    "# Append column by setting \"axis=1\".\n",
    "print(\n",
    "    \"After appending a column:\\n\",\n",
    "    np.append(my_array, np.array([6, 7]).reshape(2,1), axis=1),\n",
    "    \"\\n\"\n",
    ")\n",
    "\n",
    "# Append a row by setting \"axis=0\".\n",
    "print(\n",
    "    \"After appending a row:\\n\",\n",
    "    np.append(my_array, np.array([6, 7, 8]).reshape(1,3), axis=0),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "**Example 2:** insert a row in an array with **`np.insert()`**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(my_array, \"\\n\")\n",
    "\n",
    "# Insert a row in 2nd position (index = 1).\n",
    "print(\"After row insertion\\n\", np.insert(my_array, 1, [[6, 7, 8]], axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "**Example 3:** combine arrays by row or column with **`np.concatenate()`**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "array_1 = np.arange(0,6).reshape((2,3))\n",
    "array_2 = np.arange(6,12).reshape((2,3))\n",
    "\n",
    "print(\"Row concatenation:\\n\", np.concatenate((array_1, array_2), axis=0), \"\\n\")\n",
    "print(\"Column concatenation:\\n\", np.concatenate((array_1, array_2), axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## To copy or not to copy? <a id='7'></a>\n",
    "--------------------------------\n",
    "\n",
    "When operating and manipulating arrays, their data is sometimes copied into a new array and sometimes not. This is often a source of confusion for beginners.\n",
    "* If you plan to change array values but want to keep the old array untouched, then make a copy of it\n",
    "  using the **`.copy()`** method!\n",
    "\n",
    "**Example:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array = np.arange(1., 21.).reshape((4,5))\n",
    "print(my_array, \"\\n\")\n",
    "\n",
    "# When assigning \"my_array\" to a new variable, the data from \"my_array\" is NOT copied.\n",
    "# Instead, \"tmp\" simply points to the same memory location as \"my_array\".\n",
    "tmp = my_array\n",
    "print(\"Are 'tmp' and 'my_array' pointing to the same object in memory?\", tmp is my_array, \"\\n\")\n",
    "\n",
    "# As a result, a change made to \"tmp\" will also be present in \"my_array\"...\n",
    "# because they are pointing to the same object!\n",
    "tmp[1, 1] = 999\n",
    "print(\"Values of 'my_array' after we changed 'tmp':\\n\", my_array) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example 2:** even if it appears that the `my_array` and `tmp` arrays are not the same object, they do in fact access the same memory location.\n",
    "\n",
    "> This behavior may appear strange, but it is because NumPy arrays access their data by reference.\n",
    "> And thus multiple array may access the same memory space, or subset of the same memory space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array = np.arange(1., 21.).reshape((4,5))\n",
    "print(my_array, \"\\n\")\n",
    "\n",
    "tmp = my_array[1:3,1:3]\n",
    "print(\"Are 'tmp' and 'my_array' pointing to the same object in memory?\", tmp is my_array, \"\\n\")\n",
    "\n",
    "# \"tmp\" is a different object than \"my_array\", but the change in one is still present in the other!\n",
    "tmp[1,1] = 999\n",
    "print(\"tmp:\\n\", tmp, \"\\n\")\n",
    "print(\"Content of 'my_array' after changing 'tmp':\\n\", my_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "**If you plan to change array values but want to keep the old array untouched, then make a copy of it!**\n",
    "* To copy an array, use its **`.copy()`** method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array = np.arange(1., 21.).reshape((4,5))\n",
    "print(my_array, \"\\n\")\n",
    "\n",
    "# Make a copy of the original array.\n",
    "tmp = my_array.copy()\n",
    "print(\"Are 'tmp' and 'my_array' pointing to the same object in memory?\", tmp is my_array, \"\\n\")\n",
    "\n",
    "# This time, modifying the copy does not modify the original array.\n",
    "tmp[1,1] = 999\n",
    "print(\"'tmp' array after the change:\\n\", tmp, \"\\n\")\n",
    "print(\"'my_array' after having changed 'tmp':\\n\", my_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## NumPy functions <a id='8'></a>\n",
    "-------------------------\n",
    "\n",
    "NumPy provide a great number of functions. The power of NumPy functions is their **speed** and the possibility to **apply them to arrays element-wise, column-wise or row-wise**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Element-wise operations and functions\n",
    "* When used on an array, all basic math operators (**`+`**, **`-`**, **`/`**, **`*`**, **`//`**, __`**`__) are automatically applied element-wise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Create a test array.\n",
    "my_array = np.arange(1, 7).reshape((2,3))    # We use .reshape to create a 2 x 3 matrix.\n",
    "print(\"Original values:\\n\", my_array)\n",
    "\n",
    "# Apply math operations to each element of the array.\n",
    "print(\"\\nDivision:\\n\", my_array / 2)\n",
    "print(\"\\nAddition:\\n\", my_array + 10)\n",
    "print(\"\\nPower of 2:\\n\", my_array ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**When arrays are of the same size**, the requested operation is carried-out element wise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Divide arrays:\\n\", my_array / my_array)\n",
    "print(\"\\nAddition arrays:\\n\", my_array + my_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "**When arrays differ in size**, the shorter one is recycled.  \n",
    "In this example, the array `[1, 2, 3]` is added to both the 1st and the 2nd row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_array + np.array([1, 2, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Element-wise, row-wise, column-wise functions\n",
    "\n",
    "NumPy has many functions that can be applied:\n",
    "* Element wise - by default.\n",
    "* Row-wise - by passing `axis=0` as argument.\n",
    "* Column-wise - by passing `axis=1` as argument.\n",
    "\n",
    "Some examples include: **`np.sum()`**, **`np.std()`**, **`np.mean()`** and **`np.log2()`**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sum.\n",
    "print(\"Element-wise sum:\", np.sum(my_array))\n",
    "print(\"Sum of columns:\", np.sum(my_array, axis=0))\n",
    "print(\"Sum of rows:\", np.sum(my_array, axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard deviation.\n",
    "print(\"Std of the whole array:\", np.std(my_array))\n",
    "print(\"Std of columns:\", np.std(my_array, axis=0))\n",
    "print(\"Std of rows:\", np.std(my_array, axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Mean function.\n",
    "print(\"Mean of the whole array:\", np.mean(my_array))\n",
    "print(\"Mean of columns:\", np.mean(my_array, axis=0))\n",
    "print(\"Mean of rows:\", np.mean(my_array, axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Log base 2.\n",
    "print(\"Compute log2 or each element:\\n\", np.log2(my_array))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## Random numbers in NumPy <a id='9'></a>\n",
    "---------------------------------------\n",
    "\n",
    "The **`numpy.random`** module provides a large collection of distributions (uniform, normal, beta, binomial, gamma, poisson ...) to draw from.\n",
    "\n",
    "* **`np.random.rand()`**: random numbers from a uniform distribution [0,1].\n",
    "* **`np.random.randn()`**: random numbers from a normal distribution.\n",
    "* **`np.random.permutation()`**: random permutation of the array's values (reshuffling).\n",
    "* **`np.random.choice(array, size, replace=True/False)`**: random sampling from an array, with (`replace=True`) or without (`replace=False`) replacement.\n",
    "\n",
    "<br>\n",
    "\n",
    "**Examples:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Random numbers from a uniform distribution [0,1].\n",
    "print(\"single number:\", np.random.rand())\n",
    "print(\"array of random numbers:\\n\", np.random.rand(2,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random numbers from normal distribution.\n",
    "my_array = np.random.randn(100000)\n",
    "print(\"array of random numbers.\\n\", \"mean :\", np.mean(my_array) , \"\\tstandard deviation :\" , np.std(my_array) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Permutation and sampling.\n",
    "my_array = np.arange(7)\n",
    "print(\"array permutation:\\n\", np.random.permutation(my_array))\n",
    "print(\"sample from the array:\\n\", np.random.choice(my_array, size=3, replace=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## Linear algebra built-in capabilities <a id='10'></a>\n",
    "----------------------------------------------\n",
    "\n",
    "NumPy arrays can be used as matrices without any special conversion.  \n",
    "For advanced linear algebra operations there is special package: **`numpy.linalg`**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new array.\n",
    "my_array = np.arange(1., 7.).reshape((2,3))\n",
    "print(my_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transpose of a matrix.\n",
    "print(\"transposed matrix:\\n\", my_array.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### There are different types of matrix multiplication!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Matrix multiplication element-wise.\n",
    "print(my_array * my_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Matrix product.\n",
    "print(my_array.dot(my_array.T)) \n",
    "\n",
    "# Alternatively, one can also use the \"@\" operator:\n",
    "print(my_array @ my_array.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "# SciPy.stats and statistics in python <a id='11'></a>\n",
    "------------------------------------------------------\n",
    "\n",
    "**SciPy** references a comprehensive [project for scientific python programming](https://scipy.org) as well as a [library](https://docs.scipy.org/doc/scipy/reference) implementing various tools and algorithm for scientific software.\n",
    "\n",
    "Here we will give a few pointers on the **`scipy.stats`** library, which provides ways to interact with various random distribution functions, as well as implement numerous statistical tests.\n",
    "\n",
    "<br>\n",
    "\n",
    "## Manipulation of random distributions <a id='12'></a>\n",
    "--------------------------------------------------\n",
    "\n",
    "The **`scipy.stats`** module implements utilities for a large number of continuous and discrete distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_continue = [d for d in dir(stats) if isinstance(getattr(stats, d), stats.rv_continuous)]\n",
    "dist_discrete = [d for d in dir(stats) if isinstance(getattr(stats, d), stats.rv_discrete)]\n",
    "\n",
    "print('number of continuous distributions: ' , len(dist_continue))\n",
    "print('number of discrete distributions:   ' , len(dist_discrete))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Let's experiment with the **normal distribution**, or **`norm`** in `scipy.stats`\n",
    "\n",
    "A look at `help(stats.norm)` tells us that:\n",
    "\n",
    "```\n",
    " |  The location (``loc``) keyword specifies the mean.\n",
    " |  The scale (``scale``) keyword specifies the standard deviation.\n",
    "```\n",
    "\n",
    "Let's generate a normal distribution with:\n",
    "* mean = 10 -> `loc=10`.\n",
    "* stdev = 2 -> `scale=2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = stats.norm(loc=10, scale=2)\n",
    "\n",
    "print(N.stats())  # The mean and variance of a distribution can be retrieved using the .stats method.\n",
    "print(type(N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "That distribution object (`N`) we just created can now be used to interact with the distribution in many ways, as will be illustrated below.\n",
    "\n",
    "<br>\n",
    "\n",
    "### Drawing random numbers: `.rvs()` <a id='12.1'></a>\n",
    "\n",
    "* **`.rvs()`**: this method of a scipy distribution is used to draw random numbers from the distribution.\n",
    "* The **`size`** argument defines the dimensions of the returned arrays of random numbers.\n",
    "  The type returned by the method is a NumPy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.rvs(size = [5,5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Note:** as with any drawing of random numbers on a computer,\n",
    "> [one merely emulates randomness](https://en.wikipedia.org/wiki/Pseudorandom_number_generator).\n",
    "> This also means that one can make some random operation reproducible by setting the random seed to\n",
    "  a specific value.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "np.random.seed(2020) # Set the seed of the random number generator to a specific value.\n",
    "draw1 = N.rvs(size=5)\n",
    "print(\"Random numbers:\", draw1)\n",
    "\n",
    "# By setting the random seed back to 2020, we can reproduce randomness.\n",
    "np.random.seed(2020)\n",
    "draw2 = N.rvs(size=5)\n",
    "print(\"Random numbers:\", draw2)\n",
    "\n",
    "print(\"\\nAre the random draws equal?\",draw1 == draw2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Looking up quantiles and probability density functions <a id='12.2'></a>\n",
    "\n",
    "* **`.pdf()`**: Probability Density Function of distribution.\n",
    "* **`.cdf()`**: Cumulative Distribution Function of distribution.\n",
    "* **`.ppf()`**: Percent Point Function (inverse of CDF): gives the quantiles of the distribution.\n",
    "\n",
    "<br>\n",
    "\n",
    "**Example:** plotting a distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt \n",
    "\n",
    "# pdf: Probability Density Function.\n",
    "x = np.arange(0, 20, 0.1)\n",
    "plt.plot(x, N.pdf(x))\n",
    "plt.show()\n",
    "\n",
    "# cdf: Cumulative Distribution Function.\n",
    "print(\"What is the probability of drawing a number <=15.0?\", N.cdf(15.0)) \n",
    "\n",
    "# ppf: Percent Point Function (inverse of CDF): gives the quantiles of the distribution.\n",
    "P = [0.025, 0.5, 0.975]\n",
    "for threshold, quantile in zip(P, N.ppf(P)):\n",
    "    print(\"Quantile:\", threshold, \"->\", quantile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "For discrete distribution these rules change a bit: the `.pdf()` function is replaced by **`.pmf()`**.\n",
    "* **Example:** binomial distribution with 10 draws and a 0.5 probability of success."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(0,10)\n",
    "plt.scatter(x, stats.binom.pmf(x, n=10, p=0.5))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "## Statistical tests <a id='13'></a>\n",
    "-----------------------\n",
    "\n",
    "**`scipy.stats`** implements a number of statistical tests as functions.\n",
    "* Most tests return two values: the computed **test statistic** and the **p-value**.\n",
    "* We will only demonstrate a couple tests here, but you can get a more in-depth explanation and\n",
    "  demonstration of `scipy.stats` tests\n",
    "  [here](https://machinelearningmastery.com/statistical-hypothesis-tests-in-python-cheat-sheet).\n",
    "\n",
    "<br>\n",
    "\n",
    "### Student's t-test\n",
    "Imagine that we have two samples of measurement drawn from 2 sub-population:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample1 = stats.norm.rvs(size=93, loc=173, scale=20)\n",
    "sample2 = stats.norm.rvs(size=132, loc=181, scale=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test whether the samples have the same mean value, we perform a **t-test to test the equality of the means**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "statistic, pvalue = stats.ttest_ind(sample1, sample2)\n",
    "significance_threshold = 0.05\n",
    "\n",
    "if pvalue < significance_threshold:\n",
    "    print(\"We reject the null hypothesis (H0)\")\n",
    "else:\n",
    "    print(\"We do not reject the null hypothesis (H0)\")\n",
    "print(\"p-value=\",pvalue)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**`stats.ttest_ind`** has an **`equal_var`** parameter that one can set to `False` in order to perform Welsch's t-test, which is warranted when one cannot assume that the two sub-population's variances are equal.\n",
    "\n",
    "> In general, these functions have a very good documentation, detailing the tests and giving usage examples.  \n",
    "> We wholeheartedly recommend any would-be users to have a read at the `help()`.\n",
    "\n",
    "<br>\n",
    "\n",
    "### Chi-Squared Test\n",
    "\n",
    "Imagine that we count different cell types in two biopsies and report them in a list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "biopsy1 = np.array([135, 423, 24, 72])\n",
    "biopsy2 = [184, 552, 77, 101]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now want to test whether the two biopsies differ significantly in their composition using **Chi-squared** test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = [biopsy1 , biopsy2]\n",
    "\n",
    "stat, pvalue, degrees_freedom, expected_values = stats.chi2_contingency(table)\n",
    "print('Chi-square test of independence of variables')\n",
    "print('stat=%.3f, degree of  freedom=%i , p-value=%.3f' % (stat, degrees_freedom, pvalue))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Statistical modelling and regression <a id='13.1'></a>\n",
    "\n",
    "`scipy` implements methods to fit a model to some data. \n",
    "\n",
    "**`scipy.stats`** proposes a simple linear regression function between two variable,\n",
    "while **`scipy.optimize`** implements functions to fit (non-linear) models to data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = stats.uniform.rvs(size=30, loc=0, scale=100)         # generating x\n",
    "y = 1.6*x  + stats.norm.rvs(size=30, loc=0, scale=50)  # Y = 1.6 * x + some noise\n",
    "    \n",
    "# Perform the linear regression:   \n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)\n",
    "print(\"slope:\",slope,\"\\tintercept:\" , intercept)\n",
    "print(\"R-squared:\" , r_value**2)\n",
    "print(\"p-value for the slope:\", p_value)\n",
    "\n",
    "# Plot the data along with the fitted line:\n",
    "plt.plot(x, y, 'o', label='original data')\n",
    "plt.plot(x, intercept + slope*x, 'r', label='fitted line')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "\n",
    "# Example with 2 explanatory variables\n",
    "# ************************************\n",
    "\n",
    "# We define the model as a function.\n",
    "# here it is a form of exponential decay model \n",
    "# where variable1 is the time and the rate of decay changes with variable2\n",
    "# the model function takes as argument the 2 explanatory variable (grouped in a single tuple), and the 4 parameters\n",
    "def func(X , a, b, c , d):\n",
    "    x0 , x1 = X\n",
    "    return a * np.exp( -(b + d*x1 ) * x0 ) + c \n",
    "\n",
    "#realParameter values\n",
    "realParams = [ 5, 0.2, 0.5,  0.1 ]\n",
    "\n",
    "# we simulate some data, with some noise\n",
    "n=50 # number of points\n",
    "variable1 = stats.uniform.rvs(size =n , loc = 0 , scale = 10 ) # explanatory variable number 1 : some uniform variable\n",
    "variable2 = stats.bernoulli.rvs(size=n , p=0.5) # explanatory variable number 2 : can be 0 or 1\n",
    "\n",
    "y = func( (variable1 , variable2) , realParams[0], realParams[1], realParams[2],  realParams[3])\n",
    "y_noise = stats.norm.rvs(size=n , scale = 0.4) \n",
    "ydata = y+y_noise\n",
    "\n",
    "\n",
    "popt, pcov = curve_fit(func, (variable1 , variable2), ydata)\n",
    "perr = np.sqrt(np.diag(pcov))\n",
    "print('parameter estimates          :', popt)\n",
    "print('parameter standard deviation :', perr)\n",
    "print('\\nrelative estimation error    :',np.abs(popt - realParams)/realParams )\n",
    "\n",
    "\n",
    "plt.scatter(variable1, ydata, c=variable2 , label='data')\n",
    "\n",
    "x = np.linspace(min(variable1) , max(variable1) , 100)\n",
    "plt.plot(x, func( ( x , np.zeros(100) )  , *popt), '--' ,  label='fit for variable2==0')\n",
    "plt.plot(x, func( ( x , np.ones(100) )  , *popt), '--' , label='fit for variable2==1')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> You can find most of what was discussed here and more in the [official tutorial](https://docs.scipy.org/doc/scipy/reference/tutorial/stats.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "[Back to ToC](#toc)\n",
    "\n",
    "# Exercises <a id='14'></a>\n",
    "-----------------\n",
    "\n",
    "## Exercise 8.1\n",
    "\n",
    "1. Load the data in the file `data/sample_data.tsv` as a numpy array. Display it's dimensions (number of rows and columns).\n",
    "2. Log-transform the data (use base 2 log).\n",
    "3. Find the row-wise means for replicates of Sample1 and Sample2.\n",
    "4. Find the row-wise standard deviations the same way as means.\n",
    "5. Use a function *scipy.stats.ttest_ind* to calculate p-value for every row.\n",
    "6. Select p-values which are smaller than $10^{-2}$.\n",
    "7. Print how many P-values below $10^{-2}$ are found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Display the top of the input data file using the shell command \"head\".\n",
    "# Note: this command will not work on Windows.\n",
    "!head -3 data/sample_data.tsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution:\n",
    "\n",
    "* **Part 1**: load the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 1-24 solutions/solution_81.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Part 2**: log-transform the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 25-31 solutions/solution_81.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Part 3**: find the row-wise means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 32-52 solutions/solution_81.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Part 4**: find the row-wise standard deviations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 53-60 solutions/solution_81.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Part 5:** calculate p-value for every row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 61-65 solutions/solution_81.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Part 6 and 7:** select p-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load -r 66- solutions/solution_81.py"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}