2_evaluation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluation\n",
    "\n",
    "We can use the forward or backward algorithms to evaluate the probability of obtaining a sequence given a model —\n",
    "$P(x|\\theta)$\n",
    "where $x = x_1, x_2... x_N$ is a sequence and $\\theta$ is an HMM with $k$ states $\\pi_1,\\pi_2... \\pi_k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns; sns.set()\n",
    "import hmm\n",
    "\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's import the `transition`, `emission` and `init` matrices that define the HMM. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "transition = pd.read_csv(\"cpg_hmm_1/transition.csv\",index_col=0)\n",
    "emission = pd.read_csv(\"cpg_hmm_1/emission.csv\",index_col=0)\n",
    "init = pd.read_csv(\"cpg_hmm_1/init.csv\",index_col=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a DNA sequence of length 10 under the HMM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "st, seq = hmm.generate(transition, emission, init, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+++++++---\n",
      "GTGGCTGAGG\n"
     ]
    }
   ],
   "source": [
    "print(st)\n",
    "print(seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forward Algorithm\n",
    "\n",
    "The forward ($F$) matrix has to be initialised and filled up similar to the viterbi algorithm, except instead of calculating the maximum probability, we sum over the probabilities. The elements of $F$ are\n",
    "\n",
    "$$F_{i,j} = P(x_1,x_2...x_j,\\pi_i)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = np.shape(emission)[0]\n",
    "N = np.size(list(seq))\n",
    "F = np.zeros([k,N]) #forward matrix\n",
    "F[:,0] = init*emission[seq[0]] #initialize\n",
    "\n",
    "#fill up\n",
    "for j in range(1,N):\n",
    "    for i in range(0,k):\n",
    "        possibles = F[:,j-1]*(transition.iloc[:,i])*(emission[seq[j]][i])\n",
    "        F[i,j] = np.sum(possibles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk8AAACRCAYAAADabFwYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3cd3hThfv38XcntIwi0gUtQ2VZNsheLXtDgbK+spfKEBTZLZQlew8BRTYUGQVkb2hBAdlbdilNKXSvJM15/ggEaocNpOSHz/26rlza5E7O/eGcJHfOOYmFoigKQgghhBAiSyzN3YAQQgghxIdEhichhBBCCCPI8CSEEEIIYQQZnoQQQgghjCDDkxBCCCGEEWR4EkIIIYQwggxPQgghhBBGsH7fC6w3N+h9L/K9OD6sFldD4szdRrYo45YbgMiEFDN3kj0+srcCIEHz3/zJM3sbCwCStGZuJJvkfPkq9l/O91/OBpLvQ/X/S770yJ4nIYQQQggjyPAkhBBCCGEEGZ6EEEIIIYwgw5MQQgghhBFkeBJCCCGEMIIMT0IIIYQQRpDhSQghhBDCCDI8CSGEEEIYQYYnIYQQQggjyPAkhBBCCGEEGZ6EEEIIIYwgw5MQQgghhBFkeBJCCCGEMIIMT0IIIYQQRpDhSQghhBDCCDI8CSGEEEIYQYYnIYQQQggjyPAkhBBCCGEEGZ6EEEIIIYwgw5MQQgghhBFkeBJCCCGEMIK1uRswVvViH9G/VhFsrCy5FxHP9IN/k6BOSbd2VOPPuP88gc3nQwGwtbJkmNcnlHLJjYWFBTeexjL3yD3UKbr3GSGN82dOsm7lIrQaDUU++Yyvv/fFPlfuLNXEx8WyZNYknjx+gKLTUb9xS9p16QnAlQtnWbt8PlqtFtscOegzaATFS5XJ9jxBJ4+zZOFcNGo1nxUvwVi/yeTKnTtLNSkpKcyfPZ0/TgeRkqKl65e98O7Ymft3/8Z3zAjD/XU6HXf/vsO0WfPxbNCI7b9tZvOGtVhZW1OwYCHG+k0m30cfZUu+k8ePsXDeHNQaNcVLlMTPfwq5/5Evo5qUlBRmz/iR08GnSNGm8GXP3nTs1DnVfXds28rRwweZv3gZAIqisGThfPbv24udnR3lK1Tkux9GkSNHDpNlOnH8GAvmzUatVlOiREkmTJqaJlNGNSkpKcyaMY3gIH2m7r1649OpCwAPHz7Ab9wYoqOisLO3Z8q06RT75FMAhg8dzO1bN7G3twfgi6rVGDFqDElJSUydPJFrV6+i0+koW64cY8b5kTNnzg8mG4BarWbw1wPo0LETjZo0BSA2NhavujUpVuwTQ933I0dTtVr1t8r2po3r1xGweSMWFha4u7vjO3EyH3/8seH2XYE7WLt6leHv2LhYwlUqDhw+zscFCmR5OS9evGDc6B94GhqKhaUlvhP8qVCxEgCzZvzIwf37cHBwAKBIsWLMnD3vnbOB/nngO3Y0nxUvTo9efd66JjPmyrZ7VyCrf/kZCwsLctrZMXL0WDzKlE1Vs2H9Wlb+tIwCL9eVfa5c/Lp2g1HLyWx7jo6K4sepk7h79y7JyUn07T+QVq3bmiRfeu7cvsWPUycTGxuLlZUl4/38+dwj6+9PiYmJTPQdx82b19HpdHw7fAReDRpy5nQwc2ZON9QlJSfx8MEDNgZsNerx/80HtefJwc6aUY0/Y/zum3y5+i9Co5MYULtImroi+e2Y294DzxKpXxC+rOaGlaUFvddepPfaC+SwtuR/Vd3eV/vpio6KZNHMiYyYMJOFq7fh7OrGupULs1yz6delfOzoxLyfA5i+ZC37d/3GrWuX0Wg0zJk8moHDxzFnxSY6dOvDgmm+2Z4n8sULJvuNZdrMeQTs2ENBN3cWL5iT5ZrtWwN4/Ogh67cE8su6ADZvWMu1q5cp9ulnrN283XCpWr0mjZu2wLNBI0KfhLBs8Xx++mUd6wN24FqwECuWLcqWfC9evMBv/BhmzlvAjt37cHNzZ8Hc2Vmu2bplM48ePWTL9l2s27SFDetWc/XKZQCio6OYPNGP6dMmoyiK4fF27tjGiePHWL9pC5u37qCAoyOLF5jmRftVv77jRjN73kJ2/r6fQm7uzJ8zK8s1vwVs4tHDh2zdsZsNm39j/drVXLmszzT6h+/x6dSF7bv28PWgwQz/dogh2+VLF/hlzToCtgUSsC2QEaPGALDyp6WkpKSwZVsgv23fSXJSMj+v+OmDynbp4gX+18WHC3+dT7Wsy5cuUqnKF4bMAdsCTTI4Xb92lTW//sKa9ZvYFribwkWKsnjh/FQ1rdq0NSxz/ebfKFDAkVFjxxs1OAFMmzyRSpWrsH3XHqb+OJPvhw8lMTER0OeePmuOYTmmGi7u3b1Lv949OLB/7zvV/BtzZHtw/x5zZ81kyfKVBGwLpN+Arxg+dHCauksXLvD9D6MMyzd2cILMt+fxY0fh5OxCwNYdLF/5K9OnTUEVFvbO+dKTmJjIwH596Nm7LwFbd9B/4NeMHvm9UY+xbPFC7O3t2bFrLz+tWMXUSRNQhYVRvUbNVM+v4sVL0rtvf5MOTmDk8HT06FGTLtxYXxT5iJthcTyJSgIg8HIYDUs5pqlrW96VvdfDOXo7ItX1l0JiWPPHYxRAp8CdZ/E45zXdp/e3cencaT4r+TkF3QoD0KR1B04e3pvqzTOzmt7fjKDHwG8BiHwRgUajxj5XbmxsbFixeS+fFC+Foiionj4hT16HbM/zx5kgSnuUoXCRogB4d+zM/r27U+XJrOb4kUO0bNMOa2tr8uZ1oGGTZuz7fVeqZVz86xxHDx1g5Fg/QP9pSqvVEp8Qj06nIykpCVtb22zJdyY4CA+PshR52XvHTp3Z+/uuVPkyqzly6BBt2nrr8zk40KRpc37ftROAA/v24ejoyLDvfki1zBvXr+Hp1YA8efMC0KBhIw4d3G+yTKeDT1GmzOt+fTp3Yc8/MmVWc+TwIdq0e52pabMW/L57JyqVigf379G0eQsAatepR1JiIjdvXCck5DHx8fFMnuhHh3atGD92NNFRUQBUqvIF/QZ8haWlJVZWVpQqXZqnoaEfTDaADevWMmjIt5QtVz5VP5cuXiAmOooe/+uCT/u2BGwy/g0wPZ97lGHnnv3kyZOH5ORkwlUq8uXLl2H9qp9XkD9/fjr66Pd6atRqZv44lU4d2tGxXWvGjxlFXFxcmvtptVpOHD+GdwcfAEqVLk2RIkUJOnUStVrNzRvXWf3rL3Rs15rhQwe/9Xr7p00b19OmnTeNmzQzuub/ejYbW1v8/Cfj6OgE6NdlREQEGrU6Vd2lixfYs2c3Pu3bMrBfH+7cvmW4bdvWLXTu6I1P+7b079OT+/fuprusjLbn6KgozpwOZuDXgwBwdnFh3cYA8jpkz3vG6eAg3NzdqVO3HgD1PRsYhtEVPy2lU4d2+Hi34dvBXxMersowi3eHjgC4FixIjZq12b8v9eC8e1cgoaEhfDN4qMkzGDU8LViwwOQNGMMpjy3hca83qGexyeTOYY29rVWquvlH73HgxrM09z/3KIqQl4OXc54cdKhYkGP/GLDet4hnKgo4uhj+/tjRiYT4eBIT4rNUY2FhgZWVNfOnjmNYHx88ylehoLt+b5y1tQ1RL57Tv1Mz1iyfT5tO3bM9T3hYGM7Or3t1cnImPi6OhPj4LNWEq8JwdnZ94zaXNE+eBXNnMmDQUMOhQPfCRejWvTed2janRaN6XDh/lp59BmRLvrCwpzi7vNG7swtxcXHEv5EvsxqVKu1t4Sp9vo6dOjPg60HkzJl6oC9TtjzHjx0lMjISnU7H7p2BRDxLu32/daanYal6ck4vUyY1YWFPcXFxTXWbShWGKuwpjk5OWFq+fplxcnZGFRbGixcvqFajJuMn+LP5tx3Y29vjN16/56lmrdoULVoMgNDQJ6xfu9pw2OtDyAYwfdYc6tarn6YfKysr6tX34pfV61i45CfWrlnNkcOH3irbP9nY2HDk8CEae9Xl/PmztGnnnW5dZOQL1qxeZdjTB/DzyuVYWVuxacs2tmzfiaOTU5o9dABRL7fB/PnzG65zcnZGpQojPFxF1WrVGfrtcAK2BVKufHmGDv461aD6tsaM8/3XQ0gZ1fxfz1aokJthW1EUhVkzplHf0wubNz4AJiQkULTYJ/TtN4CArTto1749Xw/sR0J8POfO/smuwB2sWrOegK076Nm7b7p7roAMt+dHjx5RwNGRtatX0aNbZ7r4eHPj+nXs7OzeOV96Hj64T4ECjviNH0MXH28G9O1FSkoKuwJ3cOfObdZv2kLAtkBq16nHRN9xWcyiX1evaNRqFsybw4iRY7C2Nv0ZSkY9oik2lHdhiUW61+t0xvVVwikXk1uVZvvFp5y+H2mK1t6akkHvlpZWRtUMHTOZ/sPGMHPCCLasXUHnngMByJf/Y1YE7OPe7RtMGPEV7kU+MQxX2UGnpH/+mKWVZZZqdLq0t1m98QZ1+eIFoqOiaNKspeG6P04HcezwAQL3HSFfvo9YPH82/n5jmD1/ydvGyJCSTn//7DGzmvTyvflvk56WrdugUoUxoHcPctrb0b6DDzY2NkZ0nTklo/XxZqZMatJ7/mWUFcDSyopy5cozb8Fiw3VffTOIBvVqo1GrDW8a169dZdiQQXTu+j/q1ffMcp43mSNbZgZ89Y3h/52dneng04kjhw7i1aBhpvfLKq8GDfFq0JCtWwL4qn8fdu89mCorwNYtAXh6NsDNzd1w3Ynjx4iNjeVMcDAAGo2G/G+cL/VKRs9dK0sr3NzcWbxsheG6Hr36sHzZEp48CUm1rPftQ8mWkJCA79hRhIWFseSnlalus7e3Z9mKnw1/N2nanOXLlnL16hWCTp7g8aOH9Oj2+tzJ6OhooqOicPjH3seMtmetVsOTkBBy5c7N6vX6Q3u9unejSJEiJj/cBfq9fKdOHmfFqjWUK1eeo0cO8c3A/pSrUIFbN27Q1ac9ACk6HUlJiek+RrpZ3ngtPXhgP+7uhalUuYrJ+wcjhycvL69saSKrVLHJlHZ9faJngdw5iEnSkKTN+gnfXiUKMKzBJ8w/co9Dt8y71wmggJMLd25eNfz9POIZufPkJecbE39mNRfOBlOkWHHyF3DEzs6e2p5NOHPyCPFxsVy9eJZqtfXr7JMSpSn6aQke3v87W4cnZxdXrr08hwfgWbiKvHnzYmdnn6UaZxdXIiJe71V59kyF0xt7qQ4d2Euzlq1TvSGcPH6UOvU8yZ9f/4LYvlNXunVonS35XFwLcuWN3sPDVeTN64CdvX2WalxcC6baa/QsXJVqL1x6oqOjaNaiJX366femXbl8CffCpluHLq6uXLl8KU2/9qkyZVzj6urKszcyhb/M5OpakOcRESiKgoWFRarb/jp/jpjoaOp7NQAw1LwaPvbu+Z2pkyYyeux4mrds9UFly8yG9Wvx9GyAa8GCvAyOtc27fyp+9PAhERHPDG8Ubb3bM9nfj5iYaPLlS/3Fif179zByTOpP8zqdjpGjx1C7jv4wSkJ8PMnqZK5dvZLqk/+GgK0AxERHGw7phKtUOLs4c/vWTW7duplq74+iKNhYm27QfxsfQranoaEM+WYgxT79lJWr1qT5ckRo6BOOHT1C125fplq+tbU1KTodLVq1Ydh3Iwx5n4WHk9fBAR/vNoZ6P//JGW7Pjk76Q4Zt2ur3VhYuUoSKlSpx9crlbBmeHJ2cKFrsE8q9PKzt6dWQib7jsLSwpFefvvh07grov3QREx1NeLiKQQP7G+6/aNlyfZaIZxRw1J+6owoPp1SpUoaa/fv2GPJkB6MO2w0ZMiS7+siSsw+j+NwlD4Xy6Tes1uVcCLr7Isv3r1f8Y4Z4FuP7bdf/TwxOABWqVOf29SuEhjwC4MCu3/iiZr0s1wQfO0TAmuUoioJGrSb4+CHKVvwCSysrFs/05+bViwA8enCXJ48eUKJ09n7brlqNWly9cplHDx8AsP23zdSp75Xlmrr1vdgVuA2tVktsbAwH9++lbv0GhvteOH+WKlVTn2BbslRpgk6eIOHloc6jhw7gUTb1uSamUqNmLa5cusTDl73/tnkT9f/xoSKzmvqeXgRu36rPFxPD/r17DANERq5fvcp3Qwej0WjQarX8snI5zVu0zPQ+xmWqzeXLr/vdsnlTmp4yq6nv1YAd2/SZYmJi2Lf3dzwbNMTZxQU398Ls27sHgKBTJ7G0tKR4iRIkJMTz49TJhvOcfv3lZxo1boKVlRUH9+9j+rTJLFvx8zsNTubKlpkLf53n11X6PQjRUVFs3/YbTZo2f6eMABERzxg5YjiRkfrXwz27d/HZZ8XTDE4x0dE8evyI8hUqprq+Zq3abNqwHo1ajU6nY6LfeBbMnYNHmbKpTr61tramTt36bAnYDMDtWze5d+8uVb6ohoWlJdOnTSEk5DEAAZs2UKJEyVSHRM3h/3q26Kgoevf8Hw0aNWbGrLnpfqvUzs6OxQvmGU7uPnniOElJiZQtW44aNWuxb8/vPHsWDsCWzRvp16eHvs838nmUKZvh9uzm5k7pzz3YGbgdgOcREVy8eCFbBieA2rXrEvrkCdev6XcKnD93FiwsqPzFF2zb+pvhnLTFC+czdvQPODk5p8ri5ORMfa8GbN2iX1eqsDCCT52kbj39HmpFUTh//hzVqtfIlv4BLJT3fCyu3tygd7p/taIf0b92EWwsLXgSncTUfXcomC8nIxp+St/1l1LV/vOnCtb3rETuHNY8i0s21FwNjWXe0Xvv1BPA8WG1uBqS9iTErDj/xynWr1yEVqvBxdWNwaP8UT19wtLZk5i9fGOGNXnyOhAfF8tPc6fy6MFdLCygaq36dOoxEEtLS65dOs/qn+aRotViY2NDt76DKFuxqtH9lXHT7+2LTEj/JyH+KfjkcZYsnIdGq8HNzR3fSdMIDQlhqv941m7enmGNg0M+tFotC+fO5M8zwWg0Gtp18KFb996Gx65fozIBO35PtTdKURRWLF3EoQN7sbG1xcW1ICPH+KaqycxH9vq9HQmarD0VTp44zsJ5c9BqNLi5uzNp2nRCHj/G3288m7fuyLDmVb65s2Zw5rQ+X4eOPnT/x1eqd+7YxqED+1mw5PU3zBbNn8uRQwfRKTo8vRoyaOgwrP7lENEr9jb6PSNJ2swzLZg7W78+3AszZep0QkIeM9F3HAHbAjOsccinzzRn5nROnw5Gq9HQwaeT4WviDx8+wN9vPJGRkeTIYYvvhEmU/twDgNW//sL2rVvQ6XQUL14Sv4mTyOvgQKtmjYmNjcHJydnQX4WKlRgz3i/d3nO+3HGTUT5zZHulT88v6dylm+GcrRcvXjBpoi8P7t9Hq9XSuWs3uv0v83MRc1pnvu5eCdi0gU0bN2BtZYWjkxOjx/kSHRWVKufVK5cZNeI7du87mOq+SUlJzJk5nbNn/0SnS6FkqdL4TpiU5icdQP/GOsF3HE+ehGBhYcF3I0ZSs1ZtQH+C7i8rV6DTpeDs7MIE/ymv97JlkA2ylg9g/JhRhp8heLXn6FW29GrMmS2r+Vb8tJQlixZQvHjqwXvh0uUM/qo/i5Ytx8nJmaBTJ5k/dzYajYbcuXMzznciJV/uadm0YT0BmzdiaWlBrly5GT/Bn88+K55mWZltz09DQ5k62Z+QkMcoio5uX/YwfKHgXfJl5Py5s8yZNYPExERsbW35YdQYKlSsxLIlizi4fx9YWODq6oqf/xScnZ3T3D8hPp7JkyZw4/p1dLoU+g34ipat9HvaXrx4QSPPOpy/dM34xt6QM5Odwh/c8PR/1bsMT//XGTs8fWiMHZ4+NFkZnj5k7/IC/iHI6vD0Ifr/Yd2B5PtQZTY8fVC/8ySEEEIIYW4yPAkhhBBCGEGGJyGEEEIII8jwJIQQQghhBBmehBBCCCGMIMOTEEIIIYQRZHgSQgghhDCCDE9CCCGEEEaQ4UkIIYQQwggyPAkhhBBCGEGGJyGEEEIII8jwJIQQQghhBBmehBBCCCGMIMOTEEIIIYQRZHgSQgghhDCCDE9CCCGEEEaQ4UkIIYQQwggyPAkhhBBCGEGGJyGEEEIII8jwJIQQQghhBBmehBBCCCGMIMOTEEIIIYQRZHgSQgghhDCChaIoirmbEEIIIYT4UFi/7wXaVRz0vhf5XiReWMRtVYK528gWJZztAXgRn2LmTrJH/lxWAMSr/5ufI3LZWgCQqDFzI9nEzkb/3yStefvILjmt/9vZQPJ9qP5/yZceOWwnhBBCCGEEGZ6EEEIIIYwgw5MQQgghhBFkeBJCCCGEMIIMT0IIIYQQRpDhSQghhBDCCDI8CSGEEEIYQYYnIYQQQggjyPAkhBBCCGEEGZ6EEEIIIYwgw5MQQgghhBFkeBJCCCGEMIIMT0IIIYQQRpDhSQghhBDCCDI8CSGEEEIYQYYnIYQQQggjyPAkhBBCCGEEGZ6EEEIIIYwgw5MQQgghhBFkeBJCCCGEMIIMT0IIIYQQRpDhSQghhBDCCNbmbsBYTWt74D+4NTlsrbl65wkDJ24gNj4p3drlE//H9b+fMm/tYQAsLS2YO8qHOpU/A2D/qeuMnrv9vfWekbOnT7Lmp4VoNGqKflqcISP9sM+VO0s1yclJLJv7I3duXEOn6Cj5eVkGDhuFRq1hzNC+qR7j4b2/6fXVt7Tt9GW25gk6eZylC+ei0aj5tHgJxvpOJlfu3FmqSUlJYcGc6Zw5HUSKVkvX7r3w7tAZgPNn/2DB3BmkaFNwyJePb78fRfESpVAUheVLFnDsyEEASnuU5YfRvuS0s8uWfCdPHGPhvDloNGqKFy+Jr/8Ucv8jX0Y1KSkpzJn5I6eDTqFNSaF7z9508NHnO/vnGebNnolWqyVHzpz8MGosZcqWA2Dt6l8I3L4NKysrPvooP2P9JuLuXthkmU4cP8bCebNRa9QUL1GSCf5T02TKqCYlJYVZM6ZxOvgUKVp9po6dugDw8OEDJowfQ3RUFHb29kyeOp1in3yKoigsXjiP/fv2YmdnR/kKFfn+h9HkyJGD7t06k5SUaFjuwwf38W7vw8gx494624J5s1Gr1ZQoUZIJk9LPll7Nq2zBQS+z9eqNzxvZ/Ma9zjZlmj7bK2q1msFfD6BDx040atI01fI0ajW9evyPRo2b0KNXn7fKlZ4jhw8xbvQPBP/5V5rb9u/bw7Ili7G2ssLZxYUx4/0oWLCQUY+fmJjIRN9x3Lx5HZ1Ox7fDR+DVoCEAG9avZeVPyyhQoAAA9rly8evaDe8eCrhz+xY/Tp1MbGwsVlaWjPfz53OPMu8139PQUKZMnki4SkVKSgrDv/+BWrXrmCTf7l2BrP7lZywsLMhpZ8fI0WPxKFM23dorly/Tq3tXDh49wUcf5TdqOZltz9FRUfw4dRJ3794lOTmJvv0H0qp123fOBnD40EGWLl6ApYUlefLmZYL/FNwLp3792rh+HQGbN2JhYYG7uzu+Eyfz8ccfG7WcFy9eMG70DzwNDcXC0hLfCf5UqFiJXYE7WLt6laEuNi6WcJWKA4eP8/HL7dVoynuWs8I3b31x8xypqJ7HKB6tJyg5K3yjzPrlgLJs8/E0deXb+StH/7ipxCckK6NmbzNc33f8GuXoHzcV+0qDlNxVBivnrj5Qun6/8p16enVRFEW5FRZv9OXPm4+VL6pWU46eva7cCotXRvtNVYb+MDbLNeMmzVAGDB6m3AiNVa4/iVH6fjVYGT9lZprlzFq8UmnXsbNyLSTK6B5feR6n/dfL34/DlWrVqisXr/+tPI/TKv5TpiujxvpmuWb5L2uVHr36KKqoJOV+6HOlYeMmyokzfykPn0YqlStXUfYdPqk8j9Mq56/cVho2aqw8fZGgbA3cq7Rp114Ji0xQImI1yoCvBylzFizJUr/P47SGfHHJun+9PH4aoVSrXl25fvueEpesU6ZMm6GMHe+X5ZpfVq9TevXuq0TFq5XQ8EilceMmyplzF5XIuCSlWvXqyrmLV5W4ZJ2y58BhpVHjxkpcsk45fOyU0qRpMyXseYzhMTp17pqlfuOSdYZ8Cer0LyFhz5Vq1asrN+7cVxLUijL1R32/Wa1ZtWad0qtPXyUmQaOEPYtSGjdpovxx/pKSoFaUdt7tld+271QS1Ipy4PAxpWmz5kp8sk7ZsPk3pWWr1kpYRLSSoFaUeQsWKZOn/Jimtz37DylNmzZTVM9jMuz/lURN2ssTlb7vm3/fVxI1+r7HjffLcs2va9Ypvfv0VWITNYoqQp/tz/OXlESNPtvW7TuVRI2iHDyiz5ag1imJGkU5c/YvpVXrNkrZsmWVnbv3pulr3Hg/pVq1asqy5SvT7fvNS0bZ/nm5+fd9pUGDhkr58hXSva1ylSrK5Ws3lUSNopw6/afSrp13lh73zcvUadOV0WPGKYkaRbn38IlSq1Yt5cHjp0qiRlGGDB1m+PfI6iWzdffq8iImQalZs5Zy4PAxJVGjKHv2H1QaN2ny3vO1bNlKWb12vZKoUZS/Ll9TKlepokTHJ79zvhu37yo1a9ZSHoWqlESN/nlSt269dGtDw58rbdq2U0qUKKGEhj83Ol9m23O//gOUqT/OUBI1ivLg8VOlcpUqhuzvki8yNlEpV768cuvvB0qiRlGWr1yl9OnbL1XN+YtXlPr1PZXwFzFKokZRJk/9URk9drzR+QYNHqIsXLxUSdQoyoXL15VatWsrL2ISUtXEJKiVDh19lDXrN2Z5+0yPUYftxowZ83YTmok0rF6K89cecvfRMwCWbzlJ52ZfpKkb6FOXNTvPsPVg6k9fVlaW5LLLQQ5ba3LYWGNjY0WSWvNees/IhT/PULyUBwXdiwDQrG1Hjh/ci6IoWarxKF+JTt37YWlpiZWVFZ8UL8WzsKeplhEa8oiANSsZPm4y1tY22Zrnz9NBlPYog3vhogB4d+zM/r27U+XJrOb40UO0aN0Oa2tr8uZ1oFHjZuzfs4vHjx+SK3duvqhWA4CixT4hV67cXL18kfoNGrH8l3XY2NiSEB9P5IsXODjky5Z8p4OD8PAoS+Ei+t47durM3t93pcqXWc3Rw4do3dZbn8/BgcbNmrNn905sbGzZd+g4pUp/jqIoPAl5bMjwcYECjBnnZ9hb8rlHGZ4+De/t4RoAABhrSURBVDVhplN4eJSliKHfLulkyrjmyKFDtHkjU5OmLdizaycqlYoH9+/RtFkLAGrXqUdiYiI3b1znxvVreHo1JG/evAB4NWzMwYP7U/UVHR3FFH8/Jk2dTp48ed46W5kyr/v26dyFPelky6jmyOFDtGn3OlvTZi34ffcb2Zq/zpb0MhvAhnVrGTTkW8qWK5+mp107dxAXF0uduvXfKlN6EhMTGTNyBN//MCrd22/fvEnJkqUoXqIkAJWrfEFo6BOePAkBYMVPS+nUoR0+3m34dvDXhIer0n2cI4cP4d2hIwCuBQtSo2Zt9u/bC8ClixfYs2c3Pu3bMrBfH+7cvmWSbKeDg3Bzd6dO3XoA1PdswMzZ895rvps3bhAdE41P564AlC79Ob+u3YCFhcU757OxtcXPfzKOjk6A/vkdERGBRq1OVafT6RgzcgSDhw5L8xjbtm6hc0dvfNq3pX+fnty/dzfDfOltz9FRUZw5HczArwcB4OziwrqNAeR1cHjnfLqUFFAU4uJiAUhIiMfWNkeqms89yrBzz37y5MlDcnIy4SoV+fLpX/80ajUzf5xKpw7t6NiuNePHjCIuLi7NcrRaLSeOH8O7gw8ApUqXpkiRogSdOpmqbtXPK8ifPz8dX+7xf1tGDU83btx4p4W9KzeXjwhRRRn+fhIehUMeO/Lkypmqbtj0LWz8/Wya+6/deYbImATu7p/C/YNTufc4gj0nrmZ735l5Fh5GASdnw98FHJ1IiI8jMSE+SzWVqtag0MuhKjwslJ1b1lPLs1GqZaxdsZiW3p1xcnbN5jSgUoXh5Oxi+NvRyZn4uDgS4uOzVKNSheH8Rp9Ozi6Eq1QULlyUxMQE/jgdBMD1a1e4d+9vIiL0g7S1jQ1bNq2nbfMGREdFUs+rYfbkC3uKs8vr3p2cXYiLiyP+zXyZ1PzzNmdnF1Qq/Qu5jY0NzyMiaNqwHvNmz6RHL/1h18+Kl6DyF1UB/aGgBfNm06hxExNmCsPlHz2lzZRxjUr1FBcX11S3qVRhqMKe4ujkhKWl5Ru3OaNShVG2bDmOHztCZOQLdDodu3fuIOJZeKq+Vv28gtp16mV4+CIrwp6Gpfn3/me2zGrCwrKezcnZGVVYGADTZ82hbr36afq5c/sWG9atYfyESW+dKT2TJvrSwacTxUuWTPf2UqU/5+87t7n58jX82NEjREVFEfHsGbsCd3Dnzm3Wb9pCwLZAatepx0Tf9A+Rpv330K/PhIQEihb7hL79BhCwdQft2rfn64H9Uj3v39bDB/cpUMARv/Fj6OLjzYC+vUhJSXmv+R4+vE/BgoWYOX0a3Tp3pEe3zkQ8e4aNzbt/GC1UyM2wrSiKwqwZ06jv6YWNrW2quiUL51OmbLk0hwrPnf2TXYE7WLVmPQFbd9Czd1+GDx2cxXz67fnRo0cUcHRk7epV9OjWmS4+3ty4fh07E5z6YJ8rF+N8J9K9W2ca1q/Npo3r+Xb492nqbGxsOHL4EI296nL+/FnatPMG4OeVy7GytmLTlm1s2b4TRycn5s+Zleb+UZGR6HQ68ud/fSjT6eX6eyUy8gVrVq9ixKh33xGUpXOeRo8eDUBoaKjh/6dNm/bOCzdWRlN+SoouS/cfO6A5EZFxFGkwGrucNgTM6c/QL72Yv/aIKds0ypufgN9kaWllVM3ft64zZex3tPDuTNWadQ3XP1OFceFsMINH+pqo48zpdOmvC0sryyzVKOncZmllSa7cuZk+ZxE/LZ7PonmzqFCpMpWrVEv14tWxczc6dOrK8iULGDPiW5auXPOOadLSKen3bvXGm2hmNend9uZ9Py5QgP2HT3Dj+jUG9u3FJ59+SpGixQCIfPGCEcOHkDtPHgal8+nzbWW0PlJlyqRGp0u7fVpaWWa8ni2taNm6LSqViv69e2Bnb493B59U6zI5OZmtvwWwcfM2Y6KkoWSwLt4cejKrSS+bPnNG27BVutcDxMbGMnbUD0ybMRt7e/vM2jbK5o3rsbKypp13B8Oeln9yL1yYiZOmMtnfD7VajadXA0qULIWNjQ0njh/l6pUrdPVpD0CKTpfqnLM3pfvvYWWJvb09y1b8bLiuSdPmLF+2lKtXr1C1WvV3yqfVajl18jgrVq2hXLnyHD1yiG8G9mffoaPYvhwwsjufVqvl4oW/6NGrNyNGjubK5ct8M7Afv+3YidMbH2zfRUJCAr5jRxEWFsaSn1amuu3E8WNcuXKZpct/TnO/k8eP8fjRQ3p0e70nJTo6muioKBzypd4Dn9H2rNVqeBISQq7cuVm9fhOPHj6kV/duFClSJM25Zca6c/sWPy1dzPade3AvXJj169bw3beDCdgWmOY93atBQ7waNGTrlgC+6t+H3XsPcuL4MWJjYzkTHAyARqMhfzrnQmX8uvv6Obl1SwCeng1wc3N/p0yQxeGpXbt2APz111+G/zeHx2GRfFG2qOHvQk4OvIiOJyFJnfGd3tDGqzzDp29Bo01BE5fCul1/0K5hRbMOT47OLty+fsXw9/OIcHLnyZvqZOd/qzlxeB9L50xjwLejqN+oWarHDzp+iOp1vLC3z5XNSfRcXFy5fvWy4e9n4Sry5M2LnZ19lmqcXVwNe5Ne3ebk5IJOp8Pe3p4lK1Ybbuvs3RI398LcuX0TnU5HyVKfY2FhQet2HQjYuDab8hXk6uXXvYeHq8ib1wE7e/ss1bi4FEyVLzxchZOzC7GxsZz98wxeDfR7DUt/7kGJkiX5+85tihQtxu1btxg25Cs8vRox7PsfsMrkTdpYrq6uXL1yKdNMmdW4uroS8Sx1JmdnF1xdCxIREYGiKIYXyVe3RUdH0axFS/r0GwDAlcuXcC9cxPAYQSdPULJkKdzc3+1FzsXVlSuX0/b95vCSWY2rqyvPMsj2PINsGQkOOklsbCyjf/gOgKdPn3L6dBBxcXF8M3joW2cM3LGdpKQkfLzboNFoSE7W//+iZcsNb+xqtRr3wkVYtzEA0A8k69euplAhN1JSdPTq09dwSEqtVhMTHU14uIpBA/sblrNo2XL9v0fEMwo4OgKgCg+nVKlShIY+4djRI3Tt9vrLKIqiYG397t9JcnRyomixTyj38hCop1dDJvqOI+TxYz759NP3ks/R0Yk8efLi+XKPdtly5XBzc+P2rZsmGZ6ehoYy5JuBFPv0U1auWkPOnKmPpuzYthWVKozOHV6///bt1QP/yVNJ0elo0aoNw74bAeg/6DwLDyevgwM+3m0M9X7+kzPcnh2d9IcM27TV7+0pXKQIFStV4uqVy+88PAUHnaJCxUqGE8Q7d+nGrOnTiIqKNJzw/ujhQyIinlGpchUA2nq3Z7K/HzEx0eh0OkaOHkPtOvrDtgnx8SSrk7l29UqqPYgbArYCEBMdbTjcGK5S4ezyev3s37vnrb948k9ZOmxXtWpVqlatip2dneH/zeHw6RtULVuUTwvrN+y+Heqw+9iVf7nXaxdvPqZ940oAWFtb0rJeWf68cj9bes2qil/U4Nb1K4Q+fgjA3sDfqFa7fpZrgo4dZPn8GfjPXpJmcAK4evE85Su/v/VVtUYtrl65zONHDwDYvnUzdet5Zbmmbn0vdgduQ6vVEhsbw8H9e6nr2QALCwuGDx7Ijev6w6yHD+7D2tqaz4qX5O/bt5kyYSxJifpPk3t2B1L5i2rZkq9GzVpcuXyJRw/1vW8N2EQ9T68s19T39CJw+1Z9vpgY9u/dg6dXA6ysLJk4fiwXL+jP07v79x0e3L9PmbLlefToIQP6dKffwG/4fuRokw5O+n5rc/nSJR6+7Pe3zZuo79UgyzX1PRuw42WmmJgY9u/9HU+vhji7uODuXpj9e/cA+uHB0sKS4iVKcO3qVYYPHYRGo0Gr1fLzyp9o3qKVYXnnzv1Jteo1TJPt8uu+t2SULYOa+l4N2LHtdbZ9e3/Hs4E+m5t7Yfa9zBZ06iSWlvpsGWnStDl7Dx4hYFsgAdsCqe/pxZfde77T4ASwYfNvbAvcTcC2QBYtW06OHDkJ2BaY6k1drVbT48suhD3Vnw+5bs2vVKhUGYd8+ahZqzbbtv5mOI9k8cL5jB39A05OzoZeXz1efa8GbN2yGdAfyg0+dZK69Tyxs7Nj8YJ5XHn5oeHkieMkJSVS9uW3Rd9F7dp1CX3yhOvX9M/98+fOgoUFhdzc3lu+ChUrkSOHLceO6j9o3793l8ePH1O8RKl3zhcdFUXvnv+jQaPGzJg1N83gBDBn/kJ27Npr6BVg5arVeJQpS42atdi353eevTzsvWXzRvr16QGQKp9HmbIZbs9ubu6U/tyDnYH6b58/j4jg4sUL7zw4gf6Q6vlzZ3keEQHA0cOHKFTILdU3BSMinjFyxHAiI18AsGf3Lj77rDj58n1EzVq12bRhPRq1Gp1Ox0S/8SyYOwePMmVT5bO2tqZO3fpsCdCvv9u3bnLv3l2qvHwviImO5tHjR5SvUPGdMwFYKBkdE0rHoUOHaNjw3c4lsas46J3u36T25/gPbo2ttTX3QiLoO34NxdwKsMS3K9U7/5iq9p8/VZDfIRdzRnakQik3UnQKx/68xcg529Bqs3bYLzOJFxZxW5XwVvc9d/okq5cvRKvR4lLIjeFjJxEWGsLCGf4s+GVzhjV58jrQv0tr4uNi+fjlyYYApctU4Kvh+sOrX3dvzzffj8Oj3NtvMCWc9Z/SX8Sn/EulXvCp4yxdOA+NRkMhN3d8J00j9EkI0/zHs2bT9gxrHBzyodVqWTh3Jmf/CEaj0dC2vQ/duvcG4K/zZ5k3axpajYaPCzgyatxECr3c/bpi6UKOHNqPlZU1n3z6GcN/GEu+jz7KUr/5c+mHkXh11p4Kp04cZ+H8OWg0Gtzc3Zk0dTpPQh7j7zeeTb/tyLDmVb65s2fwx2l9vvYdfejeU/9V9fNn/2Tu7BlotVpsbW0ZNHQ4VatVZ9KE8ezZvdNw+A7A1taWNRsCstRvLlv9npHETL4bcfLEcRbOm/2y38JMnjadkMePmeg3joCtgRnWvMo0Z9Z0zrzM1KFjJ8PX7x8+fIC/33iioiLJYWvL+AmTKP25BwAL58/lyKED6BQdnl4NGTx0uGEwHPRVf5o1b0GLVm3Sb/gNdi+P9iVpM862YO5sNFp931OmTick5DETfccZ3ojSq3HI9zLbzOmcPh2MVqOhg0/abJGRkeTIYYvvG9le6dPzSzp36ZbmpwoAxo8ZxWfFi//rTxXktM442z89eRJC+zatOHPugmHPyqs9UAf372PpkkXodCkU++RTfCf489FH+dHpdCxbsoiD+/eBhQWurq74+U/B2TntHpWE+HgmT5rAjevX0elS6DfgK1q+XEdBp04yf65++8idOzfjfCdSslTmw0XOlzum/i3f+XNnmTNrBomJidja2vLDqDG4ubu/13yvfi4hMjISgEFDvjX8jMG75Fvx01KWLFpA8eKpB++FS5cz+Kv+qfYgvlLeoyTHTp02DCCbNqwnYPNGLC0tyJUrN+Mn+PPZZ8XTLCuz7flpaChTJ/sTEvIYRdHR7cse/3pSdVbX36YN69m0cR021jbkdXBg9DhfkpOSUj0HAzZtYNPGDVhbWeHo5MTocb64ubmTlJTEnJnTOXv2T3S6FEqWKo3vhElpfm4E9EPfBN9xPHkSgoWFBd+NGEnNWrUBuHrlMqNGfMfufQczbzadfOkxangyhXcdnv6vepfh6f86Y4enD42xw9OHJivD04fs34anD50xw9OHJqtvvh8qyfdhy2x4kl8YF0IIIYQwggxPQgghhBBGkOFJCCGEEMIIMjwJIYQQQhhBhichhBBCCCPI8CSEEEIIYQQZnoQQQgghjCDDkxBCCCGEEWR4EkIIIYQwggxPQgghhBBGkOFJCCGEEMIIMjwJIYQQQhhBhichhBBCCCPI8CSEEEIIYQQZnoQQQgghjCDDkxBCCCGEEWR4EkIIIYQwggxPQgghhBBGkOFJCCGEEMIIMjwJIYQQQhhBhichhBBCCCPI8CSEEEIIYQQLRVEUczchhBBCCPGhkD1PQgghhBBGkOFJCCGEEMIIMjwJIYQQQhhBhichhBBCCCPI8CSEEEIIYQQZnoQQQgghjCDDkxBCCCGEEWR4EkIIIYQwggxPQgghhBBGkOFJCCGEEMII/8nhadeuXTRv3pxGjRqxfv16c7eTLeLi4mjZsiUhISHmbsWkFi1aRIsWLWjRogUzZswwdzsmN3/+fJo3b06LFi1YtWqVudvJNtOnT2fUqFHmbsPkunfvTosWLWjTpg1t2rTh0qVL5m7JpI4cOYK3tzdNmzZl8uTJ5m7HpLZs2WJYb23atKFy5cr4+/ubuy2TCQwMNLx2Tp8+3dztmNzy5ctp0qQJrVq1YunSpeZuB5T/mLCwMMXT01OJjIxU4uPjlVatWil37twxd1smdfHiRaVly5aKh4eH8vjxY3O3YzJBQUFKp06dlOTkZEWtVivdu3dXDhw4YO62TOaPP/5QOnfurGg0GiUxMVHx9PRU7t69a+62TC44OFipVq2aMnLkSHO3YlI6nU6pVauWotFozN1Ktnj06JFSu3Zt5enTp4parVa6dOmiHDt2zNxtZYvbt28rjRo1Up4/f27uVkwiISFB+eKLL5Tnz58rGo1G6dChgxIUFGTutkwmKChIadmypRIbG6totVplwIAByv79+83a039uz1NwcDDVq1cnX7582Nvb06RJE/bt22futkwqICAAPz8/nJyczN2KSTk6OjJq1ChsbW2xsbHh008/JTQ01NxtmUzVqlVZs2YN1tbWPH/+nJSUFOzt7c3dlklFRUUxd+5cBg4caO5WTO7evXtYWFjQr18/Wrduzbp168zdkkkdPHiQ5s2b4+Ligo2NDXPnzqV8+fLmbitbTJgwgWHDhpE/f35zt2ISKSkp6HQ6EhMT0Wq1aLVacuTIYe62TOb69evUrl2b3LlzY2VlRZ06dTh06JBZe/rPDU/h4eE4Ojoa/nZyckKlUpmxI9ObMmUKVapUMXcbJle8eHEqVKgAwIMHD9izZw/16tUzc1emZWNjw4IFC2jRogU1atTA2dnZ3C2ZlK+vL8OGDSNv3rzmbsXkYmJiqFGjBosXL+bXX39l06ZNBAUFmbstk3n48CEpKSn06dOH1q1bs2HDBhwcHMzdlskFBweTlJREs2bNzN2KyeTOnZuhQ4fSrFkz6tatS6FChahUqZK52zIZDw8PTp06RVRUFMnJyRw5coSIiAiz9vSfG54URUlznYWFhRk6EW/rzp079O7dm5EjR1K0aFFzt2NyQ4YM4fTp0zx9+pSAgABzt2MyW7ZswdXVlRo1api7lWxRsWJFZsyYgb29Pfnz56dDhw4cP37c3G2ZTEpKCqdPn2bmzJkEBARw5coVtm/fbu62TG7Tpk306tXL3G2Y1M2bN9m6dStHjx7l1KlTWFpa8vPPP5u7LZOpUaMG3t7efPnll/Tt25fKlStjY2Nj1p7+c8OTs7Nzqok0PDz8P3d467/s/Pnz9OzZk++++4527dqZux2Tunv3Ljdu3ADAzs6Oxo0bc+vWLTN3ZTp79uwhKCiINm3asGDBAo4cOcLUqVPN3ZbJnDt3jtOnTxv+VhQFa2trM3ZkWgUKFKBGjRrkz5+fnDlz0qBBAy5fvmzutkxKrVZz9uxZvLy8zN2KSZ06dYoaNWrw8ccfY2tri7e3N3/++ae52zKZuLg4GjVqxK5du1i7di12dna4u7ubtaf/3PBUs2ZNTp8+zYsXL0hMTOTAgQPUrVvX3G2JLHj69CnffPMNs2bNokWLFuZux+RCQkIYN24carUatVrN4cOHqVy5srnbMplVq1axe/duAgMDGTJkCF5eXowZM8bcbZlMbGwsM2bMIDk5mbi4OLZv306jRo3M3ZbJeHp6curUKWJiYkhJSeHkyZN4eHiYuy2TunXrFkWLFv3PnWtYqlQpgoODSUhIQFEUjhw5QtmyZc3dlsmEhITwzTffoNVqiY2NZcuWLWY/7Prf+dj0krOzM8OGDaN79+5oNBo6dOhAuXLlzN2WyIKff/6Z5ORkfvzxR8N1nTt3pkuXLmbsynTq1avHpUuXaNu2LVZWVjRu3Pg/OST+V3l6ehrWn06no2vXrlSsWNHcbZlM+fLl6du3L127dkWj0VCrVi3at29v7rZM6vHjx7i4uJi7DZOrXbs2169fx9vbGxsbG8qWLUv//v3N3ZbJlCpVisaNG9O6dWtSUlLo2bOn2T94WijpnSQkhBBCCCHS9Z87bCeEEEIIkZ1keBJCCCGEMIIMT0IIIYQQRpDhSQghhBDCCDI8CSGEEEIYQYYnIYQQQggjyPAkhBBCCGEEGZ6EEEIIIYzw/wDf+o8l4+HTpQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "sns.heatmap(F,square=True,cbar=False,linewidth = 1,cmap='Blues',yticklabels=list(transition.columns),annot=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P = 1.4366902641444733e-06\n"
     ]
    }
   ],
   "source": [
    "P = np.sum(F[:,N-1])\n",
    "print(\"P =\",P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "P is the probability $P(x|\\theta)$ which can also be written as $P(x)$ since we can take for granted that all probabilities are with respect to the HMM we have been given."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Backward Algorithm\n",
    "\n",
    "We can do the same thing by filling up the backward ($B$) matrix with the backward algorithm. But now,\n",
    "\n",
    "$$B_{i,j} = P(x_{j+1},x_{j+2}...|\\pi_i)$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = np.zeros([k,N]) #forward matrix\n",
    "B[:,-1] = 1 #initialize\n",
    "\n",
    "#fill up\n",
    "for j in range(N-2,-1,-1):\n",
    "    for i in range(0,k):\n",
    "        possibles = B[:,j+1]*(transition.T.iloc[:,i])*(emission[seq[j+1]])\n",
    "        B[i,j] = np.sum(possibles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk8AAACRCAYAAADabFwYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3cd1hTh/fH8TdTwYELAgrOulfV1lZrHTgrKhUVtbZat7XW9W3rxomr7lVH60YRBME9cSG2jrp3taIgwwEOZGT9/kCjUUAigai/83oen8ckh9zzyU1uTu69iZlWq9UihBBCCCEyxdzUDQghhBBCvE9keBJCCCGEMIAMT0IIIYQQBpDhSQghhBDCADI8CSGEEEIYQIYnIYQQQggDyPAkhBBCCGEAy5xeYJIqp5eYM3JbftjZQPK9ryTf+y23JSSkfJg/x5fH2gyA+wkf5sornCf1yXnsxkMTd5I9ape2A8CmxgATd5I9Ek8tSPc22fMkhBBCCGEAGZ6EEEIIIQwgw5MQQgghhAFkeBJCCCGEMIAMT0IIIYQQBpDhSQghhBDCADI8CSGEEEIYQIYnIYQQQggDyPAkhBBCCGEAGZ6EEEIIIQwgw5MQQgghhAFkeBJCCCGEMIAMT0IIIYQQBpDhSQghhBDCADI8CSGEEEIYQIYnIYQQQggDyPAkhBBCCGEAGZ6EEEIIIQwgw5MQQgghhAFkeBJCCCGEMIAMT0IIIYQQBrA0dQOvWu+zFr8N6zEzM8PFxQWv8ZMoXLiwXs06nzX4rltL7ly5KVW6DCNHe2FXoIBBy3nw4AGjR/xK1J07mJmb4zVuAh/XqAnAtatXmDp5Eo8fP8bCwpwxYydQqXIVo+TTarV4jRrBR2XL0q17z9dunzF9Knt27cTOzg6AEqVK8dvMOQYtIzExkfFeo7l8+SIajYbBQ3/BtXETIPWx+2PJYooUKQKAbZ48rFyzLoupMt/7PydP8Nu0yahVaqxz5WLEqDFUrlLVoOWYIt+hgweYN2cmKSkplCtXnnETJ5M3b95M1ajVamZMn0LYkVDUKjVdu/fAs2NnAMLDbzJ29EgexsdjY2uL95RplCpdBoBVK5cTFBiApYUFBQsVYszYCbgUL663zCGDBmBv78DI0V7vVT6tVsvCeXPYtXMHNjY2VP+4Bj8PG0GuXLl4mpDA2DEjuX79OlqNhq892qX5WnnXsh77+y9mz5yOSqkiV+7cDBsxmqrVqgHg7+fLujWrsbC0oFgxZ8ZN9KZgwUJGy5Sew4cOMH/OLJTKFMqWLY/XBO/Xsr6pJjo6im5dOuK7MZiCBQsC8PBhPNMnT+LGjeskJyXRo08/WrV2z/Y8rzpy+CCL589BqUyhTNlyjPSaSJ5X8mVU09K1HvYODrrab7r2oHnLVty9G4v32FE8uH8PjVbLt9160sKtdY5me9XpY6H4rViEUpmCS6mP6D14NDZ5XskasoNtG9dgZmaGda7cfNfvf5QuV0mvZu2SWcTcuc3/xs/OyfaNZun4b7n4bxRz1uwzaR/v1J6nixfOs3rlclb7+BIYvJXiJUqycP5cvZpjf//Fij+XsezPVfgFBlOvfn0mjDP8jWPKpPHUrPUJm7ZsZ/LU3/h56CASExNJTEykX++efN+jF34BQfTp158Rw342Sr4b16/Tu0c3du/akW7NmdOnmDZjFn6BwfgFBhs8OAEsXjgfW1tbgrbsYMmyFUyeOI6Y6OjU+z91ip9/Ha67f2MNTpntfdTwXxk89Bf8AoPp3rM3o0cON3g5OZ3vwYMHeI0ewcw589m8bRfFnF2YO2tGpms2+vlyKzycgKCtrNuwEZ81qzh39iwAI379Gc+Ondm0ZTv9B/zE0MED0Wq1/HU0jKCAjaxZtwH/TZtp3KQpXqNH6C1zxZ/LOHXyxHuZLzgokEMHD+CzYSN+gcEUsbdnwbzU58uqlcvJlSs3gcFbWbPeD5+1qzl/7myWc2ZnVmVKCr/+PASv8ZPw37SZ3n1/YNSIXwCIiLjNgrmzWbHGh42btlC0WDF+XzDfKHkyEvfgAePGjGTG7Hls2rKTYs4uzJ8z06CarZuD6NmtC3djY/X+buzoETgoHFnvv4nfl63gtyneutdgTomLe4D3uNFMnjEH303bKFrMmUXzZ2W6Jvzmf+TLn59VvoG6f81btgJgyYI5VK5ajdUbNjF7wRJmTJnA/Xt3czTfyx7Fx7F01kQGjp7Kb39sxMGxGBtWLNSriYoIZ/0f8/h10jy8F/rg3qkH8yYN06v5+9AewkJ25mTrRlO+lIIdS36iXdOapm4FMHB42r9/f3b1AUClylXYvH0X+fLlIzk5mdiYGAq8skfp0sULfP55XRSOjgA0btKMgwdCUKakoExJ4bepk+nYvi0d2rZhzMjhPHny5LXlqFQqDh08gEd7TwAqVKxIiRIlORJ6mKNhR3B2ceHL+g0AaNio8VsNMGnxXe+De1sPmjX/Ks3bU1JSuHzpIqtWLqdD2zYMHfQTUXfu6G5ftuR3OrZvi6eHO4N/6k9sbEya9xOyby8e7TsA4FS0KHXq1mPXztSB7czpU2zfvhXPdl/Tr3dPrl29YpRsb+r9ObVGzeNHjwB4mpBArly53vl8R8NCqVKlKiVKlATAs1Nntm/bglarzVRNyL69uLf1wNLSkvx2drT4yo1tWzcTExPDzf9u0KKlGwD1vmxAUmIily9dpHCRIozyGqfbA1CpclW9x/PY339xJPQw7T07vZf5Ll64QKPGTcifPz8AjZs2Y+/uXQCo1WoSniagUqlITk5Go9FgZWWd5ZzZmdXK2po9IYeoWLESWq2WiIjbFCiQupdGo9agUqlISEhAo9GQlJiE9UvP++xyNOwIlStXpfizHB06dmLHa1nTr7kbG8P+kH3MX7RU734fPozn76Nh9PnhRwAUjo6sXreB/M/2OOeUY0fDqFi5Ci7FSwDg0aETu3ds08uXUc25M6cxN7dgQJ/v+c6zLcuXLkKtVgOgVmt48uQJWq2WpKQkLCwsMTM33b6Gc//8TelylXAslrrnuXGrdoTt36mX1dLKil6DR1GgUOpe91LlKhIfdx+VUglA5K3/2LZxDV9/Y7y9uDmpn2d9Vm/+i4A9/5i6FcDA4WnevHnZ1YeOlZUVIfv20sy1PidPHse9rYfe7VWqVuPY339x504kAMGbAlEqlcQ/jOfPP5ZiYWmBr38g/ps2Y+/g8NqnSoD4uDg0Gg2FCr3Ybe6gUBATE034zf8oUsSesWNG0tnTg769uuteUFk1crQXrdt8ne7tsbEx1P7scwYNHopfYDDVqldn0E/90Wq1bAkO4tq1q/j4+qfucfuyAeO9Rqd5P9HRUTg6OukuK55le/r0KSVLlaZX7774BQTRtl07+vfrzdOEhCxny6j3l42fOJlRI4fR1LU+3hPHM2LUGIB3Ol90VLRuWE9dniNPnjwh4aX7zajm9X4diYmJJiY6CnsHB8xf2ig7KBTEREdTtmw5Pvm0NpA6mM6dPYOmzVsAqY/19KneTJk+AwsLiyxlM1W+qtWqcWB/CHFxD9BoNGwJDuLu3dS9G9179uJOZCRNGn1JiyYNafGVG+UrVMhyzuzMCqnbrvv37tHUtT6zZ0zn+x69ACheogTduvfE3a0FjRvW48SJ4/Tq09coeTISEx2ll8MhjawZ1dg7KJg5Zz6ly3ykd7+3b92iSBF7fFavpPt3nenSsR2XL13ExsYm2zO9LCYmCoXiRe/2DgoSnjzRe71nVKNWq/j08zrMWrCURX+u4u+jR9jo6wPADz8NJvTgftybN+Kb9m3o2e9HChXSP30kJz24F0Nh+xeHFwsVcSDxaQJJT19ktVcU5ePa9YDU00PWLZ1Dzc/qY2llRVLiU5bMGEvvoWPJbZsnx/s3hiHT/Fm/7bip29Ax6JynV98Is4tr4ya4Nm5CgL8fP/TpydYde3Qb4FqffErf/j8yZOAAzM3N+LptO+zsCmBlZcWhgwd4/Pgxf4WFAaBUKilU+PUnvEarSXO5FuYWqFQqQg8fZNmK1VSrVp39IXv5sV8fdu7dj7W1cT79psfZ2YWFi5fpLnfr3pOlixcRGRnBoYP7OX/uHN94tgNArdGQlJSY5v1oNK+vJwsLc2xtbVm87E/ddc1btGTp4t85f/4ctT/7PNt6d3Z2AeD+vXtMGDuG5SvXULlKVUL27eV/QwayeduudzqfNp3ny8tDQUY1afZrbo5Gk87fvDQQPXjwgJ+HDCRv3rwMHDQEpVLJsJ+H8suwkdi/tDHNClPka93ma2JjYujdoxs2Nra06+CJlZUVAJMnTaBO3S8YOHgo9+/fo0/P7nz8cQ2aNGtuaLTXZFfW5woXKcLe/Ye5dPECvXt+T+kyZbgTGcnePbvZte8gBQsWZM6s3xgzcgTzFy3OWpg3SH87Z25QzatUKhWRkRHkyZuHFWvWc+tWOL26fYtL8RJGOzc0M7RprAsAcwvzTNW4e3TQXba2tqbTt93wX+9Dxy5dGTdqGF269cCjQydu3wrnx97fU6VqNSpVqWbcEJmkTee1ZJbGh6ekpESWzhzPg7ux/DIp9bSXP+ZMomlrT1xKluG/a5eytdf/Lwza8+Tq6ppdfQBwKzycf146h+Nrj3ZE3bnDo0cPddclJDzhk09qs2HjJtb7BdKkaeoG1c6uABqNhmEjRurOd/Hx9WfG7LlcOH8OTw933b/nnyAePXxxv7ExMSgcFdg7OFCyVGmqVasOQCPXJmg0aiJu387W7ABXr1xmy+Ygveu0Wi1Wllao1Rq69+yly7beL4BVa9YTGxujly02NgYnJyfuvnR8PiY2FoXCkTt3Ilnns+a1+7e0zPr3BjLq/bl/Tp7AqWhR3Qniro2bYGVpxX83rr/T+RydnLh398XyYmNjyJ/fDltb20zVODk5cfeV2xQKR5ycinL/3j29DyXPb4PUx7RLx/ZUrFiJ2fMWYmVtzcUL54mMjGDm9Kl4erjjv8GX3Tu3M85r1HuV72F8PF+5tWLjpi2sWbeB0mU+0h1a2bdnD+09O2Jubo69vQPNmrfg2LG/3zpfTmR9/Pgx+/bu0V1fsVJlypevwL9Xr3JgfwgNG7lSuHBhzM3N6di5C8eNlCcjjo5F08xh83LWTNS86vnQ3to99ahA8eIl+LhGTS6cP2fsCBlSODpx76XtwN3YWPLlz4+NjW2manZs3cy/Lx3Wf76tiI+L4+zpf2jTtj0ALsVLUPvzOpz+52QOpEpbYQdH4h/c112Ou3eXPHnzkzu3/t6+e7HRTBjaE3NzC0ZOW0SevPl4cDeGK+dPszNoPaN+7ELgmiVcOX+a38YMzukYHxSDhqeBAwdmVx8A3Lt3l2G/DCUu7gEA27du4aOPyurOHYDUJ3/P7t/pzmVasngRLVq6YWZmRt0v6uG7zgdlSgoajYbxY8cwb/YsKlepqntT9gsMxtLSki/rN8TfbwOQ+iZ148Z1Pvn0M+rVq8+dyEguXjgPwMkTx8HMjGLOztmaHcDM3JxpU7yJiEgd1Px811GuXHkUjo7U/aIegQEbdbkXzp/LqBG/4uCg0Mvm4KCgoWtjAvxTs8VERxMWepj6DRphY2PDwnlzdCfzHj50kKSkRKpWzfqnqYx6f65s+fL8e+0aN2/+B8DZs2dITEqkRMlS73S+OnXrcfbsGcLDbwLgv8GXhq6NM13T0LUxQYEBqFQqHj16xM4d22jUuAkKR0ecXYqzc8d2AI6EHsbc3Jyy5cpxKzycXt270eeH/vwyfKTu8Fz1j2uwe99B3ePRoWMnmrVoybgJ3u9VvgsXzjNk4ACUSiUqlYrly5bg1ir120wVK1Vi147Uc9iePn3KkdDDVKte/a3z5URWC3Nzxo4Zyalnb7D//nuN/27coGq16lSsVInDhw7oDift3bPbaHkyzvoF586e4dazHAF+vjRo5GpwzauKOTtToWIltgZvAlL3KJ85cypH9zoB1K5TlwvnznL7VjgAQQEb+LKBa6Zrbly/xrLFC1Cr1SQnJRGwYT2Nm7XArkAB7B0U7N+3G0g9zeP0PydNttcJoErNz/j38nmiI28BsG97IDXr1NerefL4Id6/9uWTLxoxYIQ31rlyA1DIXsF8n+14L/TBe6EPHt/1pXyVj/llonHO5f3/ykybU8finklSZXy7n+86fNevw9LCAnsHB0aM9uJhfDzjvUbjFxgMpP6cwQZfHzQaDTVq1mLEKC9y585NUlISs36bxvHjx9Bo1JSvUBGvcRNf+2oupL7gx3mNJjIyAjMzM/73yzDqfpF6vPjkiePMmjGdxMRErK2t+XX4SGrW+iTDvnNbvjnbc2NGDtf9VMGF8+f0sm3dEszyP5ah0ahRKBwZN8Ebp6JF0Wg0LF60gD27doKZGU5OToyd4I1CoXjt/p8mJDBp4jguXbyIRqOmd98fdF8jPhJ6mLmzZ6JUKsmbNy+jvca/8XyS3M923LwpX1q9W1haMKBfHxYsXoqDg4Ldu3aw5PdFmJlB7tw2DP35V2rW+uSdz3f40EHmzZ6JUqXE2aU43pOnERFxW2/dpVVjV6AAKpWKWb9N4+jRMFRKJe09O+q+eh8efpMJY8cQFxdHrlzWeI2bSMVKlRnvNZptWzdTsmQpXQ9W1tb4+Prr9fX7wvnExcVl+FMF72I+gPlzZ7Nv7240Gg2NXJswcPBQLCwsiIyMYMqkCURGRGBmbk7zFl/R99nJyW+bLyeynjh+jFm/TUOlUmFlbc3AwUP57PM6aLVaFi2Yx66d27G2tsbJqRijvcbpfbDISG5LSEh5u8106KGDzJ87C6VSibOLCxMnTyMy4jYTxo7Bd2NQujV2dvpf1KlZtQL7Dh3V/VRBVNQdpnqnriONRsM333Z9qy8v5LE2A+B+QiZX3ivCQg+xeP5slEoVxZxd8Jo4mcjICKZO8GKVb2C6NfntCpCUmMjMad5cOHcGlUqFa5Pm9B0wCDMzM65dvcysaZN59DAeM3NzOnTqoneYL7MK50l9ch678fANlW92+tgR/FYuRK1S4eBUjL4/jyM2KpI/53rjvdCH4PXLCVi7FJeSZfT+bviUheTL/2J9HtqzleOh+4zyUwW1S6d+ScCmxoAs31dm5eRPFSSeWpDube/c8PS+MmR4et8Y+ub0vpF877f/D/nednh612V1eHrXGXN4eheZYnjKSRkNT+/U7zwJIYQQQrzrZHgSQgghhDCADE9CCCGEEAaQ4UkIIYQQwgAyPAkhhBBCGECGJyGEEEIIA8jwJIQQQghhABmehBBCCCEMIMOTEEIIIYQBZHgSQgghhDCADE9CCCGEEAaQ4UkIIYQQwgAyPAkhhBBCGECGJyGEEEIIA8jwJIQQQghhABmehBBCCCEMIMOTEEIIIYQBZHgSQgghhDCADE9CCCGEEAaQ4UkIIYQQwgAyPAkhhBBCGECGJyGEEEIIA8jwJIQQQghhADOtVqs1dRNCCCGEEO8Ly5xeYJIqp5eYM3JbftjZQPK9r/6/5EtUmraP7GJjBQkpH+Zn3DzWZgDcvJ9k4k6yR8nCuQHYcCrSxJ1kj441igFgU2OAiTvJHomnFqR7mxy2E0IIIYQwgAxPQgghhBAGkOFJCCGEEMIAMjwJIYQQQhhAhichhBBCCAPI8CSEEEIIYQAZnoQQQgghDCDDkxBCCCGEAWR4EkIIIYQwgAxPQgghhBAGkOFJCCGEEMIAMjwJIYQQQhhAhichhBBCCAPI8CSEEEIIYQAZnoQQQgghDCDDkxBCCCGEAWR4EkIIIYQwgAxPQgghhBAGkOFJCCGEEMIAMjwJIYQQQhhAhichhBBCCAPI8CSEEEIIYQBLUzfwqhnTp7Jn107s7OwAKFGqFL/NnKNXs85nDX8sWUyRIkUAsM2Th5Vr1hm0HLVazYzpUwg7EopapaZr9x54duwMwMP4eKZOnsj169dJTk6iV59+tG7zdZazrfdZi9+G9ZiZmeHi4oLX+EkULlxYr+ba1StMnTyJx48fY2FhzpixE6hUuYpBy0lMTGS812guX76IRqNh8NBfcG3cBDDOY5eeQwcPMG/OTFJSUihXrjzjJk4mb968marJaH2Eh99k7OiRPIyPx8bWFu8p0yhVuozuPlNSUvipf1/ad+hI0+YtdNdHR0Xx7Tee+AcGU7Bgofcqm1arZeG8OezbuweAylWqMsprHDY2NrrlXbl8mf59e7HvYGiWspky366dO7CxsaH6xzX4edgIcuXKRdcunUhKTNQt9+bN//Bo78nwkaOznDOt3PPnzCRFmULZcuUZNyHt3GnVPM99NOxZ7u970OH5NuRh6jbkxvNtSO9+tDLCNsRQhw8dYP6cWSiVKZQtWx6vCd6v5XtTTXR0FN26dMR3YzAFCxbkxvV/GTnsZ93tGrWaf/+9xm+z59G4SbMcy/aqv48cYsXieSiVKZQqU44hI8eRJ0/e1+q0Wi0zvb0oUfojOnzTDYDk5CQWzJjM1UsX0Gg1VKhUlQE/jyRXrtw5HSNdV/75i72+f6BSpuBYvDTufX8ht20evZq/d27i2N7NmGFGIUVR2vT5H3ntCgIwtXdb8hcqoqv9onVHqtdrkqMZjGHp+G+5+G8Uc9bsM2kf79yepzOnTzFtxiz8AoPxCwx+bXACOHPqFD//OlxX8zZv/hv9fLkVHk5A0FbWbdiIz5pVnDt7FoAxo4bjoHDELyCIpX+sZNoUb2Kio7OU6+KF86xeuZzVPr4EBm+leImSLJw/V68mMTGRfr178n2PXvgFBNGnX39GvLSRyqzFC+dja2tL0JYdLFm2gskTx+n6N8Zjl5YHDx7gNXoEM+fMZ/O2XRRzdmHurBmZrslofYz49Wc8O3Zm05bt9B/wE0MHD0Sr1abmOX2Kbzt7cuqfk3rL2hIcRPeuXbgbG/teZtu3dw9hYUfwCwgicPM2kpISWbd2NQAqlYo1q1bSr08PEhIS3st8wUGBHDp4AJ8NG/ELDKaIvT0L5qW+1lf7+Oqen/0HDKRYMWd+/GlQlnOmlXvsmBHMmDOf4K27cHZ2Ye7s13OnV7PR35dbt8LZuGkrPr4b8Vm7inPnUnN7jRqOQuHIho1BLFm2kmlTs74NMVTcgweMGzOSGbPnsWnLToo5uzB/zkyDarZuDqJnN/3XUekyH+G7MUj37/O6X9DiKzeTDk7xcQ+Y6e3FmMkz+dN3M45Fi7F80dzX6m7dvMGwn3pzaN9uvevXr/wDtVrN76v9Wbx6IynJyfiu/jOn2n+jhEfxBC2eTqch4xg0ezUFHYqyZ/0yvZo7N65yZKsfvSfMZ8CM5RRyLEaI3woA7t25hU2efPSftkz3730bnMqXUrBjyU+0a1rT1K0ABg5PI0eOzK4+gNQ9CJcvXWTVyuV0aNuGoYN+IurOndfqzpw+xfbtW/Fs9zX9evfk2tUrutsCA/zp1MEDz3Zf06fn9/x343qaywrZtxf3th5YWlqS386OFl+5sW3rZh7Gx/PX0TD69R8AgMLRkbXr/cj/bE/Y26pUuQqbt+8iX758JCcnExsTQ4ECBfRqjoYdwdnFhS/rNwCgYaPGesPjsiW/07F9Wzw93Bn8U39iY2PSzebRvgMATkWLUqduPXbt3AFk/NhlxdGwUKpUqUqJEiUB8OzUme3btuiGnDfVpLc+YmJiuPnfDVq0dAOg3pcNSEpM5PKliwCsW7uGAQMHU7Vadd1yYmNjCAnZy4Lfl7632Zo0bcaqteuxsrYmISGBBw8eYGeX+ny5dPEiV69eYebsee9tvosXLtCocRPy588PQOOmzdi7e5deXw/j45k0YSyTpkwjX758Rsn6au7KlV9k6tCxMzvSyJ1eTcjevbh//SJ38xZubN+ymYcPU7chfX94aRuyLuvbEMPzHaFy5aoU1/XeKY186dfcjY1hf8g+5i9K/3X0z8kT7N2zm5Fe47Mzyhv9c+wo5StWoZhLCQBaeXgSsnu7XlaAzQG+NHNzp35j/UGv6sc1+eb73pibm2NhYUGZchWIjY7Ksf7f5N+zJyhapjyFnZwB+LRpG86G7tPLV7R0OQbNWUNu27woU1J4HHcPm7ypr69bVy9gZm7O8glDWfhrL/YHrEajUZsky9vq51mf1Zv/ImDPP6ZuBTBweLp06VJ29QGkvunV/uxzBg0eil9gMNWqV2fQT/31niBPnz6lZKnS9OrdF7+AINq2a0f/fr15mpDAiePH2BIcxIrVPvgFBPF9j14MHfRTmsuKjo7C0dFJd1mhcCQmJppbt25RxN6eNatW0K1LJzp7enDp4kW9wyVvy8rKipB9e2nmWp+TJ4/j3tZD7/bwm/9RpIg9Y8eMpLOnB317dUetTn2CbwkO4tq1q/j4+uMXGEy9Lxsw3ivtwxivZ1MQExOd4WOXVdFR0SgcHV9apiNPnjzR2zOSUU166yMmOgp7BwfMzV88VR0UCt2n+GkzZlG/QUO9XhwcFMyeu4AyH32U5VymzGZlZcV6n7W0aNKQuLg4XJs0BaBqtWpM9J6it7z3LV/VatU4sD+EuLgHaDQatgQHcfeu/l7C5X8uo179BlSuUtUoOV8VEx2N4xtyZ1QTE5PxNmTt6hV0+7YT3xhxG2KImOgovXXmkGa+9GvsHRTMnDOf0mXSfx3NmTmdH38a/NqhwJx2NyaaIgqF7rK9vYKnCU94+lR/2zbgfyNp8lXr1/6+1md1cS5eEoCYqDts8vPhS9em2dqzIR7ej8WusIPucv7C9iQnJpCc+FSvzsLSkkvHQ5n5oyc3L52lZsPU0xg0ajVlqtWi64ip9Bg7h3/PHOevnZtyNENWDZnmz/ptx03dhk6mznkaMWIEAHfu3NH9f8qUKUZvxtnZhYWLX+yK7Na9J0sXLyIyMgJnZxcAbG1tWbzsxe7U5i1asnTx75w/f44jhw9x+1Y43bp00t3+8OFDHsbHY/fKXh6NRv8TCYCFuTkqlZLIiAjy5M3LKp/UwxHdu3ahRIkSBp97lBbXxk1wbdyEAH8/fujTk6079ujeXFQqFaGHD7JsxWqqVavO/pC9/NivDzv37ufQwf2cP3eObzzbAaDWaEhKSkxzGWlmszDP8LGr/dnnWcql1WrSvP7lN86MatJbHxpNOn9jYfEWXcdnsWsAAA8mSURBVL4dU2br3OVbOn3ThYXz5vDzkIEsX7XWkNYzxRT5Wrf5mtiYGHr36IaNjS3tOnhiZWWlq0lOTiZgox++foGGRDFIev1ZvJQ7o5q0cptbmKNSPtuG5MnLqrWph/Z6dO1CcSNtQzJLk84608uXiZr0nDn9D/FxcXzl1urtGjQijfb1dQGZy/Gya5cvMn7EENq068TnXzQwRmtGoU3juQb6r9HnKn5aj4qf1uPEvq2snjKMQXPW8EnjF+vI0sqaum4d+GtnIHVbts+2nj90mXpmtW3blrZt21KgQAHd/7PD1SuX2bI5SO86rVaLleWLjeqdO5Gs81nzWo2lpSVqjQa31u668yV8N25ivV8A+e3s8PRw1/27cP4cTk5O3L17V3cfsbExKBSO2DukTvfuX6fuFSpeogQ1atbk/LNzGd7WrfBw/jl5Qnf5a492RN25w6NHD3XX2Ts4ULJUaao9OwTVyLUJGo2aiNu3Uas1dO/ZS5dtvV8Aq9asJzY2Ri9bbGxMarZ7L7LFxMaiUDhm+NhllaOTE/deeTzz57fD1tY2UzXprQ8np6Lcv3dPb+/j89tyiimyXbl8mUvPDk2amZnRtn0HLl288MHkexgfz1durdi4aQtr1m2gdJmPcCleQlcXevgQ5ctXwNnFJVsyAzg5OXHv3uuZbF7KnVGNUxqPiULhiMOzbUib59uQ4iX42AjbEEM5OhZNc529nC8zNenZvXMHbm3c03wDz2kOCkce3Lunu3zvbix58+Unt82bczx3YM8Ohg/qS48fBtG5W6/saPOtFSjiwJO4+7rLjx/cxSZPPqxzv9ibeT86kvDL53SXazb6ivi7MSQlPOb0od1Eh790CotWi4XFO/d9sfdKpp71tWvXpnbt2tjY2Oj+nx3MzM2ZNsWbiIjbAPj5rqNcufJ6u5VtbGxYOG+O7oTUw4cOkpSUSNWq1ahT9wt2bt+m2/3vv2E9vXumfpvi+dDhFxhM5SpVaejamKDAAFQqFY8ePWLnjm00atwEZ2cXKlaqzObg1F2a9+/d4/TpU1n+xHjv3l2G/TKUuLgHAGzfuoWPPipLgQIFdTX16tXnTmQkFy+cB+DkieNgZkYxZ2fqflGPwICNPHnyBICF8+cyasSvODgo9LI5OCho6NqYAP8NQOphh7DQw9Rv0CjDxy6r6tStx9mzZwgPvwmA/wZfGro2znRNeutD4eiIs0txdu7YDsCR0MOYm5tTtly5LPf8Lme7evUyY0eNIPHZt862BAdlee/gu5TvwoXzDBk4AKVSiUqlYvmyJbi1enE45eTxY3z2eZ1syauX6cyLTBvTy51OTcNGjQna9CL3rh3baOTahGLPtiFbXtqGnDHCNsRQdep+wbmzZ7j1rPcAP18aNHI1uCY9J08cp/Zn2buOMqtW7TpcvnCWyNvhAGwL8qfOlw0z/feHQ/awaPY0psxZjGuzltnU5dsrU+0Tbv97iftREQAc37uFCp/U1at5HHcf/3kTSXj2gfxs6D4cXEpim8+O2Ns3CfFfiUajRpmSzN+7gqhSp2FOx/igmGlfPaMuA3v37qVJk6ydoZ+kyvj2rVuCWf7HMjQaNQqFI+MmeGNhacGAfn1YsHgpDg4KjoQeZu7smSiVSvLmzctor/GUr1ABAN91PvhtWI+5uRl58uRlzLgJfPRR2deWo1KpmPXbNI4eDUOlVNLesyPduvcEIOrOHSZPmkBExG20Wg1dvutGB89Or93Hy3Jbvjmbn+86fNevw9LCAnsHB0aM9uJhfDzjvUbjFxgMpG6QZs2YTmJiItbW1vw6fCQ1a32CRqNh8aIF7Nm1E8zMcHJyYuwEbxQvHed/7mlCApMmjuPSxYtoNGp69/2BVq3dATJ87DLKBm/Od/jQQebNnolSpcTZpTjek6cREXFbL19aNXYFCmS4PsLDbzJh7Bji4uLIlcsar3ETqVipst6ye37/HZ06d9H7qQKA6pXLcyD0aIY/VZCZfKbItmjBPPbs2omFpQVlypRl+KjRejkiIyNo596av06cynC9vKv55s+dzb69u9FoNDRybcLAwUOxeHbIcsAPfWjR0k33vM1MvkTlG0vTzD1/TurrwdmlOJOmTCPi9m3Gjx2NX0BwujV2ds9yz5jGX0fDUCqVtO/w0jYk6g5Tnm9DNKnbkPZv2Iakx8YKElIyvZnWE3roIPPnznrWuwsTJ08jMuI2E8aOwXdjULo1z7+c8FzNqhXYd+goBQu++LBXt3YNNm3ekaVz7/JYmwFw837SW9/Hc8fCDrN88TxUSiVOxZz5xcub6MgIZk8dz++r/PRqZ0wao/dTBd09W5Pw5DGF7V+cV1S56scM+DlrX5IqWTj1pw42nIrM0v0AXD31F3vW/4FapaKQoigePw4nLiaK4KUz6D8t9XSXY7uDObY7GHMLC/IVLEyrHoMo6OBESnIS21bMI+LaJdRqFZU/a0CTTj0xMzPLUk8daxQDwKbGgCzny6yc/KmCxFML0r3NoOHJGN70Bvy+yszw9L7K7PD0vpJ877esDE/vg6wMT+86Yw5P7yJjDk/vIlMMTzkpo+HJ9AerhRBCCCHeIzI8CSGEEEIYQIYnIYQQQggDyPAkhBBCCGEAGZ6EEEIIIQwgw5MQQgghhAFkeBJCCCGEMIAMT0IIIYQQBpDhSQghhBDCADI8CSGEEEIYQIYnIYQQQggDyPAkhBBCCGEAGZ6EEEIIIQwgw5MQQgghhAFkeBJCCCGEMIAMT0IIIYQQBpDhSQghhBDCADI8CSGEEEIYQIYnIYQQQggDyPAkhBBCCGEAGZ6EEEIIIQwgw5MQQgghhAHMtFqt1tRNCCGEEEK8L2TPkxBCCCGEAWR4EkIIIYQwgAxPQgghhBAGkOFJCCGEEMIAMjwJIYQQQhhAhichhBBCCAPI8CSEEEIIYQAZnoQQQgghDCDDkxBCCCGEAWR4EkIIIYQwwAc5PG3ZsoWWLVvStGlTfHx8TN1Otnjy5AmtWrUiIiLC1K0Y1YIFC3Bzc8PNzY3p06ebuh2jmzt3Li1btsTNzY0VK1aYup1sM23aNIYPH27qNoyua9euuLm54e7ujru7O2fOnDF1S0YVEhKCh4cHLVq0YNKkSaZux6j8/f11683d3Z1atWoxYcIEU7dlNMHBwbpt57Rp00zdjtEtXbqU5s2b07p1a37//XdTtwPaD0x0dLS2UaNG2ri4OG1CQoK2devW2mvXrpm6LaM6ffq0tlWrVtrKlStrb9++bep2jObIkSPajh07apOTk7UpKSnarl27anfv3m3qtozm77//1nbq1EmrVCq1iYmJ2kaNGmmvX79u6raMLiwsTPvZZ59phw0bZupWjEqj0Wi/+OILrVKpNHUr2eLWrVvaevXqaaOiorQpKSnazp07aw8cOGDqtrLF1atXtU2bNtXev3/f1K0YxdOnT7Wffvqp9v79+1qlUqlt37699siRI6Zuy2iOHDmibdWqlfbx48dalUql7du3r3bXrl0m7emD2/MUFhbG559/ToECBbC1taV58+bs3LnT1G0ZlZ+fH2PHjsXBwcHUrRiVvb09w4cPx9raGisrK8qUKcOdO3dM3ZbR1K5dm9WrV2Npacn9+/dRq9XY2tqaui2jio+PZ/bs2fTr18/UrRjdjRs3MDMzo3fv3rRp04a1a9eauiWj2rNnDy1btsTR0RErKytmz55N9erVTd1Wthg3bhxDhgyhUKFCpm7FKNRqNRqNhsTERFQqFSqVily5cpm6LaO5ePEi9erVI2/evFhYWPDll1+yd+9ek/b0wQ1PsbGx2Nvb6y47ODgQExNjwo6Mz9vbm08++cTUbRhd2bJl+fjjjwG4efMm27dvp0GDBibuyrisrKyYN28ebm5u1KlTB4VCYeqWjMrLy4shQ4aQP39+U7didI8ePaJOnTosXLiQlStX4uvry5EjR0zdltGEh4ejVqvp2bMnbdq0Yd26ddjZ2Zm6LaMLCwsjKSmJr776ytStGE3evHkZNGgQX331FfXr16dYsWLUrFnT1G0ZTeXKlQkNDSU+Pp7k5GRCQkK4d++eSXv64IYnrVb72nVmZmYm6ES8rWvXrtGjRw+GDRtGyZIlTd2O0Q0cOJCjR48SFRWFn5+fqdsxGn9/f5ycnKhTp46pW8kWNWrUYPr06dja2lKoUCHat2/PwYMHTd2W0ajVao4ePcpvv/2Gn58f586dY9OmTaZuy+h8fX3p3r27qdswqsuXLxMQEMD+/fsJDQ3F3NycP//809RtGU2dOnXw8PDgu+++o1evXtSqVQsrKyuT9vTBDU8KhUJvIo2Njf3gDm99yE6ePMn333/P//73P9q2bWvqdozq+vXrXLp0CQAbGxuaNWvGlStXTNyV8Wzfvp0jR47g7u7OvHnzCAkJYfLkyaZuy2hOnDjB0aNHdZe1Wi2WlpYm7Mi4ihQpQp06dShUqBC5c+emcePGnD171tRtGVVKSgrHjx/H1dXV1K0YVWhoKHXq1KFw4cJYW1vj4eHBsWPHTN2W0Tx58oSmTZuyZcsW1qxZg42NDS4uLibt6YMbnurWrcvRo0d58OABiYmJ7N69m/r165u6LZEJUVFR/Pjjj8yYMQM3NzdTt2N0ERERjB49mpSUFFJSUti3bx+1atUydVtGs2LFCrZu3UpwcDADBw7E1dWVkSNHmroto3n8+DHTp08nOTmZJ0+esGnTJpo2bWrqtoymUaNGhIaG8ujRI9RqNYcPH6Zy5cqmbsuorly5QsmSJT+4cw0rVKhAWFgYT58+RavVEhISQtWqVU3dltFERETw448/olKpePz4Mf7+/iY/7PrhfGx6RqFQMGTIELp27YpSqaR9+/ZUq1bN1G2JTPjzzz9JTk5m6tSpuus6depE586dTdiV8TRo0IAzZ87w9ddfY2FhQbNmzT7IIfFD1ahRI93602g0fPPNN9SoUcPUbRlN9erV6dWrF9988w1KpZIvvviCdu3amboto7p9+zaOjo6mbsPo6tWrx8WLF/Hw8MDKyoqqVavSp08fU7dlNBUqVKBZs2a0adMGtVrN999/b/IPnmbatE4SEkIIIYQQafrgDtsJIYQQQmQnGZ6EEEIIIQwgw5MQQgghhAFkeBJCCCGEMIAMT0IIIYQQBpDhSQghhBDCADI8CSGEEEIYQIYnIYQQQggD/B9AgzI2Dn3LQgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "sns.heatmap(B,square=True,cbar=False,linewidth = 1,cmap='Blues',yticklabels=list(transition.columns),annot=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P =  1.4366902641444729e-06\n"
     ]
    }
   ],
   "source": [
    "P = np.sum(np.array(emission[seq[0]]*B[:,0]*init))\n",
    "print(\"P = \",P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, P is the probability $P(x)$. We can see that the forward and backward algorithms give the same answers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing two HMMs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's return to the CpG island problem and see if we can decide which model is more likely to have generated a sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "transition_1 = pd.read_csv(\"cpg_hmm_1/transition.csv\",index_col=0)\n",
    "emission_1 = pd.read_csv(\"cpg_hmm_1/emission.csv\",index_col=0)\n",
    "init_1 = pd.read_csv(\"cpg_hmm_1/init.csv\",index_col=0)\n",
    "\n",
    "transition_2 = hmm.transition_Durbin_CpG(1/20,1/20)\n",
    "emission_2 = pd.read_csv(\"cpg_hmm_2/emission.csv\",index_col=0)\n",
    "init_2 = pd.read_csv(\"cpg_hmm_2/init.csv\",index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "seq = (\"CGGCGGCGCATGCATCTCTCGAAGCAGCAGCAGCAACAGAGGCGTCT\"\n",
    "       \"CATAAGTGCATGCGCAGCAGCGGTAGCAGCGGTGGTGGCAGCAGCAGGGAGATTTCATTGCTGCAGTCGCCAAATTACATGCATCCTACGACT\")\n",
    "# a portion on the D.mel X chromosome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.4294265533073766e-84\n"
     ]
    }
   ],
   "source": [
    "P_1,_ = hmm.forward(seq,transition_1,emission_1,init_1)\n",
    "print(P_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.0719510945970887e-84\n"
     ]
    }
   ],
   "source": [
    "P_2,_ = hmm.forward(seq,transition_2,emission_2,init_2)\n",
    "print(P_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P_1>P_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, the sequence is more likely to have come from the refined HMM in this case."
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}