mlmsm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximum-likelihood MSMs\n",
    "\n",
    "For those already familiar with the interface: The corresponding [API docs](../api/generated/deeptime.markov.msm.MaximumLikelihoodMSM.rst#deeptime.markov.msm.MaximumLikelihoodMSM).\n",
    "\n",
    "This kind of estimator tries to determine the most likely transition matrix $P\\in\\mathbb{R}^{n\\times n}$ given the [transition counts](./transition-counting.ipynb) $C\\in\\mathbb{R}^{n\\times n}$ <cite data-footcite=\"schutte1999direct\">(Schuette, 1999)</cite><cite data-footcite=\"wu2020variational\">(Wu, 2020)</cite>.\n",
    "\n",
    "The likelihood that a transition matrix $P$ belongs to a trajectory $S = (s_0,\\ldots,s_T)$ is given by\n",
    "\n",
    "$$ \\mathcal{L}(P) = \\mathbb{P}(P\\mid S) = \\prod_{i,j=1}^n P_{ij}^{C_{ij}}. $$\n",
    "\n",
    "After importing the package, the estimator can be instantiated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:35.827209Z",
     "iopub.status.busy": "2020-12-02T13:40:35.825649Z",
     "iopub.status.idle": "2020-12-02T13:40:36.735316Z",
     "shell.execute_reply": "2020-12-02T13:40:36.735667Z"
    }
   },
   "outputs": [],
   "source": [
    "import deeptime.markov as markov\n",
    "\n",
    "estimator = markov.msm.MaximumLikelihoodMSM(\n",
    "    reversible=True,\n",
    "    stationary_distribution_constraint=None\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It has two main parameters:\n",
    "\n",
    "- `reversible` determines whether the estimated MSM should have a detailed balance symmetry, i.e.,\n",
    "  $$ \\pi_i P_{ij} = \\pi_j P_{ji}\\quad\\forall i,j=1,\\ldots,n, $$\n",
    "  where $\\pi$ is the stationary probability distribution vector. It is the probability vector which does not change under $P$, i.e., $P\\pi = \\pi$ and therefore is eigenvector for eigenvalue $\\lambda_0$ = 1. A sufficient condition for existence and uniqueness of $\\pi$ is the irreducibility of the markov chain, i.e., the count matrix must be connected.\n",
    "- `stationary_distribution_constraint`: In case the stationary distribution $\\pi$ is known a priori, the search space can be restricted to stochastic matrices which have $\\pi$ as stationary distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To demonstrate the estimator, we can generate a sequence of observations from a markov state model. It has three states and a probability of $97\\%$ to stay in its current state. The remaining $3\\%$ are used as jump probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:36.740413Z",
     "iopub.status.busy": "2020-12-02T13:40:36.739947Z",
     "iopub.status.idle": "2020-12-02T13:40:36.741444Z",
     "shell.execute_reply": "2020-12-02T13:40:36.741927Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "p11 = 0.97\n",
    "p22 = 0.97\n",
    "p33 = 0.97\n",
    "P = np.array([[p11, 1 - p11, 0], [.5*(1 - p22), p22, 0.5*(1-p22)], [0, 1-p33, p33]])\n",
    "true_msm = markov.msm.MarkovStateModel(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A realization of this markov state model can be generated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:36.746473Z",
     "iopub.status.busy": "2020-12-02T13:40:36.746019Z",
     "iopub.status.idle": "2020-12-02T13:40:36.752639Z",
     "shell.execute_reply": "2020-12-02T13:40:36.753095Z"
    }
   },
   "outputs": [],
   "source": [
    "trajectory = true_msm.simulate(50000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can gather statistics about it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:36.756662Z",
     "iopub.status.busy": "2020-12-02T13:40:36.756219Z",
     "iopub.status.idle": "2020-12-02T13:40:36.760045Z",
     "shell.execute_reply": "2020-12-02T13:40:36.759587Z"
    }
   },
   "outputs": [],
   "source": [
    "counts_estimator = markov.TransitionCountEstimator(\n",
    "    lagtime=1, count_mode=\"sliding\"\n",
    ")\n",
    "counts = counts_estimator.fit(trajectory).fetch_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "as well as re-estimate an MSM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:36.765428Z",
     "iopub.status.busy": "2020-12-02T13:40:36.764872Z",
     "iopub.status.idle": "2020-12-02T13:40:36.767306Z",
     "shell.execute_reply": "2020-12-02T13:40:36.766777Z"
    }
   },
   "outputs": [],
   "source": [
    "msm = estimator.fit(counts).fetch_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated transition matrix is reasonably close to the ground truth:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:36.773435Z",
     "iopub.status.busy": "2020-12-02T13:40:36.772810Z",
     "iopub.status.idle": "2020-12-02T13:40:36.776827Z",
     "shell.execute_reply": "2020-12-02T13:40:36.777530Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated transition matrix: [[0.97174192 0.02825808 0.        ]\n",
      " [0.01436794 0.97057023 0.01506183]\n",
      " [0.         0.02829029 0.97170971]]\n",
      "Estimated stationary distribution: [0.24913759 0.48999019 0.26087221]\n"
     ]
    }
   ],
   "source": [
    "print(\"Estimated transition matrix:\", msm.transition_matrix)\n",
    "print(\"Estimated stationary distribution:\", msm.stationary_distribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now since we know the stationary distribution of the ground truth, we can restrict the estimator's optimization space to it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:36.783223Z",
     "iopub.status.busy": "2020-12-02T13:40:36.782731Z",
     "iopub.status.idle": "2020-12-02T13:40:36.785670Z",
     "shell.execute_reply": "2020-12-02T13:40:36.784925Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True stationary vector: [0.25 0.5  0.25]\n"
     ]
    }
   ],
   "source": [
    "print(\"True stationary vector:\", true_msm.stationary_distribution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:36.793240Z",
     "iopub.status.busy": "2020-12-02T13:40:36.792440Z",
     "iopub.status.idle": "2020-12-02T13:40:36.796483Z",
     "shell.execute_reply": "2020-12-02T13:40:36.795707Z"
    }
   },
   "outputs": [],
   "source": [
    "estimator.stationary_distribution_constraint = true_msm.stationary_distribution\n",
    "\n",
    "msm_statdist = estimator.fit(counts).fetch_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:36.801667Z",
     "iopub.status.busy": "2020-12-02T13:40:36.800213Z",
     "iopub.status.idle": "2020-12-02T13:40:36.804934Z",
     "shell.execute_reply": "2020-12-02T13:40:36.805488Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated transition matrix: [[0.97150033 0.02849967 0.        ]\n",
      " [0.01424983 0.97116296 0.0145872 ]\n",
      " [0.         0.0291744  0.9708256 ]]\n",
      "Estimated stationary distribution: [0.25 0.5  0.25]\n"
     ]
    }
   ],
   "source": [
    "print(\"Estimated transition matrix:\", msm_statdist.transition_matrix)\n",
    "print(\"Estimated stationary distribution:\", msm_statdist.stationary_distribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the markov state model object one also has access to relaxation timescales $t_s = -\\tau / \\mathrm{ln} | \\lambda_i |, i = 2,...,k+1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:36.810012Z",
     "iopub.status.busy": "2020-12-02T13:40:36.809442Z",
     "iopub.status.idle": "2020-12-02T13:40:37.213506Z",
     "shell.execute_reply": "2020-12-02T13:40:37.213946Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtjUlEQVR4nO3deXhV1b3/8fcXCBJQBiX6Y7CAEwoZISIKRFGrgFQZRLBqRRQuihNWBK+tUGzrABXl0Vsvao1ThYqCqKhFBYE6QEBGAQdEBXwUuE0EhELg+/vjnBxDspOcxHMSknxez3OenL33Wnt/dwjne9Zee69l7o6IiEhRdao6ABEROTwpQYiISCAlCBERCaQEISIigZQgREQkUL2qDiCWmjdv7m3btq3qMEREqo1ly5Ztd/ekoG01KkG0bduWnJycqg5DRKTaMLOvStqmS0wiIhJICUJERAIpQYiISKAa1QchNdP+/fvZvHkze/furepQRKqtBg0a0Lp1axISEqKuE7cEYWYNgIXAEeHjzHT38WY2ARgObAsX/W93nxtQvxfwMFAXeMLd74tXrHJ427x5M0cddRRt27bFzKo6HJFqx93ZsWMHmzdvpl27dlHXi2cL4j/Aue6+y8wSgMVm9kZ42xR3n1xSRTOrCzwK/BLYDCw1sznu/kksA7x0zqVs+PeGYuvbN2vPzItnxvJQ8jPs3btXyUHkZzAzjjnmGLZt21Z24ULi1gfhIbvCiwnhV7RDx3YBPnf3je6+D5gOXBLrGNOOTSOhzqHNrYQ6CaQfmx7rQ8nPpOQg8vNU5P9QXDupzayuma0AvgfmuftH4U03mtkqM/ubmTULqNoK+KbQ8ubwuqBjjDCzHDPLKW92HJk6kjp26K+gjtVhZNrIcu1HRKQmimuCcPcD7p4OtAa6mFky8FfgRCAd+Bb4S0DVoFQX2Ppw92nununumUlJgQ8DliipYRKXnHRJpBWRUCeBfif1o3li83LtR2q27777jl//+teccMIJdO7cmTPPPJNZs2ZVagybNm0iOTk5cP3f//73Cu3zoYce4scff4wsH3nkkRWOT2qmSrnN1d1zgQVAL3f/Lpw4DgKPE7qcVNRm4PhCy62BrfGIrXArQq2HmmH2x1vodt+7tBv3Ot3ue5fZH2+p8L7cnX79+pGVlcXGjRtZtmwZ06dPZ/PmzcXK5ufn/5ywK6S0BFFWPEUThEhR8byLKQnY7+65ZpYInA/cb2Yt3P3bcLH+wJqA6kuBk82sHbAFGAL8Oh5xFrQiXtzwoloPNcDsj7dw58ur2bP/AABbcvdw58urAeiXEXiVslTvvvsu9evXZ+TIn744tGnThptuugmA7OxsXn/9dfbu3cvu3buZOXMmw4YNY+PGjTRs2JBp06aRmprKhAkTOPLII7n99tsBSE5O5rXXXgOgd+/edO/enffff59WrVrxyiuvkJiYyLJlyxg2bBgNGzake/fugfGNGzeOdevWkZ6eztVXX02zZs0Oiefuu+9m8uTJkWPdeOONZGZm8sMPP7B161Z69uxJ8+bNmT9/PgB33XUXr732GomJibzyyiscd9xx5f6dSc0RzxZEC2C+ma0i9IE/z91fAx4ws9Xh9T2B0QBm1tLM5gK4ez5wI/AWsA74h7uvjVegI1NH0um4Tmo91ACT3toQSQ4F9uw/wKS3it+tFo21a9fSqVOnUst88MEHPP3007z77ruMHz+ejIwMVq1axZ///Gd+85vflHmMzz77jFGjRrF27VqaNm3KSy+9BMA111zD1KlT+eCDD0qse99999GjRw9WrFjB6NGji8VTkptvvpmWLVsyf/78SHLYvXs3Xbt2ZeXKlWRlZfH444+XGbvUbHFrQbj7KiAjYP1VJZTfCvQptDwXKPZ8REwtfghadSKpXRbZvbJD675cCFuWQ/db43poiY+tuXvKtb68Ro0axeLFi6lfvz5Lly4F4Je//CVHH300AIsXL458wJ977rns2LGDvLy8UvfZrl070tPTAejcuTObNm0iLy+P3Nxczj77bACuuuoq3njjjVL28pPC8ZRH/fr16du3bySOefPmlXsfUrPU7qE2WnWCF4eGkgKEfr44NLReqqWWTRPLtb4sHTt2ZPny5ZHlRx99lHfeeeeQ+8kbNWoUee9e/F4KM6NevXocPHgwsq7wU+FHHHFE5H3dunXJz8/H3St8a2/heEo7blEJCQmRYxbEIbVb7U4Q7bJgUHYoKbz7p9DPQdmh9VItjbmwPYkJdQ9Zl5hQlzEXtq/Q/s4991z27t3LX//618i60jp2s7KyeP755wFYsGABzZs3p3HjxrRt2zaSaJYvX86XX35Z6nGbNm1KkyZNWLx4MUBkn0UdddRR7Ny5s8T9tGnThk8++YT//Oc/5OXl8c4770RdV0RjMbXLgsxrYeEDkHWHkkM1V9ARPemtDWzN3UPLpomMubB9hTqoIfTtf/bs2YwePZoHHniApKQkGjVqxP333x9YfsKECVxzzTWkpqbSsGFDnn76aQAGDhzIM888Q3p6OqeffjqnnHJKmcd+6qmnIp3UF154YWCZ1NRU6tWrR1paGkOHDqVZs0MfKzr++OO57LLLSE1N5eSTTyYj46erviNGjKB37960aNEi0g8hUpgFNYmrq8zMTC/3hEEFl5Uyr4WcJ9WCOAytW7eO0047rarDEKn2gv4vmdkyd88MKl+7LzEVJIdB2XDuXT9dbirokxARqcVqd4LYsvzQFkNBn8SW5aXVEhGpFWp3H0TQraztsnSJSUSE2t6CEBGREilBiIhIICUIEREJpAQhUgYz46qrfhohJj8/n6SkpMiwFNnZ2dx4442l7iM7OxszO+RBtVmzZmFmzJwZmr3wtddeIyMjg7S0NDp06MD//u//AqFnK8yMzz//PFJ3ypQpmBnlvq0buPvuu3n77beB2A35He2w42WVO+eccyLn1KdPH3Jzc0ssW1Wj0WZnZ7N1a/kHl37sscd45plnSi2zYsUK5s6N7whD5aEEITXL4oeK36b85cLQ+gpq1KgRa9asYc+e0HhO8+bNo1Wr8j94l5KSwgsvvBBZnj59OmlpaQDs37+fESNG8Oqrr7Jy5Uo+/vhjzjnnnEPqTp8+PbI8c+ZMOnToUKHzmThxIueffz4Quw/ZWCWIwubOnUvTpk1L3H44JogDBw4ErgcYOXJkmYM3KkGIxFOcxtfq3bs3r7/+OgAvvPACl19+ebn30aNHD5YsWcL+/fvZtWsXn3/+eWSQvp07d5Kfn88xxxwDhMZnat/+p+FB+vXrxyuvvALAxo0badKkCUETZC1ZsoQBAwYARIYN37dvH3v37uWEE04AYOjQocycOZOpU6dGhvzu2bNnZB933XUXaWlpdO3ale+++67YMd577z3S09NJT08nIyODnTt3Mm7cOBYtWkR6ejpTpkxh06ZN9OjRg06dOtGpUyfef/99gGLl9uzZw5AhQ0hNTWXw4MGRJAzQtm1btm/fzu7du7noootIS0sjOTmZGTNmlBh7gaVLl3LWWWeRlpZGly5d2LlzJ3v37uWaa64hJSWFjIyMyNPj2dnZDBgwgF69enHyySdzxx13AKEP+6FDh5KcnExKSgpTpkxh5syZ5OTkcMUVV5Cens6ePXto27YtEydOpHv37rz44os8/vjjnH766aSlpTFw4MBIEpswYQKTJ08GQi2lsWPH0qVLF0455RQWLVrEvn37uPvuu5kxYwbp6enMmDGjzL+peKvdt7lKzVN4fK0YPh0/ZMgQJk6cSN++fVm1ahXDhg1j0aJF5dqHmXH++efz1ltvkZeXx8UXXxwZk+noo4/m4osvpk2bNpx33nn07duXyy+/nDp1Qt/hGjduzPHHH8+aNWt45ZVXGDx4ME899VSxY3Tq1ImPP/4YgEWLFpGcnMzSpUvJz8/njDPOOKTszTffzIMPPsj8+fNp3jw0D0rBkN9/+tOfuOOOO3j88cf53e9+d0i9yZMn8+ijj9KtWzd27dpFgwYNuO+++w6Zd+LHH39k3rx5NGjQgM8++4zLL7+cnJycYuUefPBBGjZsyKpVq1i1alXg0OpvvvkmLVu2jCTovLw8mjRpUiz2Avv27WPw4MHMmDGD008/nR9++IHExEQefvhhAFavXs369eu54IIL+PTTT4HQN/ePP/44kphvuukmvv/+e7Zs2cKaNaEpa3Jzc2natCmPPPIIkydPJjPzp4ePGzRoEBk3a8eOHQwfPhyA3/3udzz55JOR+UMKy8/PZ8mSJcydO5c//OEPvP3220ycOJGcnBweeeSRYuWrQq1vQcRy9jE5TBQeXyvz2pg815KamsqmTZt44YUX6NOnT9kVSjBkyBCmT5/O9OnTi7VCnnjiCd555x26dOnC5MmTGTZsWGDd2bNn079//8D916tXj5NOOol169axZMkSbrvtNhYuXMiiRYvo0aNHmfEVHfJ706ZNxcp069aN2267jalTp5Kbm0u9esW/Z+7fv5/hw4eTkpLCoEGD+OSTTwKPt3DhQq688kog9DtOTU0tViYlJYW3336bsWPHsmjRIpo0aVLqOWzYsIEWLVpw+umnA6HkWq9ePRYvXhzpSzr11FNp06ZNJEGcd955NGnShAYNGtChQwe++uorTjjhBDZu3MhNN93Em2++SePGjUs85uDBgyPv16xZQ48ePUhJSeH5559n7drgqWwKWnol/Z4PB7U6QRTMPrYldw/OT7OPKUlUc18uDLUcsu4I/YzR0CkXX3wxt99+e1SXlx599NHIZZjC16u7dOnCmjVr2L59e+CAfSkpKYwePZp58+ZF5pUo8Ktf/Ypnn32WX/ziF6V+WPXo0YM33niDhIQEzj//fBYvXszixYvJyio7UUYz5Pe4ceN44okn2LNnD127dmX9+vXFykyZMoXjjjuOlStXkpOTw759+0o8ZlnDmp9yyiksW7aMlJQU7rzzTiZOnFhq+ZKGSi9t3LmgIdebNWvGypUrOeecc3j00Ue57rrrSqxfeIj1oUOH8sgjj7B69WrGjx9f4hDrBcc8nIdWr9UJItazj8lhII7jaw0bNoy7776blJSUMsuOGjWKFStWsGLFClq2bHnItnvvvZc///nPh6zbtWsXCxYsiCyvWLGCNm3aHFImMTGR+++/n7vuuqvUY2dlZfHQQw9x5plnkpSUxI4dO1i/fj0dO3YsVrYiQ35/8cUXpKSkMHbsWDIzM1m/fn2x/eTl5dGiRQvq1KnDs88+G+m8LVqu8PDoa9asYdWqVcWOt3XrVho2bMiVV17J7bffHhk2vaTYTz31VLZu3RqZ0Kmgf6fwsT799FO+/vrrQ/p5itq+fTsHDx5k4MCB3HPPPWUet8DOnTtp0aIF+/fvL3GY9pIcbkOwx3NO6gbAQuCI8HFmuvt4M5sE/ArYB3wBXOPuuQH1NwE7gQNAfkmjDf4c8Z59TKpAaeNr/cxLTa1bt+aWW24J3Jadnc3s2bMjyx9++CGtW7cOLNu7d+9i69ydBx54gP/6r/8iMTGRRo0akZ2dXazckCFDyozzjDPO4Lvvvou0GFJTUzn22GMDv1VXZMjvhx56iPnz51O3bl06dOhA7969qVOnziHDjt9www0MHDiQF198kZ49e0a+YRcdnvz666+PDI+enp5Oly5dih1v9erVjBkzhjp16pCQkBCZm6Ok2OvXr8+MGTO46aab2LNnD4mJibz99tvccMMNjBw5kpSUFOrVq0d2dvYhLYeitmzZwjXXXBOZcOnee+8FQi2EkSNHkpiYGDgd7D333MMZZ5xBmzZtSElJKdcHfs+ePbnvvvtIT0/nzjvvPOTSVVWI23DfFvprbOTuu8wsAVgM3AI0Bt5193wzux/A3ccG1N8EZLr79miPWd7hvrvd9y5bApJBq6aJ/GvcuVHvR+JLw32LxMZhM9y3h+wKLyaEX+7u/3T3ggtuHwLBX7MqQaxnHxMRqUni2gdhZnXNbAXwPTDP3T8qUmQYUNJM7A7808yWmdmIUo4xwsxyzCyn8DzB0eiX0Yp7B6TQqmkiRqjlcO+AlArPPiYiUpPE9TkIdz8ApJtZU2CWmSW7+xoAM7sLyAdK6sXp5u5bzexYYJ6ZrXf3Yj2N7j4NmAahS0zljbFfRislBBGRAJVyF1O4E3oB0AvAzK4G+gJXeAmdIO6+Nfzze2AWULz3SkRE4iZuCcLMksItB8wsETgfWG9mvYCxwMXuHjiQipk1MrOjCt4DFwBr4hWriIgUF89LTC2Ap82sLqFE9A93f83MPid06+u88G13H7r7SDNrCTzh7n2A4whdkiqI8e/u/mYcYxURkSLieRfTKnfPcPdUd09294nh9Se5+/Hunh5+jQyv3xpODrj7RndPC786uvuf4hWnSFk03Hf5lWfU1gJFHx6sSLnCA+IVPs8gs2fPLnEIkHhasGBBZPDC8sjJyeHmm28utUxubi7/8z//U9HQiqnVT1JLzXPpnEtJeTql2OvSOZdWeJ8a7rv8qipBFFb4PIMcjgmitCE3MjMzmTp1aqn7VoIQKUXasWkk1Ek4ZF1CnQTSj03/Wfut7cN9f/XVV5x33nmkpqZy3nnn8fXXXx+yrwIFLZCiw3oX9u2335KVlUV6ejrJycksWrSIcePGsWfPHtLT07niiisi59y5c2c6duzItGnTIvstWu5Pf/oT7du35/zzz2fDhp+GySkc27hx4+jQoQOpqancfvvtvP/++8yZM4cxY8aQnp7OF198cUiM3333Hf379yctLY20tLTIB/qDDz5IcnIyycnJPPTQQ0AoGZ522mkMHz6cjh07csEFF0S+TEydOjVy3CFDhrBp0yYee+wxpkyZQnp6OosWLWLo0KHcdttt9OzZk7Fjx7JkyRLOOussMjIyOOussyLntGDBgkirdcKECQwbNoxzzjmHE044IZI4xo0bxxdffEF6ejpjxowp9vdRbu5eY16dO3d2qXk++eSTqMt+v/t77/xsZ0/OTo68Oj/b2bf9uK3Cx2/UqJGvXLnSBw4c6Hv27PG0tDSfP3++X3TRRe7u/tRTT/moUaNK3UdBmdGjR/urr77qzz33nE+YMMGvvvpqf/HFF93d/dprr/WkpCQfMmSIP/fcc37gwAF3dx8/frxPmjTJ+/fv76tXr/Y//vGPnp2d7WeffbYvXbr0kOPs37/f27Zt6+7uv/3tbz0zM9MXL17sCxYs8CFDhri7H3LMNm3a+LZtP/1uAJ8zZ467u48ZM8bvueced3fv27evZ2dnu7v7k08+6ZdcckmxfRX8rtz9kN9PUZMnT/Y//vGP7u6en5/vP/zwwyF1C+zYscPd3X/88Ufv2LGjb9++vVi5nJwcT05O9t27d3teXp6feOKJPmnSpENi27Fjh59yyil+8OBBd3f/97//HRh7YZdddplPmTIlEmNubm7kWLt27fKdO3d6hw4dfPny5f7ll1963bp1/eOPP3Z390GDBvmzzz7r7u4tWrTwvXv3HnLcgn/PAldffbVfdNFFnp+f7+7ueXl5vn//fnd3nzdvng8YMKDY73T8+PF+5pln+t69e33btm1+9NFH+759+/zLL7/0jh07Bp6Te/D/JSDHS/hMVQtCapSkhklcctIlkVZEQp0E+p3Uj+aJzcuoWbraPtz3Bx98wK9//WsArrrqqsjcBxVx+umn89RTTzFhwgRWr17NUUcdFVhu6tSpkZbMN998w2effVaszKJFi+jfvz8NGzakcePGXHzxxcXKNG7cmAYNGnDdddfx8ssv07BhwzJjfPfdd7n++uuB0GirTZo0YfHixfTv359GjRpx5JFHMmDAgMicIO3atYu0Bgv/3lJTU7niiit47rnnAodFLzBo0CDq1g2N6pCXl8egQYNITk5m9OjRJQ4XftFFF3HEEUfQvHlzjj322MDJnX4uJQipcUamjqSOhf6061gdRqaNjMl+Ndz3TwrK1KtXLzKYnbuXOqx3gaysLBYuXEirVq246qqrAudpXrBgAW+//TYffPABK1euJCMjo8Rhs8saLrxevXosWbKEgQMHMnv2bHr16lVmjEG8nMOFA7z++uuMGjWKZcuW0blz5xJ/n4WHC//9739Pz549WbNmDa+++mqZw4UXPWYsKUFIjVPQijAsJq2HArV5uO+zzjor0kn+/PPP0717dyA0LeiyZcuAUJ/H/v37y9zvV199xbHHHsvw4cO59tprI8NoJyQkROrn5eXRrFkzGjZsyPr16/nwww8j9QuXy8rKYtasWezZs4edO3fy6quvFjverl27yMvLo0+fPjz00EOsWLGizBjPO++8yKixBw4c4IcffiArK4vZs2fz448/snv3bmbNmlVqq+zgwYN888039OzZkwceeIDc3Fx27dpV5u88Ly8vchNE0Ii+pYn1cOFKEFIjjUwdSafjOsWs9QBlD/fdunXryGvz5s0l7qd3797F5lH28HDf7du3Jz09nfHjx5c43HfQtJyFBQ33nZqaWupw30HzOhc2depUnnrqKVJTU3n22Wcj03cOHz6c9957jy5duvDRRx8FDutdtJN6wYIFkfmsX3rppcjvdMSIEZFLMr169SI/P5/U1FR+//vf07Vr10NiLijXqVMnBg8eTHp6OgMHDgz8wN65cyd9+/YlNTWVs88+OxLPkCFDmDRpEhkZGcU6qR9++GHmz59PSkoKnTt3Zu3atXTq1ImhQ4fSpUsXzjjjDK677joyMjJK/J0dOHCAK6+8MjIH9ujRo2natCm/+tWvmDVrVqSTuqg77riDO++8k27dukXm0YjWMcccQ7du3UhOTo5JJ3XchvuuCuUd7luqBw33LRIbh81w3yIiUr0pQYiISCAlCKkWatKlUJGqUJH/Q0oQcthr0KABO3bsUJIQqSB3Z8eOHTRo0KBc9eI6YZBILBTcFVTeGQNF5CcNGjSgdevyzfCsBCGHvYSEBNq1a1fVYYjUOrrEJCIigZQgREQkkBKEiIgEUoIQEZFAcUsQZtbAzJaY2UozW2tmfwivP9rM5pnZZ+GfzUqo38vMNpjZ52Y2Ll5xiohIsHi2IP4DnOvuaUA60MvMugLjgHfc/WTgnfDyIcysLvAo0BvoAFxuZhWbX1FERCokbgkiPFnRrvBiQvjlwCXA0+H1TwP9Aqp3AT53943uvg+YHq4nIiKVJK59EGZW18xWAN8D89z9I+A4d/8WIPzz2ICqrYBvCi1vDq8LOsYIM8sxsxw9SCUiEjtxTRDufsDd04HWQBczS46yatAUUYHjLLj7NHfPdPfMoEncRUSkYirlLiZ3zwUWAL2A78ysBUD45/cBVTYDxxdabg1sDSgnIiJxEs+7mJLMrGn4fSJwPrAemANcHS52NfBKQPWlwMlm1s7M6gNDwvVERKSSxHMsphbA0+E7kuoA/3D318zsA+AfZnYt8DUwCMDMWgJPuHsfd883sxuBt4C6wN/cfW0cYxURkSI05aiISC2mKUdFRKTclCBERCSQEoSIiARSghARkUBKECIiEkgJQkREAilBiIhIICUIEREJpAQhIiKBlCBERCSQEoSIiARSghARkUBKECIiEkgJQkREAilBiIhIICUIEREJpAQhIiKB4jblqJkdDzwD/D/gIDDN3R82sxlA+3CxpkCuu6cH1N8E7AQOAPklzXgkIiLxEc85qfOB37r7cjM7ClhmZvPcfXBBATP7C5BXyj56uvv2OMYoIiIliFuCcPdvgW/D73ea2TqgFfAJgJkZcBlwbrxiEBGRiiuzD8LMjjOzJ83sjfByBzO7tjwHMbO2QAbwUaHVPYDv3P2zEqo58E8zW2ZmI0rZ9wgzyzGznG3btpUnLBERKUU0ndTZwFtAy/Dyp8Ct0R7AzI4EXgJudfcfCm26HHihlKrd3L0T0BsYZWZZQYXcfZq7Z7p7ZlJSUrRhiYhIGaJJEM3d/R+EOppx93xCHcdlMrMEQsnheXd/udD6esAAYEZJdd19a/jn98AsoEs0xxQRkdiIJkHsNrNjCF3ywcy6UnrHMuFyBjwJrHP3B4tsPh9Y7+6bS6jbKNyxjZk1Ai4A1kQRq4iIxEg0ndS3AXOAE83sX0AScGkU9boBVwGrzWxFeN1/u/tcYAhFLi+ZWUvgCXfvAxwHzArlGOoBf3f3N6M4poiIxIi5e9mFQpeE2gMGbHD3/fEOrCIyMzM9JyenqsMQEak2zGxZSc+ZldiCMLMBJWw6xcwo3KcgIiI1T2mXmH5VyjYHlCBERGqwEhOEu19TmYGIiMjhJaonqc3sIqAj0KBgnbtPjFdQIiJS9aJ5kvoxYDBwE6FO6kFAmzjHJSIiVSya5yDOcvffAP929z8AZwLHxzcsERGpatEkiD3hnz+Gn1XYD7SLX0giInI4iKYP4jUzawpMApYTuoPpiXgGJSIiVa/MBOHu94TfvmRmrwEN3L3MoTZERKR6i6aTelS4BYG7/weoY2Y3xDswERGpWtH0QQx399yCBXf/NzA8bhGJiMhhIZoEUSc8MisAZlYXqB+/kERE5HAQTSf1W8A/ws9DODAS0MiqIiI1XDQJYiwwArie0INy/0R3MYmI1HjR3MV0EHgMeMzMjgZau3tUM8qJiEj1Fc1dTAvMrHE4OawAnjKzojPEiYhIDRNNJ3UTd/+B0BzST7l7Z0JThoqISA0WTYKoZ2YtgMuA16LdsZkdb2bzzWydma01s1vC6yeY2RYzWxF+9Smhfi8z22Bmn5vZuGiPKyIisRFNJ/VEQncy/cvdl5rZCcBnUdTLB37r7svN7ChgmZnNC2+b4u6TS6oYvpX2UeCXwGZgqZnNcfdPojiuiIjEQDSd1C8CLxZa3ggMjKLet8C34fc7zWwd0CrKuLoAn4ePhZlNBy4BlCBERCpJNJ3Up5jZO2a2Jrycama/K89BzKwtkAF8FF51o5mtMrO/mVmzgCqtgG8KLW8m+uQiIiIxEE0fxOPAnYSG+cbdVwFDoj2AmR0JvATcGu7s/itwIpBOqIXxl6BqAeu8hP2PMLMcM8vZtm1btGGJiEgZokkQDd19SZF1+dHs3MwSCCWH5939ZQB3/87dD4Sfr3ic0OWkojZz6KRErYGtQcdw92nununumUlJSdGEJSIiUYgmQWw3sxMJf4M3s0sJ9y2UJjx+05PAOnd/sND6FoWK9QfWBFRfCpxsZu3MrD6hFsucKGIVEZEYieYuplHANOBUM9sCfAlcGUW9bsBVwGozWxFe99/A5WaWTijhbAL+CyA8W90T7t7H3fPN7EZCd0/VBf7m7mujPSkREfn5zD3w0n7xgmaNgDruvjO+IVVcZmam5+TkVHUYIiLVhpktc/fMoG1ltiDCkwX9BmhL6KE5ANz95tiFKCIih5toLjHNBT4EVgMH4xuOiIgcLqJJEA3c/ba4RyIiIoeVaO5ietbMhptZCzM7uuAV98hERKRKRdOC2AdMAu7ip4fVHDghXkGJiEjViyZB3Aac5O7b4x2MiIgcPqK5xLQW+DHegYiIyOElmhbEAWCFmc0H/lOwUre5iojUbNEkiNnhl4iI1CLRzAfxdGUEIiIih5cSE4SZ/cPdLzOz1RQfatvdPS2+oYmISFUqrQVxS/jnOmBMofUGPBC3iERE5LBQYoIITxkKoVtcvyq8zcxOjWtUIiJS5Uq7xHQ9cANwgpmtKrTpKOBf8Q5MRESqVmmXmP4OvAHcC4wrtH6nu/9fXKMSEZEqV9olpjwgD7i88sIREZHDRTRPUouISC2kBCEiIoHiliDM7Hgzm29m68xsrZndEl4/yczWm9kqM5sVnrEuqP4mM1ttZivMTPOIiohUsni2IPKB37r7aUBXYJSZdQDmAcnungp8CtxZyj56unt6SfOliohI/MQtQbj7t+6+PPx+J6EH7lq5+z/dPT9c7EOgdbxiEBGRiquUPggzawtkAB8V2TSM0K20QRz4p5ktM7MRpex7hJnlmFnOtm3bYhKviIhUQoIwsyOBl4Bb3f2HQuvvInQZ6vkSqnZz905Ab0KXp7KCCrn7NHfPdPfMpKSkGEcvIlJ7xTVBmFkCoeTwvLu/XGj91UBf4Ap3LzoQIADuvjX883tgFtAlnrGKiMih4nkXkwFPAuvc/cFC63sBY4GL3T1wpjoza2RmRxW8By4A1sQrVhERKS6eLYhuwFXAueFbVVeYWR/gEULjOc0Lr3sMwMxamtnccN3jgMVmthJYArzu7m/GMVYRESkimhnlKsTdFxMaGryouQHrCi4p9Qm/3whovgkRkSqkJ6lFRCSQEoSIiARSghARkUBKECIiEkgJQkREAilBiIhIICUIEREJpAQhIiKBlCBERCSQEoSIiARSghARkUBKECIiEkgJQkREAilBiIhIICUIEREJpAQhIiKBlCBERCRQPOekPt7M5pvZOjNba2a3hNcfbWbzzOyz8M9mJdTvZWYbzOxzMxsXrzhFRCRYPFsQ+cBv3f00oCswysw6AOOAd9z9ZOCd8PIhzKwu8CjQG+gAXB6uKyIilSRuCcLdv3X35eH3O4F1QCvgEuDpcLGngX4B1bsAn7v7RnffB0wP1xMRkUpSKX0QZtYWyAA+Ao5z928hlESAYwOqtAK+KbS8ObwuaN8jzCzHzHK2bdsW07hFRGqzuCcIMzsSeAm41d1/iLZawDoPKuju09w9090zk5KSKhqmiIgUEdcEYWYJhJLD8+7+cnj1d2bWIry9BfB9QNXNwPGFllsDW+MZq4iIHCqedzEZ8CSwzt0fLLRpDnB1+P3VwCsB1ZcCJ5tZOzOrDwwJ1xMRkUoSzxZEN+Aq4FwzWxF+9QHuA35pZp8BvwwvY2YtzWwugLvnAzcCbxHq3P6Hu6+NY6wiIlJEvXjt2N0XE9yXAHBeQPmtQJ9Cy3OBufGJTkREyqInqUVEJJAShIiIBFKCEBGRQEoQIiISSAlCREQCKUGIiEggJQgREQmkBCEiIoGUIEREJJAShIiIBFKCEBGRQEoQIiISSAlCREQCKUGIiEggJQgREQmkBCEiIoGUIEREJJAShIiIBIrblKNm9jegL/C9uyeH180A2oeLNAVy3T09oO4mYCdwAMh398x4xSkiIsHiliCAbOAR4JmCFe4+uOC9mf0FyCulfk933x636EREpFRxSxDuvtDM2gZtMzMDLgPOjdfxRUTk56mqPogewHfu/lkJ2x34p5ktM7MRpe3IzEaYWY6Z5Wzbti3mgYqI1FZVlSAuB14oZXs3d+8E9AZGmVlWSQXdfZq7Z7p7ZlJSUqzjFBGptSo9QZhZPWAAMKOkMu6+Nfzze2AW0KVyohMRkQJV0YI4H1jv7puDNppZIzM7quA9cAGwphLjExER4nub6wvAOUBzM9sMjHf3J4EhFLm8ZGYtgSfcvQ9wHDAr1I9NPeDv7v5mvOIUEamOLn3hbDbs+79i69vXP5qZl78Xk2PE8y6my0tYPzRg3VagT/j9RiAtXnGJiNQEacd0ZOPWhewPfZkGIMGd9GM6xuwYepJaRKQaGtn9D9SpW/+QdXXq1mdkj4kxO4YShIhINZTUMIlLTu5PQvhjPIE69Dt5AM0Tm8fsGEoQIiLV1Mim6dTxAwDU8QOMbBbbq/NKECIi1dGXC0macwuXtMzCMPq1zKL5K7fAlwtjdgglCBGR6mjLchiUzcjuf6DTcZ1CfQ+DskPrY0QJQkSkOup+K7NzT6Tf1NW8t2AQlzy8itm5J0L3W2N2iHiO5ioiInEy++Mt3PnyavbsD/VBbMndw50vrwagX0armBxDLQgRkWpo0lsbIsmhwJ79B5j01oaYHUMJQkSkGtqau6dc6ytCCUJEpBpq2TSxXOsrQglCRKQaGnNhexIT6h6yLjGhLmMubF9CjfJTJ7WISDVU0BE96a0NbM3dQ8umiYy5sH3MOqhBCUJEpNrql9EqpgmhKF1iEhGRQEoQIiISSAlCREQCKUGIiEggJQgREQlk7l7VMcSMmW0Dvqpg9ebA9hiGUx3onGu+2na+oHMurzbunhS0oUYliJ/DzHLcPbOq46hMOuear7adL+icY0mXmEREJJAShIiIBFKC+Mm0qg6gCuica77adr6gc44Z9UGIiEggtSBERCSQEoSIiASqVQnCzHqZ2QYz+9zMxgVsNzObGt6+ysw6VUWcsRTFOV8RPtdVZva+maVVRZyxVNY5Fyp3upkdMLNLKzO+eIjmnM3sHDNbYWZrzey9yo4x1qL4225iZq+a2crwOV9TFXHGipn9zcy+N7M1JWyP/eeXu9eKF1AX+AI4AagPrAQ6FCnTB3gDMKAr8FFVx10J53wW0Cz8vndtOOdC5d4F5gKXVnXclfDv3BT4BPhFePnYqo67Es75v4H7w++TgP8D6ld17D/jnLOATsCaErbH/POrNrUgugCfu/tGd98HTAcuKVLmEuAZD/kQaGpmLSo70Bgq85zd/X13/3d48UOgdSXHGGvR/DsD3AS8BHxfmcHFSTTn/GvgZXf/GsDdq/t5R3PODhxlZgYcSShB5FdumLHj7gsJnUNJYv75VZsSRCvgm0LLm8PrylumOinv+VxL6BtIdVbmOZtZK6A/8FglxhVP0fw7nwI0M7MFZrbMzH5TadHFRzTn/AhwGrAVWA3c4u4HKye8KhHzz6/aNKOcBawreo9vNGWqk6jPx8x6EkoQ3eMaUfxFc84PAWPd/UDoy2W1F8051wM6A+cBicAHZvahu38a7+DiJJpzvhBYAZwLnAjMM7NF7v5DnGOrKjH//KpNCWIzcHyh5daEvlmUt0x1EtX5mFkq8ATQ2913VFJs8RLNOWcC08PJoTnQx8zy3X12pUQYe9H+bW93993AbjNbCKQB1TVBRHPO1wD3eegC/edm9iVwKrCkckKsdDH//KpNl5iWAiebWTszqw8MAeYUKTMH+E34boCuQJ67f1vZgcZQmedsZr8AXgauqsbfJgsr85zdvZ27t3X3tsBM4IZqnBwgur/tV4AeZlbPzBoCZwDrKjnOWIrmnL8m1GLCzI4D2gMbKzXKyhXzz69a04Jw93wzuxF4i9AdEH9z97VmNjK8/TFCd7T0AT4HfiT0DaTaivKc7waOAf4n/I0636vxSJhRnnONEs05u/s6M3sTWAUcBJ5w98DbJauDKP+d7wGyzWw1ocsvY9292g4DbmYvAOcAzc1sMzAeSID4fX5pqA0REQlUmy4xiYhIOShBiIhIICUIEREJpAQhIiKBlCBERCSQEoRIjJjZreFnDGJSTqSq6TZXkRgxs01AZln32kdbTqSqqQUhUgFm1sjMXg/PNbDGzMYDLYH5ZjY/XOavZpYTnovgD+F1NweUu8DMPjCz5Wb2opkdWVXnJVKYWhAiFWBmA4Fe7j48vNyE0JwEkZaBmR3t7v9nZnWBd4Cb3X1V4RaEmTUnNNRJb3ffbWZjgSPcfWJVnJdIYWpBiFTMauB8M7vfzHq4e15AmcvMbDnwMdAR6BBQpmt4/b/MbAVwNdAmTjGLlEutGYtJJJbc/VMz60xo7Jt7zeyfhbebWTvgduB0d/+3mWUDDQJ2ZcA8d7883jGLlJdaECIVYGYtgR/d/TlgMqGpIHcCR4WLNAZ2A3nhkUR7F6peuNyHQDczOym834ZmdkolnIJImdSCEKmYFGCSmR0E9gPXA2cCb5jZt+7e08w+BtYSGmL6X4XqTitSbijwgpkdEd7+O6rvPA1Sg6iTWkREAukSk4iIBFKCEBGRQEoQIiISSAlCREQCKUGIiEggJQgREQmkBCEiIoH+PypH79P+7/2YAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(true_msm.timescales(), 'o', label=\"Ground truth\")\n",
    "plt.plot(msm_statdist.timescales(), 'x', label=\"ML-MSM with statdist constraint\")\n",
    "plt.plot(msm.timescales(), 'v', label=\"ML-MSM without statdist constraint\")\n",
    "plt.xlabel(\"state\")\n",
    "plt.ylabel(\"timescale\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian sampling\n",
    "\n",
    "While the maximum likelihood estimation is useful, it does not say anything about uncertainties. These can be obtained by estimating a Bayesian posterior with the [Bayesian MSM](../api/generated/deeptime.markov.msm.BayesianMSM.rst#deeptime.markov.msm.BayesianMSM) <cite data-footcite=\"trendelkamp2015estimation\">(Trendelkamp-Schroer, 2015)</cite>. The posterior contains the prior around which sampling was performed as well as samples over which statistics can be collected. \n",
    "\n",
    "Similarly to the maximum-likelihood estimator it has the option to make estimation `reversible` as well as to introduce a `stationary_distribution_constraint`. Furthermore, it can be specified how many samples should be drawn as well as how many Gibbs sampling should be performed for each output sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:37.217498Z",
     "iopub.status.busy": "2020-12-02T13:40:37.217043Z",
     "iopub.status.idle": "2020-12-02T13:40:37.219512Z",
     "shell.execute_reply": "2020-12-02T13:40:37.219051Z"
    }
   },
   "outputs": [],
   "source": [
    "estimator = markov.msm.BayesianMSM(\n",
    "    n_samples=100,\n",
    "    n_steps=10,\n",
    "    reversible=True,\n",
    "    stationary_distribution_constraint=None\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To fit the model, an effective counting should be used as sliding window counting leads to wrong uncertainties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:37.223078Z",
     "iopub.status.busy": "2020-12-02T13:40:37.222620Z",
     "iopub.status.idle": "2020-12-02T13:40:37.254622Z",
     "shell.execute_reply": "2020-12-02T13:40:37.254242Z"
    }
   },
   "outputs": [],
   "source": [
    "counts_estimator = markov.TransitionCountEstimator(\n",
    "    lagtime=1, count_mode=\"effective\"\n",
    ")\n",
    "counts_effective = counts_estimator.fit(trajectory).fetch_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the estimated effective counts, we are ready to estimate the posterior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:37.259213Z",
     "iopub.status.busy": "2020-12-02T13:40:37.258798Z",
     "iopub.status.idle": "2020-12-02T13:40:37.273190Z",
     "shell.execute_reply": "2020-12-02T13:40:37.273551Z"
    }
   },
   "outputs": [],
   "source": [
    "bmsm_posterior = estimator.fit(counts_effective).fetch_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prior can be accessed via the `prior` property, the samples can be accessed via the `samples` property. Both are just instances of [MarkovStateModel](../api/generated/deeptime.markov.msm.MarkovStateModel.rst#deeptime.markov.msm.MarkovStateModel) and a list thereof, respectively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:37.276581Z",
     "iopub.status.busy": "2020-12-02T13:40:37.276168Z",
     "iopub.status.idle": "2020-12-02T13:40:37.278554Z",
     "shell.execute_reply": "2020-12-02T13:40:37.278120Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# samples: 100\n"
     ]
    }
   ],
   "source": [
    "print(\"# samples:\", len(bmsm_posterior.samples))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior instance also offers a method to gather statistics over all evaluatable methods of the markov state model. Just put the name of method or property as a string, if it is a method of something that is returned by another method, just separate with a `/`. Below a few examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:37.287601Z",
     "iopub.status.busy": "2020-12-02T13:40:37.282406Z",
     "iopub.status.idle": "2020-12-02T13:40:37.381389Z",
     "shell.execute_reply": "2020-12-02T13:40:37.381716Z"
    }
   },
   "outputs": [],
   "source": [
    "stats_P = bmsm_posterior.gather_stats(\"transition_matrix\")\n",
    "stats_timescales = bmsm_posterior.gather_stats(\"timescales\")\n",
    "stats_C = bmsm_posterior.gather_stats(\"count_model/count_matrix\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If a quantity evaluation requires arguments, they can be attached to the `gather_stats` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:37.387223Z",
     "iopub.status.busy": "2020-12-02T13:40:37.386224Z",
     "iopub.status.idle": "2020-12-02T13:40:37.400456Z",
     "shell.execute_reply": "2020-12-02T13:40:37.399968Z"
    }
   },
   "outputs": [],
   "source": [
    "stats_mfpt = bmsm_posterior.gather_stats(\"mfpt\", A=0, B=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These statistics contain mean, standard deviation, as well as confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:37.415541Z",
     "iopub.status.busy": "2020-12-02T13:40:37.414841Z",
     "iopub.status.idle": "2020-12-02T13:40:37.563832Z",
     "shell.execute_reply": "2020-12-02T13:40:37.563434Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABI0UlEQVR4nO3dd1yV1R/A8c/DBsGBeyNuRQT33ltzm5qmaGrulaY5cmW50nKUqZVaucpMy8rce++9UVFEFJANFzi/Px6vPzIH4zIufN+vF68E7vOcc6725fB9zjlfTSmFEEII82OR1h0QQgiRNBLAhRDCTEkAF0IIMyUBXAghzJQEcCGEMFNWqdlYrly5lIuLS2o2KYQQZu/kyZOPlVK5X/x6qgZwFxcXTpw4kZpNCiGE2dM07c7Lvi4pFCGEMFMSwIUQwkxJABdCCDOVqjlwIcyRwWDAx8eHyMjItO6KyODs7OwoVKgQ1tbWCXq9BHAh3sDHxwcnJydcXFzQNC2tuyMyKKUUT548wcfHh2LFiiXoGkmhCPEGkZGR5MyZU4K3SFGappEzZ85E/aYnAVyIBJDgLVJDYv+dSQAXQggzZTYB/OlTiIhI614Ikfq8vb1xc3NL627g7e3NmjVrnn9+4sQJhg8fnoY9+rcGDRpkuo2CZhPAg4LA2xv8/CAmJq17I0Tm82IAr1KlCgsXLkzDHgmzCeAADg4QEgK3b+szcikmJDKLmJgYevfujbu7O507dyY8PBzQj6eYMGECNWvWpEqVKpw6dYrmzZtTvHhxli5dCoCvry/16tXDw8MDNzc39u/fD4CjoyPjxo2jcuXKNGnShGPHjtGgQQNcXV3ZsmXLf/owfvx49u/fj4eHBwsWLGDPnj20adMGgKlTp9K7d2+aNWuGi4sLv/76Kx9++CEVKlSgRYsWGAwGAE6ePEn9+vWpXLkyzZs3x9fXF4CFCxdSrlw53N3d6datGwChoaH06dOHChUq4O7uzsaNGwEYNGgQVapUoXz58kyZMuWl79c///xDzZo1qVSpEl26dCE0NPT5GIztjBkzxiR/N2nJrJYRapoexGNj4eFDCAyEfPnAzi6teyYyi2m/X+TSg2CT3rNcgaxMeav8a19z9epVvv32W2rXrk3fvn356quvngegwoULc/jwYUaNGoWXlxcHDx4kMjKS8uXLM3DgQNasWUPz5s2ZOHEisbGxz4N/WFgYDRo0YPbs2XTo0IFJkyaxfft2Ll26RO/evWnbtu2/+jBr1izmzZvHH3/8AcCePXv+9f2bN2+ye/duLl26RM2aNdm4cSNz5syhQ4cObN26ldatWzNs2DA2b95M7ty5Wb9+PRMnTuS7775j1qxZ3L59G1tbW4KCggCYMWMG2bJl4/z58wAEBgYCMHPmTJydnYmNjaVx48acO3cOd3f35/14/Pgxn3zyCTt27CBLlizMnj2b+fPnM3ToUDZt2sSVK1fQNO15O+bMrAK4kaUlODlBVJSeVnF2hpw59a8LkREVLlyY2rVrA9CzZ08WLlz4PIAbA22FChUIDQ3FyckJJycn7OzsCAoKomrVqvTt2xeDwUD79u3x8PAAwMbGhhYtWjy/1tbWFmtraypUqIC3t3ei+9iyZcvn18fGxv7r3t7e3ly9epULFy7QtGlTAGJjY8mfPz8A7u7u9OjRg/bt29O+fXsAduzYwbp1657fP0eOHABs2LCBZcuWERMTg6+vL5cuXfpXAD9y5AiXLl16/n5FR0dTs2ZNsmbNip2dHf369aN169bPf3swZ2YZwI1sbcHGBoKD9Y88efTALiu+REp500w5pby4vCz+57a2tgBYWFg8/7Px85iYGOrVq8e+ffvYunUr7777LmPHjqVXr15YW1s/v0/8a43XJVb861+8d0xMDEopypcvz+HDh/9z7datW9m3bx9btmxhxowZXLx4EaXUf8Z9+/Zt5s2bx/Hjx8mRIwdeXl7/WTetlKJp06asXbv2P+0cO3aMnTt3sm7dOhYvXsyuXbsSPc70xKxy4C9jTKvY2sKDB3Dvnj4zFyIjuXv37vPAt3btWurUqZPga+/cuUOePHno378/7733HqdOnUpSH5ycnAgJCUnStQClS5fG39//+TgMBgMXL14kLi6Oe/fu0bBhQ+bMmUNQUBChoaE0a9aMxYsXP78+MDCQ4OBgsmTJQrZs2fDz8+Ovv/76Tzs1atTg4MGD3LhxA4Dw8HCuXbtGaGgoT58+pVWrVnzxxRecOXMmyWNJL8x6Bh6fpSVkzQqRkfpDzpw59dSKpFVERlC2bFlWrVrF+++/T8mSJRk0aFCCr92zZw9z587F2toaR0dHVq9enaQ+uLu7Y2VlRcWKFfHy8sLT0zNR19vY2PDLL78wfPhwnj59SkxMDCNHjqRUqVL07NmTp0+fopRi1KhRZM+enUmTJjFkyBDc3NywtLRkypQpdOzYEU9PT8qXL4+rq+vzNEl8uXPnZuXKlXTv3p2oZ7O5Tz75BCcnJ9q1a0dkZCRKKRYsWJCk9yE90VQqLuWoUqWKSuo6zTt39Nm2VQJ+5CgF4eH66/PmBUdHSauIpLt8+TJly5ZN626ITOJl/940TTuplKry4mvNPoXyMpoGWbLo+fH798HHR9IqQoiMJ0MGcCMrKz2tYjDoq1UeP4a4uLTulRBCmEaGDuBGdnZ6GiUgQM+Ph4WldY+EECL5MkUABz2t4ugI1tb6SpX79yE6Oq17JYQQSZdpAriRMa1iXK0SECBpFSGEecowywgTy95eD9yPH+sHZeXLp68nF0IIc5HpZuDxWVjoaRVLS7h7V98I9OzMHSHSjfRynGxi9OvXj0uXLiX5+pUrV/LgwQOT3c+U4h/ildYy7Qw8Pmtr/SMiAm7d0rfkZ8umB3ghROLExsayYsWKRF9jGW/X3cqVK3Fzc6NAgQIAib5fZiEhKh57e339uL+/vuzw2aFtQqS59HCc7J49e6hXrx4dOnSgXLlyDBw4kLhnD5AcHR35+OOPqV69OocPH/5XcYW1a9dSoUIF3NzcGDdu3PP7vXiN0S+//MKJEyfo0aMHHh4eRERE/Ot+Cel3bGwsY8eOpWrVqri7u/PNN9+89r34+++/qVSpEhUrVqRx48aAfm5KrVq18PT0pFatWly9evU/70lYWBh9+/alatWqeHp6snnzZgAuXrxItWrV8PDwwN3dnevXryfp7/2NlFKp9lG5cmWVVN7eSt25o9T9+6nz4e2t1OXLSj14oFR0dJK7LTKAS5cu/evz+vXr/+djyZIlSimlwsLCXvr977//XimllL+//3++9ya3b99WgDpw4IBSSqk+ffqouXPnKqWUKlq0qPrqq6+UUkqNHDlSVahQQQUHB6tHjx6p3LlzK6WUmjdvnvrkk0+UUkrFxMSo4OBgpZRSgPrzzz+VUkq1b99eNW3aVEVHR6szZ86oihUr/qcfu3fvVra2turmzZsqJiZGNWnSRP3888/P77V+/fp/vUfHjx9X9+/fV4ULF1aPHj1SBoNBNWzYUG3atOml17z4Hh8/fvylnyek3998842aMWOGUkqpyMhIVblyZXXr1q2XvhePHj1ShQoVUrdu3VJKKfXkyROllFJPnz5VBoNBKaXU9u3bVceOHZ+/D61bt1ZKKfXRRx+pH374QSmlVGBgoCpZsqQKDQ1VQ4cOVT/++KNSSqmoqCgVHh7+0nG+zIv/3p6N+YR6SUyVFMorGNMqYWF6EQljWkW25Iu0kF6Ok61WrRqurq4AdO/enQMHDtC5c2csLS3p1KnTf15//PhxGjRoQO7cuQHo0aMH+/bto3379q+85k0S0u9//vmHc+fO8csvvwDw9OlTrl+//tL3wvibRbFixQBwdnZ+fk3v3r25fv06mqY9L0oR3z///MOWLVuYN28eAJGRkdy9e5eaNWsyc+ZMfHx86NixIyVLlkz0OBNCAvgbODjoq1X8/PTVKnnz6qkWkXm9WMggPgcHh9d+P1euXK/9/qukl+NkX9UPOzu7f+WwjdRrzlp61TVvkpB+K6VYtGgRzZs3/8/1L74X2bNnf2k1+MmTJ9OwYUM2bdqEt7c3DRo0eOn4Nm7cSOnSpf/19bJly1K9enW2bt1K8+bNWbFiBY0aNUr0WN/kjTlwTdPsNE07pmnaWU3TLmqaNu3Z16dqmnZf07Qzzz5ambx36YSFhX7OOOiHakldTpHa0sNxsqDnhW/fvk1cXBzr169/Yz+qV6/O3r17efz4MbGxsaxdu5b69eu/sZ3kHl3bvHlzvv766+ez5mvXrhEWFvbS96JmzZrs3buX27dvAxAQEADoM/CCBQsC+kPVV7WzaNGi5z+oTp8+DcCtW7dwdXVl+PDhtG3blnPnziV5LK+TkBl4FNBIKRWqaZo1cEDTNOMhvAuUUvNSpGfpkI2N/hES8v8CElmzSlpFpLz0cJwsQM2aNRk/fjznz59//kDzdfLnz89nn31Gw4YNUUrRqlUr2rVr98Z2vLy8GDhwIPb29i8tAPEm/fr1w9vbm0qVKqGUInfu3Pz2228vfS9y587NsmXL6NixI3FxceTJk4ft27fz4Ycf0rt3b+bPn//K2fPkyZMZOXIk7u7uKKVwcXHhjz/+YP369fz4449YW1uTL18+Pv7440SPISESdZyspmkOwAFgENASCE1MAE/qcbLX/UK4fMuAR0HnBB0nmxpiY/VVKra2Upczo5PjZHV79uz5V01MkTJMfpyspmmWmqadAR4B25VSR599a6imaec0TftO07Qcr7h2gKZpJzRNO+Hv75+ogRgt2X2DIev+4aM/T+AdEJqke5iasS6nUvqSw0eP9KAuhBCpJbEz8OzAJmAY4A88BhQwA8ivlOr7uuuTOgN/GhqBa8kyRGXJS476XnRuVJf3qpfC2cH2zRenAqX0TUAgdTkzIpmBi9SUYgUdlFJBwB6ghVLKTykVq5SKA5YD1ZLc4zewt7Fk0HvDsfC/js/3Q1kxezQdFm7gu6PXiTCk/dNEqcsphEgLCVmFkvvZzBtN0+yBJsAVTdPyx3tZB+BCivQQfd1n//4fcGD/TQa+/wFRV/Zxe2l/vvr9b7qu3sOWC3eJjUu90nCvYqzLGROjn3To7y9pFSFEyknIDDw/sFvTtHPAcfQc+B/AHE3Tzj/7ekNgVAr2E4Ds2Z2ZPHku+/dfZcigMawY2J18TvZMW7OJnqt3csj70WvXnaYWOzs9jRIUpAfykBA9zSKEEKZk9kWNnz4Nokq1osTZZMGxzrvUadKeYXXLUSZPNhP0OPliYvTVKlmy6Plx2/SRtheJIDlwkZoyVVHjrFmzsfybDbjkycXjLXPYNsuLd+Z8zdRtp/ENTvvTqKQup8gIpk6d+ny7eHy//fZbko559fb2Zs2aNc8/X7lyJUOHDn3jdV988cXzg7wAWrVqRVBQUKLbTwkJHYMpmX0A1zSNBg2a888/p1mwYCXZVDh+az9i2+GjdFu9l8UHLhMcmfaHfEtdTpHSXrX9PSW9LoC/rj8vBvCEejGA//nnn2TPnj3R98kozD6AG1laWvL22705sP8a33zzM7990JumpQuw/Idvab/oZ9advkV0TNo+UZS6nCKpZsyYQZkyZWjatCndu3d/Phtu0KABEyZMoH79+nz55Zfs3LkTT09PKlSoQN++fYl6thzKxcWFx48fA3DixInn53pMnTqVvn37Pj+OdeHChc/bnDlzJqVLl6ZJkyYvPUr10KFDbNmyhbFjx+Lh4cHNmzf/0x8vL6/nB0qBfhQswPjx49m/fz8eHh4sWLAAgAcPHtCiRQtKlizJhx9++J/2Fi5cyIMHD2jYsCENGzb817i8vb0pU6YM/fr1w83NjR49erBjxw5q165NyZIlOXbsGJD4419Xr16Nu7s7FStW5N133wXg999/p3r16nh6etKkSRP8/Pz+01d/f386depE1apVqVq1KgcPHgRg7969eHh44OHhgaenZ7KOC4AMeJiVvb09bdp0BmBIlcJ8P/BrArfFMv1oe9Y16cGQRpVpUjL/Sw+vSS3GtEpEhD4bz50bsmeXAhLmosHKBm98TZtSbRhTa8zz13t5eOHl4cXj8Md03tD5X6/d47Xntfc6ceIEGzdu5PTp08TExFCpUiUqV678/PtBQUHs3buXyMhISpYsyc6dOylVqhS9evXi66+/ZuTIka+9/5UrV9i9ezchISGULl2aQYMGce7cOdatW/fKNgFq1apF27ZtadOmDZ07/39Mxv6AviX+ZWbNmvWvXZ0rV67kzJkznD59GltbW0qXLs2wYcMoXLjw82uGDx/O/Pnz2b17N7ly5frPPW/cuMHPP//MsmXLqFq1KmvWrOHAgQNs2bKFTz/9lN9++42ZM2fSqFEjvvvuO4KCgqhWrRpNmjRh6dKljBgxgh49ehAdHU1sbCwXL15k5syZHDx4kFy5cj0/I6VOnTocOXIETdNYsWIFc+bM4fPPP/9XX0aMGMGoUaOoU6cOd+/epXnz5ly+fJl58+axZMkSateuTWhoKHbJ3MKd4QJ4fDlyOLNr13nmzJnEb7+t5cy5bYw+0o0qzboyvH4FPAvmTNP+2dvrq1OMdTnz5tUfdgoR34EDB2jXrh32z47BfOutt/71/a5duwJw9epVihUrRqlSpQDo3bs3S5YseWMAb926Nba2ttja2pInTx78/PzYv38/HTp0wOFZoVjjkbUJYexPYjVu3Jhs2fTFB+XKlePOnTv/CuBvUqxYMSpUqABA+fLlady4MZqm/eeY2YQe/7pr1y46d+78/IeF8ZhZHx8funbtiq+vL9HR0c+PoY1vx44d/0otBQcHExISQu3atRk9ejQ9evSgY8eOFCpUKPFvVDwZOoADFC3qypIlaxgwYDQzZozl+J7veFihDkM2hlGnWB4G1y6Di7NTmvXPmFYxGPS0Stas+ozc2jrNuiTe4E0z5te9PpdDrkRf/6aVYlme/dR/3eusrKyeV8+JjIz81/fiH0FraWn5PHed1N9Ss8SbhcRvVylF9Gtyhq/qR0K9eJTuq46ZTejxr0qpl74Hw4YNY/To0bRt25Y9e/YwderU/7wmLi6Ow4cPP/+hazR+/Hhat27Nn3/+SY0aNdixYwdlypRJ1DjjyzS/tFesWIWff97Fzh1n+W1ENwbWLM3fq+bSZdZS5uw6z5OwyDffJAVZW/8/rXLrFgQGymoVoatTpw6///47kZGRhIaGsnXr1pe+rkyZMnh7e3Pjxg0Afvjhh+dHt7q4uHDy5EkANm7c+MY269Wrx6ZNm4iIiCAkJITff//9pa9707Gv8dvdvHnz8+Ndk3pcrCmOmU3o8a+NGzdmw4YNPHnyBHj5MbOrVq16aTvNmjVj8eLFzz8/c+YMADdv3qRChQqMGzeOKlWqcOXKlSSPBTJRAAd9RlGiRBlsrSxpVSwrlrcO8uCnD1k+fQDt5v/At0evER6dtlvzpS6neFHVqlVp27YtFStWpGPHjlSpUuV5qiE+Ozs7vv/+e7p06UKFChWwsLBg4MCBAEyZMoURI0ZQt27dBBVRqFSpEl27dsXDw4NOnTpRt27dl76uW7duzJ07F09PT27evPmf7/fv35+9e/dSrVo1jh49+nx27u7ujpWVFRUrVnz+EDMhBgwYQMuWLZ8/xEysyZMnYzAYcHd3x83NjcmTJwOwfv163Nzc8PDw4MqVK/Tq1Yvy5cszceJE6tevT8WKFRk9ejSgP/jt0qULdevWfWkuHvQHridOnMDd3Z1y5co9r0/6xRdf4ObmRsWKFbG3t6dly5ZJGoeR2W/kSY6IiHCWL/+CRYtnERERThb3prg278/ARlVoXa4QVmn8VNFg0Gfk2bJBrlySVkkr6WEjT2hoKI6OjoSHh1OvXj2WLVtGpUqV0rRPImVkqo08yWFv78Dw4RM4cvgmfbwGY+NzhvzZHZm96zy91+zn4G2/NN2ab0yrhIXpq1WCgmRLfmY1YMAAPDw8qFSpEp06dZLgLYBMPgN/UUREOHZ29uy66sOwAe2xLF6Dui27Max+BcrmzZ5yDSdAXJweyG1tpS5naksPM3CRecgMPIns7R3QNA33nDaUcnYg4J+v+PuTd+g+Yx4f/3UqTbfmS13OtJUeDkkTGV9i/51JAH+J3LnzsmnTPr7/fjMFs2fBf9NMfprak85fbWLR/ktpujXfxkZPq4SE6GmVp08lrZLS7OzsePLkiQRxkaKUUjx58iRRm3vMJoXy87H9lHV2x9khdU8ZjImJYd2679j42zo8+81l27WHOFjE0qeWG53di2Jj9eYn+ilF6nKmDoPBgI+Pz3/WTwthanZ2dhQqVAjrF1YsvCqFYhYBPNwQToF5+o6l990/4L3yw3G0SZvNN2e9fejUqhKWrtUo1cKLYU1r0qRUASzScGt+VJT+4ewMOXPqhSWEEBmHWefAHawd+LHZdqrkqc2cE5Oovs6FxWdmEWZI/QLHRbNn4e32XQg/9zdn5/di5Mfj8Fq9g1M+T1K9L0a2tnp+PDhYT6sEB0taRYjMwCxm4KAfw+rvDxcCj7Hk4lR2+/yFs10uBrmPxavcEBysU/cQkRs3rjJr1gT++utXrJ1ykrfXF9RzK83g2mUoljPttubHxuqrVRwc9NUqUkBCCPNn1ikUo+hofYt5UBCcCzzCkotT2euzjVz2eRjk/iEDKozCQkvdXyqOHz/En3/9SrFWA/nhxE2e+t6mQ92a9K9ZmlxZ0i4pHRmpv185c+qpFUmrCGG+MkQAN4ofyM8GHGLxhSnEqlh+brMLgDgVl+qBHODSjeu0aFwe6/ylydukH33btOSdSq442KTNmWFK6Q85NU2fjTs66n8WQpiXDBXAjeIH8igVRk6nLDwMv0/H3+syp+5y6hZsbLK2EsJgMLBmzQrmfj6FwCf+OJSpi0vzfgxpWZ825Qun2dZ8qcsphHkz64eYr2Jjo88sixWDvDmyEBoK/sHBFHFyxSVrcQAehj0gKjYqVfpjbW1N796DOHLoJqNGfUys90mufDWAz7YepNdP+zlwK2225r+sLmds2hYnEkKYgFnPwF8Uf0ZuaalvN3/375ZcCbzAMI8JdCvdF1vL1Jt++vn5cuDATnJ5NOHrg1e4fGQbNeo2Y0QjT8ql0dZ8pfSHnJaW+tpxKSAhRPqXIVMorxI/kB95tIPFF6Zw4tEhCmQpzAjPSbxdygsbS5sU70d8l69cpEljN6ydcpG1Tg/aduzJ4NrlKJDNIVX7YWRMqzg56QUkbFL37RBCJEKmCuBGxkAeGKg4/OgfFl2Ywmn/oxRyLMoIz0l0KdUba4vUO6P1yJF9TJs+hnNnj2OTuyg5G/bl3Q6d6FOtJFnt0iaCRkTowVzqcgqRfmXKAG4UP5Af9PubRRencNb/OEWcijHSczKdSr6LlUXqrBRRSvHHH78w89Px+Pk/Iu+AFWTLlp3eVUvQyd0F2zTYmm9Mq1hZSV1OIdKjTB3AjeIH8n0Pt7L4whRuBV/jSDdvnO1St8BxdHQ0V66cJ0uBkizZf4m/fvqCYjXfYkTrBjQtnTZb840FJKQupxDpiwTweOIH8rvh1yiftzSgGLyrOx1KvEOzogmvwG0Kly+fp1XrahhiYnGs1Joqb/VjVLPqVC788nJNKc2YVsmTR68GJGkVIdJWhlxGmFTG5YeurhoVC5YmNBTuBTzietBlAiIfA2CIMxAblzpr7cqWrcDBA9d5u3NPQk9uYd+nXek1fiyjNh7k1pOkF3BNKqnLKYR5yJQz8BcZZ+QBgXFYWCgcs1iy/MIC1lxZzqhKU3jLtUuq7ey8cuUCn8wcx9lL58nb9ysi4zRalytM/xql0mRrvtTlFCLtyQz8NYwz8uKuFmTPZkloKBSwLYmGxuBd3Wi6sSJ/3PqFOBWX4n0pU8aNH3/Yyr4dZ9jYtxntyubju2nv0XrK56w4co2w6NQtwyN1OYVIvySAxxN/Z2e7sm3Y1Owc82uvwRBn4P2dXWj2qyd/3d6UKrspc+RwJpu9DW0L25E96jE+ayfx2QfdeGv2Cjadv0NMXMr/MInPwUH/8PPTS7pFRKRq80KIl5AUymsYUytPAmL5y2cdSy5M43bwdcrn9OCDSlNpVrQtWiqsFomKimL16q/5fMF0QoKDyFKuARW7fsiwBh7Udc2bKn2ILzpaP+0wRw79tMOULDQthJBVKMliDOSPA2L4894allyczp3gm6xrtSNVD8x6+jSIxYs/Y/+xI2TrNJV7QeG4583K8PoVKJcve6r1wyg8XE+n5Mmjp1nkpEMhUoYEcBOIH8j3+G6hXakOWFho/HRlOfmzFKJR4Zap0g+lFLFxih/2n2TaoNY4Vm5L27f7MaS+OwVTeWu+1OUUIuXJQ0wTMObISxa3onP5joSFaYSGxfHdhUX8fG1VqvVD0zSsLC1oVCwnNTw8Cdz9LWvGt6fN+Kl8secCTyOiU60vlpb6eSpK6UsOHz2Skw6FSC0yA08G44z80ZNoIuKCKZgjF1cDLzD+wEBGV5pC3YJNUiU/vW/fdqZOH8vVy2exyVeS4n0+p0/NcnSumLpb85X6/8PNPHn0wC5pFSGST1IoKSj+6YeHH21nwtG++Ib5UC1vHcZUmU7tAg1TvA9xcXFs2rSGw6dPoqp157C3PzlUKMOa16ZZ6YKpujVf6nIKYVpJDuCaptkB+wBbwAr4RSk1RdM0Z2A94AJ4A28rpQJfd6+MGsCN/j8jj2Lj7RUsvfQpfuEPqJm/Ph9UnkbN/PVTrS8bdu5mtFdTHMo1pHL7gXzQpj5VUnlrflSU/p44O0tdTiGSIzkBXAOyKKVCNU2zBg4AI4COQIBSapamaeOBHEqpca+7V0YP4EbGQO73JJJfbi/nm4uf8ijiIbXyN2RM5WlUz183xfsQFBTIwoUz+fa7RcQqcKrSliZdBjGyaVWK53RK8faNpC6nEMlnkhSKpmkO6AF8ELAaaKCU8tU0LT+wRylV+nXXZ5YAbmQM5A8fR7Dh9jcsvzQL/wg/ZtRcSF+3YanSh3v3vJk1exK/bfoJ62x5KDjgG9q4udKvRilyO6bekhGpyylE0iUrgGuaZgmcBEoAS5RS4zRNC1JKZY/3mkClVI6XXDsAGABQpEiRynfu3En6KMyUMZD7Pg5nw62ldCzdlQKOBbnw+DTRcdFUylM9xftw/vwpzpw/S0Ch6vxy9jZRd87wXsdO9KxSgiw2qbcTJzJSP18lZ059I5CkVYR4M1PNwLMDm4BhwIGEBPD4MtsM/EUv1uwceqADZ/yPcqSbd6qWePv1r98Z1q8tNvlLUqT5+4zo2pm25QtjZZk6q0qlLqcQiWOSdeBKqSBgD9AC8HuWOuHZfx8lv5sZW/yzVpyc4NMqP/B1vS1YW9hgiDMwdl9/zj8+leL9aNesFQsWrCRbXBg3Vo7hw6Fd6LBgFXtvPkyVc140Tc+FW1vDvXtw/77+w00IkThvDOCapuV+NvNG0zR7oAlwBdgC9H72st7A5hTqY4ZjDORupRyp5VKF0FA453uJrd4babGpMn3/ac+FJ2dSrH1LS0vefrs3hw/e4KOPPoOHV7j4/UeM//0Ygzce5uLD1y4mMhkrK30LfmSkftJhQACk8hldQpi1hKxCcQdWAZboAX+DUmq6pmk5gQ1AEeAu0EUpFfC6e2X2FMqrGFMr9x495ccbC1l5ZT5Po4No6dKB0ZWmUi6ne4q2HxDwmBu3rvHAtiDLDl7izv4NtOnUh2GNKlMoe+rkN4yrVSwtpS6nEC+SjTxmwBjI7z4K4ofrX7DyygJCDMG0LtaZ0ZWmUMbZLcX7sHXb7wzo2xZLh2zkqPMOvd4dwHs1y5HdPnVy9MYCEk5O+moVKSAhhARws2IM5HceBbL62gJWXf2CMEMoQyqO46Nqn6V4+2fOHGfKtDGcOLYPqxwFyNeoD0Pe9aKrZ7FU25ovdTmF+D8J4GbIGMi9/QJYde1zPPNVprVrR8IMofiG+VAie5kUa1spxa5dfzFl2hieGhSO3eaQz8mBATVL07xM6mzNj4vT0ypWVvpqFYfUPWhRiHRDArgZe3H54XfXZjP7xAQOvH2dolldU7Tt2NhYHj16yIMYWxb8c5Sj6+ZSsXU/xnVsQdUiqbM1X+pyiszuVQFcaqmYAeOqlRw59ED+VsG+OFnkooiTHrzXX1tJtbx1KJathMnbtrS0JH/+guQH+pe04bD3SQ5/vp+e25vRqOtwPmhVmxK5spq83fisrfWPsDAICfl/WkW25IvMTmbgZij+jDw8Lpj6mwoTHhNG55K9GOE5KUVn5U+e+PP5/On88MNSsLDCqVoH3uk3lvdrlU2VrflxcXogt7XVf6jZ26d4k0KkOUmhZEDR0XoQv/bgId9enc3a60uJjYvh7VJeDPecSGEnlxRr+/btG3zy6UfcehxEXOPRWFhAV49ivFu5OFlsUz7HIXU5RWYiATwDMwbyqw98WXFlFutvfENsXCxdS/dlhOdECjoWSbG2DQYD/uEGZm/axuYFIynYuA+jevWmfYWiqbI1X+pyisxAAngmYAzkV+7fZ/mVz9hwYzkKRc8y7zOj1sIUrQ506tRRBg/rxT3va9gWLEeZtwYz7p1O1HPNm+JViaQup8jopCZmJmBjo89Ea5QvyOeNF7Ot9Q06F3+PqNio50E0KCpltslXqlSdA3svMmvWUuwjH3N26VD6DezGwJ8Pcd43ZbfmS11OkVnJDDwDM87IAwIUVlYa18NO0vH3OqxouomGhVukWLthYaEs/eZzLvsG8LD0WwSER1ErvwMjm1ZL8a35UpdTZESyjDATMs7Is2fXCAoC++A8dHLtS6U8NQG48Pg0eR0KkNshr0nbzZLFkQ9GTwEgPDqGT1b9xA8fDuLv6p14t89w3q9XMcW25muavuEnNhYePJC6nCJjkxl4JmKckQcGgoWFov3fntx6eo3e5QYzuOKH5LLPkyLt3r59gynTx7Dzn81YOjqTp/67DHlvIO9ULpniW/OlLqfICOQhpnjOGMhP373O0ksz2OL9E7aWdvQpN5RBFcfibJcyOyyPHz/E5KmjOX/mKLaF3ajQfwHv1ypD89IFsbRIuTyH1OUU5k4CuPgPYyA/eecqSy/N4HfvNdhbOdC3/DDedx+Ds11Ok7eplOLvv3/jygN/zmcpx+WHAeSLfsRH3dpRrUhuk7cXX0yMnh93cJC6nMK8SAAXr2QM5Ce8L/PVpen8eWc9DtZZeK/8CIZ5fISDdco8eIxTiqmLvuDb2aNxKFOXBt1GMq5DkxTfmi91OYW5kQAu3sgYyI95X+Sri9M4H3CMfW9fxdbSljgVh4Vm+lWnoaEhLF4ym6XfzCfGYMDRsyWdvEYzoml18jil3D55qcspzIkEcJFgxkB+71EI2eyc0KwjabO5KgPdx9KlVK8UadPPz5dZcz7m5/XfYVugNIV7f043z5Tfmh8To+fHnZwgd2595Y4Q6Y1s5BEJZlx+WKGUE1mzwsOgpxTOUoK8DgUAfTNQSHSwSdvMmzc/Cz5fzu7dF/h63lc0KJGPlQfP03jUBDacuklMbMoUy5S6nMKcyQxcvFH85YeWljD77Ch+ub6a9yt8QN/yw3C0cUqRdmctnMui2R9inasoJdsM4iOvd2lYIn+Kbc2XupwivZIZuEgy44y8WDF9ttq8QA88ctZk9omJ1FhXjCVnZhNmCDV5u+OGjWHp0vU422pcWjme/r1b0m3B95z3fW3t7CTTND1oW1rCvXtw/77+sFOI9Epm4CLRjDPyvTeOsfjCFPb5/o2zXS6GVBxHr7KDTL5qJTo6mtU/fsPceVOxyO1Kjk5TaVA8H4Nql6FwCm7Nl7qcIr2Qh5jC5IyBfNf1wyy5OJUDvv+Qyz4Pg93H0avcQOytTFvEMiQkmEcBT9j7MIbvdx3iybHN9Og7msGNq5LDIWUWdUtdTpEeSApFmJwxtdKxak02tt/GmsYHKJWtAtOPfsDNp9dM3p6TU1aKFy1G3+ol6V0gguDjm1k6ogWN3h/CigMXiDSY/ghCCwt956alJdy9C76+klYR6YfMwIXJGGfkx70vUdq5HPb2MOPoGIpnL02PMv1N3t6VKxeYPG0Mh/Ztw9IpN0Wb92P8+4NoUaZQim3NDw/XZ+VSl1OkJpmBixRnnJE39ShH1qwQGGzg7KNT3Aq6Dujb6KNjo03WXpkybvy89m82bNhFsYIFMNy/xMwd5+izbj9H7vibrJ34HBz0Dz8/uHPn/0fXCpEWZAYuUoxefFnhHxCDnbU1R59s48P9/RnuOZGupfpgY2m6XTNxcXFERkZwyOcpc9dt5Ma276nddSSTu7enZO6U2ZovdTlFapEZuEh1NjaQN69GqeLWZM0KljFZyWNfkPEHBlJ3Qyl+urIcQ5xpEsoWFhY4OGShSakCDPbIh8WT2+yc7UXrHu0Zt/Yv/EJMP1W2sdGXVYaE6JuAnj7V15ILkVpkBi5SjXFG/vvlbSy6MIVzT45R2NGFkZUm06nku1hbmG7L/NOnQcz/cibff7eQuDhFjtpdGTx0Iu9WKY5jCmzNl7qcIiXJMkKRbhgD+W+X/mTxhSlcCDhJUSdXRlSaTKcSPbGyMF0u4v79u0z/dAKPrHPh59qUbLZW9KpclM6epbC2NP0voFFR+oezs55WkZMOhSlIABfpjjGQ/3rxDxZfmMKlwNPMrL0Er3KDTd6WUopr/sGM++JLTv32FcWb92PiwEE0LlXA5FvzpS6nMDUJ4CLdMhj0wssbL2yhTsEm5HTKwo67fxBieEr74t1NeoztkSP7GPPRUG5fO49N3uJ4dhrO1D49cS/gbLI2jGJj9SNrpS6nSC4J4CLdMxj0A7MCA2HYwY7cD7vFto6nTH4OeVxcHL9uWsO0T8YT8Og+jh4t6Th4OoNqlaFIDkeTtgVSl1MknwRwYTYMBngSEMd1Xz/yO+YnyiIQr3/a8F75EbRx7WyygB4ZGcny7xZxNzYLp61ciYgIp0mRrIxoVRdnE2/Nl7qcIjkkgAuzY5yRH7l9jtGHunEz+DJlcrgxuvJUWrp0MOnMPCA8ivcnfsjhX5fiXK0Dg4d8iFftithZm3a6LHU5RVLIOnBhdqyt9SDXspI7+7ufZ26Nn4iONTBgR2ea/1qJv71/w1QTEGcHW+aPHkWzFh15cmg9M/s1pP6gYfx25iaxcaab5FhZ6Q81DQbw9obHj/VcuRBJITNwYTYMBnj8JJYfz65lycVp3Am5gVtOTz6oPJWmRd4y2WqS8+dPMW7yKM4e34dDufrU8JrOkNplqF40t0lXrEhdTpFQkkIRGYbBAP5PYlh95ie+ujide6G38Co3hJm1F5usDaUUe/Zs41aYxl+Pbbj7wIcSVqFM7fMOpXJnM1k7IHU5xZslOYBrmlYYWA3kA+KAZUqpLzVNmwr0B4ynBk1QSv35untJABemZDDAo8cGVp35gdI53KhRuBq+YT5cCTxPw0ItTDZbNsTG4TVyIHt+XY598aq06TOO8Z1bkc/J3iT3NzIWkMidG7JnlwIS4v+SE8DzA/mVUqc0TXMCTgLtgbeBUKXUvIR2QgK4SAnxlx8uOPcR3135nMPdbpM/S0GTtREREcHXyxawePEsoiLCyOrehN6DJjC0WS2Tbs2XupziZZL8EFMp5auUOvXszyHAZcB0/2cIkUzGh52urjCpzjS+b7CTbFpBlIJPjn7I/vs7k/2w097entEjJnDi6G169BpC6IVdfL/kE7qs2s36M7cxxJqmlL3U5RSJkagcuKZpLsA+wA0YDXgBwcAJ4AOlVOBLrhkADAAoUqRI5Tt37iS700K8jnFGfuvhYzps8+RhuA/V89VlTOXp1CrQwCRt3L17m7vBkWy4/pRDJ49j9/g6Hw//gKZli5j0QafU5RRggoeYmqY5AnuBmUqpXzVNyws8BhQwAz3N0vd195AUikhNBgP4+key/OQKll/+DL+IB9TM34AxladRI389k7ShlGLo+OH89uNirHLkp0LbwXwyZBAeBXOa5P4gdTlFMgO4pmnWwB/ANqXU/Jd83wX4Qynl9rr7SAAXacEYyL85uYzllz7DP/IhtQs0YkzlaVTLVyfZ91dKsWPnn0z4eDQP7lzDpkBpGvYYy5TeXSlqwq35BoM+I8+WDXLl0lNHInNIcg5c038f/Ba4HD94P3u4adQBuGCKjgphatbWUKSAHVNbDOf4u7cY7zmfKwEX6PB7XQbs6JLs+2uaRtMmrTmy/xKfzV6GbWQgJw/voueP+5i7+zwB4VEmGIU+jqxZ9bXjt2/r9UelgETmlpBVKHWA/cB59GWEABOA7oAHegrFG3hfKeX7unvJDFykBwYD3H8UztcnvsbSQmNIpdGA4mLAWdxyeiT7/hER4TwJi2Td+Qes2fQzUbeO8/6wSQxoVN1kW/Pj4vRAbmOjp1XsTbuiUaQzspFHiBfEX364+8FmBu1rz5qW26hfqJnJ2vhs/kyWLJgCltbkq92FcaM+omOlMlhamOZBp9TlzBzkLBQhXhB/+WGzko2YWGkhntkboRRs897M+cenkt3GR6Mnsn/fFWrXb47vnh8Y/U4tWnz4MYe9H5nkHBepy5m5yQxciGeMM/KAAEXrv8pyK/gqzYu2Y3TlqSZJrZw8eYSxk0YSU7IeMSXqU6lAdobWLU+ZvNmTfW+QupwZmaRQhEgggwHuPHzKwmNfsvLqfEIMT2nl0pHRladS1rlCsu6tlCI6JpbNF+8xd9FcnpzfS4veHzL53S7kz2qa9YHx63I6O0taJSOQFIoQCWRtDSUKZ+Pzth9z8l1vBpf/mH33d9BkozsDdnThSkDSF1xpmoattRVvexRjVLPqWIc9YvOnXjRo34zpG/4iJCr52y5tbfWDsYKD9SNrg4MlrZJRyQxciDcwGOCWbwBfHl3A6mtfEh4TyluubzOv3gqyWCdvnXdYWCjzF83m2+XzMURHkb9JP8YOH0dH96JYWyZ/fiV1OTMGSaEIkUwGA9z0fcKCI59zIfAY61tvx8JC42lUENlssyfr3v7+fkz5dBKB+Ty5Y1+UPDaK/jXL0Mq9uEm25ktdTvMmAVwIEzEY9AedQUEaQQZ/mm4pwYRqs+lVbqBJ7n/kjj8jPxzCgzM7KduyH5+O+oDKRfMk+75Sl9N8SQ5cCBOxtoa8eTVcXSFHNiu6FB+AZ44GKAV3gm9x++mNZN2/RtHcfPnBCFyKFOPCz3Po9FYV3pk5n9tPQpJ1X+NJh7a28OAB+PjoM3NhvmQGLkQyxd8Q9MHhbvx99xc6l+zFCM9JFM3qmuT7KqXY8sevfDx9LI8f3CZ7rW54DZ7Ae9VL4eyQ/GR2ZKTe95w59Y1AklZJvySFIkQKMxjgio8v847MZt2NpcSpWN4u5cVwz4kUdnJJxn0NfLvqG3yyuLD3iQVa8CPali/MkFb1sbdO3hrB+HU5jWkVkf5IABcilRgMcOnefeYemcWGm8sARddSfRnuOYGCjkWSde+7gaF09+qI96nd5KnWhrGjp9C1ZsVkb82XupzpmwRwIVKZwQAX7/kw9/Bn/HxrOQDdS/fjg8pTyWWf9IeSfn6+TJzxEX/99gOatR0ujXowY+wEGpQumuwVK1KXM32SAC5EGjEY4Pzdu8w5/Clb765lb5cr5MuSH6VUsgLutWuX+GDSKE4d/IdstbrRqPsIhtYpS5k82ZLVX6nLmf5IABcijRkMcNcvmJiwrFhYKPrvbUH1/PUY4TkxWfc9eGgPN2Kzsf7yYx7eOItHbls+GdCXAtmSF3mNBSScnPRDv6SARNqRZYRCpDFrayheKCuurmCbJYpcNoXIouVEKYiOjcY/3C9J961dqwG963qyoXcDctzYzrYvhlOvZS0mrvqZ4Mikb803FpCIjIRbt/RVNnGmqd0sTEQCuBCpzFghaG33bxlYZSChofDTxVXUWFeMGUfH8iTCP0n3dbS15q91vzNu0ueogHusnPA2Ndo1Y8lfe4mOiU1yf+3t9TSKv79+tkp4eJJvJUxMArgQaST+eeRNijekeeHOLDs/nxrrivHpsfEERD5O9D1tbGwYPmg0Z45707PfGEKvHubrdWt558e9bL/2gLgkpkwtLPQlhpaWcPcu+PrqKRaRtiQHLkQ6YTDA8VtX+ezwdLbeWYuDdRb6lB/G+xU+wNkuaVXufX3vcysclh29xZl9W3BWoXw2diI1ShRKVl/Dw/V0inG1imzJT1nyEFMIM2EwwJGbl5h1eDp/3d1AFmtH3nMbwYAKo8lumyNJ94yNU/QY0IP9f63F0ik31ToNZs6oUbjmSvqKFanLmXrkIaYQZsLaGuqWKcdvPdexp+t56uRrwZenP2HvvX+SfE9LC411K9bw09rt5Mubj8Mrp9G4qScjvl7Fk7DIJN3TwkJfoaJpcOcO+Pnpa8hF6pEZuBDpnMEAh29eJDdlsbayYM3NhYQYnjLCcyIWWuLnYHFxcaz95SdmfjoBu/rvkaNMDbp7FuOdSsVxsEn61vzwcH0NeZ48+uoVSauYjszAhTBT1tZQr0x5ShS3IGtWOOd3llMPj6M9+9/XEJe4p4kWFhb0ePtdzhy7yeapY6hRNDfzP59K9fYt+W7nIWKSuFbQwUGvw/nwoT4jj0zaxF4kggRwIcyEcdXK2u7fsqr1L4SGwmW/G1RdU5hFZz4jzBCaqPvZ2NhQOIcjM1tV5i23ogRd3MfkPg2o07MH285eJSm/nVta6mkVpfQlh5JWSVkSwIUwM9bWUCi/Da6u4JAljnI5qjDr+ARqrCvGV2fnEG4IS/Q9Z0+ZxeGD16jXtAP39q2nX8cqdJoym8t+QUnqo7EuZ0gI3L4tdTlTiuTAhTBzBgPsuHKUzw5PZb/v3+S0y83giuPoXW4Q9laJr3R//sIZRk8cRUyVbkRlK0T9wtkY2tCDgtmTdtas1OVMPllGKEQGZzDAP5cP8+mhKRzy205u+7wMrjiOd8sOxN4q8Wv8wqIM/HjqFgumDCT6iQ8d3vuIKX17ktUuaWfNSl3OpJOHmEJkcNbW0Nq9Jnve+4ff2x+gRFY3ph0Zzcyj45J0vyy21rxfszQfefXGJiaMdTPfo2qz2szb+GeStubb2uq7OYOC9LRKSIikVZJLZuBCZFAGA/x5cR9OcUUpkrUoN8JOcfLRId4p0x9by8TlMSIjI/l8yecs/3o2hohQinf+kEmDh9GkVAEskrBeMDZWX3bo4KA/mJW0yutJCkWITMpYs3PSno/4+eYyDne7TVbbrEm6V2BgAB/Pnoq/SyPuRltTxCKUIU2qUbdMsSTdT+pyJowEcCEyuehoxUUfHxwMhUGLZeD+NrR06cDbpbywsUxcXjtOKbZduc/I91oS/uguHm368fn4CZTKl/gzW6Qu55tJDlyITM7GRsPTtTCurmCw8ScwIpBxB96n7oZSrLmyIlEbgiw0jZZlC7H2q+9xLenGqZ/n07SxG31nzsE/NCJR/dI0PWhbW4OPD9y/rz/sFG8mAVyITMbaGtxc8nFswGHWtPqTHDZ5GLu/P/U2lGb91e+JiUv4zptKHlXY8+cBvvnud3Jkzcq2r8bRYvQEVhy5Rnh04nbwWFn9v4DE7dsQECAFJN5EUihCZHLR0YoNZ7by2ZEpXAo8RdGsxRnl+TEdSryDlUXCz0aJjY3l+7WruOFYhn13A7H2u0LXGu4MaN0Yq0RWR5a6nP8mOXAhxGtFRyvWnf6dWUemcDnoDPUKNmVtq6SdgHjBN4C329fk6f0b5KvcgikTZvJWtYqJLuIsdTl1kgMXQryWjY1Gr+ptOTPoFN83/5WuxYYQHg4Rhki23NxAbFzC13675XfmwF8HadmlH36n/2Hw2zVo1LcPx2/eTVSfpC7n68kMXAjxUsblh9+eXMmEY334tc1+quevk+j73Lp1g5GTP+Dkni3kajuOdu26MrBmaQpkS9w2/7g4Pa1iZaUXkHBI/CkBZktSKEKIJImKjmPL+Z24OzbF0hJ+uDGfgo5FaFWsY6LOIz959hQnwhxYd+Y2Qae3UbNUUT4bOphs9onbxWMw6DPyrFkhV67MkVZJcgDXNK0wsBrIB8QBy5RSX2qa5gysB1wAb+BtpVTg6+4lAVwI82UwwOMnsTReX4XLQWco6+zO6EpTaOHSPlGB/FFIBC3b1OThjbPYFyxDn2FTGdO9M7ZWidvFExGh7+jMDHU5k5MDjwE+UEqVBWoAQzRNKweMB3YqpUoCO599LoTIoKytIX8+S04NPMGSxj8Sboig/45OtNhUmW3emxN8fngeJ3uO7TrJuKmLIfQJX43vRuUW9Vi9Yz9xicgI2NvraZRHj/SzxyMSt/w8Q3hjAFdK+SqlTj37cwhwGSgItANWPXvZKqB9CvVRCJGO2NlaMrhODy4PucSiRqsIiQqh7/b2tNxUhe13/khQILe0tGR4/yGcP+lNz4ETCLl9li93HKff+oOcvPc4wX3J7HU5E5UD1zTNBdgHuAF3lVLZ430vUCn12pLZkkIRIuOJiIph+dEfmXtsOj5ht/HIXY2NbfZiZ2WX4HsEhwRz6H4w3xy+xuXfFlI0b24+nzwdt8L5E9WXjFqXM9nLCDVNcwQ2AiOVUsGJuG6ApmknNE074e/vn9DLhBBmwt7WiuH1vLg69CrzG6ygYo5axEXboRRcenIuQTPyrE5ZaVGmEGt61MXVPpYLW1fQqnE5ekycxIPAkAT3JbPV5UzQDFzTNGvgD2CbUmr+s69dBRoopXw1TcsP7FFKlX7dfWQGLkTGZ1x+ePLOFVr/WY7pNRfS121oou6x7/ABxk4aic+Vk1g7F6T3hwv4sGsHstgkfGdoVJT+kSOHftqhVcIvTXeSPAPX9K1T3wKXjcH7mS1A72d/7g1sNkVHhRDmzVh8uV4FV2bX+4om+boRHg7HHx7iwP1dCZqR16tZhyM7jjNv8QayOjry151Quq7ew4aT14hJ4E6ezFCXMyHLCOsA+4Hz6MsIASYAR4ENQBHgLtBFKRXwunvJDFyIzMc4I+++pR277m+her56jK0ynZr56yfoeqUUl/yCWHzgCv8sGEIWRycmfvQpb9evleCt+eZel1M28ggh0lRIRCSLDi3nixOf4R/pS638DRlTeRrV89dN0PUxMTGMmT6RX39YRGxMNC612zPn48+oXa5kgvtgrnU5JYALIdKF4PAIvjz4DQtPzeJxpB91CjRmTOVpVM1XO0HX+/r5Murjcezf+hOatS0thn/BtL49KZjArfnGkw417f8FJNL7ahUJ4EKIdOVpeDhfHlzKlydnERDlT/2CzZhbbwUFHQsn6Przly8yfubHhFXuQZy1HU3y2zCkWS1yOiUskJtTXU4J4EKIdCkoLIwFB75i7eXVbGxxFGcnB0IMT8lqky1B1/uHRrLs0CWWj2mHhWZBt4ET+bhfX+ysE7bsJDJST6vkypV+63JKABdCpGvR0YqgIA3/JzG0+rMMLYu15+Ma8xJ0rVKKVT+vZdZnHxHy6C5OLu4MH/sJA9u1wSIB+ZH0XpdTzgMXQqRrNjYaefJAEZcYepTrSxXnJoSHQ2BkIBcen37ttZqm4fX2O5w/doOBY+cQ9fgeM4e0pd2MLzmRgK355lqXU2bgQoh0ybj8cPq+qSy5OI3mRdvzQeWplM9Z8Y3XBocEM33RAm7mrYVfWDTFo+4xskMLKpdwSVDbERH6mSrGkw4TWRHO5CSFIoQwS/4hQczd9yVLzy4gxPCUVi6dGF15CmWdK7zx2qiYWNafvMbkd2sTGxNF9bf6sWDSFIrkdn7jtempLqekUIQQZim3U3bmtJ7CjaG3GV15Mnt9/qHJRnfe3/E2VwMuvvZaWytLelUvy2+/7qZEhRoc+WUhdeuWYtDMTwmOiHrttZqmB21LS7h3T0+rGAymHFnyyQxcCGFW/IIDmLV3PsvOfUlETBhtXbsytsoMimUr8cZrt+7czsSpo/G/dYESfRcwvEt72pUvgpXlm+eyxrRKnjyQLVvqplUkhSKEyFB8nz7ms72f8+35RaxstJ1aRWoC6o3b65VSrPvrTw5EOHP6fgDW13bh1bIpXq2av/HatKrLKQFcCJEhPQ55SlxENgIDYcbJIcRZRDO37rIEBfLdV+7wXseqRAc/pmClpnw6ZQ5Nqni8sc3UrsspOXAhRIaUyykbefKAqyvkdMyGvZadiAgNpeBh2INXXqdpGo3KunDqyHVadBvKg/P76N2xKs169+Si973XtmltrZ90GB6un3QYGJg2Jx3KDFwIkaEYlx9uv7qf3rsb0qWkFyMrTaKwk8trr7t19w4jJo3h9J7NFHpvCd0a1sKrakmy29u89rq4OH0TkI2NnlaxtzfhYJ6RGbgQIlMwnkde160EfdyG8OuNH6mzviRj9w3AJ+TOK69zLVKU31f/zM7912lfpzq/nPWmTrdOjJwzj/CoV+/qScu6nDIDF0JkaN4B95mx+zN+uLQcUHQr9R7DPCe88dCsS/ce0rVrIwLuXMYhX3HeHzmN0T3feePW/JSoyykzcCFEpuTiXJBvOy3myuAb9Cjbj3XXvqX2+hJMODCEB6E+r7yuXOF8nD1wgbEzlqGiw1kwvicVG1bjj6OnXtteatbllAAuhMgUXHMW5vvOX3Fp0HW6lfHipyvLWH1hxWsfPlpYWDCyb3/OH79F98FTCQ/wZcaea4zefIyrfoGvvM7SUk+rKAXe3npaJTbW9GOSFIoQIlO67u+NFpmD2PBs7PHdwvHHuxhX9ROyWL/6KMKIqGg2XfRh5bGr3Ph2FGU9arBw6mcUL5j/ldcopdfjLFBA3wCUFJJCEUKIeErmdqFE4Wy4usLtiHPsvbed2Cg7lAJD3Mv3zNvb2vBOJVdWdalOqZJlOfPnahrUK4XXhPEEBoe+9BpN0zf+pAQJ4EKITM3aGma1nMTpgadwzm6FX2AotdeVZPqRsTyJ8H/pNflzOvP3mk2s2XSYAiXc2b5qNpVqluCrrTuJiY176TUpQQK4EEIAjva25MkDeQqGUbNgXZZfmE+NdcWYeXQ8AZEvP1O8frXqHN12kNlfbyJX0bL8cC2Enj/tY8vx86RGeloCuBBCxFMoR15+7v4DZ/pfpLlLW74+N4caa4sx69hEAiMDXnpNz7btObZ1F/Pa10QZIhnaswGeTWqxed+BFO2rBHAhhHiJCvnL8GuPNZzqf57GRVux6OynVF/rwpzjkwmK+u8KFE3TqF0sL9+9U5/27wwm4M5FBnevR/0ubTl97WqK9FFWoQghRAKcun+eKbum8cetjdTJ34x1rbe9dpOOr78/I6Z+xKHfV+NQoho//biRdlXyJqltOY1QCCFM4ITPWQKfxlDEqjJBhkds9F7Ge27DcbLJ+tLXn792lQ2n7jCtR0MK5U3asYWvCuAptLhFCCEypiqFKkIh/dCsJYf+YP6pqTQt8Dbl8mZ96Yy8QqnSuBYsjZOd6fsiOXAhhEgCa2sYWb8vVwbfpGKhUoSGwqhd/Vl0ZhZhhpevCTc1CeBCCJEMJXIVJU8eKFQ0moAYX2Yd/4jqa4vx1dm5hBvCUrRtCeBCCGECWexs2O71B/t7H8YjT2VmHvuQGutcWXr2cyJiwlOkTQngQghhQnVcarCr79/s7XWQ8rncmXFsDGMO9k6RtuQhphBCpIB6xWqx973t7Ll1gNiILClSN1MCuBBCpKAGrnVS7N6SQhFCCDMlAVwIIcyUBHAhhDBTEsCFEMJMSQAXQggzJQFcCCHMlARwIYQwUxLAhRDCTKXqeeCapvkDd5J4eS7g5YXpMi4Zc+YgY84ckjPmokqp3C9+MVUDeHJomnbiZQeaZ2Qy5sxBxpw5pMSYJYUihBBmSgK4EEKYKXMK4MvSugNpQMacOciYMweTj9lscuBCCCH+zZxm4EIIIeKRAC6EEGYq3QVwTdNaaJp2VdO0G5qmjX/J9zVN0xY++/45TdMqpUU/TSkBY+7xbKznNE07pGlaxbTopym9aczxXldV07RYTdM6p2b/TC0h49U0rYGmaWc0Tbuoadre1O6jqSXg33U2TdN+1zTt7LMx90mLfpqSpmnfaZr2SNO0C6/4vmnjl1Iq3XwAlsBNwBWwAc4C5V54TSvgL0ADagBH07rfqTDmWkCOZ39umRnGHO91u4A/gc5p3e8U/jvODlwCijz7PE9a9zsVxjwBmP3sz7mBAMAmrfuezHHXAyoBF17xfZPGr/Q2A68G3FBK3VJKRQPrgHYvvKYdsFrpjgDZNU3Ln9odNaE3jlkpdUgpFfjs0yNAoVTuo6kl5O8ZYBiwEXiUmp1LAQkZ7zvAr0qpuwBKqcwwZgU4aZqmAY7oATwmdbtpWkqpfejjeBWTxq/0FsALAvfife7z7GuJfY05Sex43kP/CW7O3jhmTdMKAh2ApanYr5SSkL/jUkAOTdP2aJp2UtO0XqnWu5SRkDEvBsoCD4DzwAilVFzqdC/NmDR+pbeixtpLvvbiOseEvMacJHg8mqY1RA/gKVclNXUkZMxfAOOUUrH6BM2sJWS8VkBloDFgDxzWNO2IUupaSncuhSRkzM2BM0AjoDiwXdO0/Uqp4BTuW1oyafxKbwHcBygc7/NC6D+dE/sac5Kg8Wia5g6sAFoqpZ6kUt9SSkLGXAVY9yx45wJaaZoWo5T6LVV6aFoJ/Xf9WCkVBoRpmrYPqAiYawBPyJj7ALOUnhy+oWnabaAMcCx1upgmTBq/0lsK5ThQUtO0Ypqm2QDdgC0vvGYL0OvZ09wawFOllG9qd9SE3jhmTdOKAL8C75rxjCy+N45ZKVVMKeWilHIBfgEGm2nwhoT9u94M1NU0zUrTNAegOnA5lftpSgkZ81303zjQNC0vUBq4laq9TH0mjV/pagaulIrRNG0osA39KfZ3SqmLmqYNfPb9pegrEloBN4Bw9J/iZiuBY/4YyAl89WxGGqPM+CS3BI45w0jIeJVSlzVN+xs4B8QBK5RSL12KZg4S+Hc8A1ipadp59NTCOKWUWR8xq2naWqABkEvTNB9gCmANKRO/ZCu9EEKYqfSWQhFCCJFAEsCFEMJMSQAXQggzJQFcCCHMlARwIYQwUxLAhRDCTEkAF0IIM/U/YMaevrUwSosAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot([0, 1], stats_timescales.mean, label=\"bmsm timescales\")\n",
    "plt.fill_between([0, 1], stats_timescales.L, stats_timescales.R, color='b', alpha=.1)\n",
    "plt.plot([0, 1], bmsm_posterior.prior.timescales(), \"k--\", label=\"bmsm prior timescales\")\n",
    "plt.plot([0, 1], true_msm.timescales(), \"g-.\", label=\"ground truth timescales\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To, e.g., visualize samples differently or evaluate other kinds of moments, one can get access to the evaluated quantities directly by setting `store_samples` to `True`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:40:37.568607Z",
     "iopub.status.busy": "2020-12-02T13:40:37.567912Z",
     "iopub.status.idle": "2020-12-02T13:40:37.584612Z",
     "shell.execute_reply": "2020-12-02T13:40:37.584188Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# samples: 100\n"
     ]
    }
   ],
   "source": [
    "stats_mfpt = bmsm_posterior.gather_stats(\"mfpt\", store_samples=True, A=0, B=1)\n",
    "print(\"# samples:\", len(stats_mfpt.samples))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. rubric:: References\n",
    "\n",
    ".. footbibliography::"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}