{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# NumPy introduction and statistics in python with sciPy\n", "\n", "\n", "**www.numpy.org**\n", "\n", "NumPy is the fundamental package for scientific computing with Python.\n", "\n", "#### Highlights\n", "\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 capabilities\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Why NumPy?\n", "\n", "* NumPy arrays is faster than standard python lists\n", "* NumPy functions allows to write many operations with much less code\n", "* NumPy functions are faster than naive Python implementation\n", "* Great collection of mathematical functions available (numpy, scipy, sympy)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How to start?" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# loading numpy module\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The heart of NumPy: array\n", "\n", "NumPy's main object is the **homogeneous** multidimensional array.\n", "\n", "Arrays are very efficient for operations with large numerical data and in general outperform standard Python lists." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "my_array:\n", " [[1. 2. 3.]\n", " [4. 5. 6.]]\n", "my_array dimentions: (2, 3)\n", "my_array number of elements: 6\n", "my_array type of elements float64\n" ] } ], "source": [ "my_array = np.array([[1.,2.,3.],[4.,5.,6.]])\n", "print(\"my_array:\\n\", my_array)\n", "print(\"my_array dimentions:\", my_array.shape)\n", "print(\"my_array number of elements:\", my_array.size)\n", "print(\"my_array type of elements\", my_array.dtype)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One great functionnality 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": 7, "metadata": {}, "outputs": [], "source": [ "my_result = my_array * 3 # multiply by 3 all elements in the array\n", "\n", "#compare with a native python equivalent :\n", "my_list=[[1.,2.,3.],[4.,5.,6.]]\n", "my_result2 = [[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_result2[i][j] = my_list[i][j] * 3\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "native timing : 0.514258861541748\n", "numpy timing : 0.00655055046081543\n", "numpy acceleration factor : 78.50620564149227\n" ] } ], "source": [ "from time import time\n", "native_data = [x for x in range(10**7)]\n", "numpy_data = np.array(native_data)\n", "\n", "t0 = time()\n", "numpy_data *= 3\n", "t1 = time()\n", "numpyTime = t1-t0\n", "\n", "t0 = time()\n", "native_data = [ x*3 for x in native_data ]\n", "t1 = time()\n", "nativeTime = t1-t0\n", "\n", "print(\"native timing :\",nativeTime)\n", "print(\"numpy timing :\",numpyTime)\n", "print(\"numpy acceleration factor :\" , nativeTime / numpyTime )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, let's familiarize ourselves with the numpy array." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How to create NumPy array?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Creating arrays from lists" ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "my_array: [1 2 3]\n" ] } ], "source": [ "# one dimentional array\n", "my_array = np.array([1,2,3])\n", "print(\"my_array:\", type(my_array), my_array)" ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "my_array: \n", "[[1 2 3]\n", " [4 5 6]]\n" ] } ], "source": [ "# two dimentional array\n", "my_array = np.array([[1,2,3],[4,5,6]])\n", "print(\"my_array:\", type(my_array))\n", "print(my_array)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Creating arrays with functions" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0. 0. 0. 0. 0.]\n", " [ 0. 0. 0. 0. 0.]\n", " [ 0. 0. 0. 0. 0.]]\n" ] } ], "source": [ "# array filled with zeroes\n", "my_array = np.zeros((3,5)) # ( number of rows , number of columns )\n", "print(my_array)" ] }, { "cell_type": "code", "execution_count": 82, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 1.]\n", " [ 1. 1.]\n", " [ 1. 1.]\n", " [ 1. 1.]]\n" ] } ], "source": [ "# array filled with ones\n", "my_array = np.ones((4,2))\n", "print(my_array)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[42 42 42]\n", " [42 42 42]]\n" ] } ], "source": [ "# array filled with desired number\n", "my_array = np.full((2,3), 42)\n", "print(my_array)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1. 0. 0. 0.]\n", " [0. 1. 0. 0.]\n", " [0. 0. 1. 0.]\n", " [0. 0. 0. 1.]]\n" ] } ], "source": [ "## identity matrix\n", "my_array = np.eye(4,4,0) \n", "# the first 2 arguments give the matrix dimensions and the 3rd argument specify where the main diagnoal will be\n", "print(my_array)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.2066392 0.5245414 ]\n", " [0.28290859 0.2211281 ]]\n" ] } ], "source": [ "my_array = np.random.rand(2,2) # random numbers from a uniform distribution between 0.0 and 1.0\n", "print(my_array)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Arange function\n", "\n", "Generate one-dimentional array of evenly spaced numbers." ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 1 2 3 4 5]\n" ] } ], "source": [ "my_array = np.arange(6)\n", "print(my_array)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.1 2.2 3.3 4.4 5.5]\n" ] } ], "source": [ "# support float start/end points, as well as float steps\n", "my_array = np.arange(1.1, 6 , 1.1) #start , stop , step\n", "print(my_array)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Array reshaping\n", "\n", "**numpy.reshape** changes the shape of an array without changing its data" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.1 1.8 2.5 3.2 3.9 4.6 5.3 6. ]\n" ] } ], "source": [ "# 1D array\n", "my_array = np.arange(1.1, 6, 0.7)\n", "print(my_array)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1.1 1.8 2.5 3.2]\n", " [3.9 4.6 5.3 6. ]]\n" ] } ], "source": [ "# 2D array\n", "print(np.reshape(my_array, (2,4)))" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[[1.1 1.8]\n", " [2.5 3.2]]\n", "\n", " [[3.9 4.6]\n", " [5.3 6. ]]]\n" ] } ], "source": [ "# 3D array\n", "print(np.reshape(my_array, (2,2,2), order=\"C\"))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reading arrays from files" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**numpy.loadtxt**(fname, dtype=<type float>, comments='#', delimiter=None, converters=None, skiprows=0, usecols=None, unpack=False, ndmin=0)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "a b c\r\n", "0.286368595465 0.0284349169133 0.561298539899\r\n", "0.662679670119 0.718228561506 0.79446312338" ] } ], "source": [ "!cat test.tab" ] }, { "cell_type": "code", "execution_count": 140, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[['a' 'b' 'c']\n", " ['0.286368595465' '0.0284349169133' '0.561298539899']\n", " ['0.662679670119' '0.718228561506' '0.79446312338']]\n" ] }, { "data": { "text/plain": [ "array([['a', 'b', 'c'],\n", " ['0.286368595465', '0.0284349169133', '0.561298539899'],\n", " ['0.662679670119', '0.718228561506', '0.79446312338']],\n", " dtype='start:stop:step notation can also be used when accessing array elements" ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 2. 3. 4. 5.]\n", " [ 6. 7. 8. 9. 10.]\n", " [ 11. 12. 13. 14. 15.]\n", " [ 16. 17. 18. 19. 20.]]\n", "subset:\n", " [[ 1. 3. 5.]\n", " [ 6. 8. 10.]\n", " [ 11. 13. 15.]\n", " [ 16. 18. 20.]]\n" ] } ], "source": [ "my_array = np.arange(1., 21.).reshape((4,5))\n", "print(my_array)\n", "# accessing a subset with step argument\n", "print(\"subset:\\n\", my_array[:,0::2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparison operations" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "my_array:\n", " [[1. 2. 3.]\n", " [4. 5. 6.]\n", " [7. 8. 9.]]\n", "my_array > 5:\n", " [[False False False]\n", " [False False True]\n", " [ True True True]]\n" ] } ], "source": [ "# comparison operators return array of boolean values\n", "my_array = np.arange(1., 10.).reshape((3,3))\n", "print(\"my_array:\\n\", my_array)\n", "print(\"my_array > 5:\\n\", my_array > 5)" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "All values in my_array are greater than 5: False\n", "There is at least one value in my_array greater than 5: True\n" ] } ], "source": [ "# evaluation of boolean arrays\n", "results = my_array > 5\n", "print(\"All values in my_array are greater than 5:\", results.all()) \n", "print(\"There is at least one value in my_array greater than 5:\",\n", " results.any()) " ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Values which are greater than 5: [6. 7. 8. 9.]\n" ] } ], "source": [ "# extracting values with boolen arrays\n", "print(\"Values which are greater than 5:\", my_array[my_array > 5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Iterating through the numpy arrays" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1. 2. 3.]\n", "[4. 5. 6.]\n", "[7. 8. 9.]\n" ] } ], "source": [ "# standard for loop iterates over rows\n", "for x in my_array:\n", " print(x)" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0 2.0 3.0 \n", "4.0 5.0 6.0 \n", "7.0 8.0 9.0 \n" ] } ], "source": [ "# iterating over all elements in standard way\n", "for x in my_array:\n", " for y in x:\n", " print(y, \" \",end=\"\")\n", " print()" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 " ] } ], "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": [ "## How to change array values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Assignment" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1 2]\n", " [3 4]]\n", "After changing one element:\n", " [[1 5]\n", " [3 4]]\n" ] } ], "source": [ "# changing one element\n", "my_array = np.arange(1,5).reshape((2,2))\n", "print(my_array)\n", "my_array[0,1] = 5\n", "print(\"After changing one element:\\n\", my_array)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1 2]\n", " [3 4]]\n", "After changing one row:\n", " [[5 6]\n", " [3 4]]\n", "After changing one column:\n", " [[5 7]\n", " [3 8]]\n" ] } ], "source": [ "# changing rows and columns\n", "my_array = np.arange(1,5).reshape((2,2))\n", "print(my_array)\n", "my_array[0,:] = np.array([5,6])\n", "print(\"After changing one row:\\n\", my_array)\n", "my_array[:,1] = np.array([7,8])\n", "print(\"After changing one column:\\n\", my_array)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding rows and columns to an existing array" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0 1 2]\n", " [3 4 5]]\n", "after column adding\n", " [[0 1 2 6]\n", " [3 4 5 7]]\n" ] } ], "source": [ "# appending\n", "my_array = np.arange(0,6).reshape((2,3))\n", "print(my_array)\n", "# append column\n", "my_array = np.append(my_array, np.array([6, 7]).reshape(2,1), axis=1) # 0 : row , 1: column\n", "print(\"after column adding\\n\", my_array)" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0 1 2]\n", " [3 4 5]]\n", "after row insertion\n", " [[0 1 2]\n", " [6 7 8]\n", " [3 4 5]]\n" ] } ], "source": [ "# insertion\n", "my_array = np.arange(0,6).reshape((2,3))\n", "print(my_array)\n", "# insert row\n", "my_array = np.insert(my_array, 1, [[6, 7, 8]], axis=0)\n", "print(\"after row insertion\\n\", my_array)\n" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "row concatenation:\n", " [[ 0 1 2]\n", " [ 3 4 5]\n", " [ 6 7 8]\n", " [ 9 10 11]]\n", "column concatenation:\n", " [[ 0 1 2 6 7 8]\n", " [ 3 4 5 9 10 11]]\n" ] } ], "source": [ "# concatenation\n", "my_array = np.arange(0,6).reshape((2,3))\n", "my_array2 = np.arange(6,12).reshape((2,3))\n", "print(\"row concatenation:\\n\", np.concatenate((my_array, my_array2),\n", " axis=0))\n", "print(\"column concatenation:\\n\", np.concatenate((my_array, my_array2),\n", " axis=1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## To copy or not to copy?\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." ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 2. 3. 4. 5.]\n", " [ 6. 7. 8. 9. 10.]\n", " [11. 12. 13. 14. 15.]\n", " [16. 17. 18. 19. 20.]]\n", "tmp is my_array True\n", "my array after we changed tmp:\n", " [[ 1. 2. 3. 4. 5.]\n", " [ 6. 999. 8. 9. 10.]\n", " [ 11. 12. 13. 14. 15.]\n", " [ 16. 17. 18. 19. 20.]]\n" ] } ], "source": [ "my_array = np.arange(1., 21.).reshape((4,5))\n", "print(my_array)\n", "tmp = my_array\n", "print(\"tmp is my_array\", tmp is my_array)\n", "tmp[1,1] = 999\n", "print(\"my array after we changed tmp:\\n\", my_array) \n", "# the change is tmp is present in my_array, because they are the same object" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 2. 3. 4. 5.]\n", " [ 6. 7. 8. 9. 10.]\n", " [11. 12. 13. 14. 15.]\n", " [16. 17. 18. 19. 20.]]\n", "tmp is my_array False\n", "tmp:\n", " [[ 7. 8.]\n", " [12. 13.]]\n", "my_array after we changed tmp:\n", " [[ 1. 2. 3. 4. 5.]\n", " [ 6. 7. 8. 9. 10.]\n", " [ 11. 12. 999. 14. 15.]\n", " [ 16. 17. 18. 19. 20.]]\n" ] } ], "source": [ "my_array = np.arange(1., 21.).reshape((4,5))\n", "print(my_array)\n", "tmp = my_array[1:3,1:3]\n", "print(\"tmp is my_array\", tmp is my_array)\n", "print(\"tmp:\\n\", tmp)\n", "tmp[1,1] = 999\n", "print(\"my_array after we changed tmp:\\n\", my_array)\n", "# tmp is not my_array, but the change in one is reported in the other" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> This behavior may appear strange, but it is because numpy arrays access their data by reference. And thus multiple array may access the same memory space, or subset of the same memory space \n", "> See " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### If you plan to change array values but want to keep the old array untouched, then make copy of it!" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 2. 3. 4. 5.]\n", " [ 6. 7. 8. 9. 10.]\n", " [ 11. 12. 13. 14. 15.]\n", " [ 16. 17. 18. 19. 20.]]\n", "tmp is my_array False\n", "my array after we changed tmp:\n", " [[ 1. 2. 3. 4. 5.]\n", " [ 6. 7. 8. 9. 10.]\n", " [ 11. 12. 13. 14. 15.]\n", " [ 16. 17. 18. 19. 20.]]\n", "tmp array after change:\n", " [[ 1. 2. 3. 4. 5.]\n", " [ 6. 999. 8. 9. 10.]\n", " [ 11. 12. 13. 14. 15.]\n", " [ 16. 17. 18. 19. 20.]]\n" ] } ], "source": [ "my_array = np.arange(1., 21.).reshape((4,5))\n", "print(my_array)\n", "tmp = my_array.copy()\n", "print(\"tmp is my_array\", tmp is my_array)\n", "tmp[1,1] = 999\n", "print(\"my array after we changed tmp:\\n\", my_array)\n", "print(\"tmp array after change:\\n\", tmp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Functions\n", "\n", "Numpy provide and great number of functions. The power of numpy functions are speed and possibility to apply it arrays element-wise, column-wise or row-wise." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Element-wise functions\n", "This includes all basic math operators like +, -, /, \\*, //, \\*\\*" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 2. 3.]\n", " [ 4. 5. 6.]]\n", "division\n", "[[ 0.5 1. 1.5]\n", " [ 2. 2.5 3. ]]\n", "sum\n", "[[ 11. 12. 13.]\n", " [ 14. 15. 16.]]\n", "power of 2\n", "[[ 1. 4. 9.]\n", " [ 16. 25. 36.]]\n", "log2\n", "[[ 0. 1. 1.5849625 ]\n", " [ 2. 2.32192809 2.5849625 ]]\n" ] } ], "source": [ "# element-wise operation with scalars\n", "my_array = np.arange(1., 7.).reshape((2,3))\n", "print(my_array)\n", "print(\"division\")\n", "print(my_array/2)\n", "print(\"sum\")\n", "print(my_array + 10)\n", "# element-wise functions\n", "print(\"power of 2\")\n", "print(my_array**2)\n", "print(\"log2\")\n", "print(np.log2(my_array))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sum, division ... of arrays" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 2. 3.]\n", " [ 4. 5. 6.]]\n", "division\n", "[[ 1. 1. 1.]\n", " [ 1. 1. 1.]]\n", "sum\n", "[[ 2. 4. 6.]\n", " [ 8. 10. 12.]]\n" ] } ], "source": [ "# element-wise operation with arrays\n", "my_array = np.arange(1., 7.).reshape((2,3))\n", "print(my_array)\n", "print(\"division\")\n", "print(my_array / my_array)\n", "print(\"sum\")\n", "print(my_array + my_array)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Element-wise, row-wise, column-wise functions" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "element-wise sum: 21.0\n", "sum of columns: [ 5. 7. 9.]\n", "sum of rows: [ 6. 15.]\n" ] } ], "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": 53, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "std of the whole array: 1.70782512766\n", "std of columns: [ 1.5 1.5 1.5]\n", "std of rows: [ 0.81649658 0.81649658]\n" ] } ], "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": 54, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mean of the whole array: 3.5\n", "mean of columns: [ 2.5 3.5 4.5]\n", "mean of rows: [ 2. 5.]\n" ] } ], "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": "markdown", "metadata": {}, "source": [ "## Random numbers in NumPy\n", "\n", "The numpy.random mudule provides a large collection of distributions (uniform, normal, beta, binomial, gamma, poisson ...) to draw from." ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "single number: 0.8840433667904156\n", "array of random numbers:\n", " [[ 0.22295178 0.52276298 0.54877142]\n", " [ 0.49866561 0.8327065 0.85795209]]\n" ] } ], "source": [ "# random numbers from 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": 71, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array of random numbers.\n", " mean : -0.0011793176474030839 \tstandard deviation : 1.000647567911789\n" ] } ], "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": 73, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array permutation:\n", " [0 6 4 5 3 2 1]\n", "sample from the array:\n", " [6 0 1]\n" ] } ], "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,\n", " size=3,\n", " replace=True))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear algebra built-in capabilities\n", "NumPy arrays could be used as matrices without any special conversion.\n", "For advanced linear algebra operations there is spacial package for it (**numpy.linalg**)" ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1. 2. 3.]\n", " [4. 5. 6.]]\n" ] } ], "source": [ "my_array = np.arange(1., 7.).reshape((2,3))\n", "print(my_array)" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "transposed matrix:\n", " [[ 1. 4.]\n", " [ 2. 5.]\n", " [ 3. 6.]]\n" ] } ], "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": 60, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 4. 9.]\n", " [ 16. 25. 36.]]\n" ] } ], "source": [ "# matrix multiplication element-wise\n", "print(my_array * my_array)" ] }, { "cell_type": "code", "execution_count": 77, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[14. 32.]\n", " [32. 77.]]\n", "[[14. 32.]\n", " [32. 77.]]\n" ] } ], "source": [ "# matrix product\n", "print(my_array.dot(my_array.T)) \n", "# one can also use the @ operator :\n", "print(my_array @ my_array.T) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# SciPy.stats and statistics in python\n", "\n", "\n", "SciPy references a comprehensive [project for scientific python programming](https://scipy.org) regrouping as well as a [library](https://docs.scipy.org/doc/scipy/reference/) (which is part of the project) 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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## manipulation of random distributions\n", "\n", "scipy.stats implements utilisties for a large number of continuous and discrete distributions :" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "number of continuous distributions: 98\n", "number of discrete distributions: 14\n", "['alpha', 'anglit', 'arcsine', 'argus', 'beta', 'betaprime', 'bradford', 'burr', 'burr12', 'cauchy', 'chi', 'chi2', 'cosine', 'crystalball', 'dgamma', 'dweibull', 'erlang', 'expon', 'exponnorm', 'exponpow', 'exponweib', 'f', 'fatiguelife', 'fisk', 'foldcauchy', 'foldnorm', 'frechet_l', 'frechet_r', 'gamma', 'gausshyper', 'genexpon', 'genextreme', 'gengamma', 'genhalflogistic', 'genlogistic', 'gennorm', 'genpareto', 'gilbrat', 'gompertz', 'gumbel_l', 'gumbel_r', 'halfcauchy', 'halfgennorm', 'halflogistic', 'halfnorm', 'hypsecant', 'invgamma', 'invgauss', 'invweibull', 'johnsonsb', 'johnsonsu', 'kappa3', 'kappa4', 'ksone', 'kstwobign', 'laplace', 'levy', 'levy_l', 'levy_stable', 'loggamma', 'logistic', 'loglaplace', 'lognorm', 'lomax', 'maxwell', 'mielke', 'moyal', 'nakagami', 'ncf', 'nct', 'ncx2', 'norm', 'norminvgauss', 'pareto', 'pearson3', 'powerlaw', 'powerlognorm', 'powernorm', 'rayleigh', 'rdist', 'recipinvgauss', 'reciprocal', 'rice', 'semicircular', 'skewnorm', 't', 'trapz', 'triang', 'truncexpon', 'truncnorm', 'tukeylambda', 'uniform', 'vonmises', 'vonmises_line', 'wald', 'weibull_max', 'weibull_min', 'wrapcauchy']\n" ] } ], "source": [ "from scipy import stats\n", "\n", "dist_continu = [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", "print('number of continuous distributions: %d' % len(dist_continu))\n", "print('number of discrete distributions: %d' % len(dist_discrete))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "let's experiment with the normal distribution, or `norm` in `scipy.stats`\n", "\n", "\n", "A look at `help(stats.norm)` tells us that \n", "```\n", " | The location (``loc``) keyword specifies the mean.\n", " | The scale (``scale``) keyword specifies the standard deviation.\n", "```\n" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array(10.), array(4.))" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## we can generate a specific normal distribution :\n", "N = stats.norm(loc = 10 , scale = 2)\n", "\n", "# the mean and variance of a distribution can be retrieved using the .stats method :\n", "print(N.stats())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That object can then be used to interact with the distribution in many ways.\n", "\n", "### drawing some random numbers : rvs" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 7.45179805, 9.87769114, 10.12902769, 10.8202259 , 8.85423502],\n", " [ 8.39733275, 12.62407038, 12.54939775, 7.57128479, 10.62743881],\n", " [ 7.11035717, 9.2620774 , 8.46154685, 10.78523221, 10.11458767],\n", " [14.17995768, 10.08394262, 9.90331856, 8.97369216, 9.83082144],\n", " [ 7.56909984, 7.17413853, 7.0261789 , 10.76444972, 11.875346 ]])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# draw some random number in this distribution : rvs\n", "# the size argument is 1 or several integers and defines the dimensions of the returned arrays of random numbers\n", "N.rvs(size = [5,5]) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "as with any drawing of random variable on a computer, [one merely emulates randomness](https://en.wikipedia.org/wiki/Pseudorandom_number_generator). This also means that one can make some random operation reproducible by setting up the random seed.\n", "\n" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "are the ramdom draws equal? [ True True True True True]\n" ] } ], "source": [ "import numpy as np\n", "np.random.seed(2020) # we set the random seed\n", "draw1 = N.rvs(size=5)\n", "np.random.seed(2020) # we set the random seed back to 2020\n", "draw2 = N.rvs(size=5)\n", "print(\"Are the ramdom draws equal?\",draw1 == draw2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### looking up the quantiles and probability density functions\n" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3xU9Zn48c8zmVwIkEBuXAOBEC5RFDUgF1EKYtG2YlutuLZia9fart3udnd/62537a7b3dd2L7XbXV/bsmqrVuutdaUrFrGKyFWCIFeBJNzCJVcICeQ+z++POaPjMCETksmZzDzv12syZ875njPPnMzMM+d7vuf7FVXFGGNM4vG4HYAxxhh3WAIwxpgEZQnAGGMSlCUAY4xJUJYAjDEmQXndDqAncnJytKCgwO0wjDFmQNm2bVutquaGzh9QCaCgoIDS0lK3wzDGmAFFRI6Em29VQMYYk6AsARhjTIKyBGCMMQnKEoAxxiQoSwDGGJOgIkoAIrJERPaLSJmIPBRm+XdFZK+I7BSR34vI+KBly0XkoHNbHjT/GhHZ5WzzJyIiffOSjDHGRKLbBCAiScBjwM1AMXCXiBSHFNsOlKjqFcDLwL8462YB3weuBWYB3xeR4c46/w3cDxQ5tyW9fjXGGGMiFsl1ALOAMlWtABCR54GlwN5AAVV9O6j8ZuDLzvSngTWqWu+suwZYIiJrgQxV3eTMfxq4DXi9V6/GmH7W0eljU0UdB6qaEGDKyKFcOyELb5LVrprYF0kCGAMcC3pcif8XfVfu4+Mv8nDrjnFulWHmX0BE7sd/pMC4ceMiCNeY/vHqjuP82xv7OVbf/In547PT+fObpvC5K0e7FJkxkYkkAYSrmw87ioyIfBkoAW7oZt2It6mqK4AVACUlJTZ6jXFde6ePv//tHn65+ShXjM3ke7dMo6QgC4D3DtXz2NtlfPtX29l25DTf+8w0ku1owMSoSBJAJZAf9HgscCK0kIjcCHwPuEFVW4PWXRCy7lpn/tjutmlMrPH5lD9/6QNe3XGCB24o5C8+PYUkz8e/Z26ZPoqbikfwz69/yOPrD9HY0sG/3XEF1sbBxKJIfppsBYpEZIKIpADLgJXBBUTkKuBnwK2qWh20aDVwk4gMd07+3gSsVtWTQKOIzHZa/9wDvNoHr8eYqPrh7z7k1R0n+H9LpvDQzVM/8eUf4E3y8DefLeZPb5zMr9+v5N/fOOBCpMZ0r9sjAFXtEJEH8X+ZJwFPquoeEXkEKFXVlcC/AkOAl5xfOkdV9VZVrReRf8CfRAAeCZwQBr4J/AIYhP+cgZ0ANjHt9/uq+Nm6Cr48exzfvKGw2/J/vGgSJ840819vlzFrQhbXT76gM0ZjXCUDaVD4kpIStd5AjRtqm1r59KPryB2ayqsPziPVmxTRei3tnXzuP9fT0NzO7/7kerIGp0Q5UmMuJCLbVLUkdL6dnTImAv/02j4aWzr4j2VXRfzlD5CWnMSPl82g/lwbP3z9wyhGaEzPWQIwphvbjpzmN9uP8/X5E5gycmiP179sdCZfnVfAi9uOsbPyTBQiNObSWAIw5iJUlb//7R5GZKTyR5+adMnb+faiIrIHp/D3v93LQKp2NfHNEoAxF7FmbxU7Kxv4s5umMDj10gfQy0hL5ruLp7DtyGnWHqjpwwiNuXSWAIzpgqry4zcPUpCdzheuCnuheo/cfs1YxgwbxI/XHLCjABMTLAEY04U1e6vYe/Is315Y1Cd9+6R4PXx74SQ+qGxg7X47CjDuswRgTBdWrKsgP2sQS2f0XZ8+X7xmLKMz01ixrqLPtmnMpbIEYEwYO46dofTIab46d0Kf9uyZnORh+dwCNlXUsedEQ59t15hLYQnAmDCeWH+IoalevjQzv/vCPbRs1jjSU5J4Yv2hPt+2MT1hCcCYENVnW1i16yR3zsxnSC9a/nQlc1AyXyrJ57cfnKCuqbX7FYyJEksAxoR4aVslnT7l7tnjuy98ie6+dhztncqv36/svrAxUWIJwJggPp/yq/eOMrcwmwk5g6P2PEUjhjKzYDi/eu+YNQk1rrEEYEyQ9WW1VJ5u5q5Z0R997q5Z4zhUe47NFfXdFzYmCiwBGBPk5W2VDE9P5qbLRkT9uW6ZPoqhaV5eLD3WfWFjosASgDGOc60drNlbxS3TR/Wox89LlZacxGevGMXvdp/iXGtH1J/PmFARJQARWSIi+0WkTEQeCrP8ehF5X0Q6ROT2oPmfEpEdQbcWEbnNWfYLETkUtGxG370sY3ruzX1VNLd3snRG77t9iNTnrxpLc3snq/ec6rfnNCag2wQgIknAY8DNQDFwl4gUhxQ7CtwLPBc8U1XfVtUZqjoDWAicB94IKvIXgeWquuPSX4YxvffqjhOMzkyjZPzwfnvOkvHDyc8axCvbj/fbcxoTEMkRwCygTFUrVLUNeB5YGlxAVQ+r6k7Ad5Ht3A68rqrnLzlaY6Kk/lwb6w7U8LkZo/GEGec3Wjwe4fMzxrC+rJbqxpZ+e15jILIEMAYIPktV6czrqWXAr0Lm/aOI7BSRR0UkNdxKInK/iJSKSGlNjXWgZaJj1a6TdPiUpVf2X/VPwGeuGI0qvLGnqt+f2yS2SBJAuJ9DPWq4LCKjgOn4B5YP+CtgKjATyAL+Mty6qrpCVUtUtSQ31wbVNtGxcscJJuUNYdqono/41VuTRwxhYu5gXt99st+f2yS2SBJAJRDcIcpY4EQPn+dLwCuq2h6Yoaon1a8V+Dn+qiZj+t3xM828d7iepVeORqT/qn8CRIRbLh/F5op66xrC9KtIEsBWoEhEJohICv6qnJU9fJ67CKn+cY4KEP8n7jZgdw+3aUyf+O0H/t8zt/Zht889dfP0kXT6lDV7rRrI9J9uE4CqdgAP4q++2Qe8qKp7ROQREbkVQERmikglcAfwMxHZE1hfRArwH0G8E7LpZ0VkF7ALyAF+0PuXY0zP/W73KaaPyWR8dvS6fuhO8agMxmens2q3NQc1/Seirg5VdRWwKmTew0HTW/FXDYVb9zBhThqr6sKeBGpMNFSfbWHHsTP82eLJrsYhItx8+Sgef7eCM+fbGJae4mo8JjHYlcAmob25rxqAxf3Q9UN3bpk+kg6rBjL9yBKASWhr9p4iP2sQU0b0f+ufUNPHZDJm2CBet2og008sAZiE1dTawYayOm4qHulK659QIsIt00fy7sEazra0d7+CMb1kCcAkrHUHamjr9LG42P3qn4CbLhtJe6ey/mCt26GYBGAJwCSsNXurGJae3K99/3TnqvxhZA5K5q0Pq90OxSQASwAmIbV3+njrw2oWTs3DmxQ7HwNvkocbJueydn81Pp+NFGaiK3be+cb0o62H62lobuem4pFuh3KBRdPyqG1qY+fxBrdDMXHOEoBJSGv2VpHq9XD95By3Q7nADZNz8Qi8tc+ag5rosgRgEtLa/TXMKcwmPSWiayH71bD0FK4eN5y39tt5ABNdlgBMwjlSd45DtedYMDl2e5f91NQ8dh8/S9VZGyPARI8lAJNw1h3wjytxw5Q8lyPp2qJp/tjettZAJoosAZiEs3Z/DeOy0inITnc7lC5NGTGU0Zlp1hzURJUlAJNQWjs62Vhex4IpuTFx9W9XRIRPTc1jfVktLe2dbodj4pQlAJNQSg+fprm9kxtiuP4/YOHUPM63dbLtyGm3QzFxyhKASShr91eTkuRh9sRst0Pp1uyJ2SQnCesO2ljYJjosAZiE8s6BGmZOGM7g1Nhr/hlqcKqXa8YP590D1i+QiY6IEoCILBGR/SJSJiIPhVl+vYi8LyIdInJ7yLJOEdnh3FYGzZ8gIltE5KCIvOAMN2lM1Jw408yBqiYWTI7d1j+h5hflsvfkWWoabaxg0/e6TQAikgQ8BtwMFAN3iUhxSLGjwL3Ac2E20ayqM5zbrUHzfwg8qqpFwGngvkuI35iIvfNR88/Yr/8PuL7IH+uGMjsKMH0vkiOAWUCZqlaoahvwPLA0uICqHlbVnYAvkid1BoJfCLzszHoK/8DwxkTNugM1jMpMoyhviNuhROyy0RkMT0+28wAmKiJJAGOAY0GPKwkzxu9FpIlIqYhsFpHAl3w2cMYZcP6i2xSR+531S2tq7ENgLk2nT9lYXsd1k3JiuvlnKI9HuK4ol3cP1qJqvYOavhVJAgj3aenJO3GcqpYAfwD8WEQKe7JNVV2hqiWqWpKbO3AO3U1s2XfyLA3N7cydFPutf0LNL8qhprGV/VWNbodi4kwkCaASyA96PBY4EekTqOoJ574CWAtcBdQCw0Qk0BSjR9s0pqc2lvvr0OcWxl7vn92ZX+SPOdCFhTF9JZIEsBUoclrtpADLgJXdrAOAiAwXkVRnOgeYB+xV/7Hs20CgxdBy4NWeBm9MpDaU1VGYO5gRGWluh9JjozIHUZQ3hHdtmEjTx7pNAE49/YPAamAf8KKq7hGRR0TkVgARmSkilcAdwM9EZI+z+jSgVEQ+wP+F/8+qutdZ9pfAd0WkDP85gSf68oUZE9DW4WPr4XrmTRp4v/4D5hflsuVQvXULYfpURFfDqOoqYFXIvIeDprfir8YJXW8jML2LbVbgb2FkTFR9UHmG822dzC0cePX/AfMn5/DkhkNsPVzP/CI7F2b6hl0JbOLexrI6RBgQ3T90ZVZBFl6PsLG8zu1QTByxBGDi3sbyWi4bncGw9IF7sfngVC9XjRvGRrsgzPQhSwAmrjW3dbL96JkB2fon1JzCHHYdb6Chud3tUEycsARg4lrpkXraOn0Duv4/YF5hNj6FLRVWDWT6hiUAE9c2lNXh9QgzC7LcDqXXZowbRlqyx84DmD5jCcDEtU3ltczIHzYgun/uTqo3iZkFWR9d1GZMb1kCMHGrobmdXccbmDuA2/+HmluYw4GqJqobW9wOxcQBSwAmbm2pqMOnxEX9f8A8py+jTVYNZPqAJQATtzaW15GW7OGqccPcDqXPXDY6k4w0LxvLLAGY3rMEYOLWxvJaZhZkkepNcjuUPpPkEWZPzGZjhZ0HML1nCcDEpZrGVg5UNTEnjqp/AuYWZnOsvplj9efdDsUMcJYATFza5LSVnxcHF4CFCnRqZ62BTG9ZAjBxaWNZLUPTvFw+JtPtUPrcpLwh5A5NZYOdBzC9ZAnAxKWN5XXMnphNkmfgDP8YKRFhbmE2G8vrbJhI0yuWAEzcOVZ/nqP15+Oq+WeouYXZ1Da1crC6ye1QzABmCcDEnUAb+YE8AEx3Ap3bWe+gpjciSgAiskRE9otImYg8FGb59SLyvoh0iMjtQfNniMgmEdkjIjtF5M6gZb8QkUMissO5zeibl2QS3YbyWnKGpFCUN8TtUKImPyud/KxBbLALwkwvdNtBiogkAY8Bi/EPEL9VRFYGDe0IcBS4F/jzkNXPA/eo6kERGQ1sE5HVqnrGWf4Xqvpyb1+EMQGqysbyOuYU5iASf/X/weYV5vDarpN0+jQuz3WY6IvkCGAWUKaqFaraBjwPLA0uoKqHVXUn4AuZf0BVDzrTJ4BqwMazM1FTXtNETWMr8+K4/j9gTmE2jS0d7DnR4HYoZoCKJAGMAY4FPa505vWIiMwCUoDyoNn/6FQNPSoiqV2sd7+IlIpIaU1NTU+f1iSYQNPIeBgApjuBi9yse2hzqSJJAOGOLXvU9kxERgHPAF9V1cBRwl8BU4GZQBbwl+HWVdUVqlqiqiW5uXbwYC5uY3ktY4cPYlx2utuhRF3e0DSK8oZYAjCXLJIEUAnkBz0eC5yI9AlEJAN4DfgbVd0cmK+qJ9WvFfg5/qomYy5Zp0/ZVF4X180/Q80tzGbroXraOnzdFzYmRCQJYCtQJCITRCQFWAasjGTjTvlXgKdV9aWQZaOcewFuA3b3JHBjQu09cZazLR0JUf0TMKcwh+b2Tj6oPNN9YWNCdJsAVLUDeBBYDewDXlTVPSLyiIjcCiAiM0WkErgD+JmI7HFW/xJwPXBvmOaez4rILmAXkAP8oE9fmUk4gb5xEukIYPbELESw7qHNJYlonDxVXQWsCpn3cND0VvxVQ6Hr/RL4ZRfbXNijSI3pxobyOiblDSEvI83tUPrNsPQULhudwcbyWr5zY5Hb4ZgBxq4ENnGhrcPH1kP1CdH8M9Tcwhy2Hz1Dc1un26GYAcYSgIkLO46dobm9kzkJVP8fMKcwm7ZOH9uOnHY7FDPAWAIwcWFjeS0i/jrxRDOzIAuvR2x8ANNjlgBMXNhYXsflozMZlp7idij9bkiqlyvzh9n1AKbHLAGYAe98Wwfbj55OqNY/oeYWZrPreAONLe1uh2IGEEsAZsArPXya9k5lbhx3/9ydOYXZdPqUrYfr3Q7FDCCWAMyAt6G8Fq9HmFkw3O1QXHP1uOGkeD12PYDpEUsAZsDbVF7HVeOGkZ4S0WUtcSktOYmS8cPtPIDpEUsAZkBrON/O7uMNCdX9Q1fmFmaz9+RZTp9rczsUM0BYAjAD2uZDdfg0sbp/6ErgGojNFXYUYCJjCcAMaJvK60hL9nDVuMSt/w+4Ymwmg1OSrBrIRMwSgBnQNpTVMrMgixSvvZWTkzzMmpBlF4SZiNmnxgxY1Y0tHKxusvr/IHMLcyivOUfV2Ra3QzEDgCUAM2Btcqo65k2y+v+AwDCRm6wayETAEoAZsDaW1ZGR5uWy0ZluhxIzikdlkDko2aqBTEQsAZgBa2NFLbMnZpPkCTdsdWLyeIQ5E7PtRLCJSEQJQESWiMh+ESkTkYfCLL9eRN4XkQ4RuT1k2XIROejclgfNv0ZEdjnb/IkzNKQxETlSd45j9c3MS+DuH7oyd1I2laebOVZ/3u1QTIzrNgGISBLwGHAzUAzcJSLFIcWOAvcCz4WsmwV8H7gW/6Dv3xeRQHu9/wbuB4qc25JLfhUm4awv81dxXFdkCSBU4JoIqwYy3YnkCGAWUKaqFaraBjwPLA0uoKqHVXUn4AtZ99PAGlWtV9XTwBpgiTMgfIaqblJVBZ7GPzC8MRFZf7CWUZlpTMwZ7HYoMacwdwi5Q1OtGsh0K5IEMAY4FvS40pkXia7WHeNMd7tNEblfREpFpLSmpibCpzXxrNOnbCyv47pJOVjN4YVEhLmF/vMA/t9XxoQXSQII9wmL9F3V1boRb1NVV6hqiaqW5ObmRvi0Jp7tPt5AQ3O7Vf9cxNzCbGoaWymvaXI7FBPDIkkAlUB+0OOxwIkIt9/VupXO9KVs0yS4QP2/nQDuWuDiuA3WPbS5iEgSwFagSEQmiEgKsAxYGeH2VwM3ichw5+TvTcBqVT0JNIrIbKf1zz3Aq5cQv0lA6w/WMm1UBjlDUt0OJWblZ6UzLiuddw/aiWDTtW4TgKp2AA/i/zLfB7yoqntE5BERuRVARGaKSCVwB/AzEdnjrFsP/AP+JLIVeMSZB/BN4HGgDCgHXu/TV2biUnNbJ9uOnOY6u/q3W/OLcthUXkt7Z2jbDGP8IhpBQ1VXAatC5j0cNL2VT1bpBJd7EngyzPxS4PKeBGvMe4fraev0cV2RnQ/qzvyiXJ7dcpTtR88wa0KW2+GYGGRXApsBZf3BGlKSPMwqsC+07syd5L9Ket0Baz1nwrMEYAaU9WV1XDN+OINSktwOJeZlpCVzVf4w3j1oCcCEZwnADBg1ja3sO3nWmn/2wPyiXHYeb7BhIk1YlgDMgBHo2uA6a/4ZsfmTc1CFDdYthAnDEoAZMNYfrCVzUDKXj7HunyN1xZhMMtK8vHvAEoC5kCUAMyCoKuvLapk3ybp/7glvkod5k3J492CNdQthLmAJwAwIFbXnONnQYlf/XoL5RbmcaGixbiHMBSwBmAFhvXNF6/xJ1v6/p+Y7J83XWTWQCWEJwAwI7xyoYXx2OuOy090OZcDJz0pnQs5gaw5qLmAJwMS8lvZONpbXsmCy/fq/VPOLcthcUU9rR6fboZgYYgnAxLz3DtXT0u5jwZQ8t0MZsOYX5dLc7u9HyZgASwAm5q3dX0OK18PsidYB3KWaU5iN1yPWO6j5BEsAJuatPVDN7InZ1v1DLwxJ9XL1+OHWL5D5BEsAJqYdqz9PRc05q//vA9cX5bDnxFmqG1vcDsXECEsAJqat3V8NwIIplgB6K3AOZe1+OwowfhElABFZIiL7RaRMRB4KszxVRF5wlm8RkQJn/t0isiPo5hORGc6ytc42A8vsDJ+5wNr9NYxzmjGa3rlsdAYjMlJ5+8Nqt0MxMaLbBCAiScBjwM1AMXCXiBSHFLsPOK2qk4BHgR8CqOqzqjpDVWcAXwEOq+qOoPXuDixXVXtXmk/wN/+sY8GUXPwjh5reEBEWTs3j3YO1tHXYKGEmsiOAWUCZqlaoahvwPLA0pMxS4Cln+mVgkVz4ib0L+FVvgjWJZevheprbO636pw8tnDqCptYOth6u776wiXuRJIAxwLGgx5XOvLBlnDGEG4DQNnt3cmEC+LlT/fO3YRKGSXCB5p9zJlr/P31l3qRsUrwe3rJqIENkCSDcF3Not4IXLSMi1wLnVXV30PK7VXU6MN+5fSXsk4vcLyKlIlJaU2MnrxKFqvLmvirmFlrzz76UnuJlzsRsSwAGiCwBVAL5QY/HAie6KiMiXiATCD7GXEbIr39VPe7cNwLP4a9quoCqrlDVElUtyc21qoBEUVbdxJG689w4bYTbocSdhVPzOFR7jgrrHTThRZIAtgJFIjJBRFLwf5mvDCmzEljuTN8OvKVO5+Mi4gHuwH/uAGeeV0RynOlk4LPAboxxrNlXBWAJIAoWTvU3uLOjANNtAnDq9B8EVgP7gBdVdY+IPCIitzrFngCyRaQM+C4Q3FT0eqBSVSuC5qUCq0VkJ7ADOA78T69fjYkbb+6t4oqxmYzMTHM7lLiTn5VOUd4Q3t5vCSDReSMppKqrgFUh8x4Omm7B/ys/3Lprgdkh884B1/QwVpMgahpb2X7sDH9642S3Q4lbC6fl8cS7hzjb0k5GWrLb4RiX2JXAJua89WEVqlb9E02Lp42gw6d2VXCCswRgYs6avdWMGTaIaaOGuh1K3Lp63HByhqSyevcpt0MxLrIEYGJKc1sn68tquHFanl39G0Uej7C4eARr91fT0m6DxCQqSwAmpmwoq6Wl3cfi4pFuhxL3llw+knNtnWwoszECEpUlABNT1uytYmiql1kTstwOJe7NmZjN0DQvq/dYNVCisgRgYkZ7p4/Ve0+xcFoeKV57a0ZbitfDwql5rNlbRUendQ6XiOxTZmLGpvI6zpxv55bpo9wOJWEsuWwkp8+3s/WwjRWciCwBmJjx+u6TDE5J4gYb/avf3DAll1Svx6qBEpQlABMTOjp9rN5TxaJpI0hLts7f+kt6ipf5Rbm8secUTu8tJoFYAjAxYXNFPfXn2qz6xwVLLh/JiYYWdhw743Yopp9ZAjAx4bVd/uofG/yl/y0uHkFKkofffnDS7VBMP7MEYFznr/45xUKr/nFF5qBkFkzJ5bc7T9Dps2qgRGIJwLhuyyF/9c9nptvFX265dcZoahpb2VJR53Yoph9ZAjCue23XSdJTklgwJc/tUBLWoqkjGJySxMoPQsd6MvHMEoBxVWtHJ6/tPMniYqv+cdOglCQWF4/g9d2naOuwi8IShSUA46q3P6yhobmdL1w91u1QEt6tM0bT0NzOugPWRXSiiCgBiMgSEdkvImUi8lCY5aki8oKzfIuIFDjzC0SkWUR2OLefBq1zjYjsctb5iVjXjwnpN+9Xkjs0lXmF2W6HkvDmF+UyPD3ZqoESSLcJQESSgMeAm4Fi4C4RKQ4pdh9wWlUnAY8CPwxaVq6qM5zbA0Hz/xu4Hyhybksu/WWYgej0uTbe3l/NbTNG402yg1G3JSd5uHn6KNbsreJ8W4fb4Zh+EMmnbhZQpqoVqtqGf3D3pSFllgJPOdMvA4su9oteREYBGaq6yRk8/mngth5Hbwa0/9t5gvZOteqfGLL0ytE0t3fyOxsoJiFEkgDGAMeCHlc688KWcQaRbwACx/QTRGS7iLwjIvODyld2s00AROR+ESkVkdKaGqubjCe/fv84U0cOZdqoDLdDMY5ZE7IoyE7nha3Hui9sBrxIEkC4X/KhV4t0VeYkME5VrwK+CzwnIhkRbtM/U3WFqpaoaklurl0lGi/Ka5rYcewMX7Rf/zFFRLijJJ8th+o5XHvO7XBMlEWSACqB/KDHY4HQs0QflRERL5AJ1Ktqq6rWAajqNqAcmOyUD/7kh9umiWP/u/04HoGlM0a7HYoJ8cWrx+IReLHUjgLiXSQJYCtQJCITRCQFWAasDCmzEljuTN8OvKWqKiK5zklkRGQi/pO9Fap6EmgUkdnOuYJ7gFf74PWYAaCj08fL2yqZX5RLXkaa2+GYECMz01gwJY+Xt1XaQDFxrtsE4NTpPwisBvYBL6rqHhF5RERudYo9AWSLSBn+qp5AU9HrgZ0i8gH+k8MPqGq9s+ybwONAGf4jg9f76DWZGPfWh9WcbGjhD64d53YopgtfKsmnurGVdQftvFs880ZSSFVXAatC5j0cNN0C3BFmvV8Dv+5im6XA5T0J1sSHX245ysiMNBZNta4fYtWiaXnkDEnhha3HWDh1hNvhmCixxtemXx2tO8+6AzUsm5Vvbf9jWHKShy9cPZbf76umprHV7XBMlNgn0PSrZ987QpJHWDbTqn9i3Z0z8+nwKc+/d9TtUEyUWAIw/aa1o5OXSiu5cVoeIzPt5G+sK8wdwvyiHJ7ZfMQ6iItTlgBMv/nd7lPUn2vjy7PHux2KidDXrptAdWMrr++20cLikSUA02+e2XSEgux05hXmuB2KidANRblMzBnMk+sP2aDxccgSgOkX247UU3rkNPfMKcDjsY5fBwqPR7h3XgEfVDbw/lEbND7eWAIw/eKn71QwLD2ZZbPyuy9sYsoXrx7L0DQvP99wyO1QTB+zBGCirqy6iTV7q7hnTgHpKRFdemJiyOBUL3eW5PP67lOcbGh2OxzThywBmKhbsa6ctGQPy+fYyd+BavncAgBWrKtwNxDTpywBmKg61dDCK9uP86WSfLKHpLodjrlE+VnpfP6qMTy35SjVjS1uh2P6iCUAE1U/33CITp/y9esmuh2K6aUHPzWJ9k4fK96xo4B4YQnARE1tUyu/3HyEz1wxmv/lD8gAABAlSURBVHHZ6W6HY3qpIGcwt80Ywy+3HLHuIeKEJQATNY+9XUZLh48/ubHI7VBMH3lw4STaOnw8/q4dBcQDSwAmKo7Vn+fZzUe545qxFOYOcTsc00cm5g7h1itH8/SmI9Q12VHAQGcJwETFj988CALfsV//cefBhUW0dnTyn2+VuR2K6SVLAKbP7T/VyG+2V3Lv3AJGZQ5yOxzTxyblDWHZrHE8s/kIZdWNbodjeiGiBCAiS0Rkv4iUichDYZanisgLzvItIlLgzF8sIttEZJdzvzBonbXONnc4NxsdJE782xv7GZLi5Zs3FLodiomS7y6eTHpyEv/42j63QzG90G0CcMb0fQy4GSgG7hKR4pBi9wGnVXUS8CjwQ2d+LfA5VZ2Of8zgZ0LWu1tVZzi36l68DhMj1u6vZs3eKh5YUMjwwSluh2OiJGdIKn+8qIi399ewdr99dAeqSI4AZgFlqlqhqm3A88DSkDJLgaec6ZeBRSIiqrpdVU848/cAaSJiVwPFqZb2Th5+dQ8Tcwfz9fkT3A7HRNnyuQUUZKfzg9f20W6Dxw9IkSSAMcCxoMeVzrywZZxB5BuA7JAyXwS2q2pw04GfO9U/fysiYbuIFJH7RaRUREpramyA6lj22NtlHK0/zw9uu5xUb5Lb4ZgoS/F6+OtbplFW3cRTGw+7HY65BJEkgHBfzKEdg1+0jIhchr9a6BtBy+92qobmO7evhHtyVV2hqiWqWpKbmxtBuMYNZdVN/PSdcr5w1RjmWn//CWNx8QgWTc3j397Yz6Hac26HY3ookgRQCQT34TsWONFVGRHxAplAvfN4LPAKcI+qlgdWUNXjzn0j8Bz+qiYzAPl8yt/87y4GJSfx15+Z5nY4ph+JCP/0hemkJHn4i5c+oNNng8YMJJEkgK1AkYhMEJEUYBmwMqTMSvwneQFuB95SVRWRYcBrwF+p6oZAYRHxikiOM50MfBbY3buXYtzyxPpDbK6o569vmUaOdfiWcEZkpPH9z11G6ZHTNmbAANNtAnDq9B8EVgP7gBdVdY+IPCIitzrFngCyRaQM+C4QaCr6IDAJ+NuQ5p6pwGoR2QnsAI4D/9OXL8z0j12VDfzL6g/59GUjuHOmDfaSqL5w9RgWTc3jX1fvp6Kmye1wTIRkII3zWVJSoqWlpW6HYRxNrR189ifv0trh4/XvzGdYujX7TGRVZ1v49I/XMTIjjVe+NY9BKdYQIFaIyDZVLQmdb1cCm0uiqnz/1T0crT/Pj++cYV/+hhEZafz4zhnsr2rkod/stEHkBwBLAOaSPLnhML9+v5JvLyzi2omhLX5NolowJY8/WzyZV3ec4OcbDrsdjumGJQDTY6v3nOIHr+3l5stH8p1F1tmb+aRvLZjE4uIR/OOqfWwsr3U7HHMRlgBMj+w4dobvPL+dK8cO49E7Z+DxhL1+zyQwj0f40ZeuZELOYL7x9DZ2VTa4HZLpgiUAE7Gy6ka+/tRWcoem8vjyEtKS7SSfCW9oWjLP3DeLjEHJ3PPkFg5UWa+hscgSgInI3hNnufNnmwHh5/fOsvb+plujMgfx3B9eS3KShy8/voXDdqVwzLEEYLq149gZ7vqfzaR4Pbz4jdlMyrMRvkxkxmcP5pdfv5b2Th+3/3QTOyvPuB2SCWIJwFzU7/dV8eXHt5AxyMuL35jDRBve0fTQ5BFDeemBOaQle7jzZ5t5c2+V2yEZhyUAE5bPp/xozQHue6qUgpx0XvzGHPKz0t0OywxQk/KG8ptvzaVoxBDuf6aUx9+twGf9BrnOEoC5QG1TK/c9tZWf/P4gt18zlpcfmGtDO5peyxuaxvP3z2Zx8Qh+8No+7v3FVqrPtrgdVkKzBGA+oqq8VHqMG3/0DhvK6viH2y7nX2+/wlr7mD6TnuLlp1++hn+47XK2VNSx5D/eZdWuk3bVsEu8bgdgYsO+k2f5wWt72VBWR8n44fzzF6czKW+o22GZOCQifGX2eOZMzOKPf7WDbz37PnMmZvPw54qZNirD7fASinUGl+DKqht59M2DvLbzJENTvfy/m6dy96xxdoGX6RcdnT6ee+8oP1pzgLPN7Xzx6rF844ZCa2nWx7rqDM4SQALq9CnrDtTw7JYjvPVhNWnJSXx1XgF/OH+idepmXNFwvp3/+P1Bnt1yhNYOHzdOy+Nr8yYwe2K2/RjpA5YAEpyq8kFlA6v3nGLljhMcP9NMzpBU7pw5lq/Nm0C2XdhlYkBdUytPbzrC05sOc/p8O6Mz01h61Rg+d8Vopo0aShdDh5tuWAJIQMfqz7PlUD1bKupYX1bLyYYWvB5hTmE2y2aOY3HxCFK81g7AxJ7mtk7e2HuK/91+nHUHa+n0KSMyUrm+KJfrinK4etxwxg4fZAkhQr1KACKyBPgPIAl4XFX/OWR5KvA0cA1QB9ypqoedZX8F3Ad0An+sqqsj2WY4lgAu5PMp1Y2tHK0/z9H685RVN7H35Fn2njhLbVMrAJmDkpk9MYubikeyaFqeVfOYAaW2qZW39lXzzoEa3j1Yw9mWDgByhqRwxdhhTMobwsScwUzMHUJh7mCyBqdYYghxyQlARJKAA8Bi/IO/bwXuUtW9QWW+BVyhqg+IyDLg86p6p4gUA7/CP+D7aOBNYLKz2kW3GU48JQBVpcOndHQqbZ0+2gO3DqW5vZOm1g7OObem1g7Ot3Vy5nw7dedaqWtqo6apldrGVirPNNPW4ftou8lJQlHeUC4bncHlYzKZNSGLKSOGWj2qiQsdnT4+PNXI9mNn2HH0DHtONHCo9hytQZ+BQclJ5GWkMmJoGrkZqeQNTSUjLZmhad6P7oc692nJSaR4PaR4PaQ69ylJ/ul4SiJdJYBImoHOAspUtcLZ0PPAUiD4y3op8HfO9MvAf4l/7y0FnlfVVuCQM2bwLKdcd9vsM3/9yi62VNQBoM6fQNpTVRQI5EFFP57Wj8sE1lX1l+Gj6eCyQesGbfujsqp0+pR250v/UgxLTyZ7cArZQ1KZOmooNxaPID8rnfFZ6YzLSmf0sEFWrWPiljfJw+VjMrl8TCZfmT0e8DdqOHGmmfKaJsprznHyTDNVja1Un21h74mzvNPYSlNrR4+fKzlJEBE8AkkieEQQ8Xd3nSQfL/OIkORxljnzgoUmkgvSinS9LHjdJ5fPZFx2316NH0kCGAMcC3pcCVzbVRlV7RCRBiDbmb85ZN0xznR32wRARO4H7gcYN25cBOGGeQHDBjF1ZMZHe1f82/1oZ4t8PC+w/OOy8tHyj8v65wUeB5Z+vJ2u1hW8HiHZ6yE5yUNy0HRKkuBN8uD1COkpXganJjEk1Ut6ipchqf7HGYOSSU6yL3djgiV5hPysdPKz0lkwJXyZTp/S1NLB2ZZ2Gls6aHTuWzt8tHV20tbh8087960d/iNyn/p/1Pl8ik/Bpxp0C8wPWub7+EcffPwj8qPHIXEF18BcUBcTMiMaP+wiSQDhjoNCY+2qTFfzw72SsHVRqroCWAH+KqCuw+zaH31q0qWsZoyJE0keITM9mcz0ZLdDiSmRpJRKID/o8VjgRFdlRMQLZAL1F1k3km0aY4yJokgSwFagSEQmiEgKsAxYGVJmJbDcmb4deEv9xzYrgWUikioiE4Ai4L0It2mMMSaKuq0Ccur0HwRW42+y+aSq7hGRR4BSVV0JPAE845zkrcf/hY5T7kX8J3c7gD9S1U6AcNvs+5dnjDGmK3YhmDHGxLmumoFakxJjjElQlgCMMSZBWQIwxpgEZQnAGGMS1IA6CSwiNcCRS1w9B6jtw3D6SqzGBbEbm8XVMxZXz8VqbJca13hVzQ2dOaASQG+ISGm4s+Bui9W4IHZjs7h6xuLquViNra/jsiogY4xJUJYAjDEmQSVSAljhdgBdiNW4IHZjs7h6xuLquViNrU/jSphzAMYYYz4pkY4AjDHGBLEEYIwxCSruEoCILBGR/SJSJiIPhVmeKiIvOMu3iEhBP8SULyJvi8g+EdkjIt8JU2aBiDSIyA7n9nC043Ke97CI7HKe84Ke9sTvJ87+2ikiV/dTXFOC9sUOETkrIn8SUqZf9pmIPCki1SKyO2heloisEZGDzv3wLtZd7pQ5KCLLw5Xp47j+VUQ+dP5Xr4jIsC7Wvej/PQpx/Z2IHA/6X93SxboX/fxGKbYXguI6LCI7ulg3mvss7HdE1N9nqho3N/xdS5cDE4EU4AOgOKTMt4CfOtPLgBf6Ia5RwNXO9FDgQJi4FgD/58I+OwzkXGT5LcDr+Ed3mw1scen/egr/xSz9vs+A64Grgd1B8/4FeMiZfgj4YZj1soAK5364Mz08ynHdBHid6R+GiyuS/3sU4vo74M8j+D9f9PMbjdhClv878LAL+yzsd0S032fxdgTw0QD2qtoGBAabD7YUeMqZfhlYJCLhhq7sM6p6UlXfd6YbgX18PDZyrFsKPK1+m4FhIjKqn2NYBJSr6qVeBd4rqroO/zgXwYLfR08Bt4VZ9dPAGlWtV9XTwBpgSTTjUtU3VDUwAvpm/KPt9asu9lckIvn8Ri0253vgS8Cv+vI5I3GR74iovs/iLQGEG8A+9Iv2EwPYA4EB7PuFU+V0FbAlzOI5IvKBiLwuIpf1U0gKvCEi20Tk/jDLI9mn0baMrj+UbuwzgBGqehL8H14gL0wZt/fd1/AfvYXT3f89Gh50qqae7KIqw+39NR+oUtWDXSzvl30W8h0R1fdZvCWA3gxgH3UiMgT4NfAnqno2ZPH7+Ks4rgT+E/jf/ogJmKeqVwM3A38kIteHLHdtfwGIf8jQW4GXwix2a59Fys332vfwj8L3bBdFuvu/97X/BgqBGcBJ/FUtoVx9rwF3cfFf/1HfZ918R3S5Wph5Ee23eEsAvRnAPqpEJBn/P/ZZVf1N6HJVPauqTc70KiBZRHKiHZeqnnDuq4FX8B+GB4tkn0bTzcD7qloVusCtfeaoClSFOffVYcq4su+ck4CfBe5Wp5I4VAT/9z6lqlWq2qmqPuB/ung+195rznfBF4AXuioT7X3WxXdEVN9n8ZYAejOAfdQ4dYtPAPtU9UddlBkZOBchIrPw/2/qohzXYBEZGpjGfwJxd0ixlcA94jcbaAgckvaTLn+VubHPggS/j5YDr4Ypsxq4SUSGO1UeNznzokZElgB/Cdyqque7KBPJ/72v4wo+b/T5Lp4vks9vtNwIfKiqleEWRnufXeQ7Irrvs2ic0Xbzhr/VygH8rQm+58x7BP8HAiANf3VCGfAeMLEfYroO/yHZTmCHc7sFeAB4wCnzILAHf8uHzcDcfohrovN8HzjPHdhfwXEJ8JizP3cBJf34v0zH/4WeGTSv3/cZ/gR0EmjH/2vrPvznjX4PHHTus5yyJcDjQet+zXmvlQFf7Ye4yvDXBwfeZ4EWb6OBVRf7v0c5rmec989O/F9qo0Ljch5f8PmNdmzO/F8E3ldBZftzn3X1HRHV95l1BWGMMQkq3qqAjDHGRMgSgDHGJChLAMYYk6AsARhjTIKyBGCMMQnKEoAxxiQoSwDGGJOg/j8XPZF79OLN6AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "what is the probability of drawing a number <=15.0 ? 0.9937903346742238\n", "quantiles: [0.025, 0.5, 0.975] -> [ 6.08007203 10. 13.91992797]\n" ] } ], "source": [ "# pdf: Probability Density Function\n", "# I know this is not the plotting lesson, but here is a small recipe to plot the distribution\n", "import matplotlib.pyplot as plt \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", "\n", "# ppf: Percent Point Function (Inverse of CDF) , gives the quantiles of the distribution\n", "P = [0.025,0.5,0.975]\n", "Q = N.ppf(P)\n", "print( 'quantiles:', P , '->' , Q )\n" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAARwElEQVR4nO3df2xd533f8fdnlJywDVq5sf6RZEXKpmoxli4abpVuxlJg+SEFGyytSBBlyOAOAYxu9dYtmwZrBRrA/aNZNQzdH8YWo8lQrD9c1xEEYXPHZbW7f4ZkosIkmuxxVdTUJpkhKhxlw8LFkvLdH7yyrxjKPIxIHum57xdA6J7nec7ll0f3fHj4nIeXqSokSe36M30XIEnaWAa9JDXOoJekxhn0ktQ4g16SGrel7wKWu++++2rPnj19lyFJd5Vz5879aVVtX6nvjgv6PXv2MD093XcZknRXSfInt+pz6kaSGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjesU9EkOJ5lNcjHJYyv0fyLJC0m+muQPkrxtpO96ki8PP86sZ/GSpNWt+u6VSSaAJ4D3A3PA2SRnquqFkWEzwKCqvpPk7wK/Cnxk2LdYVe9a57olSR11uaI/CFysqktV9SrwFHBkdEBVPV9V3xlufgHYtb5lSpJ+UF3ej34n8PLI9hzw7jcY/3Hg90e235xkGrgGfKqqTi/fIckjwCMAu3fv7lCSdPc5PTPPyalZFq4ssmPbJMcP7efogZ19l6Ux0CXos0JbrTgw+RgwAH56pHl3VS0keTvwXJLzVfW1m56s6kngSYDBYLDic0t3s9Mz85w4dZ7Fq9cBmL+yyIlT5wEMe224LlM3c8D9I9u7gIXlg5K8D/hF4KGq+u6N9qpaGP57CfhD4MBt1CvdlU5Ozb4W8jcsXr3OyanZnirSOOkS9GeBfUn2JrkHOAbctHomyQHg0yyF/DdH2u9N8qbh4/uAB4HRm7jSWFi4srimdmk9rRr0VXUNeBSYAl4Enq6qC0keT/LQcNhJ4C3A7y1bRvkOYDrJV4DnWZqjN+g1dnZsm1xTu7SeOv1x8Kp6Fnh2WdsvjTx+3y32+6/AO2+nQKkFxw/tv2mOHmBy6wTHD+3vsSqNi05BL+n23Ljh6qob9cGglzbJ0QM7DXb1wve6kaTGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjdvSdwHSRjs9M8/JqVkWriyyY9skxw/t5+iBnX2X1QuPxXgy6NW00zPznDh1nsWr1wGYv7LIiVPnAcYu4DwW48upGzXt5NTsa8F2w+LV65ycmu2pov54LMZXp6BPcjjJbJKLSR5bof8TSV5I8tUkf5DkbSN9Dyf5o+HHw+tZvLSahSuLa2pvmcdifK0a9EkmgCeADwIPAB9N8sCyYTPAoKp+AngG+NXhvj8GfBJ4N3AQ+GSSe9evfOmN7dg2uab2lnksxleXK/qDwMWqulRVrwJPAUdGB1TV81X1neHmF4Bdw8eHgM9X1StV9S3g88Dh9SldWt3xQ/uZ3DpxU9vk1gmOH9rfU0X98ViMry43Y3cCL49sz7F0hX4rHwd+/w329a6PNs2Nm4yuNPFYjLMuQZ8V2mrFgcnHgAHw02vZN8kjwCMAu3fv7lCS1N3RAzsNsyGPxXjqMnUzB9w/sr0LWFg+KMn7gF8EHqqq765l36p6sqoGVTXYvn1719olSR10CfqzwL4ke5PcAxwDzowOSHIA+DRLIf/Nka4p4ANJ7h3ehP3AsE2StElWnbqpqmtJHmUpoCeAz1bVhSSPA9NVdQY4CbwF+L0kAC9V1UNV9UqSX2bpmwXA41X1yoZ8JZKkFaVqxen23gwGg5qenu67DEm6qyQ5V1WDlfr8zVhJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1rlPQJzmcZDbJxSSPrdD/niRfSnItyYeW9V1P8uXhx5n1KlyS1M2W1QYkmQCeAN4PzAFnk5ypqhdGhr0E/CzwT1Z4isWqetc61CpJ+gGsGvTAQeBiVV0CSPIUcAR4Leir6uvDvu9tQI2SpNvQZepmJ/DyyPbcsK2rNyeZTvKFJEdXGpDkkeGY6cuXL6/hqSVJq+kS9FmhrdbwOXZX1QD4W8CvJfmz3/dkVU9W1aCqBtu3b1/DU0uSVtMl6OeA+0e2dwELXT9BVS0M/70E/CFwYA31SZJuU5egPwvsS7I3yT3AMaDT6pkk9yZ50/DxfcCDjMztS5I23qpBX1XXgEeBKeBF4OmqupDk8SQPAST5ySRzwIeBTye5MNz9HcB0kq8AzwOfWrZaR5K0wVK1lun2jTcYDGp6errvMiTprpLk3PB+6PfxN2MlqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklq3Ja+C1C7Ts/Mc3JqloUri+zYNsnxQ/s5emAtf1derfK1sbkMem2I0zPznDh1nsWr1wGYv7LIiVPnATyhx5yvjc3n1I02xMmp2ddO5BsWr17n5NRsTxXpTuFrY/MZ9NoQC1cW19Su8eFrY/MZ9NoQO7ZNrqld48PXxuYz6LUhjh/az+TWiZvaJrdOcPzQ/p4q0p3C18bm82asNsSNm2qurNByvjY2X6qq7xpuMhgManp6uu8yJOmukuRcVQ1W6nPqRpIaZ9BLUuMMeklqnEEvSY3rFPRJDieZTXIxyWMr9L8nyZeSXEvyoWV9Dyf5o+HHw+tVuCSpm1WDPskE8ATwQeAB4KNJHlg27CXgZ4HfXrbvjwGfBN4NHAQ+meTe2y9bktRVlyv6g8DFqrpUVa8CTwFHRgdU1der6qvA95btewj4fFW9UlXfAj4PHF6HuiVJHXUJ+p3AyyPbc8O2Ljrtm+SRJNNJpi9fvtzxqSVJXXQJ+qzQ1vW3rDrtW1VPVtWgqgbbt2/v+NSSpC66BP0ccP/I9i5goePz386+kqR10CXozwL7kuxNcg9wDDjT8fmngA8kuXd4E/YDwzZJ0iZZNeir6hrwKEsB/SLwdFVdSPJ4kocAkvxkkjngw8Cnk1wY7vsK8MssfbM4Czw+bJMkbRLf1EySGuCbmknSGDPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNa5T0Cc5nGQ2ycUkj63Q/6Ykvzvs/2KSPcP2PUkWk3x5+PFv1rd8SdJqtqw2IMkE8ATwfmAOOJvkTFW9MDLs48C3qurPJTkG/HPgI8O+r1XVu9a5bklSR12u6A8CF6vqUlW9CjwFHFk25gjwG8PHzwDvTZL1K1OS9IPqEvQ7gZdHtueGbSuOqaprwLeBtw779iaZSfJfkvzVlT5BkkeSTCeZvnz58pq+AEnSG+sS9CtdmVfHMd8AdlfVAeATwG8n+ZHvG1j1ZFUNqmqwffv2DiVJkrpadY6epSv4+0e2dwELtxgzl2QL8KPAK1VVwHcBqupckq8BPw5M327hurXTM/OcnJpl4coiO7ZNcvzQfo4eWP5DmKRxOVe6XNGfBfYl2ZvkHuAYcGbZmDPAw8PHHwKeq6pKsn14M5ckbwf2AZfWp3St5PTMPCdOnWf+yiIFzF9Z5MSp85yeme+7NOmOMk7nyqpBP5xzfxSYAl4Enq6qC0keT/LQcNhngLcmucjSFM2NJZjvAb6a5Css3aT9uap6Zb2/CL3u5NQsi1ev39S2ePU6J6dme6pIujON07nSZeqGqnoWeHZZ2y+NPP5/wIdX2O9zwOdus0atwcKVxTW1S+NqnM4VfzO2MTu2Ta6pXRpX43SuGPSNOX5oP5NbJ25qm9w6wfFD+3uqSLozjdO50mnqRnePGysGxmElgXQ7xulcydIKyDvHYDCo6WlXX0rSWiQ5V1WDlfqcupGkxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TG+acE19npmfmx+NNkktbPRueGQb+OTs/Mc+LUeRavXgdg/soiJ06dBzDsJa1oM3LDqZt1dHJq9rX/rBsWr17n5NRsTxVJutNtRm4Y9Oto4crimtolaTNyw6BfRzu2Ta6pXZI2IzcM+nV0/NB+JrdO3NQ2uXWC44f291SRpDvdZuRGp6BPcjjJbJKLSR5bof9NSX532P/FJHtG+k4M22eTHFq3ypc5PTPPg596jr2P/Qce/NRznJ6Z36hPdUtHD+zkV37mnezcNkmAndsm+ZWfeac3YiXd0mbkRqrqjQckE8D/BN4PzAFngY9W1QsjY/4e8BNV9XNJjgF/s6o+kuQB4HeAg8AO4D8DP15V15d/nhsGg0FNT0+v6YtYftcalr4jGrKSxkWSc1U1WKmvyxX9QeBiVV2qqleBp4Ajy8YcAX5j+PgZ4L1JMmx/qqq+W1V/DFwcPt+6crWLJN1al6DfCbw8sj03bFtxTFVdA74NvLXjvrfN1S6SdGtdgj4rtC2f77nVmC77kuSRJNNJpi9fvtyhpJu52kWSbq1L0M8B949s7wIWbjUmyRbgR4FXOu5LVT1ZVYOqGmzfvr179UOudpGkW+sS9GeBfUn2JrkHOAacWTbmDPDw8PGHgOdq6S7vGeDYcFXOXmAf8N/Wp/TXudpFkm5t1fe6qaprSR4FpoAJ4LNVdSHJ48B0VZ0BPgP8uyQXWbqSPzbc90KSp4EXgGvAz7/RipvbcfTAToNdklaw6vLKzfaDLK+UpHF3u8srJUl3MYNekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcXfc2xQnuQz8yW08xX3An65TOXc7j8XNPB6v81jcrIXj8baqWvFP9N1xQX+7kkzf6j2Zx43H4mYej9d5LG7W+vFw6kaSGmfQS1LjWgz6J/su4A7isbiZx+N1HoubNX08mpujlyTdrMUreknSCINekhrXTNAnOZxkNsnFJI/1XU+fktyf5PkkLya5kOQX+q6pb0kmkswk+fd919K3JNuSPJPkfwxfI3+575r6kuQfDc+R/57kd5K8ue+aNkITQZ9kAngC+CDwAPDRJA/0W1WvrgH/uKreAfwU8PNjfjwAfgF4se8i7hD/CviPVfXngb/ImB6XJDuBfwAMquovABPAsX6r2hhNBD1wELhYVZeq6lXgKeBIzzX1pqq+UVVfGj7+PyydyDv7rao/SXYBfx349b5r6VuSHwHeA3wGoKperaor/VbVqy3AZJItwA8BCz3XsyFaCfqdwMsj23OMcbCNSrIHOAB8sd9KevVrwD8Fvtd3IXeAtwOXgX87nMr69SQ/3HdRfaiqeeBfAC8B3wC+XVX/qd+qNkYrQZ8V2sZ+3WiStwCfA/5hVf3vvuvpQ5K/AXyzqs71XcsdYgvwl4B/XVUHgP8LjOU9rST3svST/15gB/DDST7Wb1Ubo5WgnwPuH9neRaM/gnWVZCtLIf9bVXWq73p69CDwUJKvszSl99eS/Ga/JfVqDpirqhs/4T3DUvCPo/cBf1xVl6vqKnAK+Cs917QhWgn6s8C+JHuT3MPSDZUzPdfUmyRhaQ72xar6l33X06eqOlFVu6pqD0uvi+eqqsmrti6q6n8BLyfZP2x6L/BCjyX16SXgp5L80PCceS+N3pje0ncB66GqriV5FJhi6c75Z6vqQs9l9elB4G8D55N8edj2z6rq2R5r0p3j7wO/NbwougT8nZ7r6UVVfTHJM8CXWFqpNkOjb4XgWyBIUuNambqRJN2CQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIa9/8BRCDF1LgPBoMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# For discrete distribution these rules change a bit , the pdf function is replaced by pmf:\n", "X = np.arange(0,10)\n", "plt.scatter( X , stats.binom.pmf( X , n = 10 , p = 0.5 ) ) # binomial distribution with 10 draws and a 0.5 probability of success\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## statistical tests\n", "\n", "`scipy.stats` implements a number of statistical tests as functions.\n", "\n", "Most return two values : the computed test statistic and the p-value.\n", "\n", "We will only demonstrate a couple tests here.\n", "You can get a more in-depth explaination and demonstration of scipy.stats tests [there](https://machinelearningmastery.com/statistical-hypothesis-tests-in-python-cheat-sheet/)\n", "\n" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "We reject the hypothesis of equality of means(H0). p-value : 0.0011955595116142342\n" ] } ], "source": [ "#Imagine we have two samples of measurement drawn from 2 sub-population :\n", "sample1 = stats.norm.rvs(size = 93 , loc = 173 , scale = 20)\n", "sample2 = stats.norm.rvs(size = 132 , loc = 181 , scale = 20)\n", "\n", "# we perform a t-test to test the equality of the means\n", "statistic , pValue = stats.ttest_ind(sample1 , sample2) \n", "significanceThreshold = 0.05\n", "if pValue < significanceThreshold:\n", " print( \"We reject the hypothesis of equality of means(H0). p-value :\" , pValue )\n", "else:\n", " print( \"We do not reject the hypothesis of equality of means(H0). p-value :\" , pValue )\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`stats.ttest_ind` has a `equal_var` parameter that one can set to `False` in order to perform Welsch's t-test, which is warranted when one cannot assume the two sub-population's variances to be equal.\n", "\n", "> In general, these functions have a very good documentation, detailing the tests and giving usage examples. We heartily recommend any would-be users to have a read at the `help()`." ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Chi-square test of independence of variables\n", "stat=14.555, degree of freedom=3 , p-value=0.002\n" ] } ], "source": [ "# Example of the Chi-Squared Test\n", "# imagine you count different cell types in two biopsies and report them in a list :\n", "biopsy1 = np.array([135 , 423 , 24 , 72])\n", "biopsy2 = [184 , 552 , 77 , 101]\n", "\n", "table = [biopsy1 , biopsy2]\n", "\n", "stat, pValue, degreeOfFreedom, expectedValues = 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, degreeOfFreedom, pValue))\n", "# here the two biopsies seem to differ significatively in their composition\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### statiscial modelling and regression\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": 75, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "slope: 1.618178 intercept: -13.018001\n", "R-squared: 0.526027\n", "p-value for the slope: 0.000006\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df3QU9b3/8ee7GCH+KAFFCwkC3uuNIj8SjFSlV6RcRJEi0t4Wj9Zo24ui1tp7L4rtbbW1HrCgUnvVFq1gW9vqxRipilSQHmtVNIgNCPIVFWsCBUSDKEFD/Hz/2An54e4mm93ZnZl9Pc7JSXZmsvPZzPLis+/PzGfMOYeIiETTZ3LdABER8Y9CXkQkwhTyIiIRppAXEYkwhbyISIQdlOsGtHXkkUe6wYMH57oZIiKhsmbNmnecc/3irQtUyA8ePJiamppcN0NEJFTM7K1E61SuERGJMIW8iEiEKeRFRCIsUDX5eJqamqirq2Pfvn25bkpe69WrFyUlJRQUFOS6KSKSgsCHfF1dHYcffjiDBw/GzHLdnLzknGPXrl3U1dUxZMiQXDdHRFIQ+JDft2+fAj7HzIwjjjiCnTt35ropIpFTvbaeecs3sbWhkQFFhcyaWMrU8uKMPX/gQx5QwAeAjoFI5lWvree6qnU0NjUDUN/QyHVV6wAyFvQaeBURyZF5yzcdCPgWjU3NzFu+KWP7UMhn0KRJk2hoaEi6zQ9/+ENWrFjRref/85//zOTJkzvd7owzzuj0orIFCxawd+/ebrVDRDJja0NjSsu7IxTlmlT4Xd+KxzmHc47HH3+8021//OMf+9qWrlqwYAEXXnghhxxySK6bIpK3BhQVUh8n0AcUFWZsH5HqybfUt+obGnG01req19an9by33norw4YNY9iwYSxYsACALVu2cMIJJ3D55ZczatQo3n77bQYPHsw777wDwI033sjxxx/PhAkTOP/885k/fz4AF198MUuWLAFi0zhcf/31jBo1iuHDh/Pqq68C8MILL3DaaadRXl7OaaedxqZNyT+6NTY2Mn36dEaMGMHXvvY1Ghtb3zQzZ86koqKCE088keuvvx6A22+/na1btzJu3DjGjRuXcDsR8desiaUUFvRot6ywoAezJpZmbB+R6sknq291tze/Zs0aFi1axOrVq3HO8fnPf56xY8fSp08fNm3axKJFi7jzzjvb/U5NTQ0PPfQQa9euZf/+/YwaNYqTTjop7vMfeeSRvPTSS9x5553Mnz+fe+65h+OPP56nn36agw46iBUrVvC9732Phx56KGEb77rrLg455BBqa2upra1l1KhRB9bddNNN9O3bl+bmZsaPH09tbS1XXXUVt956K6tWreLII49MuN2IESO69TcTka5pyaW8P7umq/yobz3zzDOcd955HHrooQBMmzaNv/zlL0yZMoVBgwZxyimnxP2dc889l8LC2EeuL33pSwmff9q0aQCcdNJJVFVVAbB7924qKyt57bXXMDOampqStvHpp5/mqquuAmDEiBHtwvnBBx9k4cKF7N+/n23btrFhw4a44d3V7UQks6aWF/taUo5UuSZRHSud+layG523BH8qv9NRz549AejRowf79+8H4Ac/+AHjxo1j/fr1/PGPf+zS1b7xTnF88803mT9/PitXrqS2tpZzzjkn7nN1dTsRCZ8uh7yZDTSzVWa20cxeMbPveMv7mtmTZvaa972Pt9zM7HYz22xmtWY2Kvke0udHfev000+nurqavXv38uGHH/Lwww/zr//6r0l/5wtf+MKBcP7ggw947LHHUtrn7t27KS6O/c++ePHiLrXx/vvvB2D9+vXU1tYC8P7773PooYfSu3dvtm/fzrJlyw78zuGHH86ePXs63U5Ewi2Vcs1+4L+ccy+Z2eHAGjN7ErgYWOmcm2tms4HZwLXA2cBx3tfngbu8777xo741atQoLr74YkaPHg3At771LcrLy9myZUvC3zn55JOZMmUKI0eOZNCgQVRUVNC7d+8u7/Oaa66hsrKSW2+9lS9+8Yudbj9z5kwuueQSRowYQVlZ2YG2jhw5kvLyck488USOPfZYxowZc+B3ZsyYwdlnn03//v1ZtWpVwu1EJNwsldJCu180ewT4X+/rDOfcNjPrD/zZOVdqZr/0fv69t/2mlu0SPWdFRYXreH73xo0bOeGEE7rVxlz64IMPOOyww9i7dy+nn346CxcubDcgGkZhPRYiUWdma5xzFfHWdWvg1cwGA+XAauDoluD2gv4ob7Ni4O02v1bnLWsX8mY2A5gBcMwxx3SnOYE0Y8YMNmzYwL59+6isrAx9wItIOKUc8mZ2GPAQcLVz7v0kc5rEW/Gpjw3OuYXAQoj15FNtT1D97ne/y3UTRERSO7vGzAqIBfz9zrkqb/F2r0yD932Ht7wOGNjm10uArek1V0REUpHK2TUG/ArY6Jy7tc2qpUCl93Ml8Eib5Rd5Z9mcAuxOVo8XEZHMS6VcMwb4OrDOzF72ln0PmAs8aGbfBP4O/Lu37nFgErAZ2AtckpEWi4hIl3U55J1zzxC/zg4wPs72Driim+0SEZEMiNQVr365/fbbOeGEE7jgggtYunQpc+fOBaC6upoNGzYc2G7x4sVs3ZrasMOWLVsYNmxY0uU1NTUHpi0QEUlFpOau8cudd97JsmXLDtzfdMqUKUAs5CdPnszQoUOBWMgPGzaMAQMGZHT/FRUVVFTEPQVWRCQp9eQ7cdlll/HGG28wZcoUbrvtNhYvXsyVV17Js88+y9KlS5k1axZlZWXcfPPN1NTUcMEFF1BWVkZjYyNr1qxh7NixnHTSSUycOJFt22LjzmvWrGHkyJGceuqp3HHHHZ22oe3NQm644Qa+8Y1vcMYZZ3Dsscdy++23H9jut7/9LaNHj6asrIxLL72U5ubmRE8pInkiXD35q6+Gl1/ufLtUlJWBN0d8PL/4xS944oknDkzL2zKXzGmnncaUKVOYPHkyX/nKVwBYtmwZ8+fPp6KigqamJr797W/zyCOP0K9fPx544AG+//3vc++993LJJZfw85//nLFjxzJr1qyUm/zqq6+yatUq9uzZQ2lpKTNnzmTz5s088MAD/PWvf6WgoIDLL7+c+++/n4suuqhbfxYRiYZwhXyIbNq0ifXr1zNhwgQAmpub6d+/P7t376ahoYGxY8cC8PWvfz3lCcHOOeccevbsSc+ePTnqqKPYvn07K1euZM2aNZx88slA7EYiRx11VCfPJCJRF66QT9LjDhrnHCeeeCLPPfdcu+UNDQ1xpwVORcv0xNA6RbFzjsrKSubMmZPWc4tItKgmn4a20/V2fFxaWsrOnTsPhHxTUxOvvPIKRUVF9O7dm2eeeQbgwBTB6Ro/fjxLlixhx47YBcfvvvsub731VkaeW0TCSyGfhunTpzNv3jzKy8t5/fXXufjii7nssssoKyujubmZJUuWcO211zJy5EjKysp49tlnAVi0aBFXXHEFp5566oG7R6Vr6NCh/OQnP+HMM89kxIgRTJgw4cBAr4jkr25PNeyHKE01HEU6FiLBlGyqYfXkRUQiTCEvIhJhoQj5IJWU8pWOgUg4Bf4Uyl69erFr1y6OOOKItE89lO5xzrFr1y569eqV66aIhFL12vqM3ns6FYEP+ZKSEurq6ti5c2eum5LXevXqRUlJSa6bIRI61Wvrua5qHY1NsWlG6hsaua5qHUBWgj7wIV9QUHBgYjARkbCZt3zTgYBv0djUzLzlm7IS8qGoyYuIhNXWhsaUlmeaQl5ExEcDiuJf8JhoeaYp5EVEfDRrYimFBT3aLSss6MGsiaVZ2X/ga/IiImHWUnfX2TUiIhE1tbw4a6Hekco1IiIRppAXEYkwlWtERHyUy6tdQSEvIuKbXF/tCirXiIj4JtnVrtmikBcR8Umur3YFhbyIiG9yfbUrKORFRHyT66tdIYWQN7N7zWyHma1vs+wGM6s3s5e9r0lt1l1nZpvNbJOZTcx0w0VEgm5qeTFzpg2nuKgQA4qLCpkzbXhgz65ZDPwv8OsOy29zzs1vu8DMhgLTgROBAcAKM/sX51wzIiJ5JJdXu0IKPXnn3NPAu13c/FzgD865j5xzbwKbgdHdaJ+IiKQhEzX5K82s1ivn9PGWFQNvt9mmzlv2KWY2w8xqzKxGd38S+bTqtfWMmfsUQ2Y/xpi5T1G9tj7XTYqMQPxt//EPuOgiWLfOl6dPN+TvAv4JKAO2Abd4y+PdjDXunaCdcwudcxXOuYp+/fql2RyRaGm5mKa+oRFH68U0Cvr05fRv+8knMHcumEH//vCb38A11/iyq7SueHXObW/52czuBh71HtYBA9tsWgJsTWdfIvko17eOi7Lu/m07m6Yg6foXX4Tx42HPnvZPesst8N3vZuy1tZVWT97M+rd5eB7QcubNUmC6mfU0syHAccAL6exLJB8F4WKaqOrO37az3n+89b+8e1msx24Go0e3Bvz48bB9OzgH//mfsfU+6HJP3sx+D5wBHGlmdcD1wBlmVkasFLMFuBTAOfeKmT0IbAD2A1fozBqR1A0oKqQ+Tuhk82KaTMn1RF0ddedv21nvv+36LTdPjv8ky5fDmWd2v+Ep6nLIO+fOj7P4V0m2vwm4qTuNEpGYWRNL201wBdm/mCYTgjBRV0fd+dt21vu/4VffY8Lm1fF/+YMP4NBDu9/gbtIVryIBFoSLaTIhCBN1ddSdv228Xn6/D97lzZsng9mnAv7ek6Yw+NpHGTNnZU4CHjTVsEjg5fpimkwI6thCqn/btr3/hOUYYPC1jx74OdefvBTykhNBq8+Kv6IytjB1xe+Y+pMEpzo+9RSMG0f12nqKA/TeVshL1gWxPiv+CvXYQlMTHHxw4vWu/SVAQfvkpZq8ZF0Q67Pir1COLbSc9hgv4BsbY+Hu4l7jGSjqyUvWBbU+K/4KWg83rj/9CSYmmDR39myYMye77ckAhbxkXVTqsxIhSS5EGnLto6219Sw2KVMU8pJ1oa7PSnT06QMNDXFXrfy/lVz5t48jMW6kmrxkXSjrsxIN9fWttfZ4Ae/V2X+4mciMG6knLzkRivqsREeyeWE++eRT66M0bqSevERGIOYGl+CYNau1197RLbe0nh0TZ30QbsCdKerJSyRk6tx7XaQVcs3NcFCSWOviKY9RGjdST14iIRPn3usGHSHW0mOPF/C7dqV8TnuUxo3Uk5dIyEQNVTfoCJnnn4dTT42/bsQI+Nvf0nr6qIwbKeQlJUEtZ2Ti3PsoDbZFWrJB1BBcgZptKtdIlwW5nDFrYimFBT3aLUu1hhqlwbbIOeusxIOoq1eHZoqBXFDIS5cFec6ZTNRQM/EfhWTQe++1Bvvy5Z9e3xLso0dnv20honKNdFnQyxnp1lBbfjeI5ai8kqwc09wMn1HfNBUKeemyfJhzJiqDbaGzcCFcemn8dfPmwX//d3bbEyEKeemyKJ07LAHgXPJeuWrsGaGQly5TOUMyIlk5Zvt2OOqo7LUlDyjkJSUqZ0i3rFsXO3c9nooKePHF7LYnjyjkRcQ/Oqc95zRMLSKZdeGFic9pX7FC57RnmXryIpK+Dz+Eww5LvF6hnjMKeRHpvmTlmI8/hoKC7LVF4lK5RkRS88ADicsx117bWo5RwAdCl3vyZnYvMBnY4Zwb5i3rCzwADAa2AF91zr1nZgb8DJgE7AUuds69lNmmi0hWaRA1lFLpyS8GzuqwbDaw0jl3HLDSewxwNnCc9zUDuCu9ZopITvTtm7jXvmWLBlFDoMsh75x7Gni3w+Jzgfu8n+8DprZZ/msX8zxQZGb9022siGTB66+3Bvt777VfN3Bga7APGpSb9klK0h14Pdo5tw3AObfNzFouVSsG3m6zXZ23bFvHJzCzGcR6+xxzzDFpNkdEuk3lmEjya+A13rsl7rvEObfQOVfhnKvo16+fT80RkbiuvjpxOaaqSuWYCEi3J7/dzPp7vfj+wA5veR0wsM12JcDWNPclIpnw8cfQs2fi9Qr1SEm3J78UqPR+rgQeabP8Ios5BdjdUtYRkRxp6bHHC/i9e9Vrj6guh7yZ/R54Dig1szoz+yYwF5hgZq8BE7zHAI8DbwCbgbuByzPaahHpmmXLEpdjvvWt1mAvjM49AaS9LpdrnHPnJ1g1Ps62Driiu40SkTRpEFU8uuJVJCqGDk3ca9+4UeWYPKW5a0TCbOtWKE4wv/9BB0FTU3bbI4GjkBcJo2TlmE8+Sb5e8orKNSJhceONicsxixa1lmMU8NKGevIiQdbcHCu7JKIau3RCPXmRIGrpsccL+N27NYgqXaaQFwmKJ59MXI45+eTWYP/sZ7PfNgktlWtEck3ntIuP1JMXyYV//ufEvfZnnlE5RjJGPXmRbNmxA44+OvF6hbr4QCEv4rdk5ZjmZviMPlCLf/TuEvHDTTclLsf84Aet5RgFvPhMPXmRTOkstFWOkRxQN0IkXS099ngBv22bBlElpxTyIt2xenXickxJSWuwf+5z2W+bSBsq14ikQue0S8ioJ++D6rX1jJn7FENmP8aYuU9RvbY+102SdIwdm7jX/vjjKsdIoKknn2HVa+u5rmodjU3NANQ3NHJd1ToAppYnmPdbguf996F378TrFeoSEurJZ9i85ZsOBHyLxqZm5i3flKMWSUpaeuzxAr6pSb12CR2FfIZtbWhMabkEwIIFicsxM2e2BnuyKX9DQGXE/BTud20ADSgqpD5OoA8oKsxBa/xVvbaeecs3sbWhkQFFhcyaWBquklQeDaKqjJi/1JPPsFkTSyks6NFuWWFBD2ZNLM1Ri/zREhr1DY04WkMj8L3Dlh57vIB//fXIlmNURsy+oHxyUshn2NTyYuZMG05xUSEGFBcVMmfa8Mj1lkIVGrW1iYMdWoP92GOz264sUhkxu4LUCVK5xgdTy4sjF+odhSI08qgc05l8KiMGQbJOULazQT156ZZE4ZDz0PjSlxL32u+7L7LlmM7kSxkxKILUCVJPXrpl1sTSdgN5kMPQ2LcPCpP855KHod5RS+8x1APlIRKkT04KeemWQIRGsnLM3r3Jgz8P5UMZMSiC1AlSyEu35SQ0qqrgy1+Ov66yEhYvzmpzROIJRCfIk5GQN7MtwB6gGdjvnKsws77AA8BgYAvwVefce5nYn+QhDaJKyATlk1MmB17HOefKnHMV3uPZwErn3HHASu+xSNf17Zt4EHXjxrwdRBVJhZ9n15wL3Of9fB8w1cd9SVTU17cG+3txPvi1BPvxx2e/bSIhlKmQd8CfzGyNmc3wlh3tnNsG4H0/Kt4vmtkMM6sxs5qdO3dmqDkSOi3BXlLy6XUtwa5eu0jKMhXyY5xzo4CzgSvM7PSu/qJzbqFzrsI5V9GvX78MNUdC4X/+J3E55p57FOwiGZCRgVfn3Fbv+w4zexgYDWw3s/7OuW1m1h/YkYl9Scjt3w8FBYnXK9RFMirtnryZHWpmh7f8DJwJrAeWApXeZpXAI+nuS0KspcceL+A/+CCveu1BmbhK8kMmevJHAw9b7CP3QcDvnHNPmNmLwINm9k3g78C/Z2BfgRT6KXf9smIFTJgQf92kSfDYY9ltTwBoyl/JtrRD3jn3BjAyzvJdwPh0nz/o9I82Dp3TnlCQJq6S/KAJytIUqil3/TR8eOJB1JqavCrHJBOkiaskP2hagzTl9T/aHTvg6KMTr1eof0qQJq6S/KCefJoCO+Wun1p67PEC/pNP1GtPQlP+SrYp5NOUN/9of/nLxOWY+fNbgz1ZPV4S3jkM0Bk34guVa9IUpNnmMs45+EySfoB6693SceIqDd6LnxTyGRCU2eYyJllv/N13oU+f7LUlD+iMG/GTyjUS87e/JS7HTJ3aWo5RwGdcXg/ei+/Uk893Oqc953TGjfhJPfl8VFmZuNe+erXOjsmyvBm8l5xQTz5f7NkDn/1s4vUK9ZyJ9OC95JxCPuqSlWP274cePRKvl6yJ3OC9BIbKNVH06KNdO6ddAS8SeerJR4XOaReRONSTD7vTT4/12OMF/LZtGkQVyXPqyYfRW2/B4MHx1114IfzmN1ltjkgu6D4OXZPXIR+6N4nOaRcBNBVEKvK2XNPyJqlvaMTR+iYJ3MRQP/1p4kHU558PdTlGt8GT7tJ9HLoub3vygZ4v5KOPoFev+OsOPji2PuTUE5N0aCqIrsvbnnwg3yQtPfZ4Af/xx7EeewQCHtQTk/Tk5X0cuilvQz4wb5K//CVxOeZnP2stxxQUZLddPgvkf7ISGpoKouvytlwza2Jpu3IBZPlNkueDqJqUS9KhqSC6Lm9DPidvki9/Gaqq4q+rr4cBA/zbd8Dk/D9ZCT1NBdE1eRvykKU3yfbt8LnPxV83aRI89pi/+w8o9cREsiOvQ95XeV6O6Qr1xET8l7cDr76oqko8iLpyZajPaReRcFJPPl2ffJJ8NkeFuojkkHry3fXd78Z67PECft8+9dpFJBB8D3kzO8vMNpnZZjOb7ff+fPXOO63lmAUL2q979NHWYO/ZMzftExHpwNdyjZn1AO4AJgB1wItmttQ5t8HP/XaU9kRkiU59HDAgduqjiEhA+V2THw1sds69AWBmfwDOBbIW8t2eI+Xpp2Hs2Pjrdu9Ofr9UEZGA8LtcUwy83eZxnbfsADObYWY1Zlazc+fOjDcgpTlSPvqotRzTMeCrqlrLMQp4EQkJv0M+3sni7UYjnXMLnXMVzrmKfv36ZbwBXZoj5Uc/ij8x2Mknx86ecQ7OOy/jbRMR8Zvf5Zo6YGCbxyXAVp/32U6iOVJO3v9u4guWtmyBQYP8bZiISBb43ZN/ETjOzIaY2cHAdGCpz/tsp91sdc7x5D0z2XLzZB685aL2G958c2s5RgEvIhHha0/eObffzK4ElgM9gHudc6/4uc+OppYXc/jmVxn/1bPjb/Dxx5GbxldEpIXvV7w65x4HHvd7P5+ybx9cdRXcfTfjO6577jk45ZSsNylVobsHrYgETvSmNXj4YZg27dPLH3oo/vKA0u3xRCQTojGtgXNw5pmxgdS2Qf4f/wF798bWhyjgQbfHE5HMiEZPfu1aePLJ2M/FxfDEEzBsWG7blCbdHk9EMiEaIV9eDn//O5SUJJ/HPUTy8fZ4GoMQybxolGvMYODAyAQ85N+NilvGIOobGnG0jkFUr9XcQCLpiEZPPoK6e3u8sPaGk41BhKH9IkGlkA+wVG+PF+YzcjQGIeKPaJRrBAj3GTmJxhqiPAYhkg0K+QgJc28438YgRLJFIR8hYe4NTy0vZs604RQXFWJAcVEhc6YND3yZSSToVJOPkFkTS9vV5CFcveFUxyBEpHMK+Qjp7hk5IhJdCvmIUW9YRNpSTV5EJMIU8iIiEaaQFxGJMIW8iEiEKeRFRCJMIS8iEmEKeRGRCFPIi4hEmEJeRCTCFPIiIhGmaQ0iIqx3hBIRfynkIyDMd4QSEX+pXBMBYb4jlIj4SyEfAWG+I5SI+EshHwFhviOUiPgrrZA3sxvMrN7MXva+JrVZd52ZbTazTWY2Mf2mSiK6P6qIJJKJgdfbnHPz2y4ws6HAdOBEYACwwsz+xTnXHO8JJD26I5SIJOLX2TXnAn9wzn0EvGlmm4HRwHM+7S/v6Y5QIhJPJmryV5pZrZnda2Z9vGXFwNtttqnzln2Kmc0wsxozq9m5c2cGmiMiIi067cmb2Qrgc3FWfR+4C7gRcN73W4BvABZnexfv+Z1zC4GFABUVFXG3SUYXAYmIJNZpyDvn/q0rT2RmdwOPeg/rgIFtVpcAW1NuXSd0EZCISHLpnl3Tv83D84D13s9Lgelm1tPMhgDHAS+ks694cnERUPXaesbMfYohsx9jzNynqF5b79u+RETSle7A60/NrIxYKWYLcCmAc+4VM3sQ2ADsB67w48yabF8EpE8OIhI2aYW8c+7rSdbdBNyUzvN3ZkBRIfVxAt2vi4CSfXJQyItIEIX6itdsXwSk6QNEJGxCHfJTy4uZM204xUWFGFBcVMicacN961Vr+gARCZvQTzWczYuAZk0sbVeTB00fICLBFvqQzyZNHyAiYaOQT5GmDxCRMAl1TV5ERJJTyIuIRJhCXkQkwhTyIiIRppAXEYkwcy7l2X19Y2Y7gbe6uPmRwDs+Nieo9Lrzi153/unOax/knOsXb0WgQj4VZlbjnKvIdTuyTa87v+h1559Mv3aVa0REIkwhLyISYWEO+YW5bkCO6HXnF73u/JPR1x7amryIiHQuzD15ERHphEJeRCTCQhnyZnaWmW0ys81mNjvX7fGLmQ00s1VmttHMXjGz73jL+5rZk2b2mve9T67b6gcz62Fma83sUe/xEDNb7b3uB8zs4Fy3MdPMrMjMlpjZq95xPzUfjreZfdd7j683s9+bWa8oHm8zu9fMdpjZ+jbL4h5fi7ndy7laMxvVnX2GLuTNrAdwB3A2MBQ438yG5rZVvtkP/Jdz7gTgFOAK77XOBlY6544DVnqPo+g7wMY2j28GbvNe93vAN3PSKn/9DHjCOXc8MJLY64/08TazYuAqoMI5NwzoAUwnmsd7MXBWh2WJju/ZwHHe1wzgru7sMHQhD4wGNjvn3nDOfQz8ATg3x23yhXNum3PuJe/nPcT+wRcTe733eZvdB0zNTQv9Y2YlwDnAPd5jA74ILPE2idzrNrPPAqcDvwJwzn3snGsgD443sXtbFJrZQcAhwDYieLydc08D73ZYnOj4ngv82sU8DxSZWf9U9xnGkC8G3m7zuM5bFmlmNhgoB1YDRzvntkHsPwLgqNy1zDcLgGuAT7zHRwANzrn93uMoHvdjgZ3AIq9MdY+ZHUrEj7dzrh6YD/ydWLjvBtYQ/ePdItHxzUjWhTHkLc6ySJ8HamaHAQ8BVzvn3s91e/xmZpOBHc65NW0Xx9k0asf9IGAUcJdzrhz4kIiVZuLxatDnAkOAAcChxEoVHUXteHcmI+/5MIZ8HTCwzeMSYGuO2uI7MysgFvD3O+eqvMXbWz62ed935Kp9PhkDTDGzLcTKcV8k1rMv8j7OQzSPex1Q55xb7T1eQiz0o368/w140zm30znXBFQBpxH9490i0fHNSNaFMeRfBI7zRt4PJjZAszTHbfKFV4f+FbDROXdrm1VLgUrv50rgkWy3ze8cM74AAADwSURBVE/OueuccyXOucHEju9TzrkLgFXAV7zNovi6/wG8bWal3qLxwAYifryJlWlOMbNDvPd8y+uO9PFuI9HxXQpc5J1lcwqwu6WskxLnXOi+gEnA/wNeB76f6/b4+Dq/QOzjWS3wsvc1iVh9eiXwmve9b67b6uPf4AzgUe/nY4EXgM3A/wE9c90+H15vGVDjHfNqoE8+HG/gR8CrwHrgN0DPKB5v4PfExh2aiPXUv5no+BIr19zh5dw6YmcfpbxPTWsgIhJhYSzXiIhIFynkRUQiTCEvIhJhCnkRkQhTyIuIRJhCXkQkwhTyIiIR9v8BbPFPehbGb/gAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "\n", "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", " \n", "slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)\n", "print(\"slope: %f intercept: %f\" % (slope, intercept))\n", "print(\"R-squared: %f\" % r_value**2)\n", "print(\"p-value for the slope: %f\" % p_value)\n", "\n", "# Plot the data along with the fitted line:\n", " \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", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 133, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "parameter estimates : [5.33171291 0.18878757 0.21927505 0.08021771]\n", "parameter standard deviation : [0.26920328 0.02583608 0.29269083 0.02316083]\n", "\n", "relative estimation error : [0.06634258 0.05606215 0.56144989 0.1978229 ]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3RUxR7A8e/sZjebSjoQQggl0jsBpEgHUUCkiIooFhAVu9grivrkqTyxoFJERUSko9J7J/QuLUAIEFJI32TLvD8uBEJ6spsQmM85OYRbf8mB386dO/MbIaVEURRFufnoyjsARVEUxTlUglcURblJqQSvKIpyk1IJXlEU5SalEryiKMpNyqW8A7hWQECADAsLK+8wFEVRKowdO3bESSkD89p3QyX4sLAwIiMjyzsMRVGUCkMIcSq/faqLRlEU5SalEryiKMpNSiV4RVGUm9QN1QevKLc6i8VCdHQ0ZrO5vENRbjAmk4mQkBAMBkORz1EJXlFuINHR0Xh5eREWFoYQorzDUW4QUkri4+OJjo6mZs2aRT7vpumikZlrsSeMwB5/H/bUKUh7enmHpCjFZjab8ff3V8ldyUEIgb+/f7Gf7G6KFrw99WtI+xFkhrbBchhpngv+cxDCVL7BKUoxqeSu5KUk/y4qfAte2uIh9furyR0AM1ijkenzyi0uRVGU8lbhEzyW3SDyeumQAZmryzwcRanovvrqK+rXr8/QoUNZuHAhn376KQDz58/n4MGDeZ5z8eJF2rRpQ/PmzVm/fn1Zhptt0qRJ/PzzzwUeM3z4cP78889c29esWUOfPn0KPHf58uW0bNmSxo0b07JlS1atWlWqeK83ffp0wsPDCQ8PZ/r06Q65ZsXvotH5YLULpu/qxENNNuLqYr2yA/R5zt5VFKUA3377Lf/880/2y7x+/foBWoLv06cPDRo0yHXOypUrqVevXrESk81mQ6/XOyRmq9XKqFGjHHKt/AQEBLBo0SKCg4PZv38/vXr14uzZsw65dkJCAh988AGRkZEIIWjZsiX9+vXD19e3VNet+C14Q3N2xDRg3Lr+vLb8fq4uUGVEuA8tz8gUpcIZNWoUJ06coF+/fnz55Zf89NNPjB49mk2bNrFw4ULGjBlDs2bNOH78ePY5u3fv5tVXX+Xvv/+mWbNmZGRkMHPmTBo3bkyjRo147bXXso/19PTk3XffpU2bNmzevDl7+6FDh2jdunX236OiomjSpAkAY8eOJSIigkaNGjFy5EiurELXuXNn3nzzTTp16sT//vc/3n//ff773/8C8OOPPxIREUHTpk0ZOHAg6elXB12sWLGCjh07ctttt7F48eJcv4O0tDQee+wxIiIiaN68OQsWLACgefPmBAcHA9CwYUPMZjOZmZml/p0DLF26lB49euDn54evry89evRgyZIlpb5uhW/BC6GjbY1nWFT9ZQYdeYWJfn14ru0q8HofYWhAYmwSmxdGYrPauL1vSwKq+Zd3yIpSZEO+35xrW58mVRl2exgZWTaGT9uWa/+gliEMblWdhLQsnvp1R459s568vcD7TZo0iSVLlrB69WoCAgL46aefAGjXrh39+vWjT58+DBo0KMc5zZo1Y+zYsURGRvL1118TExPDa6+9xo4dO/D19aVnz57Mnz+f/v37k5aWRqNGjRg7dmyOa9SvX5+srCxOnDhBrVq1mDVrFvfddx8Ao0eP5t133wVg2LBhLF68mL59+wJw6dIl1q5dC8D777+ffb0BAwYwYsQIAN5++22mTJnCs88+C2gfHmvXruX48eN06dKFY8eO5Yhl3LhxdO3alalTp3Lp0iVat25N9+7d8fDwyD5mzpw5NG/eHFdXVwCGDBnCkSNHcv0+X3rpJR5++GFmzJjB+PHjc+2vU6cOf/75J2fPnqV69erZ20NCQhzydFDhEzyASMmi8cX9fBWylCc396Vm6HPc0zyMVTPX8/kTk9AJgZSSSS/9xIjPHqL/6LvKO2RFuWlt376dzp07ExiodZEOHTqUdevW0b9/f/R6PQMHDszzvPvuu48//viD119/nVmzZjFr1iwAVq9ezWeffUZ6ejoJCQk0bNgwO8EPGTIkz2vt37+ft99+m0uXLpGamkqvXr1y3Een0xEeHk6tWrU4fPhwjnOXLVvGwoULs58GzGYzp0+fpn79+gAcOHCA1157jWXLlmWfcyXW/AwdOpShQ/PvUchrbWxHjKa6KRI8Ne+ARoPoeWgubUULXvo5i6l3jcMcl4zdZs9x6I+v/kqrns0IuS24nIJVlKIrqMXtZtQXuN/Pw1hoi90Z8kpWV5hMpnz73YcMGcLgwYMZMGAAQgjCw8Mxm808/fTTREZGUr16dd5///0cY8GvbVVfa/jw4cyfP5+mTZvy008/sWbNmux91yfO6/8upWTOnDnUrVs313Wjo6O59957+fnnn6ldu3aO2EvTgg8JCckRY3R0NJ07d87zZyuOit8Hf5ns+SGZmXZejPkvwbOWkH7hUq7kDmCz2lkza2M5RKgoFZuXlxcpKSmFHtemTRvWrl1LXFwcNpuNmTNn0qlTp0LPq127Nnq9ng8//DC7ZX4lmQcEBJCamprnCJi8pKSkULVqVSwWCzNmzMixb/bs2djtdo4fP86JEydyJfJevXoxceLE7A+qXbt2AVp30N13380nn3xC+/btc5wza9Ysdu/enevr4YcfBrQWfF77r/w8vXr1YtmyZSQmJpKYmMiyZctyPHWU1E2T4A/sS+bXDeG0qXmBLt7aI1dK/ZrYTK45jpN2O1aLrTxCVJQK7f7772f8+PE0b948x0vW61WtWpVPPvmELl260LRpU1q0aME999xTpHsMGTKEX3/9Nbv/3cfHhxEjRtC4cWP69+9PREREka7z4Ycf0qZNG3r06EG9evVy7Ktbty6dOnWid+/eTJo0CZMp52TId955B4vFQpMmTWjUqBHvvPMOAF9//TXHjh3jww8/pFmzZjRr1ozY2NgixVMYPz8/3nnnHSIiIoiIiODdd9/Fz8+v1NcVBT1OlbVWrVrJki74Mf/rf5j6+nQGtvqXBdtCSTD4cGrEAFzPx1Ft1jJ0Ni2pu7oZmbDhI+o0L3o9B0UpK4cOHcru61WU6+X170MIsUNK2Sqv42+aFnzVmkGgN/DrutqkmA0YUlKp/Pd6zCGVudCnI+gEru5G+j1zp0ruiqLcEm6Kl6zSnkirHiF4+3mRmZ5F1UopvHbvPr5cnIrY5MO5ds3wr1uFsYObcVvL2oVfUFEU5SZQoVvw0haDPf5+ZGwHREI3pm46zp0PVyHDaiLYN4PXBp/gt2+G8sjtNdiCiYwgNQZeUZRbR4VtwUtpRcbfD/ZYQBst4yJO8vzHFxk14S/E3n+ovfxFyFrLu30fpmfDKjSt7lO+QSuKopShituCz1wPMoUryT2btOIqlmBq9yjUaA/L30WffpH2dQIA2HYyga0n4ss+XkVRlDJWcRO87SxIax47zGA7BUJAnwlgSYdNEwGw2yUfLDrAEz9Hcvh8ctnGqyiKUsacmuCFEFFCiH1CiN1CiJKNf8yPoTF5hi/cEcaW2veBt8HDC6CrNo5VpxN8P6wl7kY9j0zdRnSiWvVJUa6nygXnzdnlgu+88058fHwKjaM4yqIPvouUMs7hVzU0AWNzyNoJXJm6bARdZTBdMwOsRjvtT3My6PSE+How/bHWDJ60mYenbuPPUe3w8zA6PDxFqahUueC8ObNcMMCYMWNIT0/n+++/d9g1K2wXjRAC4fsDeI4CXYiW2N2HIvxnI8R1CducDN/eDqvGAVCvijdTHongbGIG0zdFlX3winKDUuWCy6dcMEC3bt3w8vJy2PXA+S14CSwTQkjgeynlD9cfIIQYCYwECA0NLdbFhTAiPJ8Gz6cLPtDkDbf1hK3fQaMByGpNiKi6nbmPnqdesAdS1kGICvtZp9zMpt2de1vD/tB6BGSlw4zBufc3exCaD4W0ePjj4Zz7Hv2rwNupcsHlUy7YWZyd4NtLKWOEEEHAciHEYSnlumsPuJz0fwCtVIHTIun+Pvy7FLlgFHKACXRp1PfIhGQj0ecaMP3ga7xxVxP0OrXgsaKUhioXnFth5YKdxakJXkoZc/nPWCHEPKA1sK7gs5zEVAn6TED8NhgiPZERlz+NpZW1x0xM3hBNilnw6cDGalV75cZRUIvb6F7wfg//QlvszqDKBV9107bghRAegE5KmXL5+57A2EJOcyoZ3gnC3eBiFkh3bSgl8GCTDZxPC+TrrXdQyd3AG73rqSSvKNcpTrng559/nri4OHx9fZk5c2Z290hBilou+PouorxcXy64WrVq2ftmz57NI488wsmTJ7PLBW/ZsiV7/5VywRMnTkQIwa5du2jevHmh5YILUl4teGd2PFcGNggh9gDbgL+klKVfZLBUJLKTD7K3T3Zyv+LF21cxrG0Nflh3gm9WH8vnfEW5dalywc4rFwzQsWNHBg8ezMqVKwkJCWHp0qWlvuZNUy64qOzxQ8GyA1IscMECdUyAAdwGgtcHvDJ7D0cupDDnqXaYDI4ZwqUoRaXKBSsFKW654Apbi6akRKVPkPH3ISLPwNFUpJ8nBIYgvF5C6ASfDWpChsWGyaBHSqm6ahRFqbBuubGBwiUUEbQaenwMJg/EWiPCdx5CpxUic9Hr8DIZMFtsjPg5kkV7Yso5YkVRlJK55RI8gBBuiIBHEH0mIS6cRGz8OtcxUkJyhpUXZu1m6YHz5RCloihK6dySCT5bg37QaBCs/Q+c25tjl5tRz9RHI2hcrRKjf9vJ6iOOe5miKIpSFm7tBA9w13hoPRJ8w3Lt8nR1YfpjrbmtshdP/rKDjcccX1JHURTFWVSCd/eDOz/RyhnY7bl2V3Iz8MvjbWgR6kMlN0M5BKgoilIyKsFfkRgF398BURtz7fLzMDJzRFsaVasEwNlLGWUcnKKUn2uLeOWloDLCSvlSCf4K9wDISoH5oyAz92y9K8Mlp2+KoscXa9kelVDWESrKDUkl+BuXSvBXuHrCvd9DUjT881q+h/VuVIUq3iaGT92WI8mnXkpj+S9r+XvySuLOqiUBlbIxf9dZ2n+6ipqv/0X7T1cxf5dj6pOPGzeOunXr0r179+waK3mV4M2rjHBBpXqVsqUS/LVC20KHl2D3DDgwP89DgrxNzBzZlsreJh65nOQ3L4rk/mojmfjMZL57cRqPhD/LnAm560wriiPN33WWN+bu4+ylDCRa1+Ebc/eVOsnv2LGD33//nV27djF37ly2b98OaCV4t2/fzp49e6hfvz5TpkzJLiM8fvx4du/eTe3atfM8TikfKsFfr/PrENwCIqdqg+HzUNnbxO8j21KlkolHp23jg+HfkJmRRUaqGXNaJllmC9PemsnJ/afLOHjlVjJ+6REyLLYc2zIsNsYvzV3VsDjWr1/Pvffei7u7O97e3tkrOu3fv5+OHTvSuHFjZsyYwYEDB/I8v6jHKc6nEvz19AZ4YCYM/TNXQbJrBXmb+H1EW4YFGTBaci/+bcmysnJG+axNqdwaYvJ52Z/f9uLIq0TH8OHD+frrr9m3bx/vvfdejrK9JTlOcT6V4PPiVQVcjJBxCY78k+9hQd4mGrtq9aPTalYjrUbV7H3SLrFkWsoiWuUWFezjVqztRXXHHXcwb948MjIySElJYdGiRUDuErxXXF9GOL/jlLKnEnxBVn0Esx6CszvyPaR172bYbHYS2jfj3MDupNXS6k67uhvpOKBNWUWq3ILG9KqL23UVT90Mesb0yr1QRXG0aNGCIUOG0KxZMwYOHEjHjh2B/EvwXl9GuKBSvUrZuuXKBRdLRiJ810Hrthm1HlzzXhD3zy8XM/mjOZzs14XMAB9Cl27i3tvDeGHSk6oapVIsxS0XPH/XWcYvPULMpQyCfdwY06su/ZtXK/xEpUJS5YIdyc0XBk6Gn+6Cv16BAd/nedigF/vQskcT/vp1A7/Z7UTf3ZFa9zVTyV1xuv7Nq6mEruRLddEUpsbt0Ok12Ps77P0j38NqNgpl9KcPsuzDPrSu6ceO04llGKSiKEpuqgVfFHeMgaxUCOtY6KGeri789GhrjHrtszMpw6Jq2CjFohaaUfJSku50leCLQqeHnh9p39ttYLeCi2ueh26cv43Jr8/g3MkL+NQJ5t8B3Xm4Ux2e7xZe6H9aaT2GTP4AsiJBuILbQITXGIQwFXjetex2O2lJ6bh7u+W7er1y4zKZTMTHx+Pv76+SvJJNSkl8fHyu9WMLoxJ8cdisMGMg+IRCv4m5dm+Yt5VPH/qKzIwsAOIOR2Pfe4wJFkmK2cpbd9VHp8v7P620XUTGDwGZCkiQ6ZD+B9IahfAr2kzAhd8u4ad3Z5GRasbVzcgDb9zLfWPuUYmiAgkJCSE6OpqLFy+WdyjKDcZkMhESElKsc1SCLw69izbLdcMXUKMDNB2SY/eUN3/LTu4AQkoCFq3HJGAKcCndwn8GNsZFn/vVh0yfATITuPYxLBOytiOtxxAudQoMben01fzw6q9kpmcCYM2y8uvYP9Eb9Ax6sW9Jf2KljBkMBmrWrFneYSg3CfWStbi6vAWh7WDxCxB7KMeuc8cv5DpcAJ4L1/NCJz1zdkbz1apjeV/XcgDIyr1duID1eKFh/fLB7OzkfoU5PZPfxs0tUd+doigVn0rwxaV3gUFTwegBs4blKC0cFOqf5yn+QRaea/4aX/ZP4omO+bTODA0BY+7t0goutQsNKz4m71E7qYmp2Ky2PPcpinJzUwm+JLyraknelgXJMdmbh3/4AK7uOZO0q5uNYa+cBzK4J+w/eLlaMVtsvDF3H7HJV2t0CPcHQRjR2vzZZ4OxZaHdMwDV6wbnuT2wegAuBtUTpyi3IpXgS6rmHTA6EgKvTgvv+kAHnv9uJP5V7YDEL8jCqLEx3PXQ5brxQg/WY/x7IYUFu88ycNImTsalXd4VhPCbBYYIQAfCDdwGIXy/LVI4I8c/jKvbdR8u7kZGfjbMET+toigVkCpVUFo2Cyx/Fxrco9WTB+wJT2DLWEfuUYquiMDlCH0V9py5xKM/aXW2pw6PoFl1n+yjSjoOeufKfUx98zdOHz5LcK0ghn/4AG37tCzpT6YoSgVQUKkCleBLy5wEP3SGrHR4ci14VUFmbkUmjgCuLZNqBGMEOr9p2VtOXEzlkWnbiEvJ4oeHW9IxPLCso1cUpYIrKMGrLprSMlWCIb9CZjLMHg42C8K1DXi/B8IThAdacm+L8PlfjlNrBXoy56l2tArzpVopS7wqiqJcz+kteCGEHogEzkop+xR0bIVswV+x70+Y8zi0Hgl3jQdAyiywngKdL0IfUOglpJQsO3iBng0qq8lJiqIUSXm34J8HDhV6VEXXeBDcPhp2/gKJpwAQwogwhBcpuQOsOBTLk7/s4OXZe8iy2p0ZraIotwCnJnghRAhwNzDZmfe5YXT/QKsb71ujZKfXD+KlHrcxd+dZHv1pG0kZakUoRVFKztkt+AnAq8Ct0RzVu0BAuPb97pmQfK5YpwsheK5bOJ8Pbsq2kwkM+m4TZxLSnRCooii3AqcleCFEHyBWSpn/enfacSOFEJFCiMibpsBSynn462WYNRQsxV8AeWDLEH5+rA0JaVlExac5IUBFUW4FTnvJKoT4BBgGWAET4A3MlVI+lN85Ffol6/UOLdYSfOPBMOBHKMFL07RMKx6u2izUk3Fp1AzwcHSUiqJUcOXyklVK+YaUMkRKGQbcD6wqKLnfdOr3ga7vwL7ZWvXJEriS3Dcfj6fb52v4auVRVThMUZQiU+Pgnanjy1oLftVHEHe0xJdpUcOH/s2q8cXyf3lh1m7MFlU8TFGUwpVJFSop5RpgTVnc64YihLYwSJP7r758LQFXFz2f39eU2kGejF96hKj4dH4c1pIg7+Kt7qIoyq1FteCdzeAG4d21709tgqSzJbqMEIJnutRh0kMtOXohheWHcteeVxRFuZaqI1tWMlNh1kPgHQyPLgFXzxJd5s5GVWgc0ongSlrrPTbFTJCXaskripKbasGXFVdPuPcHuHAA/nxMW9+1hKr5uCGEICoujW7/Xcun/xzGZlcvXxVFyUkl+LIU3l2rU3N0KfzzKpRyREywjxv9mgUzae1xHp++Xc18VRQlB5Xgy1rEE9D+eYicAkf+LtWljC46xt3bmI/6N2LD0Tju+XoDR86nFH6ioii3BJXgy0O392HAZLitt0Mu91DbGswc2Za0LBsztp5yyDUVRan41EvW8qDTQZPB2veJUdq6rjXaleqSEWF+/PVsByq5GwA4l5RBoKcrLvrSfYZLmQUYVPliRamAVAu+iKT1GPaUCdiT/4PM2uW4Cy96HmbcB+f2lvpSQd4mXF30mC02HvhhC8OmbCMuNbNE17JnLMMe2wV5oTEyNgJ76g9qFq2iVDAqwReBPe0XZNwASPse0qciE4djT3rPMQnvnm/B5A0zBmmteQcwGfQ806UOO08n0uerDew8nVis82Xmekh6BexnAQkyGVK/QaZ945D4FEUpGyrBF0LaYiHlP2jrq9rQEl4GZMwHiwNa8pWqwUNzwZoJvwyAVMdU1Bzcqjpzn26HwUVw36TNTN1wssgfSDJ1AjnXkwXIgLQpSKlG6ihKRaESfGEy14LQ57HDjDT/45h7BNWDB//Q+uJXfuCYawINgyuxeHRHOtcN4u9954o+Vt56Ou/t0gr2JIfFpyiKc6mXrIURLiDzesEoQBgdd5/QNvDwAqjcwHHXBCq5G/jx4ZakZFpx0euIT83kQnImDYK98z/JpQ5Y8ijjL4yg83FofIqiOI9qwRfGtSt5L0hlRJj6OvZeoW3A1Usra7B2PNgc0x0ihMDbpI2uGff3Ifp/u5Hftp7Ot8tGeL2EVsL/2o1u4DkaIVSbQFEqCpXgCyF0laDS54Ar4I6W+FzB8zmEoZ5zbnp0Gaz+COY/BXbHlgZ+8676tKnpx5vz9vH877tJMef+EBHGCITv9+BSHzCCLhi83kG4D3doLIqiOJfTVnQqiRt5RSdpTwTzKiALXDsj9FWde8P1X2j98S2HQ58JJVoRKj92u+TbNcf4Yvm/hPq5M/mRVtQJ8nLY9RVFKTsFreiknreLSOh8wX1g2d2w40uQlQrrPwejJ/T8yGFJXqcTjO4aTpta/nz89yH8PFwdcl1FUW4sqovmRtb1HWgzCvb+AWmOX5A8IsyPuU+1w8/DiMVm55N/DpV4YpSiKDceleBvZEJAr0/gyXXgGVTq6pN530J7Kth3NolpG6Po/b/1rD/q+A8TRVHKnkrwNzqdDryrasl92duw4Uun3KZFqC8LnmlPJTcDw6Zs46PFB8m0qrVfFaUiUwm+opB2SDkPK96HjV855Rb1q3qzaHQHhrWtweQNJ3lp1h6n3EdRlLKhXrJWFDo93Ps9SBssf0fb1v45h9/Gzajnw/6N6Fw3kMqXF/XOtNow6HTodKqipKJUJCrBVyR6F62OPGhJXgho96xTbtWtfuXs799feJDoxHTGD2pKlUpq/VdFqShUF01FcyXJNxkCPqFlcstG1byJjEqk14R1LNwTUyb3VBSl9NREp5tB7CEIrOfQyVDXO3ExlZdn72HX6Uvc3aQqH93TCF8PB9biURSlRAqa6KRa8BXd2R3wXXtY9aFThlFeUSvQk9lP3s6YXnXZeCyO1Eyr0+6lKIpjqARf0VVtDi2GaTNel7zh1CTvotfxTJc6rH+1C9X93JFSMnn9CZLSVY14RbkRqQRf0el0Wq2atk/D1u9g4bMOL1B2Pa/LlSkPxCTzyT+H6TlhLSsPXXDqPRVFKT6V4G8GQkCvj+GOV2HXL3BwfpnctlG1Ssx/uj0+bkYenx7Ji7N2cyk9CwC7Pa8Sy4qilCX1kvVmc2IN1Ozk1Beu18uy2vl69TG+XX2MysJO0PdzSE1IpVp4VZ76cjht7mpRZrEoyq2mXF6yCiFMQohtQog9QogDQgjHrUWn5K9WZy25X/wXZj4IGcVbcLskjC46XupxG4/rUjH8uZrUhFSkTsepmEt8eN/n7Fl7oMDzbTYbx/dEcebIWccsZK4oCuDcLppMoKuUsinQDLhTCNHWifdTrpVwHI4th2l3ayUOnCwr08Km/y3EcFRbzzWxdSNOPd6f2Do1+Om9P/I9b+eKvdwfPJIXO77DUy1f5bH6L3DqULTT4y3I8T1RrJ29udzjUJTSctpMVqk1xVIv/9Vw+Us1z8pK3d4wdDb8PhQm94BhcyEg3CGXllJC5lJk2q8gU8HUm0sJd+YYwON5JIr0sGBi72zPivMXOX4xldqBnjmuE3smjvfu/Qxz2tUSxWePxvBKl/eZeWYSLoaynWidnpLBW3d/zNGdJ9G76LBZbDTp3JD3547B6Goo01gUxRGc+pJVCKEXQuwGYoHlUsqteRwzUggRKYSIvHhRlal1qFqdkQ/PRlpTkVN6aBOiHECmfIpMeg0s28B6EFK/xt80Cle3qxnemJhMtd+XELRkI+YAX3pPWM8fkWdyXGfpT6uxWXKO+JESsjKy2L5kt0NiLY5vnp/Kke3HyEzPJD05g8yMLPas2c/0d38v81gUxRGcmuCllDYpZTMgBGgthGiUxzE/SClbSSlbBQYGOjOcW4qUEnvy50iXUch7PKB6FnbdVOz2rNJd13Ye0meAzLhmaybCHsMr31bH5H51dSgBBB07xdS76tCrURXqVtaWBbTbtQ+Ci2fisWTlnjBls9lJPH+pVHHm59iuk/z42i/88OrPHNl+LHu73W5n9cwNWK6bwJWVYeGfKaucEouiOFuZPANLKS8JIdYAdwL7y+KetzqZ/gtk/AyYwRtkV0+wLoAzi7BfaIuImIIQ+uJf2LILhAHk9R8UGbTpnsLjnw5l5ifzSLqYTI2GIYz6/BGad21Ix2uOfGv+PrKsks6dGrF65kbMaeacsUtJo471ix9bIX4ZO5tZn83HYrYggYXfLuWe0b0Z8elDSLvEasl7/kBmulrlSqmYnJbghRCBgOVycncDugP/cdb9lOukTb6ula0Re5MRO+YhY2PgriXaRKni0Pnns0MP+ir0H92b/qN753u6lBI/DyM/rDvBCqOe6t1aYly+nawM7QPD5OFKhwFtCK1XrWbqW+cAACAASURBVHhxFSL66Dl+/3QeWears24z07NYMPEfug/tSM3GNajXug6HthzNcZ7QCVr0aOrQWBSlrDizi6YqsFoIsRfYjtYHv9iJ91OuZc+7i0O29EA2cENEbkX++TBYcn8IFMjQCkQlcv/TMSDcHyj0dCEEY3rV45/nO1K/qjf764WT/OIDVG7fgPptb+O5b0YwZtozxYupCLYsisxzCKYly8rGBdsBeGHSk7h7u2G4/ELVaDLg6ePBU1884vB4FKUsFNqCF0KMBmZIKYs1oFpKuRdoXtLAlFIyNAJLHpPGdALZ0QvpbUJsWQxJfeCBmdqar0UghA78fkYmjgLbWRA6QAfenyBc6hQ5vDpBXswc0ZYFu2P4csW/fLbodYJ93Ip8fnG5GF0QeTyt6PQ6DEYtoddqUoMpByeweNIyTuw5Rb02dbhrRHd8Ais5LS5FcaZCZ7IKIT4C7gd2AlOBpdJJs1HUTFbHkZZ9yPiHgHxa6MIDEf8uYuU4GL4YvKoU/x7W4yDTwaUeQpR8GKHNLtHrBFJKRs/cxe21/HmgdSh6B64gFReTwCN1RufoogGtlT7l4ASqhBXtA05RbjSlmskqpXwbCAemAMOBo0KIj4UQtR0apeJQwtAY4T8bjO3QxrNcyw08nkI0uBee3qIld5sVTucaxVrwPVxqa/cpRXIHshN5aqaV+NRM3p6/n74TN7A9KqFU171WQLAfL01+CqPJgMnDFZOHK0aTgdETH1fJXblpFbkWjRCiKfAo2kiY1UBbtH71Vx0VjGrBO4fM2o1M+RQsh0AfAB6jEG6DENfWq9k0EZa9Az0+gHbPlWktmxyxSslf+84x7q9DnEsy07dpMO/3bYC/p2vhJxdBUlwymxdq/fFt+7TEt7KPQ66rKOWloBZ8UfrgnwMeAeKAycAYKaVFCKEDjgIOS/CKcwhjM4R/IZN1Wj2uLR6y/F24cBD6TgCD8/rE8yOEoE+TYLrWC+L7tSdYsPssJkMJhnNeQ0rJ3nUHWfHzWux2SdcHO9Cie5OcH3CKchMqSh/8WGCKlPJUHvvqSykdMz0S1YIvd1LCuvGwehwEN4chM6CSY4crFpfFZseg15FltTN82jYGtQyhf7Nq6IrRPz/ppZ/468cVZKZnIqU2FLPrAx148YdRToxcUcpGafvg380ruV/e57DkrtwAhIBOr8L9v0FSNKTFlndEGPTaP9G41ExSM6289Mce+n+7ka0n4ot0ftSBMyz+fjnmtMzsWjnmtExW/rY+x0xWRbkZqQU/lNzq3Q3P79Va8QAn1zt1KcCiCPZxY/7T7fl8cFNikzMZ8sMWnvwlkqSMgpcL3L5kNzZb7sVHLGYL2/7Z5axwFeWGoBK8kjeju/bn8VUwvQ/MexKy0ss1JJ1OMLBlCKtf6cwrPW8jMd2Cl6v2GsmSRxIHMLkb0bvk/meuN+gxeZicGq+ilDeV4JWC1ewMXd6CvX/AlB4Qf7y8I8LNqGd013BmjWyLTie4lJ5F5/FrmLjyKOnXFS/rOKhtnkWqhU5H5yHtyihiRSkfKsErBdPptH75oX9C8ln4vhMcujEqTlwZBWO22GlUzZvPl//LHZ+t4Zctp7Jb9D6BlXjr9xcxebji7u2Gu5cbru5GXp32DIEh+dXVUZSbg1qTVSm6S2dg9nC4/WloNLC8o8llx6lE/vPPYbZFJRDm787CZzvgbdImYWWkZrBj+V6kXdKiRxM8vN3LOVpFcYyCRtGoBK8Uj90Gusvj0g8ugKpNwTesXEO6lpSS1Udi2XIigTfv0koOHz6fTN3KXmrcu3JTKtVEJ0XJ4Upyz0qHv14BaybcMxEa3FO+cV0mhKBrvcp0rVcZgKi4NO7+agNNQirxSs+6tK8TUOg1pJTsXr2fzYsi8fB2o9tDnQgJr+rs0BXF4VQLXim5xCiY/SjE7ISWj0Kvj6+OvrlBWGx2/twRzVcrj3IuyUzbWn683LMuEWF+eR5vt9sZd/8Etv2zE3NaJnoXPXqDnhcmjaTHsE5lHL2iFK5UE50UJV++YfDYUq12zY5pMLkbWMyFnlaWDHodD7QOZfUrnXmvbwOOxabx4I9buJiS9ypNW//amZ3cAWxWG1kZWUwY9QNpSWllGbqilJrqolFKx8UIPT+EWp3h/F4wOG9seWJsEunJ6VStVRldMVeiMhn0PNq+JvdHhLLjVCKBXlrxsv8uPUKXekG0rOELwOrfN2Yn92u5GPTsXLmfjgPalP4HUZQyohK84hh1umlfAFEbYcOXcM/XxaozL23nkcnjIHOttu6rqT/C62WS4i18/OD/2L/hMHq9DjcvEy/9+BRt+7QsdphuRj0dwrV++Ispmfy27TRfrz5Gx/AAnusWjtFkQIi8J+4ajOq/i1KxqD54xfH2/gELn9OqUfb5Ehr2L/QUaU9FxvUEewJwZVaqEQwNGH1nGCf3ns6xKLaruysTN4+jZuMapQo1LdPKr1tO8cO6E8SnZdHE34T5iz+QF3LWonf3dmP2+ckYTcZS3U9RHE31wStlq8l98OQ68AmF2Y/A3JGQkfcasVfIjIVgT+NqcgfIwp55GDfXf3MkdwBLpoW5//u71KF6uLrwZKfabHitK+/0aYDZxYUBj3fGaDKgC/DG5GnC5OHKB/NeVcldqXDUM6fiHIG3wRMrYP3nsPYzqN4aIp7I/3jrPvJaXlBKSc0GZvZuyrlqlN1m59zJCw4L182o5/EONXmsfRhCCPo+1oWB03Zi0Ate6F2fxhGle1JQlPKgWvCK8+gN0Pl1eGojtHxM23ZmO5iTcx/rUhfI/YJWp9dz+t/cSwIaTQZadGvs4ICvlj+oElaZV/o3weDpxsvzDtD18zX8tvU05uueJBTlRqYSvOJ8QfW1mjaWDPj9Afj2dji2Ischwu1eEK7kXD/WgHCpTq0WAzB5XF2yz8Wgx9PXkz6jejotZL1OMKhlCMtf6sSkh1pQyc3Am/P2sexg0Z8aMlIzmDNhMS93fo+xgz9n77qDTotXUfKiXrIq+bo6qmU1CD2Y7kJ4vYHQeZf8otGRMP8piPsXmj4IvcaBuzbpSFpPIpPeAUskoAdTD4T3+yAqseKXdcyZsJjUxDRu79uKB98aUKbrqUop2XwinogwPwx6HdM3RXEmIZ3HOtQk2Cf30oYZaWZGt36dC1EXyczIArQXw499/AADnru7zOJWbn6qFo1SbFJmIC/2AHs8cKVbwgAutRD+C9CW5C0hi1lbGnDDl+Dur3XheAZdc28roCvdPa6RlpyOJdOCT2Alh1wPYOyig0zfHIUA+jSpyhMda9Go2tXrz/vqL6a88Vt2cr/C1c3IrHM/qmJnisOoUTRK8WX8DfZUriZ3AAvYzkDWltJd22CCbu/AyDXQ8pGryf3ygiJCuDgkuSdeuMTrvT5kUNDjPFB9FI/Vf55DW4+W+roA7/ZtwLpXuzC8XRjLD16gz8QN/GfJ4ez9G+dvz5XcAVyMLhx2UAyKUhiV4JU8SethII8VnKQVrA5ay7RqE+j6tvb9xSPwZUPY8h3YrAWfVwRSSl7p+j67Vx/AmmXFmmXlzJEYXusxlrizRVvPtTDVfNx4u08DNr/Zjbfuqk+XutoHVXRiOmdrhSLzmBhlt9nx8vN0yP0VpTAqwSt5Ei7hQO6+ZYQLuNR0/A0NblCtBSx5HSZ3hbM7SnW5AxsPc/FMPDZrzlEvVouVv35ckc9ZJeNtMjDijlq0rqm9S1iy/zybAypz4qn7uNglAkslLaELncCvqg/hLWo59P6Kkh+V4JW8me4GnTs5/4m4gK4KGNs7/n4+odqqUYOnQ2os/NgN/n61xIt9n4+6mOd2S6aV6CMxpYm0UE90rMXcp9vR1MfIpRb1iRo5kAsDu1G1dhU+WfK2qkuvlBmnJXghRHUhxGohxCEhxAEhxPPOupfieELngfD/E4wdAD1gAFNPhP9Mh738zH1ToZU1eGYbtBml9dVfSYbFTPThLWthy2Mhbld3Vxp3bOCIaAvUItSXee/dzZoXOzK4ViW63dmEnw7/j6o1K7N4bwwpZovTY1AUp42iEUJUBapKKXcKIbyAHUB/KWW+g4HVKJobk5R2QJR9y1NKLcGfWAMrP4Ten0FI0QuMvT9wPJFLdme/7NQb9PhWrsTUgxNw88yj+8mJLl1MYv5Xf7Nu0zFWRbTApBcMaFWdh9rUoEFwKYadKre8clnRSUp5Djh3+fsUIcQhoBqgZntUME5rsRd+Y+1PSwYkndH65ps9BN3fyzGsMj9v//4if36xmL++X0ZmRhbt7ongkbH3l3lyjz+XyKjmY0hL0oZrhh6JIblFfWbbJb9tPU3LGr58NqgJtQPVy1fFscpkHLwQIgxYBzSSUiZft28kMBIgNDS05alTp5wej1IBZaZoNW22fAcuJi3Jtx5R3lEVyVfP/Mjfk1diu67MgXuwH3fPGMOCPeeYObIt3iYDkVEJ+LgbqRNUfsleSsmKX9ex4JslZKSa6TT4dga+2EeN3b9BletEJyGEJ7AWGCelnFvQsaqLRilU3DFY/g7U6gJtRoLdrrX0b+AXl8NqPcP5qNhc200erny34zNCbgvO3nbP1xvYE51E65p+PNC6Or0bVcVk0JdluHw56ntWzVifvfCJ0WSgco1Avtv5Ga5uroWcrZS1cpvoJIQwAHOAGYUld0UpkoA68MDMq633Xb/A1F5aEbMblHdA3q1xm9WOp69Hjm2TH4ng9d71iE028+KsPUSMW8G0jSfLIkwAzp28wIqf1+ZY1SrLbOFidDyrfttQZnEojuHMUTQCmAIcklJ+4az7KLeoKy12V09t8e8p3eGPRyD+eLmGlZeBL/bNUSwNwMWop2mnBrnKJwR6uTKqU21WvdyZmSPa0q1eEH4eWh36hLQspm+K4lJ67hmyjnJoy1H0eTwxmNMy2bF8r9PuqziHM+vBtweGAfuEELsvb3tTSln6VRoU5YpGAyG8F2yaqH0dXgxd34EOL5R3ZNm63N+eUwfOMPuLRRhdDVizrNRpUZM3fst/5LBOJ7i9tj+31/bP3rbqcCzvLTzAuL8O0aNhZQa1DKFjnQBc9I5rp/lVybuAm4tBT1BogMPuo5QNVWxMuXmkXIB1n0F4T7itl/Zi1m4Dt7KrOlmQ5IQUTu49jX+wb45+9+I4EJPE7MhoFuw+S2K6hSreJla83AlPV8e01ex2Ow/XHs3FM3HY7Vdzg6u7Kz/s+S/BtYu+xq5SNlQ1SeXWtPJD2D5Za823HglGj8LPqSCyrHZWHY5l/9kkXulVF4BP/j5EkLeJfk2DCfQq+cvQ81GxvD9wPGcOnUXnosfkZuTV6aOJuLO5o8JXHEgleOXWdG4vrPoQji4Dj0Do8BK0ekybIXuTsdrsDJq0md1nLqHXCTqGB9C/WTV6NqyMex5Fz4rifFQs5rRMqtcLRq8v25E8StGpBK/c2k5vhdUfwcl12iIj935X3hE5zdELKczbdZb5u84Sk2RmTK+6PNOlDpbLZRsMDuyvV24MKsErCsDJ9doCI5UbaKNtji7X6tEbynZma1mw2yXboxIIC/CgsreJv/ae4+35++jduCr9mgYTEeaHXnfjzh1Qiq5cShUoyg2nZser3x+cDyvHwvrPod1oaPW4NuTyJqHTCdrUujoCp5qvGx3DA5m38yy/bT1NkJcrdzepymt31ivziVRK2VEteKXYpD0NmTEXMteCvirC/UGEoX55h1V8URu1pQNPrAY3X2j/PHR4sbyjcqr0LCsrD8WyeG8Mp+LT+ef5jgghWLQnhmq+bjQL8UGnWvYVimrBKw4j7SnI+AFguwCYAR0yYwGy0qfo3O4q7/CKJ6y99hUdqbXkk6Kv7kuNLVJBs4rG3ehC36bB9G0ajM0uEUJgt0s+WHSQuNRMqlYy0btRVe5qXIUWob4q2VdwqgWvFIs99WtI/R7IzLlDeCGCtqBVp6ig7DbQ6eH0Fvjpbmh8H7R/DoIq4NNJMSVlWFh56AJ/7zvHun/jyLLZebpzbV69sx52u8QmpXpBe4NSLXjFcczLyZXcAbCD9QgYGpV1RI6ju9wXXSkEIp6AnT/Dnt+gTg+4/Rmo1fmGLmpWGpXcDAxoEcKAFiGkmC2sOhxLvSpanfrtUQmM/GUH3etXplfDynQMD8TNqPrtKwLVgleKxR7/EFi25bHHhAhYhHCp4dT7S8teZNoMsF8E1y4It4EInZPK2KYnaBOltv2gte5fOlhhR9xI6wlk2jRtwXRjc4T7Iwh95SKdezAmmcnrT7Di0AWSzVZMBh13hAcy7t7GpZpQpTiGGiapOIw0L0NeGgNkXLNVBy63oQtY6NR729PnQPIHQBZgB0zaS17/OQjd1REwUmpPGEI4KPlYzBB3BKo21coTT+sNtTppk6a8Sj91X0o72BNA54kQjp+EJbO2IxOeQPu92QADCJP2e3MJy3mseSUy9SuwnQWXugivVxBGbQarxWZn28kElh04z9aTCSx6tgMGvY4ZW0+RarbSrX5lagd6qDVni+F8VCzbl+zG5O5Ku3ta4VGp+LOtVYJXHEZKiUz5AtKngTACdtAFIvymI/Qlq69StPtmIGPbgsy4bo8JPJ9H5/k40haDTHoDsi4/YRgjEN4fI1xCHBdIWjwseBr+XQI6g7aGbOuREBJRou4be8bfkPIh2FO1DW73IrzfRgijw0K2X7wTbCeu26oD127ofL+5elz6PEh+D+3l+RUmhN80hDH/pRKf+nUH/+w/D0CYvzvd6lemV8MqtK7p57CfoSSyzFmsmrmRbX/vJCDEjz5P9iS0XrVyjelaP3/wB7P+Mx+hE+h0Oux2yXtzXiGiV7NiXUcleMXhpC0OLHtAFwCGJk5vtcmsSGTiSJCpuXcamiH8fkVe7Ar2OLTWPYAOdP6IwFWOa81fEX8ctv0Iu2dAZjIMmw+1uxTrEjJzi/YzXZdQMd2FzudTh4Qp7anI2NaANfdO4Ymu8k7tOCmRF9uBPT73cYaW6PxnFnifmEsZrDwcy4qDF9h8Ip5eDasw8QGt5b9oTwxtavoR5F12JSIyUjN4rt1bnD+plVvQu+hwMbjw+q/P0eHeNmUWR34ObvmXV7uPJTM95/ssk4crf5yfjJtH0X9X6iWr4nBCHwD6bmV4Q2+Qtnz2+YB5Gcg0riZ3tO9lmrbPra9j4/GvDb0/ha5vwYF5EHZ5EtXG/0FyDLR8FILqFXgJmfYtOZM72t/Ni5H2NxE6ByzGLYzku+yDuGZil0wBe3Lex1mPFHqbYB83hrWtwbC2NUjPspKcoX2gnIpP49mZuwBoGOxN57qBdK4bRPPqPg4tc3y9Bd8sJebYBbLMWu18m9WOzZrFfx//lrZ9WuJiKN/Ut/znNWRl5K7rr9PpiFy6h44DHPMhpMY9KRWDSzjog8n9T9YN4TEMbGfy6L5B22Zz4jq/rl7Q4mHQX04YyTGwfQp82wam3gl7ftcWDc+L9Uze24VBe4nsAEIYwXQXcH2XjwncH77mQHftvnnRFe89g7vRhSqVtBZoqJ87/zzfkVfvrIuH0YVJa08weNJm/r7cpZOYlkXMpXx+P6Wwbvam7OR+LbvNzvE9Zb/uc9zZeM4cOYvdrjVALJlW8uo9kYA1K4+nrRJSLXilQhBCgO8PyMTHLic/HUgLeD6JcO2IlFYQbpdb8dee6A4uBbekHar3f+COMVrXzY6fYN6TcGIN3Dsp97HGpmA+R86nDgAJese9NxDe7yHtcZC1XWvRy0xw64PwePTqMcIF6T4c0qaR8wW6G8Lr2SLdJ/5cItFHYgiuU4XAEP/L1xXUr+pN/arePN25DkkZFjYei6N9bW3xkDk7o/nor0PUDvSgY3ggd9wWQJua/niUsr69m1feo53sNplrdS1nijsbz9jBn3N8dxQ6vQ43TxNjpj1D5yHtWfvHphxLIwLYLFZa9mzisPurPnilQpFSgmUvyEQwNEXofC9vtyHj+4E1CrBcPtoA+lBt+KYoXcKw2WxELtnNsd1RVKtThXb9W2N0LWRSl90OpzaCqRJUbQIXj8CcJ6D5MGg8CGlMRMYPvPzkceX/oRt4jkbnOaJU8eZFWk9rTzoudfIcIimlHZk6AdKna91hwg08X0bncX+B17VZbXwxYhKrf9+I0WQgK9NC2z4tef2X5wr9HZ2OT2fZwfOsPxrH1pPxmC12TAYdu9/ticmgJzbZjJ+HsdjdOev+3Mz4R7/JkUCFEITUDWbKgS/LZKSPlJLHGrxAzLHz2G1XP8Rd3V35bsd/mDFuLhvnbSUzPROdXo/eoGf0V4/R+/HidX2ql6zKLUHak7URPua/AAmmuxFeLyF0lQo9tyBpSWm80PEdLpy6SGZaJq7urrh7u/HVpnEEhQYW/UKnt8I/Y+DcHtAb4bY7kY07I4O2gG0P6AMRHqMQjn5fUExSWrQ+eVEJIQqf0PTzB3/wx2cLyLymT9noZuTuEd15esKjBZyZk9liY8epRI5fTOXh28MAuP+HzRw4m0ybWn60qx1Auzr+3BbkVWgJBSkl378ynYXfLsNgdEECnj7ujF/5HtXqVC1yTKVxcPMRXu/1ERmpOd+z6F309Hu6F099OZx96w+xacF23Dxd6Tb0jhKt9KUSvKKUwsTRk/l78socfaM6vY6mnRvy2fJ3i3/Bc3thz0zY+wdkpcIr/2qt/PQErehZBRlHbrPa0Ol1DK78OElxKbn2u7q7sijll1K1lpceOM/afy+y8Vgcp+LTAbizYRUmDdOGbZ5JSCfE1y3fe8SeiePAxiP4Vq5Ek04N0OnK7rXj+jlb+O/j35KenPsdw+39WjF2/msOuY8aRaMopbD69425XnzZbXb2rj1AljkLo6mYY9arNtG+eoyFC/u15A7w6wDITIXGg6HxIG2kzg3o4OYj/O/pHzm593R2l0xesjIysdvtpVoNqlfDKvRqqL3kjU5MZ/PxeHzdtd93UrqFO8avJtDTlTa1/Glby482Nf1zTLYKqh5A0P3ls1h43dZ18nxh6upupEV3x/WzF0SNolGUwhTwkFuqJ2C9AYKbX7mQNrTSqwqs+QQmtoAfusBB584OLq4zR87yWo8PObHnFFLKHN0y16vVNMyhS/2F+LozuFV1ujfQ3h/o9YJx/Rtze21/tp6I5615++n+xVp+23Ya0AqoHYxJzrF4eFkKqh5Az+FdcrzUNbi64BvkQ6/hncskBtWCV5RCdLrvdpZMW52zi0YnaNyxPq5uDhqRIYS2ulTLR7SyxfvnaF8Zidr+1Fg4MB8a9HNIeYSSmv3fRbla7PJyAtW56LBb7ej0OowmA8996/gXxdfydHXhwTahPNgmFCklUfHpbD0RT/s6Wot99eFYXpi1G2+TC63C/GgV5kvrMD+ahPhgdCmbtu1z3zxB/TbhzJ/4D+nJ6XQY2Ib7xtyDm2fZ1DRSffCKUojUS2k81+4t4qLjyUg14+ZpwuRh4qvN46gS5uSa8VJqyX/PLJg3EhBQox00uAfq9wXv/F/KScsRZOo3YD2ojZzxfBphKF3XwLNt3+DwtmO5trt5udG8ayPiYxKp1bQGg1/uS/W65VsW4GJKJhuOXWTbyQS2RyVyLFabBb3p9a4E+7gRGZVAQloWLWv44u9ZcYumqZesilJKNquNrX/t5PieKIJrV6HDgNaOa70X1cUjWiv+wDy4eAiEDl7+FzwDtclU11S6lFl7kAkPo5V2tgMCcEX4fodwbV/iECaOnsxfP67AZsk5q9hoMvDLiW/wq+Jb4ms7W0JaFnvOXKJLPe1D+dmZu1i0JwaAmgEeNA/1oXWYH/e3Di3PMItNJXhFudnEHdUWJmkxTPv7jPsg6QzUvQvq3YXd9TOw7sp9nr4musClJb7tuZMXeLLZK2SkXB365+pmpNN97Rgz7ZkSX7c8mC029p9NIvJUIjtOJbLzVCKVvU38/bxWduLjvw/hohM0D/WlWXWfG7Y0skrwinKz2z5Fa9mf2gTShnTXIxu7QfPry88KROX9pVp56+S+U3z74k8c2HgEd283+j/bmwdevxe9S8VeBERKSVKGBZ/Lo3Qe+GEL26MSsF5+x1DNx40H24TyTJc6gPYBcSMsWK4SvKLcKtIT4Ogy5J4xyCA7NPMAi0SsTEJWN0INH0TtfapmexFlZNnYH5PE7tOX2B19iRahvjzeoSapmVaaj11G7UBPGlerRJPqPjQNqUTdKl64lvEHnUrwinKLsadNg5QJQAYkWBFLLiGSL/ebB9aH8O7Q6nHwq1mucVZUl9KzmLrhJHuik9gbfYnEdG1k0ft9GzC8fU1iU8wsP3iBxpXjuC0gCVf3Zgidc+rjq4lOinKLEe6PIG0XIf0X8PdEPuCONHdCxNyGOLYStkyCRoO0g89sg5hdULubNrlKte4L5eNu5KWedQGtayc6MYN9Z5NoXE2btBZ54gRvzTsJgEFnJdz/TxoGezG6Zz9q+Bd/1aaScloLXggxFegDxEopi7QSs2rBK4pjSXuqtvyevmrO+vKZqWBwB50OVnwAG77Qtleqri0uXrsLNOh/dSFypVhsF/tzJiGW/ReqceBiNQ7EhrA/NpR5I3ypEdyDXzZHMW1jFPWDvWkUXIlRnWqVuNusXLpohBB3AKnAzyrBK8oNLuEEHF8NJ1bDyXVg8NAWGRcC9s4GV09t/L2pdIXbbgXSGoWM68f1i7lICRjbovf/meUHLzA78gwHzyXjohOsGVO81cCuVS5dNFLKdUKIMGddX1EUB/KrpX1FPA42KyRHX+2qWfOx9gEgdFC1GYR1gLq9tYRfTqTMQKbPgcyVoAtAuD+EMDYtt3hykCkg9LlKXAgByCQAejSoTI/LJRfMlnxWKnOAcu+DF0KMBEYChIZWrAkGinJT0ruAb9jVvz+1GaK3ay37qA2w5TuwpGsJ3m6DFe9BSGsIvV2bdOVk0p6OjB8Etmi0VrJAmpcivd9C5z7E6fcvlEtd+mHChgAAC8FJREFUtIll13MF1+65tjpzqKVTR9FcbsEvVl00inITyUrXErxHgLb4+HftwHq5O8I/HELbQsQTENzMKbfXRgh9Se71bN0QQf9v7+6DrKrrOI6/v/fevfssuyyL8igIiJiTiRg+FgmTOvlUjo06mUPNWE2UNpCROtmUMz1oqeNok2MlFWllluiYUWqho8MISiriA2niKggLirCs+3Dvtz9+l1jYZV3ce+7Zc/fzmjmzd8+59/y+v1n47m9/55zv7wksVRNJuwci334fbL8K6CQ8SVwF6dFY0z3FWWu3B91FIyLFk60JG4S7bhZvgDfXwIbH4bUnYN194QIthKdtH7sRJhwXRvnjZkJ2kHeRvLec3smdMC3S9SxUFmfB6sFIVZ+FZ6bgu34NuU1QOQerPh9Lle4OGlCCF5HBylTCxNlhO/kbYanC3RPQ7W/D1vXw0l/D95aC0R+Ci+6CEeMLd/NUH9jdOqmG/RzIQ6p+MD0pKqs4Ehvxw1hjiCzBm9mdwBxglJm1ANe4+y+iak9EhoieqyZNPyNsu7aFefyWVbBxDdQV1oR9+Fp4+jfh4u24Y2DszFAjv3FSn/fjh9s+t/TVKKRGQ2ZGJF1KqijvorkwqnOLSMLUjITDTwtbT1PnQb4L3ngKVv4ccp1QdwgsejEcf+aPYXQ/5mhonIxv/yZ0v9D7/KmRWOPtKsGwD03RiEh8ps0LG0B3J2xeC22te46vuA5aQ7L3bB3W1IVPqoCjC3PZeYeUQWoClpmI5zZB90uQHo9lDitxZ4YeJXgRGRoy2T1LGO725cdC7fs318Cbj0LLMmxHLszwu2NLtkBtGppWkj/4NDhoPTTXQbXjFUeH+vepujh6MySo2JiIJIJ7F775+PAgEUC3Y6t3wtYcbAPb2RHed2wtflwddBi2uh6b8OVQYG30jDBVVGZ0m6SIJJ5ZBV5/Jbz7PaAdMobPbgSrAarw9k2wrRtqChd5d3bCC+vhmYV7TlLbDGffHC78trXC5udh1OHhom8Zzt8rwYtIYqRqzsPTY/G22yC3EbLHY3VfwlvPhKoUjM3ueXNTBT6/GW/LYPnrsNb/humegwprxb66Au6eH15XjoBR08L28W+FMsodOyGVgYqqkvezWJTgRSRRrPIEyM6E/NuQagoj++xJ0PFgH282qK+GERVw+IK9j035BFz8l7D8YeuL0PoSvPIvmLM4HF99Byy/GhomQNNUGDklPNg18/ODf1irRJTgRSQx3HP4jp/Art+GHZbB6xZg9YvwjoeArv18so8lCqsbQ5Kfsp9KjhNPCKP5rS+Hkgwtq6BjBxxbGPUvvxqeX1Yo1DYZGieHr0ecOWSme5TgRSQxfOfNsGsp/y9V4MCOm3BrhIYb4J3LgH2rM+ah8gNUvhx/bNj+37jjba2Qbod8Hjvkw/DuxlBpc+2fw1O7taNhxlnh/fcugLfWQuOh0HAoNEwM8/2TTznwWD4gJXgRSQT3HOxaArTvc6Qd2m4h1fwP8jXzwypWAKQBxxpuwWzw8+je9TTe/m3Y2RK+n3AKdtT1WKqxEMY70NbjKduRk2F7C2z8N6y7PzzQNfYYuPSf4fjS88OF3uYj4NM/G3R8fVGCF5Fk8PfC1pdC+YLUQVfgNZ+FjkfBaqFqXlGqN3ruDfztL4Dv2rOz41F82xeg6Z7wBG11Q9h2O2Vh2CCUVd6xCTrb9hxvnh6e3G3rq/RCcSjBi0gyWA2kRkH+rd7HKvbUoLHMJMhMKmrTvmspeOc+e7sg9wp0r4WK96mInkrDiHF77/vktUWNsc9mI29BRKQIzAzqrwT2nW6pwuqviLbx7v8A3X0cSIc1b4coJXgRSYxU9RlY461Q8RFINUH2RGzkb7DszGgbrphF718sgHcN6QqWmqIRkUSxypOxypNL2mZ36lw6224lm4VM4Y7LvFeSqp6LZYbuUqMawYuIvI9rzvs5X5k3hYf/3Mj2rWneer2CpTeMoS3/3bhD65dG8CIi/Xj1uQ08++g6OnYZP7l8z2i9sjpLzcErOH/h2TFG1z+N4EVE+vHqM6+RSvdOlR3tnbyw8uUYIho4JXgRkX6MnXoInu9dVj1bVcGko4bu/DsowYuI9Gv6cVOZOGMcmezeM9qZbIZPXTovpqgGRgleRMqG53fi+XeLek4z40fLv8PJn5lNpiJNKp1i+kencsOK7zPykMaitlVsusgqIonn3S349iugaw1geOYIrOHHWGZKUc5f11DLVb+7nFx3jlx3jmxV9v0/NARoBC8iiebeiW+7ALqeIjxt2gXdz+FbL8TzO4vaVjqTTkxyByV4EUm6jkfA24B8j50O3gHvPRBXVEOCEryIJFuuJSTzXtrx3OslD2coUYIXkWTLfAisj2kTq8Her8pjmVOCF5Fky86GzDSgsudOSI2BylPjimpIUIIXkUQzM6xxCdReAqnRocpkzYVY0+8x62Mt1mFEt0mKSOJZqgarXwT1i+IOZUiJdARvZqeb2Ytmtt7MFkfZloiI7C2yBG9maeAW4AzgSOBCMzsyqvZERGRvUY7gPwqsd/dX3L0TuAs4J8L2RESkhygT/Dig502oLYV9ezGzS81slZmt2rIlutXFRUSGmygTvPWxr1fNTXe/zd1nufus5ubmCMMRERleokzwLcCEHt+PB96MsD0REenB3HsXsi/Kic0ywEvAXOAN4EngIndf289ntgCvHUAzo4DWwcSZQMOxz6B+Dzfq98Ad6u59Tn9Edh+8u3eb2QLgb0Aa+GV/yb3wmQOaozGzVe4+axBhJs5w7DOo33HHUWrqd3FE+qCTuz8ADO9ybiIiMVGpAhGRMpX0BH9b3AHEYDj2GdTv4Ub9LoLILrKKiEi8kj6CFxGR/VCCFxEpU4lM8MOxSqWZTTCzR8xsnZmtNbPL4o6plMwsbWZPm9n9ccdSKmbWYGZ3m9kLhZ/7CXHHFDUz+0bh3/dzZnanmVXFHVNUzOyXZrbZzJ7rsW+kmf3dzF4ufG0cTBuJS/DDuEplN7DQ3WcAxwNfHSb93u0yYF3cQZTYTcCD7n4EcDRl3n8zGwd8HZjl7kcRnp+5IN6oInUHcPo++xYDD7n7NOChwvcfWOISPMO0SqW7b3T3pwqvdxD+s/cq3laOzGw88Cng9rhjKRUzOwj4GPALAHfvdPd34o2qJDJAdeFJ+BrKuLyJu68Atu2z+xxgSeH1EuDcwbSRxAQ/oCqV5czMJgHHACvjjaRkbgSuAPJxB1JChwFbgF8VpqZuN7PauIOKkru/AVwPbAA2AtvdfXm8UZXcwe6+EcKgDhg9mJMlMcEPqEpluTKzOuBPwOXu/m7c8UTNzM4ENrv76rhjKbEMMBP4mbsfA7QxyD/Xh7rCfPM5wGRgLFBrZp+LN6pkS2KCH7ZVKi2sIPwnYKm73xN3PCVyEnC2mf2XMB13qpn9Nt6QSqIFaHH33X+l3U1I+OVsHvCqu29x9y7gHuDEmGMqtbfMbAxA4evmwZwsiQn+SWCamU02syzhIsyymGOKnJkZYT52nbv/NO54SsXdv+3u4919EuFn/bC7l/2ozt03Aa+b2fTCrrnA8zGGVAobgOPNrKbw730uZX5huQ/LgEsKry8B7h3MySItNhaFD1KlskycBFwMPGtmawr7riwUdJPy9DVgaWEg8wowP+Z4IuXuK83sbuApwl1jT1PGJQvM7E5gDjDKzFqAa4AfAn8wsy8SfuGdP6g2VKpARKQ8JXGKRkREBkAJXkSkTCnBi4iUKSV4EZEypQQvIlKmlOBFRMqUEryISJlSghfZDzM7zsyeMbMqM6st1Ck/Ku64RAZKDzqJ9MPMrgWqgGpCbZgfxBySyIApwYv0o1Am4EngPeBEd8/FHJLIgGmKRqR/I4E6oJ4wkhdJDI3gRfphZssIZYonA2PcfUHMIYkMWOKqSYqUipl9Huh2998V1gJ+3MxOdfeH445NZCA0ghcRKVOagxcRKVNK8CIiZUoJXkSkTCnBi4iUKSV4EZEypQQvIlKmlOBFRMrU/wCnViBxxM26lwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "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", "\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": [ "# Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Load the data in thefile \"sample_data.tsv\" as a numpy array\n", "2. Log-transform the data\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": 62, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1.19127275e-02 7.65848100e-04 2.37491470e-05]\n" ] } ], "source": [ "# use of ttest_ind function\n", "import scipy.stats as sps\n", "\n", "# two arrays of random numbers\n", "a = np.random.randn(3,5) * 3 + 15\n", "b = np.random.randn(3,8) * 2 + 5\n", "\n", "print(sps.ttest_ind(a, b, axis=1, equal_var=False).pvalue)" ] }, { "cell_type": "code", "execution_count": 63, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sample1.1\tsample1.2\tsample1.3\tsample1.4\tsample2.1\tsample2.2\tsample2.3\tsample2.4\r\n", "6.801411216875708305e+02\t5.565263845884210241e+02\t2.525828159262368047e+02\t2.956703139364956314e+02\t9.604893343325777550e+02\t4.193181331491081210e+02\t3.126077132320110081e+02\t8.246976543530363415e+02\r\n", "5.874671746889415545e+02\t4.128408625270991479e+02\t4.378621470916139060e+02\t6.736096888155080933e+02\t7.240856598618292992e+02\t4.807175238900777003e+02\t4.731265298940686534e+02\t5.216420065104225614e+02\r\n" ] } ], "source": [ "!head -3 test_data.tab" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "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 }