IntroMC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sesión 1. Introducción a Monte Carlo\n",
    "\n",
    "## 1.2  Fundamentos del método Monte Carlo \n",
    "\n",
    "Monte Carlo es un método de integración basado en el muestreo de variables estocásitcas o aleatorias, si tenemos una integral:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "G(x)=\\int_{\\Omega} f(x)dx \n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Esta integral puede ser escrita en terminos de una función $p(x)$:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "G(x)=\\int_{\\Omega} p(x)\\frac{f(x)}{p(x)}dx \\rangle\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    " con \n",
    " \n",
    " $$\\int p(x) dx =1 $$ \n",
    " y \n",
    " $$p(x)\\geq 0, \\thinspace x\\in\\Omega$$\n",
    "\n",
    "\n",
    "Donde $p(x)$ es una distribución de probabilidad, reemplazando por $g(x)=\\frac{f(x)}{p(x)}$, sin perdida de generalidad la integral $G(x)$, el valor esperado: \n",
    "\n",
    "$$\\begin{equation}\n",
    "G(x)=\\int_{\\Omega} p(x) g(x) dx = \\langle g(x)\\rangle\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "La integral es interpretada como el valor promedio de $\\langle g(x)\\rangle$, note que este tipo de integral es similar al valor esperado de un observable fisico en mecanica estadistica.  \n",
    "\n",
    "La varianza de $g(x)$ con respecto a la districuion $p(x)$ es definida:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\sigma^2 & =  \\mathrm {Var}[g(x)]= \\int (g(x)-\\langle g(x)\\rangle)^2 p(x) dx \\\\\n",
    "&  = \\int g(x)^2 p(x) dx - 2 \\langle g(x)\\rangle \\int g(x) p(x) dx - \\langle g(x)\\rangle^2 \\int p(x) dx \\\\\n",
    "& = \\int g(x)^2 p(x) dx - 2 \\langle g(x)\\rangle \\langle g(x)\\rangle  - \\langle g(x)\\rangle^2 \\\\\n",
    "& = \\int g(x)^2 p(x) dx  - \\langle g(x)\\rangle^2\\\\\n",
    "& = \\langle g(x)^2\\rangle - \\langle g(x)\\rangle^2\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "El error estandar es definido como: \n",
    "\n",
    "$$\n",
    "\\sigma_{error}=\\frac{\\sigma}{\\sqrt{n}}\n",
    "$$\n",
    "\n",
    "\n",
    "Finalmente, la __ley de numeros grandes__ es una de las piedras angulares de los algoritmos Monte Carlo, esta dice: \n",
    "\n",
    "Si $x_{i}$ son un conjunto de numeros aleatorios independientes identicamente distribuidos de acuedo a una distrbucion $p(x)$ se cumple que: \n",
    "\n",
    "$$\\langle g(x) \\rangle =\\int_{a}^{b} g(x) p(x) \\thinspace\\simeq\\thinspace\\frac{1}{n} \\sum_{i} g(x_{i})$$\n",
    "\n",
    "La integral puede ser calculada aproximadamente como el promedio de la función sobre variables aleatorias distribuidas segun $p(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.1 Ejemplo integrando usando la ley de números grandes\n",
    "\n",
    "Evaluemos la integral \n",
    "\n",
    "$$\\int_{0}^{1} e^{x} dx = e^{x}\\bigg\\rvert_{0}^{1}= e-1 = 2.7182-1=1.7182$$\n",
    "\n",
    "Ahora podemos hacer la evaluación usan  usando la Ley de números grandes que es la piedra angular del método de Monte Carlo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        10 1.674289\n",
      "       100 1.722180\n",
      "      1000 1.723951\n",
      "     10000 1.717775\n",
      "    100000 1.716328\n",
      "   1000000 1.718227\n",
      "  10000000 1.718104\n",
      " 100000000 1.718188\n"
     ]
    }
   ],
   "source": [
    "#Cargue desde el inicio las bibliotecas necesarias para no lamarlas de nuevo a lo largo del cuaderno jupyter. \n",
    "import numpy as np\n",
    "#%matplotlib inline\n",
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "nsamples=10**np.array([1,2,3,4,5,6,7,8])\n",
    "\n",
    "\n",
    "for n in nsamples:\n",
    "    x = np.random.uniform(0, 1, n) #x viene de una distribucion uniforme de numeros independientes \n",
    "    sol = np.mean(np.exp(x))\n",
    "    print('%10d %.6f' % (n, sol))\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.2 Teorema del Límite Central y Desigualdad de Chebyshev\n",
    "\n",
    "El teorema del límite central nos dice que un conjunto de variables aleatorias $x_{i}$, independientes e idénticamente distribuidas, si la substraemos su media a $y_{i}=x_{i}-\\langle x \\rangle$, construyendo la nueva variable:  \n",
    "\n",
    "\n",
    "$$ Y=\\frac{1}{\\sqrt{n}}\\sum_{i} y_{i} = \\frac{1}{\\sqrt{n}}\\left(\\frac{1}{n}\\sum_{i} x_{i}-\\langle x \\rangle\\right)$$\n",
    "\n",
    "Por cosntruccion el valor esperado de esta nueva variable es nulo $\\langle Y \\rangle =0$, el teorema del límite central nos dice que la distribución de probabilidad de esta nueva variable $Y$ es una distribución normal o gaussiana con media 0 y varianza $\\sigma_{D}^2$:   \n",
    "\n",
    "$$p(Y_{\\lim n\\to \\infty}) =\\frac{1}{\\sqrt{2\\pi\\sigma_{D}^2}}\\exp\\left(\\frac{-Y^2}{2}\\right) $$\n",
    "\n",
    "La utilidad del límite central es dada porque no importa cual es la distribución de $x_{i}$, siempre y cuando todas sean identicamente distribuidas e independientes. Por esto podemos relacionar el teorema con la aproximación de los números largos a la integral. La ley de números largos nos dice que el valor de nuestra integral es el promedio $\\bar{x}=\\sqrt(n)(\\frac{1}{n}\\sum_{i} x_{i}$, relacionando este preomedio con la variable Y por definicion tenemos $\\bar{x}=\\langle x\\rangle+Y$m, haciendo ese cambio de variable la distribución de $\\bar{x}$, tiene media $\\langle x\\rangle$ y varianza $\\frac{\\sigma^2}{n}$:\n",
    "\n",
    "$$p(\\bar{x}) =\\sqrt{\\frac{N}{2\\pi\\sigma^2}}\\exp{\\frac{-(\\bar{x}-\\langle x\\rangle)^2}{2}} $$\n",
    "\n",
    "\n",
    "De esta manera encontramos el error de $\\bar{x}$, sin perdida de generalidad de $\\overline{g(x)}$, el cual va a tener varianza $\\sigma_{g(x)}=\\frac{\\sigma^2}{n}$. Ademas $\\sigma_{g(x)}$ esta normalmente distribuido y indiferentemente es asintoticamente estimada.\n",
    "\n",
    "Finalmente, para relacionar la integral con su distribución $p(\\bar{x})$ podemos usamos la desigualdad de Chebyshev para mostrar que el valor de $\\bar{x}$:\n",
    "\n",
    "\n",
    "$$p(|x-\\langle x \\rangle|>\\tilde{\\epsilon})\\leq \\frac{\\sigma_{D}^2}{\\epsilon^2}$$\n",
    "\n",
    "con $\\tilde{\\epsilon}=\\epsilon\\sigma_{D}$, aplicado a nuestro problema $x=\\bar{x}$ y $\\sigma_D = \\frac{\\sigma^2}{n}$ :\n",
    "\n",
    "$$p(|\\bar{x}-\\langle x \\rangle|>\\tilde{\\epsilon}\\sigma)\\leq \\frac{\\sigma^2}{ n^2\\tilde{\\epsilon}^2} $$ \n",
    "\n",
    "Podemos interpretar esta como la probabilidad de desviarse del resultado dado un valor $\\epsilon$ fijo, la probabilidad de desvise un cierto $\\epsilon$, es inversamente proporcional a la varianza entre el dado $\\epsilon$. Es decir si la varianza es del mismo orden de magnitud que $\\epsilon$ es altamente probable que que el desvio sea de esa magnitud.      \n",
    "\n",
    "\n",
    "### 1.2.3 Muestreo Directo  Monte Carlo, integración de $\\pi$\n",
    "Uno de los primeros uso de Monte Carlo conocido historicamente fue el uso de Monte Carlo en la integración de $\\pi$, la idea basica es calcular $\\pi$ como la razón del area de un cuadrado de lado $l$ y un circulo inscrito dentro del cuadrado de radio $r=\\frac{l}{2}$\n",
    "\n",
    "### 1.2.3.1 Hit or Miss (Acierta o desacierta)\n",
    "La idea es que $\\pi$ es proporcional al número de puntos debajo de una curva, para este caso, se basa en la relacion geometica de las areas de un cuadrado de lado $i$ que inscribe un circulo de radio $l/2$:\n",
    "\n",
    "$$\\frac{\\pi}{4} =\\frac{\\int_{Circulo, r=l/2} dx dy}{\\int_{Cuadrado, l} dx dy}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Dibujando el cuadrado de lado l centrado en (0,0)\n",
    "#puntos\n",
    "l =4.0; #longitud lado cuadrado\n",
    "s=np.arange(-0.5, 0.6125, 0.125);\n",
    "#print 's',s\n",
    "npt=s.shape[0]\n",
    "#print npt\n",
    "c=np.zeros(npt); #vector de ceros\n",
    "u=np.ones(npt); # vector de unos\n",
    "\n",
    "#generando los pares coordenados del perimetro del cuadrado\n",
    "x=np.concatenate((l*s,0.5*l*u,l*s[::-1],-0.5*l*u), axis=0)\n",
    "#print 'x',x\n",
    "y=np.concatenate((0.5*l*u,s[::-1]*l,-.5*l*u,s*l), axis=0);\n",
    "#print 'y',y\n",
    "plt.plot(x,y,'b-',linewidth=3);\n",
    "\n",
    "#dibuja el perimetro del cuadrado\n",
    "r=0.5*l;\n",
    "theta=np.arange(0.0,2.05*np.pi,0.05);\n",
    "rx=r*np.cos(theta);\n",
    "ry=r*np.sin(theta);\n",
    "plt.plot(rx,ry,'r-',linewidth=2);\n",
    "\n",
    "plt.axis('equal');\n",
    "plt.xlim(-l,l);\n",
    "plt.ylim(-l,l);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Luego dada de relación geometrica de las areas se puede inferir una relación de proporcionalidad, lanzando puntos aleatorios dentro del cuadrado y al contar el número de puntos dentro del circulo y todos los números dentro del cuadrado la razon de estos es proporcional al area. \n",
    "\n",
    "$$\\frac{\\pi}{4} \\propto \\frac{Area \\thinspace Circulo}{Area \\thinspace Cuadrado}\\propto \\frac{n_{circulo}}{n_{cuadrado}}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#generando una distribucion 2d de puntos uniformes\n",
    "nm=2000; #numero de muestras\n",
    "sx = np.random.uniform(-0.5*l,0.5*l,nm)\n",
    "sy = np.random.uniform(-0.5*l,0.5*l,nm)\n",
    "plt.plot(sx,sy,'g.',markersize=1.5)\n",
    "\n",
    "#Contar el numero de puntos dentro del circulo\n",
    "nmc=0;\n",
    "for i in range(0,nm):\n",
    "    rp=np.sqrt((sx[i]*sx[i])+(sy[i]*sy[i]));\n",
    "    #print rp\n",
    "    if rp<=r:\n",
    "        nmc=nmc+1;\n",
    "\n",
    "print \"Valor aproximado de pi \", (4.0*nmc)/nm\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#Implemente una funcion para calcular pi en funcion del numero de muestras\n",
    "def direct_pi(nint):\n",
    "    sx = np.random.uniform(-0.5*l,0.5*l,nint)\n",
    "    sy = np.random.uniform(-0.5*l,0.5*l,nint)\n",
    "    #Contar el numero de puntos dentro del circulo\n",
    "    nhits=0;\n",
    "    for i in range(0,nint):\n",
    "        rp=np.sqrt((sx[i]*sx[i])+(sy[i]*sy[i]));\n",
    "        if rp<=r:\n",
    "            nhits=nhits+1;\n",
    "    return (4.0*nhits)/nint;\n",
    "\n",
    "\n",
    "anp=np.arange(100,13100,100);\n",
    "vpi=np.array([ direct_pi(i) for i in anp ]);\n",
    "\n",
    "#closing last plot\n",
    "plt.close()\n",
    "\n",
    "\n",
    "#starting new plot of the density \n",
    "plt.plot(anp,vpi,\"g^-\",linewidth=2.0,markersize=8.0)    \n",
    "plt.axhline(y=np.pi,color='r',linewidth=4)\n",
    "plt.xlabel(\"Numero de muestras\")\n",
    "plt.ylabel(\"Estimado de $\\pi$ \")\n",
    "plt.title(\"Ley de numeros grandes\")\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.4. Reduciendo la Varianza, Muestreo de Importancia (Impotance Sampling)\n",
    "\n",
    "La varianza puede ser afectada por la elección de la función de distribución $p(x)$, supongamos caso mas general en que se integra $g(x)$, usando la Ley de números largos:  \n",
    "\n",
    "$$\n",
    "\\sigma^2 = \\frac{1}{N} \\sum_{i} g(x_i)^2  -\\left[\\frac{1}{N} \\sum_{i} g(x_i)\\right]^2\n",
    "$$\n",
    "\n",
    "Luego para valores grandes de $g(x)$ la varianza aumentara, la varianza aumentara mas aun en el caso que $p(x)>>0$ es alta para un $x$ donde $g(x)>>0$. \n",
    "\n",
    "Una manera de disminuir las variaciones es multiplicar por una densidad de probabilidad $h(x)\\geq p(x)g(x)$ en los valores donde mas fluctua. \n",
    "\n",
    "Si multiplicamos por $\\frac{h(x)}{h(x)}$ otra densidad de probabilidad, tal que nuestra nueva función $$\\tilde{g}(x)=\\frac{p(x)g(x)}{h(x)}$$. Entonces la varianza puede ser estimada al reemplazar $g$ con $\\tilde{g}$, escribiendo explicitamente: \n",
    "\n",
    "$$\n",
    "\\sigma^2 = \\frac{1}{N} \\sum_{i} \\left( \\frac{p(x)g(x)}{h(x)} \\right)^2  -\\left[\\frac{1}{N} \\sum_{i} \\left( \\frac{p(x)g(x)}{h(x)} \\right)\\right]^2\n",
    "$$\n",
    "\n",
    "De esta manera por defición $$\\frac{p(x)g(x)}{h(x)}<<1$$ lo cual disminuye las fluctuaciones y por tanto la varianza. \n",
    "\n",
    "## 1.3 Monte Carlo usando Cadenas de Markov (Caminos aleatorios)\n",
    "\n",
    "En el problema descrito anteriormente, podemos observar que deben existir varias formas de muestrear la distribución de tal manera que la eficiencia es aumentada esto significa que un mayor número de puntos sean muestreados dentro del circulo. \n",
    "\n",
    "Una manera de hacer esto es generando un camino aleatorio de puntos en secuencia que inicien en algun punto del cuadrado y se muevan aleatoriamente de tal manera que despues de un tiempo el número de puntods dentro del circulo sea mayor que el número de puntos fuera de este. Esta tecnica es conocida como una cadena de markov una secuencia de puntos o estados dentro del espacio de configuración muestreando una cierta distribución de probabilidad. Esta tecnica ha llegado a ser la parte central del algoritmo de Metropolis Monte Carlo, el cual es uno de los algoritmos predominantes en los que esta basado, la gran mayoria de metodos Monte Carlo.  \n",
    "\n",
    "Para desarrollar una mejor intuición acerca de las cadenas de Markov, estudiemos el siguiente algoritmo para calcular $pi$, imaginen que en vez de solo lanzar puntos dentro del cuadrado y contar los que estan dentro del circulo. Vamos a partir de un punto inicial y movernos por un vector aleatorio uniforme $\\delta=(\\delta_x,\\delta_y)$ y aceptar el nuevo punto si esta dentro del cuadrado. De lo contrario el nuevo punto de la configuración sera el punto anterior que estaba dentro del cuadraro. De nuevo se cuentan los puntos dentro del circulo y del cuadrado para estimar $\\pi$.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Cerrar Figura anterior\n",
    "plt.close()\n",
    "#x,y,rx,ry,l,r = 0.5*x,0.5*y,0.5*rx,0.5*ry, 0.5*l,0.5*r;\n",
    "\n",
    "#Dibujando el circulo y el cuadrado de nuevo\n",
    "plt.plot(x,y,'b-',linewidth=3);\n",
    "plt.plot(rx,ry,'r-',linewidth=2);\n",
    "plt.axis('equal')\n",
    "plt.xlim(-l,l)\n",
    "plt.ylim(-l,l)\n",
    "\n",
    "\n",
    "delta=0.5;\n",
    "nm=4000;\n",
    "nmc=0;\n",
    "#l=r;\n",
    "#r=1.0;\n",
    "sx,sy = 0.0*l,0.0*l;\n",
    "l2=0.5*l;\n",
    "\n",
    "dx=np.random.uniform(-0.5*delta,0.5*delta,nm);\n",
    "dy=np.random.uniform(-0.5*delta,0.5*delta,nm);\n",
    "for i in range(0,nm):\n",
    "    \n",
    "    if sx>l2:\n",
    "        print sx\n",
    "    if sy>l2:\n",
    "         print sy\n",
    "    #print dx, dy\n",
    "    xt=sx+dx[i];\n",
    "    yt=sy+dy[i];\n",
    "    sxold=sx;\n",
    "    syold=sy;\n",
    "    if np.abs(xt)<l2 and np.abs(yt)<l2:\n",
    "        sx=xt; sy=yt;   \n",
    "    rp2=((sx*sx)+(sy*sy));\n",
    "    #print rp\n",
    "    r2=r*r;\n",
    "    if rp2 < r2:\n",
    "        nmc+=1;\n",
    "        #plt.arrow(sxold, syold, sx, sy, head_width=0.005, head_length=0.1, fc='g', ec='g')\n",
    "    plt.plot(sx,sy,'yo')\n",
    "        #print nmc\n",
    "        #plt.show()\n",
    "        #print sx,sy\n",
    "          \n",
    "        \n",
    "print nmc,nm, nmc/nm\n",
    "print \"Valor aproximado de pi \", (4.0*nmc)/nm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De que factores depende el muestreo, haga una función que calcule pi, basado en el algoritmo anterior que solamente dependa del valor $\\delta$ y el número de intentos, que entregue el número de movimientos aceptados dentro del cuadrado y el valor estimado de $\\pi$.  \n",
    "\n",
    "Genere diferentes gráficas de valores en función de diferentes valores de $\\delta$ y diga cual es el valor óptimo, para calcular pi y trate de explicarlo. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "\n",
    "#define python function markov pi\n",
    "#import random\n",
    "def markov_pi(dl,nrun):\n",
    "    nmkc,naccpt,sx,sy=0,0,r,r;\n",
    "    dx=np.random.uniform(-delta,delta,nrun);\n",
    "    dy=np.random.uniform(-delta,delta,nrun);\n",
    "    for i in range(0,nrun):\n",
    "        if abs(sx+dx[i])<0.5*l and abs(sy+dy[i])<0.5*l:\n",
    "            sx=sx+dx[i]; sy=sy+dy[i];  \n",
    "            naccpt+=1\n",
    "        rp2=((sx*sx)+(sy*sy));\n",
    "        r2=r*r;\n",
    "        if rp2 < r2:\n",
    "            nmkc+=1;\n",
    "            ##plt.arrow(sxold, syold, sx, sy, head_width=0.005, head_length=0.1, fc='g', ec='g')\n",
    "            ##plt.plot(sx,sy,'yo')\n",
    "    #print nmkc,naccpt,(4.0*nmkc)/nrun;\n",
    "    return naccpt,(4.0*nmkc)/nrun; \n",
    "        \n",
    "#DL=l*np.arange(0.05,1.0,0.05);\n",
    "\n",
    "DL=np.arange(0.01*l,0.5*l,0.005*l);\n",
    "\n",
    "#for ii in DL:\n",
    "#    markov_pi(ii,4000)\n",
    "NT=2000;\n",
    "data=np.array([markov_pi(ii,NT) for ii in DL]);\n",
    "data.shape\n",
    "\n",
    "plt.plot(DL,data[:,0]/NT,'rs-', label=\"Razon de Aceptacion\")\n",
    "plt.plot(DL,data[:,1]/np.pi,'gs-',label=\"$\\pi$\")\n",
    "plt.xlabel(\"$\\delta$\")\n",
    "plt.ylabel(\"scaled values\")\n",
    "plt.legend(loc=\"best\")\n",
    "#maxElement = np.amax(data[:.0])\n",
    "\n",
    "#print max(data[:,0])/NT\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3.1 Cadenas de Markov\n",
    "\n",
    "Las cadenas de Markov consisten en una secuencia sucesiva de estados $\\left(x_{1},x_{2},\\cdots,x_{n}\\right)$ y sus respectivas probabilidades de transición que solo dependen del estado anterior es decir $\\mathbf{p}(x_{n}|x_{n-1})$, esta es llamada la matriz de probabilidades de transición y cumple los siguientes requerimientos: \n",
    "\n",
    "$$\\mathbf{p}(x_{k}|x_{k-1})\\geq 0 $$ \n",
    "\n",
    "y\n",
    "\n",
    "$$\\sum_{x_{k-1}}\\mathbf{p}(x_{k}|x_{k-1})=1$$ \n",
    "\n",
    "Si queremos usar la candena de Markov para muestrear la distribucion de probabilidad de equilibrio $P_{eq}(x)$, la transición de probabilidad debe satisfacer la condición estacionaria:\n",
    "\n",
    "$$\\sum_{i}\\mathbf{p}(x_{f}|x_{i})P_{eq}(x_{i})=\\sum_{i}\\mathbf{p}(x_{i}|x_{f})P_{eq}(x_{f})=P_{eq}(x_f)$$\n",
    "\n",
    "\n",
    "Esta condición requiere que si comienza de una condición inicial arbitraria $P_{eq}(x)$ despues de un repetido muestreo de la transición de probabilidad debe seguir muestrando la distribucion $P_{eq}(x)$. Otra condición que debe cumplir es que debe llegar a la distribucion deseada en un número finito de pasos. Esto significa que:\n",
    "\n",
    "$$\\lim_{n\\to\\infty} P(0)\\mathbf{p}^n=P_{eq}$$\n",
    "\n",
    "Es decir hay una distribucion limitante de la canena de markov dada una matriz de trnasicion y queremos que esa sea la distribucion de equilibrio. \n",
    "\n",
    "\n",
    "Por ejemplo si tenemos uma matriz de transición de probabilidad $\\mathbf{p}$ para tres estados:\n",
    "\n",
    "$$\\mathbf{p}=\\begin{bmatrix}\n",
    "0.3 & 0.25 & 0.45\\\\\n",
    "0.4 & 0.5 & 0.1\\\\\n",
    "0.6 & 0.1 & 0.3\\\\\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "Parte de la condición incial\n",
    "\n",
    "$$P(0)=\\begin{bmatrix}\n",
    "0.3 & 0.4 & 0.3\\\\\n",
    "\\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Matriz Transicion probabilidadndad \n",
    "p=np.array([[0.3, 0.25, 0.45],[0.4, 0.5,0.1],[0.6,0.1,0.3]]);\n",
    "P0=np.array([0.1,0.5,0.4]);\n",
    "\n",
    "P=P0;\n",
    "for i in range(0,15):\n",
    "    P=np.dot(P,p);\n",
    "    print P    \n",
    "    \n",
    "##cambie la condicion inicial y repital el mismo proceso    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La distribución de probabilidad $P_{eq}$, fue independiente de la condición inicial, además es un autovector de la matriz de la transición de probabilidad con autovalor 1. Los otros autovectores tienen autovalor menor que 1 y controntran la razón de transición a la probabilidad limite o de equilibrio. Esta es una consecuencia del teorema Perron-Frobenius el cual aplica a matrices no-negativas y ergodicas. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Ejercicio compruebe que Peq es un autovector de la matriz de probabilidad de transicion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3.2 Metropolis Monte Carlo\n",
    "\n",
    "Sin embargo con las dos condiciónes anteriores hay multiples elecciones para la transición de probabilidad.\n",
    "El metodo de Metropolis impone __la condición de balance detallado__: \n",
    "\n",
    "$$\\mathbf{p}(x_{f}|x_{i})P_{eq}(x_{i})=\\mathbf{p}(x_{i}|x_{f})P_{eq}(x_{f})$$\n",
    "\n",
    "Esta condición exige ergodicidad haciendo que tenga el mismo número de estados de transición cuando va del estado $x_{i}$ a $x_{f}$ y viceversa.\n",
    "\n",
    "Además, esta condición impone que la matriz limite sea $P_{eq}$, es decir con esto podemos construir una matriz transición de probailidad que converge a la distribucion de equilibrio que se requiera. \n",
    "\n",
    "Sin embargo, aun con las restricciones regulares de una cadena de Markov y la condición de balance detallado hay mucha libertad para escoger como construir $\\mathbf{p}$. Mettropolis propuso escribir la probabilidad de transición como el producto de transición de acceptacion y la probabilidad de prueba o proposicion:\n",
    "\n",
    "$$\\mathbf{p}(x_{f}|x_{i}) = p_{accept}(x_{f}|x_{i})T(x_{f}|x_{i})$$\n",
    "\n",
    "Reescrbiendo la condición de balance detallado tenemos: \n",
    "\n",
    "$$p_{accept}(x_{f}|x_{i})T(x_{f}|x_{i})P_{eq}(x_{i})=p_{accept}(x_{i}|x_{f})T(x_{i}|x_{f})P_{eq}(x_{f})$$\n",
    "\n",
    "reagrupando los terminos de la condiciónd e balance detallado tenemos:\n",
    "\n",
    "$$\\frac{p_{accept}(x_{f}|x_{i})}{p_{accept}(x_{i}|x_{f})}=\\frac{T(x_{i}|x_{f})P_{eq}(x_{f})}{T(x_{f}|x_{i})P_{eq}(x_{i})}$$\n",
    "\n",
    "Ahora si \n",
    "\n",
    "$$\\frac{T(x_{i}|x_{f})P_{eq}(x_{f})}{T(x_{f}|x_{i})P_{eq}(x_{i})}>1$$\n",
    "\n",
    "Sera aceptado con probabilidad\n",
    "\n",
    "$$p_{accept}(x_{f}|x_{i})=1$$\n",
    "\n",
    "$$p_{accept}(x_{i}|x_{f})=\\frac{T(x_{f}|x_{i})P_{eq}(x_{i})}{T(x_{i}|x_{f})P_{eq}(x_{f})}$$\n",
    "\n",
    "Lo cual satisface la condición de balance detallado, de manera similar cuando:\n",
    "\n",
    "$$\\frac{T(x_{i}|x_{f})P_{eq}(x_{f})}{T(x_{f}|x_{i})P_{eq}(x_{i})}<1$$\n",
    "\n",
    "La probabilidad de aceptación: \n",
    "\n",
    "$$p_{accept}(x_{f}|x_{i})=\\frac{T(x_{i}|x_{f})P_{eq}(x_{f})}{T(x_{f}|x_{i})P_{eq}(x_{i})}$$\n",
    "\n",
    "$$p_{accept}(x_{i}|x_{f})= 1$$\n",
    "\n",
    "Este criterio puede ser resumido como \n",
    "\n",
    "$$p_{accept}(x_{i}|x_{f})= Min\\left\\{1,\\frac{T(x_{i}|x_{f})P_{eq}(x_{f})}{T(x_{f}|x_{i})P_{eq}(x_{i})}\\right\\}$$\n",
    "\n",
    "\n",
    "Por ultimo queda a determinar la eleccion de de la probabilidad de prueba $T(x^{'}|x)$, en el algoritmo original de metropolis propone que esta sea simetrica:\n",
    "\n",
    "$$T(x_{i}|x_{f}) =T(x_{f}|x_{i})$$\n",
    "\n",
    "Luego el criterio de aceptación del algoritmo de Metropolis depende solo de $P_{eq}$ \n",
    "\n",
    "$$p_{accept}(x_{i}|x_{f})= Min\\left\\{1,\\frac{P_{eq}(x_{f})}{P_{eq}(x_{i})}\\right\\}$$\n",
    "\n",
    "\n",
    "Cuando $$T(x_{i}|x_{f}) \\neq T(x_{f}|x_{i})$$ es conocido como el __algorithmo de Metropolis-Hastings__, una adecuada eleccion de la probabilidad de prueba puede ser de gran beneficio para mejorar el muestreo de la distribucion. \n",
    "\n",
    "\n",
    "### 1.3.3 Receta Algoritmo de Metropolis\n",
    "\n",
    "\n",
    "1. Proponer un movimiento al generar una configuracion $x'$ de acuerdo a la probabilidad de proposicion T(x'|x_{n})\n",
    "\n",
    "\n",
    "2. Acepte o rechaze el movimiento, se genera un número aleatorio $\\eta\\in[0,1)$ si $p_{accept}(x'|x_{n})>1$ es aceptado  $x_{n+1}=x'$, de lo contrario es rechazado y $x_{n+1}=x_{n}$. \n",
    "\n",
    "En el caso que $x'=x+\\zeta$, con cita un número aleatorio uniformemente distribuido entre $[-a,a]$, $T(x^{'}|x)=\\frac{1}{2a}$ luego $$T(x'|x) =T(x|x')$$ que es la misma eleccion de Metropolis. El hecho que el criterio de aceptación solo quede dependiendo de la razon de la probabilidad de equilibrio entre la configuracion anterior y la propuesta, es una simplicacion substancial. De hecho ya que por elemplo en el emsamble canonico la distribucion de probabilidad de un sistema es la distribucion de Maxwell-Boltzmann: \n",
    "\n",
    "$$P_{eq}=\\frac{\\exp\\left(-\\beta E\\right)}{Z}$$ donde $$ Z = \\int \\exp\\left(-\\beta E\\right)  dp dq $$ $\\beta=\\frac{1}{k_{B}T}$\n",
    "\n",
    "La mayor parte de las veces no podemos calcular Z directamente, de lo contrario no habria necesidad de usar Monte Carlo para muestrear la distribucion y encontrar la distribucion de equilibrio. \n",
    "\n",
    "\n",
    "## 1.4 El Modelo de Ising \n",
    "\n",
    "El modelo de Ising fue introducido inicialmente para explicar la observación la observación de magnetización expontanea cuando un material ferrogmanético es enfriado por la valores debajo de la temperatura de Curie. Desde entonces historicamente ha desempeñado un rol importante en Física como un ejemplo de transiciónes de fase. Además, su modificación, extensión y generalización a diferentes situaciones ha sido particularmente prolifica en Física ayudando a entender bastas y diferentes situaciones en multiples sitemas físicos.   \n",
    "\n",
    "Basicamente es una red en la que cada sitio tiene un spin up-down $\\pm 1$ e interactua solamente con sus primeros vecinos como muestra la figura de abajo, es decir los que estan solo a un sitio de red de distancia. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "#grid LxL\n",
    "L=6;\n",
    "x=np.arange(1,L+1)\n",
    "S = [np.random.choice([1, -1]) for k in range(N)]\n",
    "\n",
    "for i in x:\n",
    "    for j in x:\n",
    "        plt.plot(i,j,'ro')\n",
    "        plt.arrow(i, j-(S[i+(j*L)]*0.25), 0.0, (S[i+(j*L)]*0.25), head_width=0.05, head_length=0.25, fc='g', ec='g')\n",
    "\n",
    "#plt.grid(b=True, which='major', color='#666666', linestyle='-')\n",
    "\n",
    "        \n",
    "plt.arrow(L/2, L/2, 0, 0.7, head_width=0.05, head_length=0.25, fc='b', ec='b')\n",
    "plt.arrow(L/2, L/2, 0.7, 0, head_width=0.05, head_length=0.25, fc='b', ec='b')\n",
    "plt.arrow(L/2, L/2, -0.7, 0, head_width=0.05, head_length=0.25, fc='b', ec='b')\n",
    "plt.arrow(L/2, L/2, 0, -0.7, head_width=0.05, head_length=0.25, fc='b', ec='b')\n",
    "#plt.show()        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El hamiltoniano es definido como:\n",
    "\n",
    "$$E=-J\\sum_{<ij>}s_{i}s_{j}-\\mathit{h} \\sum_{i}^{N} s_i $$\n",
    "\n",
    "Aqui $<ij>$ denota la suma sobre los vecinos cercanos, $J$ es el parámetro que cuantifica la energía de interación mejor conocido como la energía de intercambio. $\\mathit{h}$ es el campo magnético externo, es facil observar que el campo rompe la simetria de la enérgia de interación entre los sitios, es decir da una preferencia a minimizar la energía dependiendo de la orientacion del spin de los atomos en la red. Disminuyendo la energía en el caso que el spin es paralelo al campo magnético externo.   \n",
    "\n",
    "La diferencia entre el caso 1d y 2d, esta definido solamente por la interacción con sus vecinos, en el caso 1d cada sitio de red interactuaria solamente con dos vecinos uno de cada lado. En el caso 2d interactua adicionalmente con los vencinos de arriba y abajo. Lars Onsager resolvio el modelo de Ising en una red cuadrdada de forma análitica con campo magnetico externo cero. Osanger mostro que la hay una transición de fase en la magnetización con el cambio de temperatura critica $T_c$  \n",
    "\n",
    "$$\n",
    "\\sinh\\frac{2J}{k_BT_c}=1\n",
    "$$\n",
    "\n",
    "\n",
    "Para esta parte del tutorial queremos mostrar como usar el algoritmo de Metropolis para mostrar como ocurre la transición de fase con la temperatura, además de observar como Cambia la energía con la temperatura el valor de $T_c$. Para esto debemos construir el algoritmo Metropolis en el ensamble canonico, para la red cuadrada del modelo de Ising. \n",
    "\n",
    "__Algoritmo Metropolis Monte Carlo del Modelo Ising __\n",
    "\n",
    "0. Inicialize una red cuadrada de spins aleatorialmente. \n",
    "1. Cambie el spin de un sitio de la red Escogido aleatoriamente.\n",
    "2. Calcule el cambio de la energía en la red al cambiar el spin.\n",
    "3. Use la probabilidad del sistema canonico con el criterio de Metropolis para aceptar o rechazar el movimiento.\n",
    "4. Repita 1-3.\n",
    "\n",
    "Antes de continuar es posible observar que los movimientos son locales es decir solo cambian la energía de intercambio con sus vecinos cercanos, la energía de un solo sitio en la red es:\n",
    "\n",
    "$$\n",
    "E_i = s_i (J\\sum_{j}s_j+h)\n",
    "$$\n",
    "\n",
    "Al cambiar el spin de un sitio $i$ de $\\pm 1$ a $\\mp 1$, el cambio de energía es\n",
    "\n",
    "$$\n",
    "\\Delta E = 2s_i (J\\sum_{j}s_j+h)\n",
    "$$\n",
    "\n",
    "__Observables__\n",
    "\n",
    "En el ensamble canonico el valor esperado de un observable $O$ a una temperatura T, esta dado por:\n",
    "\n",
    "$$\n",
    "\\langle O\\rangle=\\frac{\\sum_{i\\in\\mathcal{S}}O_ie^{-\\frac{E_\\gamma}{k_bT}}}{\\sum_{\\gamma\\in\\mathcal{S}}e^{-\\frac{E_\\gamma}{k_bT}}}\n",
    "$$\n",
    "\n",
    "Donde se suma sobre todas las configuraciones posibles $\\gamma$, lo cual crece factorialmente con el tamaño de la red, solo para el caso 2D existe la solución de Osanger en el limite termodinamico. El valor esperado del observable es estimado como el promedio sobre las configuraciones despues de un tiempo de termalizacion. \n",
    "\n",
    "$$\n",
    "\\langle O\\rangle\\approx\\frac{1}{\\#\\, of\\, Monte\\, Carlo\\, Steps}\\sum{O_i}\n",
    "$$\n",
    "\n",
    "La magnetizacion es calculada promediando el spin de los sitios de la red, para las configuracion generadas, si $N=L\\times L$, el numero de sitios en la red cuadrada: \n",
    "\n",
    "\n",
    "$$\n",
    "M  = \\frac{1}{N}\\sum_{i=1}^{N}s_i\n",
    "$$\n",
    "\n",
    "La energía por spin:\n",
    "$$\n",
    "\\epsilon_i = \\frac{\\frac{-J}{2}\\sum_{(i,j)}s_is_j-h\\sum_i s_i}{N}\n",
    "$$\n",
    "\n",
    "El factor $\\frac{1}{2}$ es para evitar soble conteo de la energía sobre pares de espin. \n",
    "\n",
    "Desarrolle los siguientes pasos para terminar con el tutorial.\n",
    "\n",
    "1. En lo siguiente defina una red cuadrada de lado $L=16$, con una matriz que contiene los indices de los vecinos. \n",
    "\n",
    "2. Escriba el algoritmo metropolis para una temperatura fija T=0.2, grafique la configuracion inicial y final de spin asigne dos colores para $\\pm1$, escribalo como una matriz.   \n",
    "\n",
    "3. Use el algorimo metropolis dentro un loop para 20 temperaturas, para canda temperatura calcule el valor esperado de la energía y de la magnetizacion. \n",
    "\n",
    "4. Grafique la energía y magnetizacion con el cambio de temperatura, observe que pasa. \n",
    "\n",
    "Para esta parte del tutorial vamos a tomar unidades normalizadas $J=1$, $k_B=1$ y campo externo $h=0$. Ejecute  en una celula diferente cada paso. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "L = 16; \n",
    "N = L * L;  #tamanho de la red\n",
    "\n",
    "#iniciliza la red de spins\n",
    "S = [np.random.choice([1, -1]) for k in range(N)];\n",
    "\n",
    "\n",
    "#define los indices de los vecinos tomando condiciones periodicas\n",
    "nbr = {i : ((i // L) * L + (i + 1) % L, (i + L) % N,\n",
    "            (i // L) * L + (i - 1) % L, (i - L) % N)\\\n",
    "                                   for i in range(N)}\n",
    "\n",
    "##guarda la configuracion inicial co\n",
    "SQ0=np.reshape(S,(L,L));\n",
    "\n",
    "\n",
    "T=0.2;\n",
    "npasos = 10000 * N;\n",
    "beta = 1.0/T;\n",
    "for paso in range(npasos):\n",
    "        k = np.random.randint(0, N - 1)\n",
    "        delta_E = 2.0 * S[k] * sum(S[nn] for nn in nbr[k])\n",
    "        if np.random.uniform(0.0, 1.0, 1) < np.exp(-beta * delta_E):\n",
    "            S[k] *= -1\n",
    "\n",
    "\n",
    "            \n",
    "        \n",
    "SQF=np.reshape(S,(L,L));\n",
    "        \n",
    "plt.subplot(211)\n",
    "plt.imshow(SQ0,cmap='Wistia',  interpolation='nearest');\n",
    "plt.title(\"Configuracion Inicial\")\n",
    "plt.subplot(212)\n",
    "plt.imshow(SQF,cmap='Wistia',  interpolation='nearest')\n",
    "plt.title(\"Configuracion Final\")\n",
    "plt.savefig('blkwht.png')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "\n",
    "def energy(S, N, nbr):\n",
    "    E=0.0;\n",
    "    for k in range(N):\n",
    "        E -=  S[k] * sum(S[nn] for nn in nbr[k])\n",
    "    return 0.5*E;\n",
    "    \n",
    "\n",
    "#número de pasos\n",
    "npasos = 10000 * N\n",
    "#temperatura\n",
    "list_T = [1.0 + 0.2 * i for i in range(20)]\n",
    "list_av_m = []\n",
    "list_av_E = []\n",
    "\n",
    "#calcula la magenization incial\n",
    "M = sum(S)\n",
    "#calcula la energia inicial\n",
    "\n",
    "\n",
    "##haga la curva con diferentes temperaturas    \n",
    "for T in list_T:\n",
    "    print 'T =', T\n",
    "    beta = 1.0/T;\n",
    "    M_tot = 0.0;\n",
    "    E_tot = 0.0;\n",
    "    n_medidas = 0;\n",
    "    for paso in range(npasos):\n",
    "        k = np.random.randint(0, N - 1)\n",
    "        delta_E = 2.0 * S[k] * sum(S[nn] for nn in nbr[k])\n",
    "        if np.random.uniform(0.0, 1.0, 1) < np.exp(-beta * delta_E):\n",
    "            S[k] *= -1\n",
    "            M += 2 * S[k]\n",
    "        if paso % N == 0 and paso > npasos / 2:\n",
    "            M_tot += abs(M)\n",
    "            E_tot += energy(S, N, nbr);\n",
    "            n_medidas += 1\n",
    "    list_av_m.append(abs(M_tot) / float(n_medidas * N))\n",
    "    list_av_E.append(E_tot / float(n_medidas * N))\n",
    "\n",
    "\n",
    "#print list_av_E\n",
    "#print list_av_m\n",
    "\n",
    "plt.subplot(211)\n",
    "#plt.savefig('blkwht.png')    \n",
    "plt.title('$%i\\\\times%i$ lattice' % (L, L))\n",
    "plt.xlabel('$T$', fontsize=16)\n",
    "plt.ylabel('$<|M|>/N$', fontsize=16)\n",
    "plt.plot(list_T, list_av_m, 'bo-', clip_on=False)\n",
    "plt.ylim(0.0, 1.0)\n",
    "plt.subplot(212)\n",
    "#plt.savefig('blkwht.png')    \n",
    "plt.title('$%i\\\\times%i$ lattice' % (L, L))\n",
    "plt.xlabel('$T$', fontsize=16)\n",
    "plt.ylabel('$<|E|>/N$', fontsize=16)\n",
    "plt.plot(list_T, list_av_E, 'go-', clip_on=False)\n",
    "#plt.ylim(0.0, 1.0)\n",
    "\n",
    "\n",
    "#plt.show()\n",
    "#plt.savefig('plot_local_av_magnetization_L%i.png' % L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Referencias\n",
    "1. Statistical Mechanics: Algorithms and Computations, W. Krauth, Oxford University Press, 2006\n",
    "2. Quantum Monte Carlo Approaches for Correlated Systems, Federico Becca, Sandro Sorella, Cambridge University Press, 2017\n",
    "3. Monte Carlo Methods, Malvin H. Kalos, Paula A. Whitlock, John Wiley & Sons, 2009\n",
    "4. Exactly Solved Models in Statistical Mechanics, Dover books on physics, Rodney J. Baxter, Courier Corporation, 2007\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}