notebooks/T5 - 2 - Logistic Regression - Implementación.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementación el método de la máxima verosimilitud para la regresión logística"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definir la función de entorno L(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$L(\\beta)=\\sum_{i=1}^n P_i^{y_i}(1-Pi)^{y_i}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, Math, Latex\n",
    "display(Math(r'L(\\beta)=\\sum_{i=1}^n P_i^{y_i}(1-Pi)^{y_i}'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood(y, pi):\n",
    "    import numpy as np\n",
    "    total_sum = 1\n",
    "    sum_in = list(range(1, len(y)+1))\n",
    "    for i in range(len(y)):\n",
    "        sum_in[i] = np.where(y[i]==1, pi[i], 1-pi[i])\n",
    "        total_sum = total_sum * sum_in[i]\n",
    "    return total_sum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calcular las probabilidades para cada observación"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$P_i = P(x_i) = \\frac{1}{1+e^{-\\sum_{j=0}^k\\beta_j\\cdot x_{ij}}} $$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(Math(r'P_i = P(x_i) = \\frac{1}{1+e^{-\\sum_{j=0}^k\\beta_j\\cdot x_{ij}}} '))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logitprobs(X,beta):\n",
    "    import numpy as np\n",
    "    n_rows = np.shape(X)[0]\n",
    "    n_cols = np.shape(X)[1]\n",
    "    pi=list(range(1,n_rows+1))\n",
    "    expon=list(range(1,n_rows+1))\n",
    "    for i in range(n_rows):\n",
    "        expon[i] = 0\n",
    "        for j in range(n_cols):\n",
    "            ex=X[i][j] * beta[j]\n",
    "            expon[i] = ex + expon[i]\n",
    "        with np.errstate(divide=\"ignore\", invalid=\"ignore\"):\n",
    "            pi[i]=1/(1+np.exp(-expon[i]))\n",
    "    return pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calcular la matriz diagonal W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$W= diag(P_i \\cdot (1-P_i))_{i=1}^n$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(Math(r'W= diag(P_i \\cdot (1-P_i))_{i=1}^n'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def findW(pi):\n",
    "    import numpy as np\n",
    "    n = len(pi)\n",
    "    W = np.zeros(n*n).reshape(n,n)\n",
    "    for i in range(n):\n",
    "        print(i)\n",
    "        W[i,i]=pi[i]*(1-pi[i])\n",
    "        W[i,i].astype(float)\n",
    "    return W"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Obtener la solución de la función logística"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\beta_{n+1} = \\beta_n -\\frac{f(\\beta_n)}{f'(\\beta_n)}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$f(\\beta) = X(Y-P)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$f'(\\beta) = XWX^T$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(Math(r\"\\beta_{n+1} = \\beta_n -\\frac{f(\\beta_n)}{f'(\\beta_n)}\"))\n",
    "display(Math(r\"f(\\beta) = X(Y-P)\"))\n",
    "display(Math(r\"f'(\\beta) = XWX^T\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logistics(X, Y, limit):\n",
    "    import numpy as np\n",
    "    from numpy import linalg\n",
    "    nrow = np.shape(X)[0]\n",
    "    bias = np.ones(nrow).reshape(nrow,1)\n",
    "    X_new = np.append(X, bias, axis = 1)\n",
    "    ncol = np.shape(X_new)[1]\n",
    "    beta = np.zeros(ncol).reshape(ncol,1)\n",
    "    root_dif = np.array(range(1,ncol+1)).reshape(ncol,1)\n",
    "    iter_i = 10000\n",
    "    while(iter_i>limit):\n",
    "        print(\"Iter:i\"+str(iter_i) + \", limit:\" + str(limit))\n",
    "        pi = logitprobs(X_new, beta)\n",
    "        print(\"Pi:\"+str(pi))\n",
    "        W = findW(pi)\n",
    "        print(\"W:\"+str(W))\n",
    "        num = (np.transpose(np.matrix(X_new))*np.matrix(Y - np.transpose(pi)).transpose())\n",
    "        den = (np.matrix(np.transpose(X_new))*np.matrix(W)*np.matrix(X_new))\n",
    "        root_dif = np.array(linalg.inv(den)*num)\n",
    "        beta = beta + root_dif\n",
    "        print(\"Beta: \"+str(beta))\n",
    "        iter_i = np.sum(root_dif*root_dif)\n",
    "        ll = likelihood(Y, pi)\n",
    "    return beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comprobación experimental"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array(range(10)).reshape(10,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0],\n",
       "       [1],\n",
       "       [2],\n",
       "       [3],\n",
       "       [4],\n",
       "       [5],\n",
       "       [6],\n",
       "       [7],\n",
       "       [8],\n",
       "       [9]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = [0,0,0,0,1,0,1,0,1,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "bias = np.ones(10).reshape(10,1)\n",
    "X_new = np.append(X,bias,axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 1.],\n",
       "       [1., 1.],\n",
       "       [2., 1.],\n",
       "       [3., 1.],\n",
       "       [4., 1.],\n",
       "       [5., 1.],\n",
       "       [6., 1.],\n",
       "       [7., 1.],\n",
       "       [8., 1.],\n",
       "       [9., 1.]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter:i10000, limit:1e-05\n",
      "Pi:[array([0.5]), array([0.5]), array([0.5]), array([0.5]), array([0.5]), array([0.5]), array([0.5]), array([0.5]), array([0.5]), array([0.5])]\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "W:[[0.25 0.   0.   0.   0.   0.   0.   0.   0.   0.  ]\n",
      " [0.   0.25 0.   0.   0.   0.   0.   0.   0.   0.  ]\n",
      " [0.   0.   0.25 0.   0.   0.   0.   0.   0.   0.  ]\n",
      " [0.   0.   0.   0.25 0.   0.   0.   0.   0.   0.  ]\n",
      " [0.   0.   0.   0.   0.25 0.   0.   0.   0.   0.  ]\n",
      " [0.   0.   0.   0.   0.   0.25 0.   0.   0.   0.  ]\n",
      " [0.   0.   0.   0.   0.   0.   0.25 0.   0.   0.  ]\n",
      " [0.   0.   0.   0.   0.   0.   0.   0.25 0.   0.  ]\n",
      " [0.   0.   0.   0.   0.   0.   0.   0.   0.25 0.  ]\n",
      " [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.25]]\n",
      "Beta: [[ 0.43636364]\n",
      " [-2.36363636]]\n",
      "Iter:i5.777190082644626, limit:1e-05\n",
      "Pi:[array([0.08598797]), array([0.12705276]), array([0.18378532]), array([0.2583532]), array([0.35019508]), array([0.45467026]), array([0.56329497]), array([0.66616913]), array([0.75533524]), array([0.82687453])]\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "W:[[0.07859404 0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.11091035 0.         0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.15000827 0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.19160683 0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.22755849 0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.24794521\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.24599375 0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.22238782 0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.18480392 0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.14315304]]\n",
      "Beta: [[ 0.60426056]\n",
      " [-3.34641372]]\n",
      "Iter:i0.9940407075349087, limit:1e-05\n",
      "Pi:[array([0.0340128]), array([0.06053134]), array([0.10546805]), array([0.1774629]), array([0.28305225]), array([0.41943069]), array([0.56933774]), array([0.7075284]), array([0.81572841]), array([0.89011647])]\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "W:[[0.03285593 0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.0568673  0.         0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.09434454 0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.14596982 0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.20293367 0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.24350859\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.24519228 0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.20693196 0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.15031557 0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.09780914]]\n",
      "Beta: [[ 0.65761412]\n",
      " [-3.66759924]]\n",
      "Iter:i0.10600674406802064, limit:1e-05\n",
      "Pi:[array([0.02490177]), array([0.04697681]), array([0.0868775]), array([0.15515129]), array([0.26170168]), array([0.40624059]), array([0.56907679]), array([0.71823018]), array([0.83108181]), array([0.90473054])]\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "W:[[0.02428167 0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.04476999 0.         0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.0793298  0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.13107937 0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.19321391 0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.24120917\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.2452284  0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.20237559 0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.14038483 0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.08619319]]\n",
      "Beta: [[ 0.66217766]\n",
      " [-3.6953843 ]]\n",
      "Iter:i0.000792835124600839, limit:1e-05\n",
      "Pi:[array([0.02423594]), array([0.04594805]), array([0.08540873]), array([0.15331276]), array([0.25986436]), array([0.40504298]), array([0.56897776]), array([0.71907124]), array([0.83230289]), array([0.90586963])]\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "W:[[0.02364856 0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.04383683 0.         0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.07811408 0.         0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.12980796 0.         0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.19233487 0.\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.24098316\n",
      "  0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.24524207 0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.20200779 0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.13957479 0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.08526985]]\n",
      "Beta: [[ 0.66220827]\n",
      " [-3.69557172]]\n"
     ]
    }
   ],
   "source": [
    "a = logistics(X,Y,0.00001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "ll = likelihood(Y, logitprobs(X,a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.32622426e-06])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ll"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = 0.66220827 * X -3.69557172"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Con el paquete statsmodel de python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "logit_model = sm.Logit(Y,X_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.431012\n",
      "         Iterations 6\n"
     ]
    }
   ],
   "source": [
    "result = logit_model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                        Results: Logit\n",
      "==============================================================\n",
      "Model:              Logit            No. Iterations:   6.0000 \n",
      "Dependent Variable: y                Pseudo R-squared: 0.360  \n",
      "Date:               2018-04-05 18:54 AIC:              12.6202\n",
      "No. Observations:   10               BIC:              13.2254\n",
      "Df Model:           1                Log-Likelihood:   -4.3101\n",
      "Df Residuals:       8                LL-Null:          -6.7301\n",
      "Converged:          1.0000           Scale:            1.0000 \n",
      "----------------------------------------------------------------\n",
      "         Coef.    Std.Err.      z      P>|z|     [0.025   0.975]\n",
      "----------------------------------------------------------------\n",
      "x1       0.6622     0.4001    1.6551   0.0979   -0.1220   1.4464\n",
      "const   -3.6956     2.2889   -1.6145   0.1064   -8.1818   0.7906\n",
      "==============================================================\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(result.summary2())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}