{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#Quickstart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Install\n",
    "```\n",
    "# install via pip\n",
    "# requires python 3; tested for python 3.4 and up\n",
    "pip install pytrails \n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's assume, we have observed a sequence of rainy and sunny days.\n",
    "To explain what we have observed, we have two hypotheses:\n",
    "\n",
    "* `sticky`: when it is sunny it stays sunny, and when it rains it stays rainy\n",
    "* `random`: the weather is totally random\n",
    "\n",
    "We want to test, which one of these hypotheses is more plausible given our observations.\n",
    "To this end we use [HypTrails as introduced by Singer et al. in 2015](http://www.www2015.it/documents/proceedings/proceedings/p1003.pdf).\n",
    "\n",
    "In particular, we \n",
    "\n",
    "* first we calculate transition counts between states\n",
    "* define hypotheses as transition probability matrices\n",
    "* calculate the evidence for each hypothesis using a set of increasing concentration factors\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd4VFX+x/H3NwUCqAgEQWpQAaUH\nIhYEREBBVMT9oSCiIOquBdRdsawoy4rrFuxie0Rsq6tiV1BB4UFWXQ2shSKKAitSDKEYQErI+f1x\nEpIAKZCZuVM+r+eZJzP3zsz9TtDPnJx77jnmnENEROJfUtAFiIhIZCjwRUQShAJfRCRBKPBFRBKE\nAl9EJEEo8EVEEoQCX0QkQSjwRUQShAJfRCRBpARdQEnp6ekuIyMj6DJERGLK/Pnz1zvn6lf0vKgK\n/IyMDLKzs4MuQ0QkppjZyso8T106IiIJQoEvIpIgFPgiIglCgS8ikiAU+CIiCUKBLyISkIYNwWzf\nW8OG4TmeAl9EJCDr1h3Y9qqKqnH4IpK4Gjbcf9A1aABr10a+nvI4B/n5sHMn7Njhf6akQHq63//V\nV7BtW+n9DRtCly5+/xNPwK+/Rr5uBb6IRIXyWru5ubBrV3FXx8qVkJNTHKY7dkByMvTt6/e/9x78\n8EPpwK1dG0aP9vvvvhuWLCm9PyMD7r3X7x82DBYtKv3+J5wAr7zi97dsCd9/X7rOgQPh9df9/b59\n4eefS+8fNgyee87fHzNGgS8ihaKxteuc/2kGeXk+0H791d+2bfM/Tz0V0tIgOxs+/rh4e9FzJk3y\n+6dMgRdfLP3a8qSnQ61asGWLf3zLLfDCC6Wf07AhrFnj70+eDG+9VXr/MccUB/68efDZZ1C9OlSr\n5n/WrFn6eBkZxfuqVYM2bYr3jxkDv/xSev/RRxfvf+YZKCgo/f71S0x8sGyZ316/wskQQstc0b9i\nFMjKynKaWkHEh2pZSv4v65xvfRaF5uGH+2DcuBG++GLfwD3zTGja1Hc5PP30voH997/DscfCa6/B\nH/+47+u/+gratoX774frrtu3thUroHlzuPNOGDeueHtamg/U776DunXhwQfh+eehRg1/q1kTpk0r\n+zM/8IB/3mWX+cfZ2f6Lr2Tg1qwJHTv6/Tk5sHt36f0pKeX/XoNQ2X/nit/H5jvnsip6nlr4IgHY\ntcuH3+bNpW9du0KnThW/PjsbevSA7dtLB8OLL8L558P8+cXdGyW9844P/BUr4LHHSgdujRrFLe26\ndaFDh9L7atb028G/99NPF28vek6DBn7/tdfC737nt1evDkl7DQ8ZPbq4tV2kvPDb+7lZFURbpFvO\nB6tBg7L/kgsHtfBFKqmgwPfnpqX5x/Pm7RvYHTtC//4+iAcP3nf/9dfD7bf71umRR+57jLvugptv\nrrjlt2qVb2XvHdi9e/uuhY0b4csv9w3sevV8AEejULV2E1FUtPDNbDRwNbAbeMc5d2M4jyfxKRT9\n2fn5PnDz84tbTzNn+teXDOSjj4bf/tbvHzDAnxws2peXBxddBM8+6/f37euDvaQrrvCBX60arF4N\nhx0GLVr4E4aHHVbcMq1XD/71L7+9aF/t2n57ZTRpAv/4R9n769Tx/emxJNKt3UQUtsA3s17AQKCj\nc26HmR0RrmNJfCtv9Mb995cO7Lp1fT80+K6Nolb4tm1+W7dufhv4PujFi4vfLy0NzjqrOPDr1PHb\nSoZy587Fz3/3Xd9yLrm/Rg2/LynJd6uUJTUVLrjgwH8X8Szahl7Go3C28K8E/uqc2wHgnPu5gudL\ngtq6FX76yXdTrF/vgxrgb38rHgZXlqIThzVr+tBt3754X7t2xS3nolvJ9XVefdUP5SsK6727OoqG\n0JWlZ89KfbyDotauhEM4A78V0N3M7gS2Azc45z4P4/EkyhQU+AD/6afiQP/pJ7jtNt/lceedvlti\n8+bi15jBoEG+BVxQ4FvZ5Vm/3od1auq++26/vfzXtm594J8pUtTalXCo0klbM5sF7G/Wh1uBO4HZ\nwBjgeOBF4Ci31wHN7ArgCoBmzZp1WbmyUgu3SBTYuNFfnFIU6EWhfu+90KiR71q56abSr0lK8iNE\nmjb1Q/8+/BAaNy59a9my9KgOncwTKV9ETto65/qUU8CVwKuFAf+ZmRUA6UDOXu/xOPA4+FE6ValH\nqs45f0tK8hexzJhRunX+00/w6KP+qsMZM/zVg0Vq1vSBnZvrA79vX9+n3bixP8nYuLHvkkgp/K9u\n0CB/E5HICGeXzutAL2C2mbUCqgHrw3i8hFCVESu7d/vn1KjhT26uWeNPeu7d5fLEEzB0qB8nPmqU\nf239+sUt8KLWd69e/sRlUaDXrl26NZ6Z6W9Vpf5skdAIZ+A/CTxpZguBncAle3fnyIErb8TKt9/6\nwD7iCH815MaNPrCLAn3NGt8v/re/wY03+iGF99zjW+ONG/twPuss36UCcPzxfj6SRo32P3b7yCP3\nP5Y81NSfLRIaYQt859xO4KJwvb/sq+gk5OjR/lL0WrX8l0CjRn4ekKKWeLdu/nkZGT70974KskiN\nGn4MuYjEB02tEOWc87P6zZrlb+V55hkf6q1a+cfVqsHChWU/v2ixBRFJDAr8KJSbW3zF5cCBxbP+\nHXNM+a8bPjy8dYlIbNOKV1Fg0yY/j/Y11/iZChs1Kr4ydPhwfxJ1+XJ/ElVE5GCphR+Aon7zatVg\n6lQ/5WtBge9z79HDz8eSn++fO3hw6ddqxIqIHCwFfgQUFPi5yYv64efN83OBn3uuH88+bhz06ePv\nV6tW/ntpxIqIHCwFfhgULUqRluZnTGzfHjZs8PvatvUt+KOO8o/btIEJE4KrVUQShwI/RH7+2U8T\n8MEHvhV/6qm+u+bII+HCC33rvXfvyIxbFxHZHwX+QcrPLz1FQNHixbVrw2mn+S4a8MMeH3wwmBpF\nREpS4FdSfj58/nlxP/yyZfDjj/7k66mn+qtS+/Tx86Wn6LcqIlFI0VSGokkgzGDKFL80XV6ef5yZ\n6Vc++vVXP7Lm2muDrVVEpDIU+CWsWuX74Iv64V9+2U9D0LKl74fv08dPGFbZZehERKJJQge+c77F\n/v33ftKwb77x2+vX9ydYi5ar69HD30REYllCBf6OHfDxx771/sEH/uTqX/7iJxRr1Qouv9y34tu1\nK3tCMRGRWBXzgd8w+WfWFey7PnqDpJ9Zu7t4+5Ah8Oabvt89ORm6di1e37R6dXjjjQgVLCISkJgP\n/P2FfdH2U07xV7WCb8UXteB79vTroIqIJJKYD/zyHHMM7NzppyuYNCnoakREghXXgf/UU0FXICIS\nPeL71OTu3UFXICISNeI78O+4I+gKRESiRswHfoOkn/e/3dbBVVf5B0uWwP/+F8GqRESiT8wH/trd\nR+Ac+9zWFjSAIwpH8Fx9tV9KasKE4qWkREQSTMwHfqU89RScfTb86U9w3HHw0kvFk+WIiCSIxAj8\nZs3gxRdhzhyoUwcuuACeeSboqkREIioxAr9Iz54wf75v8V9wgd/2+eewfn2gZYmIREJiBT74eRUu\nucSvP7h7t59zoWVLv0pJ0crhIiJxKPECv6TkZD/BTlYWjBkDnTr5mdVEROJQYgc++FXF338fXnvN\nj+Dp2xfmzg26KhGRkFPgg58U/9xzYfFiv/J49+5++9y5sGVLsLWJiIRI2ALfzF40sy8KbyvM7Itw\nHStk0tJgxAj/BbB5s18VpXVr+Oc/NYxTRGJe2ALfOXeBc66Tc64T8ArwariOFRa1a8O778KRR/oF\nbE85xY/wERGJUWHv0jEzA84HXgj3sULu5JPhs8/8KubLlvlVU5YvD7oqEZGDEok+/O7AOufcdxE4\nVuglJcGll8K338Kzz0KLFn77rFl+sn0RkRhRpcA3s1lmtnA/t4ElnjaUclr3ZnaFmWWbWXZOTk5V\nygmv2rXhwgv9/e+/hzPOgA4dYMaMYOsSEakkc2E8GWlmKcBPQBfn3KqKnp+VleWys7PDVk9IvfMO\nXH89fPcdDBgA997rL+ASEYkwM5vvnMuq6Hnh7tLpA3xTmbCPOQMGwMKF8I9/+OGbWVmQlxd0VSIi\nZQp34A8hFk/WVla1anDDDb5//5ln4NBD/fDNGTOgoCDo6kRESglr4DvnRjjnHg3nMaJCw4YwsPC0\nxZw5cOaZcNJJ8J//BFqWiEhJutI21Hr29K39H3+EE0+Eiy+G1auDrkpERIEfcklJMHw4LF0Kt9zi\n5+Hv2VNdPCISOAV+uBx6KPzlL35+nkcf9V8Eu3bBe+9pmgYRCYQCP9yOPhp69/b3n30W+vXztyVL\ngq1LRBKOAj+Shg+H++7zJ3M7dPDj+DdtCroqEUkQCvxISk2Fa6/1F2uNHAn33w+DBwddlYgkiJSg\nC0hI9evD44/DlVf6ZRYBNm6ERYv8rJwiImGgFn6QMjP9FboAd9/tF1658EJYFX8XJotI8BT40eKW\nW+C22/xSi61bw8SJ8OuvQVclInFEgR8tatWCP//Zj97p39+H/7XXBl2ViMQR9eFHm4wMmDYNZs+G\nZs38thUr/MRs7dsHWZmIxDi18KNVr15+DD/AuHHQqRNccw1s2BBsXSISsxT4seD++/2InkcegfR0\nv8j63reGDYOuUkSinAI/FtSrBw89BF98Ufa0DOvWRbYmEYk5CvxYoj58EakCBb6ISIJQ4MeTa66B\nrVuDrkJEopQCP548/LCflO2HH4KuRESikAI/1jRoUPb2OXP8dA1Nm0a0JBGJDQr8WLN2rR+ps/dt\n7Vro0cNftJWa6idj691b6+qKyB4K/Hi1cqWfhvnkk/08PTt2BF2RiARMgR+vOnWCr7+GESPgr3/1\ns3IuWBB0VSISIAV+PKtdG6ZMgbffhtxcmDAh6IpEJECaPC0RDBjgF1fZudM/XrkSfvlFF3KJJBi1\n8BNFnTrFI3zGjvVdPHfdBfn5wdYlIhGjwE9EkyfDOefAH//ol1RcujToikQkAhT4iah+fXjpJXjh\nBT+Sp1MnmDUr6KpEJMwU+InKDIYM8X37w4dD165+e0FBsHWJSNiELfDNrJOZfWpmX5hZtpl1Ddex\npAoaNoTHH4fDDvNj9U86yU/FrOAXiTvhbOH/HZjgnOsE3F74WKLZli3+5O7o0dC3rx/NIyJxI5yB\n74DDCu/XBlaH8VgSCvXqwYwZvsX/2Wd+2OYTT5S96IqIxBRzYfqf2cyOA94DDP/FcrJzrtwmY1ZW\nlsvOzg5LPXKAVqyASy+FTZv8fDypqUFXJCJlMLP5zrmsip5XpQuvzGwWsL/FVG8FegPXO+deMbPz\ngSlAn/28xxXAFQDNmjWrSjkSShkZfuRObm7xZGzvvutP9JoFXZ2IHIRwtvA3A4c755yZGbDZOXdY\nea9RCz+KTZwIt90GAwfCY4+VPU2ziERcZVv44ezDXw30LLx/GvBdGI8l4XbLLTBpkm/lt23rx/GL\nSEwJZ+BfDtxtZl8Cf6Gw20ZiVHIy/OEP8N//wlFHwQUX+C8AEYkZYZs8zTk3D+gSrveXgBx3HHz8\nMdxzD1x4od+2fTukpQVbl4hUSFfayoFLSYEbb4RGjfwFWv37w8UX+xO7IhK1FPhSNQUFfmnF55+H\ndu38OH4RiUoKfKmalBS/sMp//uOv0j3zTLjsMsjLC7oyEdmLAl9Co0sXmD8fbroJZs/W1bkiUUiB\nL6FTvbpfP/frr4snY5s4EbZuDboyEUGBL+FQs6b/OWsW3H47dOwI8+YFW5OIKPAljAYMgDlzfPdO\njx5+HP+vvwZdlUjCUuBLePXoAV9+CVde6cfujxwZdEUiCStsF16J7HHIIX4d3UGD/IIr4EfxVKvm\n+/1FJCLUwpfI6dPHj9UHuOoqyMqCBQuCrUkkgSjwJRhDhvipl084Af70J9i1K+iKROKeAl+CMWAA\nLFzoJ2GbMMEH/7ffBl2VSFxT4Etw6taF556DV1+Fbdt8X7+IhI0CX4I3aBAsWlQ8Gdt118E33wRd\nlUjcUeBLdEhO9j+//x6efRYyM/0wzt27g61LJI4o8CW6tGzpW/t9+/oLtU49FZYtC7oqkbigwJfo\n07AhvPEGPPWUn5fn2GP9wul734rG9ItIpSjwJTqZwSWX+JE8ZXXrrFsX2ZpEYpwCX6JbkyZBVyAS\nNxT4Ets0775IpSnwJbbdfHPQFYjEDAW+xLYLL/Q/f/7ZT9UgImVS4Ev0a9Cg7O0dO/r7N9zgh3Q+\n+KDm5REpgwJfot/atb6vfu/b2rXFz7npJr+u7pgx0KkTzJwZXL0iUUqBL/GhbVt4/30/fn/HDjj9\ndHjggaCrEokqCnyJH2Zwzjn+St2//x3OP99vX7ECfvkl0NJEooECX+JP9eowdmzxlbgXXwytWsHU\nqX5yNpEEpcCX+Hf33dCiBVx6KXTtCh9/HHRFIoEIW+CbWUcz+8TMvjazt8zssHAdS6Rcxx/vQ/65\n52DNGujWzc/BL5JgwtnCfwK42TnXHngNGBvGY4mUzwyGDYOlS+Guu6B/f7998WL49ddgaxOJkHAG\nfitgbuH9mcBvwngskco55BB/dW6NGn68/llnwXHHwcsva5oGiXvhDPxFwMDC+4OBpmE8lsiBS031\nJ3IPP9yP6Dn1VPjii6CrEgmbKgW+mc0ys4X7uQ0ELgWuMrP5wKHAzjLe4wozyzaz7JycnKqUI3Lg\nevaE+fPhscd8907nzvDpp0FXJRIW5iLwZ6yZtQKec851Le95WVlZLjs7O+z1iOzXpk3w5JN+Td2k\nJP9F0L49VKsWdGUi5TKz+c65rIqeF85ROkcU/kwCxgGPhutYIiFx+OHw+9/7sN+4EU47zQf+9OlB\nVyYSEuHswx9qZt8C3wCrgalhPJZIaB1+ODz/vD+RO2CAvy1dGnRVIlUSkS6dylKXjkSdnTv9DJx/\n/rMfvvndd9C8edBViZQSeJeOSFyoVg3+8Acf9JMnF4f9v/9d9lq7IlFKgS9SGUccAZdf7u9/9x30\n6AFZWTB3bvmvE4kiCnyRA3XMMfDCC7Bhgx/WecEFsHJl0FWJVEiBL3KgzPyFWkuWwIQJ8NZbftGV\nvLygKxMplwJf5GDVrAm33+5H7zz2GBx6qN/+wQeapkGikgJfpKqaNi1ebGXOHOjTB045BTTiTKKM\nAl8klLp3hylTYNkyP/f+pZeWXntXJEAKfJFQSk72If/dd3DDDX4O/lNO0RBOiQoKfJFwOOwwv67u\nokXw8MP+iyA/3y+0rv59CYgCXyScWraE00/3959/Hs44A/r18zNzikSYAl8kUoYOhQcegM8+gw4d\nYMwYP5ZfJEIU+CKRkpoKo0f7/v0rrvBTNZx3XtBVSQJJCboAkYSTnu779X/3O9i+3W/bvNnPv3/a\nacHWJnFNLXyRoHTo4Idugu/q6d3bt/h/+CHYuiRuqYUvEg3GjoWUFLjzTr+oemoqbN267/MaNNC4\nfjloauGLRIO0NLjlFvj2WxgyZP9hD7BuXWTrkriiwBeJJo0awdNPB12FxCkFvkis0YVbcpAU+CKx\nJjMT3ngj6CokBinwRWLN9u0we7a/7xwUFARbj8QMBb5INGrQoOztixbBxIn+8Zw50K6dn7ZBE7RJ\nBRT4ItFo7Vrfet/7tnatn4jtkEOKn5ucDMOGQZs28MwzfpI2kf1Q4IvEsl694Msv4ZVXoEYNuOQS\nv8C6TuzKfijwRWJdUpK/QnfBAnj9dbj6ar/ubn4+PPss7NwZdIUSJRT4IvEiKQkGDvTdOwDTp8PF\nF8Mxx8Ajj8COHcHWJ4FT4IvEq7PPhnffhSZN4Kqr4Oij4aGHYNeuoCuTgCjwReKVmV9w5d//hpkz\noUUL39JPTvb71c+fcBT4IvHODPr0gblz/S0pCfLy/HDOu+8ue94eiTsKfJFEYQb16vn7ublw5JF+\nofUWLfz6u1u2BFufhF2VAt/MBpvZIjMrMLOsvfbdYmbLzGypmZ1RtTJFJKQyMmDWLJg3Dzp3hptu\n8ts09XJcq2oLfyFwHjC35EYzawMMAdoC/YCHzSy5iscSkVDr1s2f2P30U7jySmjY0G9/5x3YtCnY\n2iTkqhT4zrklzrml+9k1EPiXc26Hc245sAzoWpVjiUgYnXAC3HGHv79+vR/Xn5EB48fDxo2Bliah\nE64+/MbAjyUeryrctg8zu8LMss0sOycnJ0zliEilpaf7Fn/v3vDnP0Pz5jBuHGzYEHRlUkUVBr6Z\nzTKzhfu5DQxFAc65x51zWc65rPr164fiLUWkqjIz/XQNX34J/fr5k7pFLX0N54xZFa5p65zrcxDv\n+xPQtMTjJoXbRCSWdOgAL70Ea9b4UT3gl2Bs2tSvw1vWrJ4SlcLVpfMmMMTMqptZC6Al8FmYjiUi\n4VYU9vn5UL063Huv7+O//nr/ZSAxoarDMgeZ2SrgJOAdM3sPwDm3CHgJWAy8C1ztnNNk3SKxLiXF\nT8H8zTe+pf/gg34c//TpQVcmlWAuivrjsrKyXHZ2dtBliEhl/fADTJoEd94Jder4GTvT06FZs6Ar\nSyhmNt85l1XR83SlrYgcvKOOgocf9mEPfpK2Y46B3/4WVqwItDTZlwJfRELn5Zfh8svhqaegZUsY\nNQq+/z7oqqSQAl9EQqdpU5g82Xf1XHWVX2v3ww+DrkoKKfBFJPQaN4b77/fBf8klftvDD8NFF8GS\nJcHWlsAU+CISPkceCdWq+ftbtsBrr0HbtjB0KCxaFGxtCUiBLyKRceON/kTuTTfB22/7+fj/9Keg\nq0ooCnwRiZz69eGuu3zwjxvnZ+sEyMnxwznN9r0VzeApVVbh1ApB27VrF6tWrWL79u1BlxLV0tLS\naNKkCampqUGXIlKxevWKZ+cE39+fm7v/565bF5maEkDUX3i1fPlyDj30UOrVq4eZBVRZdHPOkZub\nS15eHi1atAi6HJEDt3kzHH542fujKKeiUdxceLV9+3aFfQXMjHr16umvIIldtWuXv3/btsjUEeei\nPvABhX0l6Hckca1BA7j6ag3prKKYCPxodN9997GtRKvjzDPPZFM5S8KNGDGCadOmRaI0kfjzm9/A\nE09AmzZw+unw1luwW/MxHqj4CvyGDSN2ln/vwJ8+fTqHl9cHKSLlK2tu/QYN/FQNP/4IEyfC4sVw\nwQW+3x/Uv38A4ivwyzqbX8Wz/Fu3bmXAgAF07NiRdu3aMWHCBFavXk2vXr3o1asXABkZGaxfvx6A\nZ555hg4dOtCxY0eGDx++z/vddtttjBgxgg8++IBzzz13z/aZM2cyaNCgKtUqErPWrvXhvfdt7Vq/\n/4gj4NZbYfly+OgjqFvX7+/Vyy/AvnhxsPXHgKgflrmPU0/dd9v55/t5Oyqyfj383/+V3jZnToUv\ne/fdd2nUqBHvvPMOAJs3b2bq1KnMnj2b9PT0Us9dtGgREydO5OOPPyY9PZ0Ne60DOnbsWPLy8pg6\ndSoAV199NTk5OdSvX5+pU6dy6aWXVvw5RBJZaip06eLv79gBRx/t/wJ49FG/Du/o0XDWWZCcHGiZ\n0Si+Wvhh0r59e2bOnMlNN93ERx99RO1yRhR8+OGHDB48eM8XQd26dffsu+OOO9i8eTOPPvooZoaZ\nMXz4cJ577jk2bdrEJ598Qv/+/cP+eUTiRloaTJniu3vuuguWLoVzz/WTtsk+Yq+FX4kWeZnS0w/q\n9a1atWLBggVMnz6dcePG0bt374M6/PHHH8/8+fPZsGHDni+CkSNHcvbZZ5OWlsbgwYNJSYm9fxKR\nwKWnw803ww03wBtvwJln+u2PPw7Z2b7V3759sDVGAbXwK2H16tXUrFmTiy66iLFjx7JgwQIOPfRQ\n8vLy9nnuaaedxssvv0xu4VWDJbt0+vXrx80338yAAQP2vLZRo0Y0atSIiRMnMnLkyMh8IJF4lZLi\nR/TUqOEfr1kDzz7rF2Pv1ctP3pafH2yNAYqv5mSDBvs/QVvW2f9K+vrrrxk7dixJSUmkpqbyyCOP\n8Mknn9CvXz8aNWrE7Nmz9zy3bdu23HrrrfTs2ZPk5GQyMzN56qmn9uwfPHgweXl5nHPOOUyfPp0a\nNWowbNgwcnJyOO6446pUp4jsZfx4uOYa3+0zeTKcd54/5/fii0FXFoion1phyZIlcR+E11xzDZmZ\nmYwaNapK75MIvyuRg5af78fv160LPXv60T+33ea/EDp2DLq6KombqRXiXZcuXfjqq6+46KKLgi5F\nJL6lpMCgQT7swfft//Of0KmT3zZtWtx398RXl04Mmj9/ftAliCSms86CVavgySfhoYdg8GBo3tyP\n569ZM+jqwkItfBFJXHXr+pE9338Pr78OI0YUh/3998N//xtoeaGmwBcRSU6GgQOLV+DasMEv0NK5\nM3TvDi+/DLt2BVpiKCjwRUT2Vrcu/O9/cPfdsHq1H9nTogX8+99BV1YlCnwRkf2pUwd+/3v49lt4\n800/lr9VK7/v009hwYJg6zsICvwAlJxoTUSiXHIynH02TJ/u1+QFP5yzSxe/Ju+LL8ZMd09cBX4k\nZkd2zlFQUBC6NxSR2DNtGtx7r7/Qc8gQyMjw8/VHuSoFvpkNNrNFZlZgZlklttczs9lmtsXMHqp6\nmZUTptmRWbFiBa1bt+biiy+mXbt2jBo1iqysLNq2bcv48eP3PC8jI4Px48fTuXNn2rdvzzfffANA\nbm4up59+Om3btuWyyy6j5MVu99xzD+3ataNdu3bcd999e4537LHHMmLECFq1asWwYcOYNWsW3bp1\no2XLlnz22WdV+0AiUjW1a8N11/nunrff9vP0FLXyt271Y/yjkXPuoG/AcUBrYA6QVWJ7LeAU4HfA\nQ5V9vy5duri9LV68uNTjnj33vU2e7PftfzJtf3POuZycfV9bGcuXL3dm5j755BPnnHO5ubnOOefy\n8/Ndz5493Zdffumcc6558+bugQcecM45N3nyZDdq1CjnnHOjR492EyZMcM459/bbbzvA5eTkuOzs\nbNeuXTu3ZcsWl5eX59q0aeMWLFjgli9f7pKTk91XX33ldu/e7Tp37uxGjhzpCgoK3Ouvv+4GDhy4\n3zr3/l2JSAQVFPifjz7qQ+fEE517/nnnduwI+6GBbFeJjK1SC985t8Q5t3Q/27c65+YBcbOqdvPm\nzTnxxBMBeOmll+jcuTOZmZn4YAFDAAAGBElEQVQsWrSIxSUWXjjvvPMAfwXtihUrAJg7d+6eK2kH\nDBhAnTp1AJg3bx6DBg2iVq1aHHLIIZx33nl89NFHALRo0YL27duTlJRE27Zt6d27N2ZG+/bt97yv\niESRonWlhw71Y/hzc+HCC/3FXBMmRMVVvIFfaWtmVwBXADRr1qzC5wcwOzIAtWrVAmD58uVMmjSJ\nzz//nDp16jBixAi2by/+XqtevToAycnJ5FfhH7jofQCSkpL2PE5KSqrS+4pImB12GIwZ4+foee89\nePBBeP99P5EbwMqV/ksgABW28M1slpkt3M9tYCgKcM497pzLcs5l1S86Ax7FfvnlF2rVqkXt2rVZ\nt24dM2bMqPA1PXr04PnCBRlmzJjBxo0bAejevTuvv/4627ZtY+vWrbz22mt07949rPWLSIQkJUH/\n/n50zwcf+G3r10Pr1nDCCX4enwiuww2VaOE75/qE5chhEKbZkUvp2LEjmZmZHHvssTRt2pRu3bpV\n+Jrx48czdOhQ2rZty8knn7znL5nOnTszYsQIunbtCsBll11GZmamumxE4k1amv9ZowZMmuTn7ilv\nwsSqjjQpQ0imRzazOcANzrnsvbaPwJ/MvaYy75Oo0yOHin5XIjGioABmzoR+/cp+zgFkc0SmRzaz\nQWa2CjgJeMfM3iuxbwVwDzDCzFaZWZuqHEtEJG4kJcEZZ0T8sFU6aeucew14rYx9GVV5bxERCa24\nutJWRETKFhOBH4rzDPFOvyORGFTWiJJQjjQpIfBx+BVJS0sjNzeXevXqYUUXNkgpzjlyc3NJKxoJ\nICKxYe3aiB4u6gO/SZMmrFq1ipycnKBLiWppaWk0adIk6DJEJIpFfeCnpqbSokWLoMsQEYl5MdGH\nLyIiVafAFxFJEAp8EZEEEZKpFULFzHKAlVV4i3QgkdYOTLTPC/rMiUKf+cA0d85VOPtkVAV+VZlZ\ndmXmk4gXifZ5QZ85Uegzh4e6dEREEoQCX0QkQcRb4D8edAERlmifF/SZE4U+cxjEVR++iIiULd5a\n+CIiUoa4CHwz62dmS81smZndHHQ94WZmTc1stpktNrNFZnZt0DVFipklm9l/zeztoGuJBDM73Mym\nmdk3ZrbEzE4KuqZwM7PrC/+7XmhmL5hZ3M0KaGZPmtnPZrawxLa6ZjbTzL4r/Fkn1MeN+cA3s2Rg\nMtAfaAMMTYDVtfKBPzjn2gAnAlcnwGcuci2wJOgiIuh+4F3n3LFAR+L8s5tZY2AMfmnUdkAyMCTY\nqsLiKWDv9Q1vBj5wzrUEPih8HFIxH/hAV2CZc+4H59xO4F/AwIBrCivn3Brn3ILC+3n4EGgcbFXh\nZ2ZNgAHAE0HXEglmVhvoAUwBcM7tdM5tCraqiEgBaphZClATWB1wPSHnnJsLbNhr80Dg6cL7TwPn\nhvq48RD4jYEfSzxeRQKEXxEzywAygf8EW0lE3AfcCBQEXUiEtABygKmF3VhPmFmtoIsKJ+fcT8Ak\n4H/AGmCzc+79YKuKmAbOuTWF99cCIV8FJR4CP2GZ2SHAK8B1zrlfgq4nnMzsLOBn59z8oGuJoBSg\nM/CIcy4T2EoY/syPJoX91gPxX3aNgFpmdlGwVUWe88MnQz6EMh4C/yegaYnHTQq3xTUzS8WH/T+d\nc68GXU8EdAPOMbMV+G6708zsuWBLCrtVwCrnXNFfb9PwXwDxrA+w3DmX45zbBbwKnBxwTZGyzsyO\nBCj8+XOoDxAPgf850NLMWphZNfwJnjcDrimszK/1OAVY4py7J+h6IsE5d4tzrolzLgP/b/yhcy6u\nW37OubXAj2bWunBTb2BxgCVFwv+AE82sZuF/572J8xPVJbwJXFJ4/xLgjVAfIOpXvKqIcy7fzK4B\n3sOf0X/SObco4LLCrRswHPjazL4o3PZH59z0AGuS8BgN/LOwMfMDMDLgesLKOfcfM5sGLMCPRvsv\ncXjVrZm9AJwKpJvZKmA88FfgJTMbhZ81+PyQH1dX2oqIJIZ46NIREZFKUOCLiCQIBb6ISIJQ4IuI\nJAgFvohIglDgi4gkCAW+iEiCUOCLiCSI/wckbrPsCb6u6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcdaf801128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.sparse import csr_matrix \n",
    "from pytrails.hyptrails import MarkovChain\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# sequence of observations (0 = rainy, 1 = sunny)\n",
    "observations = [0, 1, 0, 0, 1, 0, 0, 1, 1]\n",
    "\n",
    "########\n",
    "# calculate transition counts between states\n",
    "########\n",
    "transition_counts = np.zeros((2, 2))\n",
    "for i in range(len(observations) - 1):\n",
    "    transition_counts[observations[i], observations[i+1]] += 1\n",
    "transition_counts = csr_matrix(transition_counts)\n",
    "\n",
    "########\n",
    "# define hypotheses as transition probability matrices\n",
    "########\n",
    "hyp_sticky = csr_matrix([\n",
    "    [1.0, 0.0],\n",
    "    [0.0, 1.0]])\n",
    "    \n",
    "hyp_random = csr_matrix([\n",
    "    [0.5, 0.5],\n",
    "    [0.5, 0.5]])\n",
    "\n",
    "########\n",
    "# calculate the evidence for each hypothesis using a set of increasing concentration factors:\n",
    "########\n",
    "\n",
    "concentration_factors = np.array([0, 1, 2, 3, 4, 5])\n",
    "\n",
    "# scale concentration factors by number of states\n",
    "concentration_factors = concentration_factors * transition_counts.shape[0]\n",
    "\n",
    "# calculate evidences\n",
    "evidence_sticky = [MarkovChain.marginal_likelihood(transition_counts, hyp_sticky * c) for c in concentration_factors]\n",
    "evidence_random = [MarkovChain.marginal_likelihood(transition_counts, hyp_random * c) for c in concentration_factors]\n",
    "\n",
    "########\n",
    "# plot the results:\n",
    "########\n",
    "\n",
    "plt.plot(concentration_factors, evidence_sticky, 'rs--', label=\"sticky\")\n",
    "plt.plot(concentration_factors, evidence_random, 'bs--', label=\"random\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When comparing our results, we see that, our `random` hypothesis results in\n",
    "higher evidence values (i.e., marginal likelihood). \n",
    "**That is, assuming random weather explains our observations \n",
    "better than believing in sticky conditions.**"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}