01_Aleatoreidad.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si estas en Google Colab, ejecutá esta celda de abajo. Si clonaste el repositorio con git o lo descargaste y corres localmente, ignorá esta primer celda."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!mkdir library\n",
    "!wget https://raw.githubusercontent.com/exord/UNSAM_IA/master/library/utils.py\n",
    "!mv utils.py library/utils.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Por si alguien corre en python2\n",
    "from __future__ import division\n",
    "\n",
    "# Preparamos todo para correr\n",
    "import numpy as np\n",
    "from math import sqrt\n",
    "from matplotlib import pylab as plt\n",
    "\n",
    "# Importamos un módulo nuestro\n",
    "from library import utils"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *Aleatoreidad y probabilidad*\n",
    "La aleatoredidad es un concepto central de la teoría de probabilidad. En general, la aleatoredidad viene asociada con la incetidumbre, que puede surgir por los errores de los datos o por los tamaños finitos de las muestras.\n",
    "\n",
    "En la computadora, la aleatoreidad se simula usando como base un [generador de números pseudo aleatorios](https://en.wikipedia.org/wiki/Pseudorandom_number_generator) (GNPA). Un GNPA genera una serie de números cuyas propiedades son similares a las de una serie de números aleatorios.\n",
    "\n",
    "En <tt>python</tt> muchas de las funciones relacionadas están implementadas en el módulo <tt>random</tt>, pero una implementación más práctica aparece en el paquete <tt>numpy</tt>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Importamos el paquete random de numpy\n",
    "from numpy import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generamos un número al azar en el intervalo [0.0, 1.0)\n",
    "x = random.rand()\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fíjense que si corren la celda de arriba muchas veces, el resultado obtenido es diferente, pero siempre un número entre 0 y 1. Pero en realidad, los GNPA utilizan algoritmos, cuyas propiedades estás perfectamente definidas por el valor inicial (la semilla, o <em>seed</em>). Por lo tanto, no son realmente aleatorios. En los viejos tiempos, en <tt>python</tt>, el <em>seed</em> podía definirse simplemente con la función homónima."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "random.seed(1234)\n",
    "print(random.rand())\n",
    "random.seed(1234)\n",
    "print(random.rand())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una característica práctica de la implementación en <tt>numpy</tt> es que se puede generar <tt>arrays</tt> de números aleatorios, pasándole el tamaño en cada dimensión. El jueves vamos a ver más detalle de las características de los <tt>arrays</tt> de <tt>numpy</tt>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "random.seed(1234)\n",
    "print('Un array de 5x1')\n",
    "print(random.rand(5,))\n",
    "print('Un array de 3x2')\n",
    "print(random.rand(3,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Y por supuesto, si vuelvo a la misma semilla, los resultado son idénticos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "random.seed(1234)\n",
    "print('Un array de 5x1')\n",
    "print(random.rand(5,))\n",
    "print('Un array de 3x2')\n",
    "print(random.rand(3,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como se puede ver arriba, para una semilla dada, la serie de números generada es idéntica. En la actualidad, fijar el <em>seed</em> de la manera que aparece arriba se considera mala práctica, y la sugerencia es hace lo siguiente:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numpy.random import RandomState, SeedSequence\n",
    "rs = RandomState(random.MT19937(SeedSequence(1234)))\n",
    "print(rs.rand(5,))\n",
    "# Later, you want to restart the stream\n",
    "rs = RandomState(random.MT19937(SeedSequence(1234)))\n",
    "print(rs.rand(5,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora vamos a intentar visualizar una muestra de una distribución uniforme. Generamos una muestra de 100 elementos de vectores en dos dimenciones, en la región -1 < x < 1; -1 < y < 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Genaramos la muestra\n",
    "sample = random.rand(100, 2)\n",
    "\n",
    "# Como queremos que la muestra esté entre -1 y 1 en ambas dimensiones, tenemos que multiplicar por 2 y restar 1.\n",
    "sample = sample * 2 - 1\n",
    "\n",
    "# Imprimimos la forma de la muestra y la representamos en un plano.\n",
    "print(sample.shape)\n",
    "l0, = plt.plot(sample[:, 0], sample[:, 1], '.')\n",
    "\n",
    "# Para que sea más bonito, ponemos los ejes del mismo tamaño y le agregamos lineas verticales y horizontales\n",
    "ax = plt.gca()\n",
    "ax.set_aspect('equal')\n",
    "ax.axvline(0, color='0.5', zorder=-1)\n",
    "ax.axhline(0, color='0.5', zorder=-1)\n",
    "ax.set_xlim(-1, 1)\n",
    "ax.set_ylim(-1, 1)\n",
    "ax.set_xlabel(u'$\\\\theta_0$', fontsize=16)\n",
    "ax.set_ylabel(u'$\\\\theta_1$', fontsize=16)\n",
    "#ax.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos a contar ahora cuantos de los elementos de la muestra caen en el cuadrante superior derecho, y veamos a qué fracción del total corresponde. Pero antes, <b>¿qué valor esperamos? ¿por qué?</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Escribimos la condición de pertenencia al cuadrante de interés\n",
    "cond = (sample[:, 0] > 0) * (sample[:, 1] > 0)\n",
    "\n",
    "# Veamos la pinta que tiene cond\n",
    "print(cond)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Sumar los elementos de cond nos dice en cuantos casos la condición se satisface\n",
    "n_ad = np.sum(cond)\n",
    "print('Hay {} elementos en el cuadrante superior derecho.'.format(n_ad))\n",
    "\n",
    "# Calculo la fracción a la que correponde\n",
    "f_ad = n_ad / len(sample)\n",
    "print('La fracción de elementos en el cuadrante superior derecho es {}.'.format(f_ad))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como es de esperar, la fracción se aproxima a 1/4. Pero, <b>¿cuán cerca de 0.25 tenemos que tiene que estar el número obtenido para que nos quedemos contentos? ¿y cómo cambia esto con el tamaño de la muestra?</b>\n",
    "\n",
    "## Incertezas\n",
    "\n",
    "Para responder a esas preguntas, por un lado tenemos que definir algún concepto de incerteza de f_ur. Ya veremos que la variable f_ur sigue una distribución Poissoneana, cuya desviación estándar (ya lo definiremos bien, pero podemos pensarlo como un error), es $\\sqrt{n_{UR}}/N$, donde $N$ es el tamaño de la muestra. Entonces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sd_f = sqrt(n_ad)/len(sample)\n",
    "\n",
    "print('El error de f_ur es {}'.format(sd_f))\n",
    "print('La fracción es entonces {:.3f} \\pm {:.3f}'.format(f_ad, sd_f))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vemos, entonces que la fracción es perfectamente compatible con el valor esperado 1/4. \n",
    "\n",
    "***\n",
    "### Preguntas\n",
    "\n",
    "<ul>\n",
    "    <li>Si tuviéramos solo esta medición (es decir, la fracción de elementos de la muestra que cae en el cuadrante superior derecho), pero si no conociéramos la forma en la que la muestra fue generada, ¿podríamos decir que la distribución es uniforme en la región de $\\mathbb{R}ˆ2$ tal que $-1 < x < 1$ y $-1 < y < 1$? </li>\n",
    "    <li> A la inversa, si el número fuera, por ejemplo $0.245 \\pm 0.001$, ¿qué podríamos decir del proceso de generación de la muestra?\n",
    "    </ul>    \n",
    "    \n",
    "***\n",
    "\n",
    "## Tamaño de la muestra\n",
    "\n",
    "Estudiemos ahora cómo cambia nuestra precisión con el tamaño de la muestra. Para eso, escribamos una pequeña función que genere una muestra, y que calcule la fracción de elementos que caen en el cuadrante superior derecho y su incerteza. La función tiene un solo argumento: el tamaño de la muestra.\n",
    "\n",
    "Para escribirla, hay que ir recolectando las líneas relevantes de los celdas de arriba (Consejo: usar <tt>random.rand</tt> para generar la muestra). Borrar el comando <tt>pass</tt> y a codear!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def fraccion_ad(n, seed=None):\n",
    "    \"\"\"\n",
    "    :param int n: Tamaño de la muestra\n",
    "    \n",
    "    :returns: fracción de la muestra en el cuadrante superior derecho e incerteza de esta fracción\n",
    "    :rtype: float, float\n",
    "    \"\"\"\n",
    "    # Fijar la semilla\n",
    "    if seed is not None:\n",
    "        np.random.seed(seed)\n",
    "        \n",
    "    # Crear la muestra.    \n",
    "    pass\n",
    "\n",
    "    return\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora vamos a testear la función usando el paquete utils que está en el repositorio. Después los invitamos a chusmearlo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "utils.test_fraccion(fraccion_ad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora calculemos la fracción para diferentes tamaños de la muestra, y guardemos todo en unas listas, para plotearlas después. <em>NB: el jueves veremos con más detalle cómo funcionan las listas, arreglos, etc.</em>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Crea dos listas vacías\n",
    "f = []\n",
    "sd = []\n",
    "\n",
    "# Crea un arreglo con los tamaños de las muestras que vamos a usar. Generamos valores desde 2^4 hasta 2^13.\n",
    "# Vamos a espacirlo en log, para tener más resolución en los tamaños pequeños.\n",
    "sizes = 2**np.arange(4, 14)\n",
    "\n",
    "# Iteramos en los tamaños (usamos siempre el mismo seed?)\n",
    "for size in sizes:\n",
    "    output = fraccion_ad(size)\n",
    "    f.append(output[0])\n",
    "    sd.append(output[1])\n",
    "    \n",
    "# Veamos que obtuvimos\n",
    "for i in range(len(f)):\n",
    "    print('{}, {}'.format(f[i], sd[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Y ahora el momento del plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.errorbar(sizes, f, sd, fmt='.')\n",
    "\n",
    "# Para que sea más legible, hacemos el eje x logarítmico\n",
    "ax = plt.gca()\n",
    "ax.set_xscale('log')\n",
    "ax.axhline(0.25, color='0.5', zorder=-1)\n",
    "\n",
    "ax.set_xlabel('Tamaño de la muestra', fontsize=16)\n",
    "ax.set_ylabel('Fracción en AD', fontsize=16)\n",
    "ax.tick_params(axis='both', which='major', labelsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "### Preguntas\n",
    "\n",
    "<ul>\n",
    "    <li> Sí repetimos el experimento, ¿cómo va a cambiar el plot de arriba?</li>\n",
    "    <li> ¿En qué fracción de los casos veremos puntos cuyas barras de error no tocan 0.25?</li>\n",
    "    <li> ¿Cuál es la <em>probabilidad</em> de que un elemento de la muestra caiga en el cuadrante superior derecho?</li>\n",
    "</ul>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eficiencia\n",
    "\n",
    "Claro que al correr la función muchas veces, cometimos el error de generar muchas muestras. Una forma más práctica sería generar una sola muestra grande, e ir obteniendo sub-muestras en las cuales hacer el cálculo (ver [Ejercicios](Ejercicios_01_Aleatoreidad.ipynb)). Veamos cuánto tiempo nos tarda la generación de este plot con la versión actual. Para eso, podemos usar muchas mediciones diferentes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def corrida():\n",
    "    f = []\n",
    "    sd = []\n",
    "\n",
    "    # Crea un arreglo con los tamaños de las muestras que vamos a usar. Generamos valores desde 2^4 hasta 2^13.\n",
    "    # Vamos a espacirlo en log, para tener más resolución en los tamaños pequeños.\n",
    "    sizes = 2**np.arange(4, 14)\n",
    "\n",
    "    # Iteramos en los tamaños (usamos siempre el mismo seed?)\n",
    "    for size in sizes:\n",
    "        output = fraccion_ad(size)\n",
    "        f.append(output[0])\n",
    "        sd.append(output[1])\n",
    "        \n",
    "    return\n",
    "\n",
    "nruns = 1000\n",
    "# Podemos usar el módulo timeit\n",
    "import timeit\n",
    "dt = timeit.timeit(corrida, number=nruns)\n",
    "print('Tiempo de ejecución medida con timeit para {} corridas: {:.3f} segundos'.format(nruns, dt))\n",
    "\n",
    "# O contar el tiempo a mano con el módulo time\n",
    "import time\n",
    "ti = time.process_time()\n",
    "for i in range(nruns):\n",
    "    corrida()\n",
    "print('Tiempo de ejecución medida con time para {} corridas: {:.3f} segundos'.format(nruns, (time.process_time() - ti)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}