logi_note.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "!pip install seaborn emcee --user"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "import emcee\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.set_style('white')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Crecimiento Exponencial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\frac{dG}{dt} = rG$,\n",
    "\n",
    "donde G representa el numero de personas a un tiempo data y r es la tasa de crecimiento. La solucion general a esta ecuación diferencial es:\n",
    "\n",
    "$G = Ae^{rt}$,\n",
    "\n",
    "donde A es una normalización."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Crecimiento Logístico"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\frac{dG}{dt} = rG\\left(1-\\frac{G}{K}\\right)$,\n",
    "\n",
    "donde el parametro nuevo, K, es capacidad máxima de población que un sistema ecológico puede soportar.La solución a esta ecuación es:\n",
    "\n",
    "$G = \\frac{K}{1+Ae^{-rt}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exp_growth(pars,time):\n",
    "    r,logA = pars\n",
    "    A = 10**logA\n",
    "    return A*np.exp(r*time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logistic(pars,time):\n",
    "    r,logK,logA = pars\n",
    "    A = 10**logA\n",
    "    K = 10**logK\n",
    "    val = K / (1. + A*np.exp(-r*time))\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Growth_Models(object):\n",
    "    def __init__(self,time,model):\n",
    "        self.model = model\n",
    "        self.time = time\n",
    "\n",
    "    def __call__(self,pars):\n",
    "        if self.model == 'exp':\n",
    "            return exp_growth(pars,self.time)\n",
    "        if self.model == 'logistic':\n",
    "            return logistic(pars,self.time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lnhood(pars,data,err,time,model):\n",
    "    init_fun = Growth_Models(time,model)\n",
    "    model_val = init_fun(pars)\n",
    "\n",
    "    p = ((data - model_val)**2 / err**2) + np.log(2. * np.pi * err**2)\n",
    "    return -0.5 * np.sum(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def priors(pars,plist,model):\n",
    "    if model == 'exp':\n",
    "        r,A = pars\n",
    "        if plist[0]<r<plist[1] and plist[2]<A<plist[3]:\n",
    "            return 0.0\n",
    "        return -np.inf\n",
    "    if model == 'logistic':\n",
    "        r,K,A = pars\n",
    "        if plist[0]<r<plist[1] and plist[2]<K<plist[3] and\\\n",
    "            plist[4]<A<plist[5]:\n",
    "            return 0.0\n",
    "        return -np.inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lnpost(pars,data,err,time,plist,model):\n",
    "    p = priors(pars,plist,model)\n",
    "    if not np.isfinite(p):\n",
    "        return -np.inf\n",
    "    return p + lnhood(pars,data,err,time,model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Graficamos los datos primero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt('Datos_Poblacionales_Gdl.dat').T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(data[0],data[1],yerr=data[2],fmt='o',color='orangered');\n",
    "plt.xlabel('Year');\n",
    "plt.ylabel('Number of people');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = np.arange(0,80,np.diff(data[0])[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print len(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximizamos el Likelihood para inicializar nuestros caminadores en una zona de alta probabilidad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "empezamos con el modelo exponencial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fun = lambda *args: -lnhood(*args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize as op"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = op.minimize(fun,[0.5,6.,3.0],args=(data[1],data[2],time,'logistic'),method='TNC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print result.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(data[0],logistic(result.x,time),'g',label='Inference');\n",
    "plt.errorbar(data[0],data[1],yerr=data[2],fmt='o',color='orangered',label='data');\n",
    "plt.xlabel('Year');\n",
    "plt.ylabel('Number of people');\n",
    "plt.legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ahora sigue el MCMC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior_list = [0.,1.,0,8.0,0.,8.]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nwalkers = 60\n",
    "ndim = (len(prior_list)/2) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.zeros((ndim,nwalkers))\n",
    "h = 1e-2\n",
    "\n",
    "pos_i=[]\n",
    "\n",
    "for i in range(ndim):\n",
    "    z[i,:] = result.x[i] + h*np.random.randn(nwalkers)\n",
    "\n",
    "for i in range(nwalkers):\n",
    "    pos_i.append(np.array([z[0,i],z[1,i],z[2,i]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b_steps, steps = 150,500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampler = emcee.EnsembleSampler(nwalkers,ndim,lnpost, \n",
    "                                args=(data[1],data[2],time,prior_list,'logistic'),\n",
    "                                threads = 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos,prob,state=sampler.run_mcmc(pos_i, b_steps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print sampler.acceptance_fraction.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampler.reset()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_,_,_=sampler.run_mcmc(pos, steps, rstate0=state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print sampler.acceptance_fraction.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Graficamos las cadenas para checar convergencia "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(sampler.flatchain[:,0]);\n",
    "plt.figure()\n",
    "plt.plot(sampler.flatchain[:,1]);\n",
    "plt.figure()\n",
    "plt.plot(sampler.flatchain[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chains_df = pd.DataFrame(sampler.flatchain,\n",
    "                         columns=['r','log(K)','log(A)'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = sb.PairGrid(chains_df)\n",
    "c.map_lower(sb.kdeplot, cmap=\"Blues_d\")\n",
    "c.map_diag(plt.hist,bins=50,histtype='step')\n",
    "for i,j in zip(*np.triu_indices_from(c.axes, 1)):\n",
    "    c.axes[i,j].set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "models_exp=[]\n",
    "for i in range(len(sampler.flatchain[:,0])):\n",
    "    models_exp.append(logistic([sampler.flatchain[i,0],\n",
    "                                sampler.flatchain[i,1],\n",
    "                               sampler.flatchain[i,2]],time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print np.shape(models_exp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_exp = np.mean(models_exp,axis=0)\n",
    "std_exp = np.std(models_exp,axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(data[0],data[1],yerr=data[2],fmt='o',color='orangered',label = 'data');\n",
    "plt.plot(data[0],mean_exp,'g',label='inference');\n",
    "plt.fill_between(data[0],mean_exp-std_exp,mean_exp+std_exp,alpha=0.2,color='g');\n",
    "plt.xlabel('Year',fontsize=14);\n",
    "plt.ylabel('Number of people',fontsize=14);\n",
    "plt.legend(loc='best');"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}