05-Scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SciPy - Librería de algorítmos científicos para Python\n",
    "\n",
    "Versión original en inglés de J.R. Johansson (robert@riken.jp) http://dml.riken.jp/~rob/\n",
    "\n",
    "Traducido/Adaptado por [G.F. Rubilar](http://google.com/+GuillermoRubilar) y Stefan Vogt.\n",
    "\n",
    "La última versión de estos [notebooks de IPython](http://ipython.org/notebook.html) está disponible en [http://github.com/gfrubi/clases-python-cientifico](http://github.com/gfrubi/clases-python-cientifico).\n",
    "\n",
    "La última versión del original (en inglés) de estos [notebooks de IPython](http://ipython.org/notebook.html) está disponible en [http://github.com/jrjohansson/scientific-python-lectures](http://github.com/jrjohansson/scientific-python-lectures).\n",
    "\n",
    "Los otros notebooks de esta serie están listados en [http://jrjohansson.github.com](http://jrjohansson.github.com)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# ¿qué hace esta línea?. La respuesta en la clase 4\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introducción\n",
    "\n",
    "El paquete SciPy agrega características a los algorítmos de bajo nivel de NumPy para arreglos multidimensionales, y provee un gran número de algorítmos de alto nivel de uso científico. Algunos de los tópicos que cubre SciPy son:\n",
    "\n",
    "* Funciones especiales ([scipy.special](http://docs.scipy.org/doc/scipy/reference/special.html))\n",
    "* Integración ([scipy.integrate](http://docs.scipy.org/doc/scipy/reference/integrate.html))\n",
    "* Optimización ([scipy.optimize](http://docs.scipy.org/doc/scipy/reference/optimize.html))\n",
    "* Interpolación ([scipy.interpolate](http://docs.scipy.org/doc/scipy/reference/interpolate.html))\n",
    "* Transformada de Fourier ([scipy.fftpack](http://docs.scipy.org/doc/scipy/reference/fftpack.html))\n",
    "* Procesamiento de señales ([scipy.signal](http://docs.scipy.org/doc/scipy/reference/signal.html))\n",
    "* Álgebra lineal ([scipy.linalg](http://docs.scipy.org/doc/scipy/reference/linalg.html))\n",
    "* Problemas de Eigenvalores de matrices dispersas ([scipy.sparse](http://docs.scipy.org/doc/scipy/reference/sparse.html))\n",
    "* Estadística ([scipy.stats](http://docs.scipy.org/doc/scipy/reference/stats.html))\n",
    "* Procesamiento de imágenes multi-dimensional ([scipy.ndimage](http://docs.scipy.org/doc/scipy/reference/ndimage.html))\n",
    "* Entrada/Salida desde/hacia archivos ([scipy.io](http://docs.scipy.org/doc/scipy/reference/io.html))\n",
    "\n",
    "Cada uno de estos submódulos provee un muchas funciones y clases que pueden ser usadas para resolver problemas en sus respectivos tópicos.\n",
    "\n",
    "En esta clases veremos cómo usar algunos de estos subpaquetes.\n",
    "\n",
    "Para acceder al paquete SciPy en un programa Python, comenzamos importando todo desde el módulo `scipy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si sólo necesitamos usar una parte de las rutinas de SciPy podemos incluir selectivamente sólo aquellos módulos en los que estamos interesados. Por ejemplo, para incluir el paquete de álgebra lineal bajo el nombre `la`, podemos incluir:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy.linalg as la"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Funciones Especiales\n",
    "\n",
    "En muchos problemas de física y química computacional son importantes varias funciones matemáticas especiales. SciPy provee implementaciones de muchas de estas funciones especiales. Para más detalles, ver la lista de funciones en la documentación http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special. \n",
    "\n",
    "Para demostrar el uso típico de estas funciones especiales nos concentraremos en la función de Bessel de primera especie:\n",
    "\n",
    "$J_\\alpha(x) = \\sum_{m=0}^\\infty \\frac{(-1)^m}{m! \\, \\Gamma(m+\\alpha+1)} {\\left(\\frac{x}{2}\\right)}^{2m+\\alpha}$\n",
    "\n",
    "y de segunda especie:\n",
    "\n",
    "$J_\\alpha(x) = \\sum_{m=0}^\\infty \\frac{(-1)^m}{m! \\, \\Gamma(m+\\alpha+1)} {\\left(\\frac{x}{2}\\right)}^{2m+\\alpha}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.special import jn, yn, jn_zeros, yn_zeros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 0    # orden de la función\n",
    "x = 0.0\n",
    "\n",
    "# Función de Bessel de primera especie\n",
    "print(\"J_%d(%f) = %f\" % (n, x, jn(n, x)))\n",
    "\n",
    "x = 1.0\n",
    "# Función de Bessel de segunda especie\n",
    "print(\"Y_%d(%f) = %f\" % (n, x, yn(n, x)))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# El módulo scipy.special incluye muchas funciones de Bessel \n",
    "# Aquí usaremos las funciones jn e yn, que son las funciones de Bessel \n",
    "# de primera y segunda especie, y de orden real. Incluimos también las  \n",
    "# funciones jn_zeros e yn_zeros que entregan los ceros de las\n",
    "# funciones jn e yn.\n",
    "#\n",
    "\n",
    "x = np.linspace(0, 10, 100)\n",
    "\n",
    "for n in range(4):\n",
    "    plt.plot(x, jn(n, x), label=r\"$J_%d(x)$\" % n);\n",
    "plt.legend();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integración\n",
    "\n",
    "### Integración numérica: cuadraturas\n",
    "\n",
    "La evaluación numérica de una función, del tipo\n",
    "\n",
    "$\\displaystyle \\int_a^b f(x) dx$\n",
    "\n",
    "es llamada *cuadratura numérica*, o simplemente *cuadratura*. SciPy suministra funciones para diferentes tipos de cuadraturas, por ejemplo las funciones `quad`, `dblquad` y `tplquad` para calcular integrales simples, dobles o triples, respectivamente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import quad, dblquad, tplquad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Las función `quad` acepta una gran cantidad de argumentos opcionales, que pueden ser usados para ajustar detalles del comportamiento de la función (ingrese `help(quad)` para más detalles).\n",
    "\n",
    "El uso básico es el siguiente:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define una función simple para ser integrada\n",
    "def f(x):\n",
    "    return x**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_inf = 0 # el límite inferior de x\n",
    "x_sup = 2 # el límite superior de x\n",
    "\n",
    "val, errabs = quad(f, x_inf, x_sup)\n",
    "\n",
    "print(\"valor de la integral =\", val, \", error absoluto =\", errabs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si necesitamos incluir argumento extras en la función integrando podemos usar el argumento `args`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# ceros de la funcion de Bessel \n",
    "n = 3 # order\n",
    "m = 2 # numero de raizes a calcular\n",
    "jn_zeros(n, m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def integrando(x, n):\n",
    "    \"\"\"\n",
    "    función de Bessel de primera especie y orden n. \n",
    "    \"\"\"\n",
    "    return jn(n, x)\n",
    "\n",
    "\n",
    "x_inf = 0 # el límite inferior de x\n",
    "x_sup = 6.38 # el límite superior de x\n",
    "\n",
    "val, errabs = quad(integrando, x_inf, x_sup, args=(3,)) # evalua la integral con n=3\n",
    "\n",
    "print(val, errabs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para funciones simples podemos usar la función lambda function (función anónima) en lugar de definir explícitamente una función para el integrando:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "val, errabs = quad(lambda x: np.exp(-x ** 2), -np.Inf, np.Inf) # Inf = infinito!\n",
    "\n",
    "print(\"Resultado numérico  =\", val, errabs)\n",
    "\n",
    "analitico = np.sqrt(np.pi)\n",
    "print(\"Resultado analitico =\", analitico)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como se muestra en este ejemplo, podemos usar 'Inf' y '-Inf' como límites de la integral.\n",
    "\n",
    "Integrales de dimensión mayor se evalúan de forma similar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def integrando(x, y):\n",
    "    return np.exp(-x**2-y**2)\n",
    "\n",
    "x_inf = 0  \n",
    "x_sup = 10\n",
    "y_inf = 0\n",
    "y_sup = 10\n",
    "\n",
    "val, errabs = dblquad(integrando, x_inf, x_sup, lambda x : y_inf, lambda x: y_sup)\n",
    "\n",
    "print(val, errabs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note como requerimos incorporar funciones lambda para los límites de la integración en y, ya que estos límites pueden en general ser funciones de x."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Álgebra lineal\n",
    "\n",
    "El módulo de álgebra lineal contiene muchas funciones relacionadas con matrices, incluyendo resolución de ecuaciones lineales, cálculo de valores propios, funciones de matrices (por ejemplo, para exponenciación matricial), varias decomposiciones diferentes (SVD, LU, cholesky), etc. \n",
    "\n",
    "Una documentación detallada está disponible aquí: http://docs.scipy.org/doc/scipy/reference/linalg.html\n",
    "\n",
    "Veremos cómo usar algunas de estas funciones:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sistemas de ecuaciones lineales\n",
    "\n",
    "Los sistemas de ecuaciones lineales de la forma\n",
    "\n",
    "$A x = b$\n",
    "\n",
    "donde $A$ es una matriz y $x,b$ son vectores, pueden ser resueltos del modo siguiente:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import linalg\n",
    "\n",
    "A = np.array([[1,2,3], [4,12,6], [7,8,9]])\n",
    "b = np.array([4,2,10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# verificamos la solución\n",
    "np.dot(A, x) - b\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos también hacer lo mismo con\n",
    "\n",
    "$A X = B$,\n",
    "\n",
    "donde ahora $A, B$ y $X$ son matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = np.random.rand(3,3)\n",
    "B = np.random.rand(3,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = linalg.solve(A, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# verificamos la solución\n",
    "np.dot(A, X) - B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Valores y vectores propios"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El problema de valores propios para la matriz $A$:\n",
    "\n",
    "$\\displaystyle A v_n = \\lambda_n v_n$,\n",
    "\n",
    "donde $v_n$ es el $n$-ésimo vector propio y $\\lambda_n$ es el $n$-ésimo valor propio.\n",
    "\n",
    "Para calcular los vectores propios de una matriz usamos `eigvals` y para calcular tanto los valores como los vectores propios, podemos usar la función `eig`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "evals = linalg.eigvals(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "evals[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "evals, evecs = linalg.eig(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "evals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "evecs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los vectores propios correspondientes al $n$-ésimo valor propio (guardado en `evals[n]`) es la $n$-ésima *columna* en `evecs`, es decir, `evecs[:,n]`. Para verificar esto, intentemos multiplicar los vectores propios con la matriz y comparar el resultado con el producto del vector propio y el valor propio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 1\n",
    "\n",
    "np.dot(A, evecs[:,n]) - evals[n] * evecs[:,n]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Operaciones matriciales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# matriz inversa\n",
    "linalg.inv(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# determinante\n",
    "linalg.det(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# norma de un matri\n",
    "linalg.norm(A, ord=2), linalg.norm(A, ord=np.Inf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimización\n",
    "\n",
    "La optimización (encontrar el máximo o el mínimo de una funciónn) constituye un campo amplio en matemáticas, y la optimización de funciones complicadas o de muchas variables puede ser complicada. Aquí sólo revisaremos algunos casos muy simples. Para una introducción detallada a la optimización con SciPy, ver (en inglés): http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html\n",
    "\n",
    "Para usar el módulo de optimización de Scipy hay que importar el módulo `optimize`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Encontrando mímimos\n",
    "\n",
    "Veamos primero cómo encontrar el mínimo de una función simple de una variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 4*x**3 + (x-2)**2 + x**4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-5, 3, 100)\n",
    "plt.plot(x, f(x));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos usar la función `fmin_bfgs` para encontrar el mínimo de la función:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_min = optimize.fmin_bfgs(f, -2) # busca un mínimo local cerca -2\n",
    "x_min "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "optimize.fmin_bfgs(f, 0.5)  # busca un mínimo local cerca 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos también usar las funciones `brent` o `fminbound`. Estas funciones tienen una sintaxis algo distinta y usan algoritmos diferentes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "optimize.brent(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "optimize.fminbound(f, -4, 2)  # busca el mínimo en el intervalo (-4,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Encontrando las raíces de una función\n",
    "\n",
    "Para encontrar las soluciones a una ecuación de la forma $f(x) = 0$ podemos usar la función `fsolve`. Ella requiere especificar un punto inicial: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "omega_c = 3.0\n",
    "def f(omega):\n",
    "    return np.tan(2*np.pi*omega) - omega_c/omega"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig, ax  = plt.subplots(figsize=(12,4))\n",
    "x = np.linspace(0, 3.0, 1000)\n",
    "y = f(x)\n",
    "mask = np.where(np.abs(y) > 50)\n",
    "x[mask] = y[mask] = np.NaN # elimina líneas verticales cuando la función cambia de signo\n",
    "ax.plot(x, y)\n",
    "ax.plot([0, 3], [0, 0], 'k')\n",
    "ax.set_ylim(-5,5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "optimize.fsolve(f, 0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "optimize.fsolve(f, 0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "optimize.fsolve(f, 1.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolación\n",
    "\n",
    "La interpolación es simple y conveniente en Scipy: La función `interp1d`, cuando se le suministran arreglos describiendo datos X e Y, retorna un objeto que se comporta como una función que puede ser llamada para un valor de x arbitrary (en el rango cubierto por X), y que retorna el correspondiente valor interpolado de y:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.interpolate import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return np.sin(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = np.arange(0, 10)  \n",
    "x = np.linspace(0, 9, 100)\n",
    "\n",
    "y_meas = f(n) + 0.1 * np.random.randn(len(n)) # simula medidas con error\n",
    "y_real = f(x)\n",
    "\n",
    "linear_interpolation = interp1d(n, y_meas)\n",
    "y_interp1 = linear_interpolation(x)\n",
    "\n",
    "cubic_interpolation = interp1d(n, y_meas, kind='cubic')\n",
    "y_interp2 = cubic_interpolation(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,4))\n",
    "ax.plot(n, y_meas, 'bs', label='datos con ruido')\n",
    "ax.plot(x, y_real, 'k', lw=2, label=u'función exacta')\n",
    "ax.plot(x, y_interp1, 'r', label=u'interpolación lineal')\n",
    "ax.plot(x, y_interp2, 'g', label=u'interpolación cúbica')\n",
    "ax.legend(loc=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estadística\n",
    "\n",
    "El módulo `scipy.stats` contiene varias distribuciones estadísticas, funciones estadísticas. Para una documentación completa de estas las características, ver (en inglés) [http://docs.scipy.org/doc/scipy/reference/stats.html](http://docs.scipy.org/doc/scipy/reference/stats.html).\n",
    "\n",
    "También existe un paquete Python muy poderoso para modelamiento estadístoco llamado statsmodels. Ver [http://statsmodels.sourceforge.net](http://statsmodels.sourceforge.net) para más detalles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "crea una variable aleatoria (discreta) con distribución poissoniana\n",
    "\n",
    "$P(k) = \\frac{\\mu^k e^{\\mu}}{k!}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "X = stats.poisson(3.5) # mu = 2.5  (goles promedio en un mundial)\n",
    "X.pmf(1)  # Probabilidad de 0 goles en un partido del mundial.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = np.arange(0,15)\n",
    "\n",
    "fig, axes = plt.subplots(3,1, sharex=True)\n",
    "\n",
    "# grafica la \"probability mass function\" (PMF)\n",
    "axes[0].plot(n, X.pmf(n),lw=1,linestyle=\"-\")\n",
    "# grafica la \"commulative distribution function\" (CDF)\n",
    "axes[1].step(n, X.cdf(n))\n",
    "\n",
    "# grafica histograma de 1000 realizaciones de la variable estocástica X rvs: random variate sample\n",
    "axes[2].hist(X.rvs(size=1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# crea una variable aleatoria (contínua) con distribución normal\n",
    "Y = stats.norm(loc=1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-5,5,100)\n",
    "\n",
    "fig, axes = plt.subplots(3,1, sharex=True)\n",
    "\n",
    "# grafica la función distribución de probabilidad (\"probability distribution function\", PDF)\n",
    "axes[0].plot(x, Y.pdf(x))\n",
    "\n",
    "# grafica función de distribución acumulada (\"commulative distributin function\", CDF)\n",
    "axes[1].plot(x, Y.cdf(x));\n",
    "\n",
    "# grafica histograma de 1000 realizaciones aleatorias de la variable estocástica Y\n",
    "axes[2].hist(Y.rvs(size=1000), bins=50);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estadística:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X.mean(), X.std(), X.var() # distribución de Poission"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Y.mean(), Y.std(), Y.var() # distribucuón normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Datos de diametro en km de cometas apocalipticos que han impactado la tierra en su historia\n",
    "\n",
    "data = [7.40, 6.90, 5.20, 6.20, 4.20, 8.10, 9.50, 7.40, 6.80, 7.50, 7.80, 6.60, 6.00, 8.60, 8.0, 7.20, 7.30, 9.00, 9.11, 9.20, 8.20, 10.01, 11.3]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mean, sigma = stats.norm.fit(data)\n",
    "print(mean)\n",
    "print(sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(3,12,20)\n",
    "\n",
    "D = stats.norm(loc=mean, scale = sigma)\n",
    "plt.hist(data,x,color=\"g\", normed=True)\n",
    "plt.plot(x,D.pdf(x),color='b');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "D.mean(), D.std(), D.var()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prob_over_10 = 1 - D.cdf(x= 10) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(prob_over_10  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Regresiónes lineales\n",
    "\n",
    "Regresiones lineales son un métdod científico simple, que permite  predecir\n",
    "una relación continua y lineal entre dos variables. En el ejemplo a continuación se quiere \n",
    "obtener la masa molecular de una gas a partir de datos reales de presión y densidad del gas. \n",
    "La masa molecular del gas se puede calcular simplemente a partir de la ecuación de los gases ideales como:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mu = \\frac{\\rho RT}{P}\n",
    "\\end{equation}\n",
    "\n",
    "donde $\\rho$ es la densidad del gas. En estricto rigor esta ecuación solamente es valida a P=0. Por \n",
    "eso que hay que utilizar una regresión lineal para extrapolar los datos a P=0.\n",
    "\n",
    "Considermos los siguientes presiones de gas en kPa: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "P = [12.223,25.20,36.97,60.37,85.23,101.3]\n",
    "P = [float(x)*1000 for x in P] # en Pa⋅⋅\n",
    "rho = [0.225, 0.456, 0.664, 1.062, 1.468, 1.734]\n",
    "rho = [float(x)*1000 for x in rho ] # en g/m^3\n",
    "P_d = []\n",
    "\n",
    "for i in range(len(P)):\n",
    "    P_d.append(float(rho[i])/float(P[i]))\n",
    "\n",
    "plt.plot(P,P_d, 'yo')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Usando numpy\n",
    "fit = np.polyfit(P,P_d,1)\n",
    "print(fit)\n",
    "# Usando scipy\n",
    "fit_scipy = stats.linregress(P,P_d)\n",
    "print(fit_scipy)\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(P,P_d)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fit_fn = np.poly1d(fit)\n",
    "plt.plot(P,P_d, 'yo', P, fit_fn(P), '--k')\n",
    "#plt.legend(\"y = \"+str(slope)+\"x +\"+str(intercept))\n",
    "R = str(round(r_value,4))\n",
    "m = str(round(slope,4))\n",
    "t = str(round(intercept,4))\n",
    "plt.legend(('Pressure', \"y = \"+m+\"x +\"+t+\"   R = \"+R), 'upper left')\n",
    "plt.autoscale(True)\n",
    "plt.grid(True)\n",
    "plt.show()\n",
    "\n",
    "R = 8.314\n",
    "T = 298.15\n",
    "\n",
    "Mu = R*T*intercept\n",
    "\n",
    "print(\"La masa molar es: \"+str(Mu)+ \" g/mol\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ecuaciones diferenencias ordinarias (EDOs)\n",
    "\n",
    "SciPy provee dos formas diferentes para resolver EDOs: Una API (Interfaz de programación de aplicaciones, del inglés \"Application programming interface\") basada en la función `odeint`, y una API orientada al objeto basada en la clases `ode`. Usualmentey `odeint` es más simplea de usar, pero la clase `ode` ofrece niveles de control más finos.\n",
    "\n",
    "Aquí usaremos las funciones `odeint`. Para mayor información sobre las clases `ode`, use `help(ode)`. Hace casi todo lo que hace `odeint`, pero de una forma más orientada al objeto.\n",
    "\n",
    "Para usar `odeint`, primero importelo desde el módulo `scipy.integrate`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint, ode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un sistema de EDOs es usualmente formulado en forma estándar antes de ser resuelto numéricamente. La forma estánder es:\n",
    "\n",
    "$y' = f(y, t)$\n",
    "\n",
    "donde \n",
    "\n",
    "$y = [y_1(t), y_2(t), ..., y_n(t)]$ \n",
    "\n",
    "y $f$ es una función que determina las derivadas de la función $y_i(t)$. Para resolver la EDO necesitamos conocer la función $f$ y una condición inicial, $y(0)$.\n",
    "\n",
    "Note que EDOs de orden superior siempre pueden ser escritas en esta forma introduciendo nuevas variables para las derivadas intermedias.\n",
    "\n",
    "Una vez definida la función `f` y el arreglo `y_0`, podemos usar la función `odeint`:\n",
    "\n",
    "    y_t = odeint(f, y_0, t)\n",
    "\n",
    "donde `t` es un arreglo con las coordenadas temporales para las que se resolverá el sistema de EDOs. El resultado `y_t` es un arreglo con una linea para cada punto de tiempo `t`, y donde cada columna corresponde a una solución `y_i(t)` para ese tiempo. \n",
    "\n",
    "Veremos cómo implementar `f` e `y_0` en código Python en los siguientes ejemplos."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo: Oscilador armónico amortiguado\n",
    "\n",
    "Problemas de EDO son importantes en Física Computacional, de modo que veremos un ejemplo adicional: el oscilador armónico amortiguado. Este problema está bastante bien descrito en wikipedia (en inglés): http://en.wikipedia.org/wiki/Damping.\n",
    "\n",
    "La ecuación de movimiento para el oscilador amortiguado es:\n",
    "\n",
    "$\\displaystyle \\frac{\\mathrm{d}^2x}{\\mathrm{d}t^2} + 2\\zeta\\omega_0\\frac{\\mathrm{d}x}{\\mathrm{d}t} + \\omega^2_0 x = 0$\n",
    "\n",
    "donde $x$ es la posición del oscilador, $\\omega_0$ la frecuencia, y $\\zeta$ es el factor de amortiguamiento. Para escribir esta EDO de segundo orden en la forma estándar, introducimos $p = \\frac{\\mathrm{d}x}{\\mathrm{d}t}$:\n",
    "\n",
    "$\\displaystyle \\frac{\\mathrm{d}p}{\\mathrm{d}t} = - 2\\zeta\\omega_0 p - \\omega^2_0 x$\n",
    "\n",
    "$\\displaystyle \\frac{\\mathrm{d}x}{\\mathrm{d}t} = p$\n",
    "\n",
    "En la implementación de este ejemplo agregaremos algunos argumentos extras a la función del lado derecho de la EDO, en lugar de usar variables globales como en el ejemplo anterior. Como consecuencia de los argumentos extra, necesitamos pasar un argumento clave `args` a la función `odeint`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dy(y, t, zeta, w0):\n",
    "    \"\"\"\n",
    "    El lado derecho de la EDO del oscilador amortiguado\n",
    "    \"\"\"\n",
    "    x, p = y[0], y[1]\n",
    "    \n",
    "    dx = p\n",
    "    dp = -2 * zeta * w0 * p - w0**2 * x\n",
    "\n",
    "    return [dx, dp]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# condición inicial: \n",
    "y0 = [1.0, 0.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# tiempos en los que se resolvera la EDO\n",
    "t = np.linspace(0, 10, 1000)\n",
    "w0 = 2*np.pi*1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# resuelve el sistema de EDOs para tres valores diferentes del factor de amortiguamiento\n",
    "\n",
    "y1 = odeint(dy, y0, t, args=(0.0, w0)) # no amortiguado\n",
    "y2 = odeint(dy, y0, t, args=(0.2, w0)) # subamortiguado\n",
    "y3 = odeint(dy, y0, t, args=(1.0, w0)) # amortiguado crítico\n",
    "y4 = odeint(dy, y0, t, args=(5.0, w0)) # sobreamortiguado\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(t, y1[:,0], 'k', label=\"no amortiguado\", linewidth=0.25)\n",
    "ax.plot(t, y2[:,0], 'r', label=\"subamortiguado\")\n",
    "ax.plot(t, y3[:,0], 'b', label=u\"amortiguado crítico\")\n",
    "ax.plot(t, y4[:,0], 'g', label=\"sobreamortiguado\")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Lectura adicional\n",
    "\n",
    "* http://www.scipy.org - La página oficial del proyecto SciPy.\n",
    "* http://docs.scipy.org/doc/scipy/reference/tutorial/index.html - Un tutorial sobre cómo comenzar a usar SciPy. \n",
    "* https://github.com/scipy/scipy/ - El códifo fuente de SciPy. \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}