Simulation - Expected number of returns of a transient state in DTMC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b421bbde",
   "metadata": {},
   "source": [
    "# <center> Author: Tianjian Che\n",
    "\n",
    "\n",
    "Starting from state $i$, denote $f^{(n)}_{ij}$ as the probability that the discrete time Markov chain (DTMC) arrives at state $j$ for the first time exactly after $n$ steps, that is:\n",
    "\n",
    "$$\n",
    "f^{(n)}_{ij} = P\\big( X_{1} \\neq j, X_{2} \\neq j, \\ldots, X_{n-1} \\neq j, X_{n} =j ~ | ~ X_{0} =i \\big)\n",
    "$$\n",
    "\n",
    "Naturally, we have:\n",
    "\n",
    "$$\n",
    "f_{ij} = \\sum\\limits^{\\infty}_{n=1} f_{ij}^{(n)}\n",
    "$$\n",
    "\n",
    "be the probability that started from state $i$ and the DTMC finally arrives at state $j$ after some steps, and we have $f_{ij} \\leq 1$ by the law of total probability."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af0f276b",
   "metadata": {},
   "source": [
    "Here it comes the definitions of **recurrent state** and **transient state**:\n",
    "\n",
    "A state $i$ is called **recurrent** if: \n",
    "\n",
    "$$\n",
    "P\\big( X_{n} = i \\mbox{ for some } n \\geq 1 ~ | ~ X_{0} = i \\big) = 1\n",
    "$$\n",
    "\n",
    "and state $i$ is called **transient** if:\n",
    "\n",
    "$$\n",
    "P\\big( X_{n} = i \\mbox{ for some } n \\geq 1 ~ | ~ X_{0} = i \\big) < 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4da5171c",
   "metadata": {},
   "source": [
    "The above two mathematical definitions cannot be further intuitive and obvious. A **recurrent state** $i$ means that, the DTMC which started from state $i$ can always return to itself ($i$) for probability $1$, it further indicates that the DTMC will visit state $i$ for infinitely many times. By definition above we know that in this case: $f_{ii} = 1$.\n",
    "\n",
    "Similarly, a **transient state** $i$ means that, started from state $i$, there exists a strictly positive probability (which is therefore $1 - f_{ii} > 0$) that the DTMC will never come back again. Hence, the DTMC will \"hit on\" this probability as long as it's strictly positive given a sufficiently long but finite time. In other words, sooner or later, the DTMC will leave from state $i$ and never come back again. In this case, the expected number of returns to state $i$ is also a finite and can be derived as following:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "S & = 1 \\cdot f_{ii} \\cdot (1-f_{ii})  + 2 \\cdot (f_{ii})^{2} \\cdot (1-f_{ii}) + 3 \\cdot (f_{ii})^{3} \\cdot (1-f_{ii}) + \\cdots\\\\\n",
    "f_{ii} S & = (f_{ii})^{2} \\cdot (1-f_{ii}) + 2 \\cdot (f_{ii})^{3} \\cdot (1-f_{ii}) + 3 \\cdot (f_{ii})^{4} \\cdot (1-f_{ii}) + \\cdots\\\\\n",
    "(1-f_{ii}) S & = f_{ii} \\cdot (1-f_{ii}) + (f_{ii}) ^{2} \\cdot (1-f_{ii}) + (f_{ii})^{3} \\cdot (1-f_{ii}) + \\cdots\\\\\n",
    "S & = \\lim\\limits_{n \\rightarrow \\infty}\\frac{f_{ii}\\big( 1 - (f_{ii})^{n} \\big)}{1 - f_{ii}} = \\frac{f_{ii}}{1 - f_{ii}}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8c6eb4c9",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-11T09:41:10.402114Z",
     "start_time": "2022-09-11T09:41:09.786101Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a2c8f6ee",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-11T09:41:14.070028Z",
     "start_time": "2022-09-11T09:41:10.404048Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(66)\n",
    "actual = []\n",
    "for p in np.linspace(0.1, 0.9, 81):  # denote as f_ii\n",
    "    count = 0\n",
    "    for n in range(10000):\n",
    "        while np.random.uniform(0, 1) <= p:\n",
    "            count += 1\n",
    "    count = count / 10000\n",
    "    actual.append(count)\n",
    "\n",
    "theoretical = list(map(lambda p: p/(1-p), np.linspace(0.1, 0.9, 81)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "234b7cf9",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-11T09:41:14.179374Z",
     "start_time": "2022-09-11T09:41:14.071030Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArkAAAHSCAYAAADohdOwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABGMElEQVR4nO3dd3yV5cH/8c+VhLBXIICACigQRAQRtxUUrYvUuuqe1VZbrdq6qo9tfZ62arWO6s9Wa22rorYOqlg3MlRwsHcYYYOQAQFkJeT+/ZGIqIwASe6Tcz7v1ysvkpz73Pnm8oDfXLnu6w5RFCFJkiQlk7S4A0iSJEnVzZIrSZKkpGPJlSRJUtKx5EqSJCnpWHIlSZKUdCy5kiRJSjoZNXHS1q1bR506daqJU0uSJEkAjBs3rjCKouxtPVYjJbdTp06MHTu2Jk4tSZIkARBCWLC9x1yuIEmSpKRjyZUkSVLSseRKkiQp6dTImtxtKS0tZfHixWzYsKG2vmRKadCgAR07dqRevXpxR5EkSYpdrZXcxYsX07RpUzp16kQIoba+bEqIooiioiIWL15M586d444jSZIUu1pbrrBhwwZatWplwa0BIQRatWrlLLkkSVKlWl2Ta8GtOY6tJEnSV1LmwrNVq1bx2GOPATBixAgGDRpUq1//H//4B0uXLt3y8ZVXXsn06dN3+TxxZJckSaprUrLk1pSysrLtPvbNkvvkk09ywAEH1GgeSZKkVJW4JTc/H3r2hIyMij/z8/fodLfddhtz586lT58+3Hzzzaxdu5azzz6bnJwcLrzwQqIoAmDcuHH079+fQw45hJNOOolly5YBMHHiRI444ggOOuggzjjjDFauXAnAgAEDuP322+nfvz8PP/zwNp//0ksvMXbsWC688EL69OnD+vXrGTBgwJa7wr311lv07duX3r17M3DgQAA+/fRTjjrqKA4++GCOOuoo8vLy9uj7lyRJSiW1trvC19xwA0ycuONjPvsM1q2reH/6dOjVCw49dPvH9+kDDz203Yfvuecepk6dysSJExkxYgSnn34606ZNo3379hx99NF89NFHHH744Vx33XW8+uqrZGdn869//Ys77riDp556iksuuYRHHnmE/v3786tf/Yq77rqLhyq/3qpVqxg5ciSlpaX0799/m89/9NFHuf/+++nXr9/XchUUFHDVVVcxatQoOnfuTHFxMQA5OTmMGjWKjIwM3nvvPW6//XZefvnlHY+ZJEmSgLhKblV8WXC39/EeOuyww+jYsSMAffr0Yf78+bRo0YKpU6dy4oknArB582b22msvSkpKWLVqFf379wfg0ksv5ZxzztlyrnPPPReAvLy8bT5/Rz7++GOOPfbYLVt/ZWVlAVBSUsKll17K7NmzCSFQWlpajd+9JElScoun5O5gxnWLnj1h5kwoL4e0NMjJgREjqi1C/fr1t7yfnp5OWVkZURTRs2dPxowZ87VjS0pKdniuxo0bA2z3+TsSRdE2d0a48847Oe644xgyZAjz589nwIABVT6nJElSqkvcNblDh1YU2/T0ij+HDt2j0zVt2pQ1a9bs8Jju3btTUFCwpaSWlpYybdo0mjdvTsuWLfnggw8AeOaZZ7bM6lbl+Tv6+kceeSQjR45k3rx5AFuWK5SUlNChQweg4qI1SZIkVV3iLlfo0gUqC2J1aNWqFUcffTQHHnggDRs2pG3btt86JjMzk5deeomf/exnlJSUUFZWxg033EDPnj355z//ydVXX826devo0qULf//733fp+ZdddhlXX301DRs2/NpMb3Z2Nk888QRnnnkm5eXltGnThnfffZdbbrmFSy+9lAceeIDjjz++2sZBkiQpFYQvdxWoTv369Yu+3DngSzNmzKBHjx7V/rX0FcdYkiSlkhDCuCiK+m3rscRdriBJkqSEt2wZ1MCc6R6z5EqSJGm3lJZCp05w551xJ/k2S64kSZJ2y8z3FrFpExzw+4uq5eZd1cmSK0mSpN0y+cpHAOgdTajY+jU3N+ZEX7HkSpIkabdMWppNfTbQnbyKexvk5cUdaQtLriRJknbLpMxD6ck0MthccfOu7t3jjrSFJXcbRowYwejRo/foHE2aNKmmNJIkSYlpUlofDmJKtd28qzol7s0gYjRixAiaNGnCUUcdFXcUSZKkhLR8OSzf0ILeXdbA3LK443xLSs3kfv/73+eQQw6hZ8+ePPHEEwC89dZb9O3bl969ezNw4EDmz5/PX/7yFx588EH69OnDBx98wGWXXcZLL7205TxfztKuXbuWgQMH0rdvX3r16sWrr74ay/clSZJU2yZP2AxA794xB9mOWGZyb7gBJk6s3nP26QMPPbTjY5566imysrJYv349hx56KKeffjpXXXUVo0aNonPnzhQXF5OVlcXVV19NkyZNuOmmmwD429/+ts3zNWjQgCFDhtCsWTMKCws54ogj+N73vkcIoXq/OUmSpAQzaXgxkM1Bx7eOO8o2pdRyhT/96U8MGTIEgEWLFvHEE09w7LHH0rlzZwCysrJ26XxRFHH77bczatQo0tLSWLJkCcuXL6ddu3bVnl2SJCmRTBr9BR3YSKtjesQdZZtiKbk7m3GtCSNGjOC9995jzJgxNGrUiAEDBtC7d2/yqrDVRUZGBuXl5UBFsd20aRMAgwcPpqCggHHjxlGvXj06derEhg0bavT7kCRJSgST8urTO0yAA06IO8o2pcya3JKSElq2bEmjRo2YOXMmH3/8MRs3bmTkyJHMmzcPgOLiYgCaNm3KmjVrtjy3U6dOjBs3DoBXX32V0tLSLeds06YN9erVY/jw4SxYsKCWvytJkqTat3EjzCjMpnf2MsjMjDvONqVMyT355JMpKyvjoIMO4s477+SII44gOzubJ554gjPPPJPevXtz7rnnApCbm8uQIUO2XHh21VVXMXLkSA477DA++eQTGjduDMCFF17I2LFj6devH4MHDyYnJyfOb1GSJKlWzJwJZVEGvXtsijvKdqXMmtz69evz5ptvbvOxU0455Wsfd+vWjcmTJ3/tcx9//PGW9++++24AWrduzZgxY7Z5zrVr1+5JXEmSpIQ1aeQqoAW9j07c+wKkTMmVJElS9Zg0YiUNqM/+A/eNO8p2pcxyBUmSJFWPSZPhQKaScUiCbpKLJVeSJEm7IIpg0qJW9G6aD82bxx1nu2q15EZRVJtfLqU4tpIkqTZ8/jkUbmpWcTvfBFZrJbdBgwYUFRVZxmpAFEUUFRXRoEGDuKNIkqQkN2nMOgB690vMrcO+VGsXnnXs2JHFixdTUFBQW18ypTRo0ICOHTvGHUOSJCW5Se+uADrR64S2cUfZoVorufXq1dty+1xJkiTVTZM+28Q+LKDlMT3jjrJDXngmSZKkKps0tzG9M2dAhw5xR9khS64kSZKqZMMGyFvVlt4dCiGEuOPskCVXkiRJVTJ9cimbyeCgXom/kYAlV5IkSVUy6e3lAPQ+tkW8QarAkitJkqQqmfThGhrxBft9d7+4o+yUJVeSJElVMnlGBr3CNNIP6B53lJ2y5EqSJGmnoggmLWtD79ZLID097jg7ZcmVJEnSTi1ZHFFc1pyDum2IO0qVWHIlSZK0U5Peq7hrbe8jGsacpGosuZIkSdqpScMKATjolMS+CcSXLLmSJEnaqckTy+lMPs2OTOzb+X7JkitJkqSdmrSgOQc1yYdGjeKOUiWWXEmSJO3Q+vUwa217encqiTtKlVlyJUmStENTPyqhnHR6H1x3qmPdSSpJkqRYTHprGQC9j28dc5Kqs+RKkiRphyZ/sp4mrKHzyYl/p7MvWXIlSZK0Q5NmNaRXvTzS2rWJO0qVWXIlSZK0XVEEkwrb03uvFXFH2SWWXEmSJG3XwlkbKClvRu8DSuOOskssuZIkSdquSa8vAqD3MU1jTrJrLLmSJEnatvx8Jv/mFQAOfOrnkJ8fc6Cqs+RKkiRp23Jz+XRtD7oyi6bzp0BubtyJqqxKJTeEcGMIYVoIYWoI4fkQQoOaDiZJkqR4lc2cw0j6cxzDobwc8vLijlRlOy25IYQOwM+AflEUHQikA+fVdDBJkiTFa+zeZ7Ca5gxkGKSlQffk2yc3A2gYQsgAGgFLay6SJEmSEsH7h90GwHFhJOTkwNChMSequoydHRBF0ZIQwv3AQmA98E4URe/UeDJJkiTFati45vRmItlFM6Fly7jj7JKqLFdoCZwOdAbaA41DCBdt47gfhRDGhhDGFhQUVH9SSZIk1Zr16+Gj+R04PmtSnSu4ULXlCicA86IoKoiiqBR4BTjqmwdFUfREFEX9oijql52dXd05JUmSVIvGfFTOxvJMBvYriTvKbqlKyV0IHBFCaBRCCMBAYEbNxpIkSVKchr1YTDplHHt63ZvFhSqU3CiKPgFeAsYDUyqf80QN55IkSVKMhr1XzmF8StOBh8UdZbdUaXeFKIp+HUVRThRFB0ZRdHEURRtrOpgkSZLiUVICn81rzcCGo6Fbt7jj7BbveCZJkqSvGTUKyqM0BvYphhDijrNbdrqFmCRJklLLsP+upwFw5Kl1cz0uWHIlSZL0DcPeKuUYPqH+sYfHHWW3uVxBkiRJWyxfDlMXNGNg2nDo1y/uOLvNkitJkqQthg+v+HNg9yXQqFG8YfaAyxUkSZK0xbB3y2nOavqekBV3lD3iTK4kSZK2GPZWKQMYQfoxR8YdZY9YciVJkgTAvHkwb2l9BjIMjjoq7jh7xJIrSZIkAN5/v+LPge2mQ8eO8YbZQ5ZcSZIkATBsGLRLW0GPY7PjjrLHLLmSJEkiiuD99zZzfPm7hKPr9lIFsORKkiQJmD4dlhekJ8V6XLDkSpIkiYqlCgADG4yG3r3jDVMNLLmSJEli2DDoUn8x+x6xF9SrF3ecPWbJlSRJSnFlZTBiRMTATW8lxVIFsORKkiSlvPHjYfXqwMDoXUuuJEmSksOX63GPYzgccUS8YapJRtwBJEmSFK9hw6BX03m06dAKWrWKO061cCZXkiQphW3YAB99lFzrccGSK0mSlNLGjIENGwLHb3zDkitJkqTkMHw4pIVyjmVUUpVc1+RKkiSlsBEjoG+rhTTfnA7du8cdp9o4kytJkpSi1q2Djz+OOK74ZVi1Cnr1gvz8uGNVC0uuJElSiho9GkpLA8eVvwdRBDNnQm5u3LGqhSVXkiQpRY0YAemUcQwfVnyivBzy8mLNVF0suZIkSSlq+HDolz6Rpqyt+ERaWtKsy7XkSpIkpaC1a+HTT+G4pmOhcWNIT4ecHBg6NO5o1cKSK0mSlIJGj4ayMhiwagjcdVfFB9OmQZcucUerFpZcSZKkFDR8OGSkl3M0H8EJJ8Qdp9q5T64kSVIKGj4cDms5hyZpjSu2DksyzuRKkiSlmDVrYOzYiAHr36iYxU1LvkqYfN+RJEmSdujDD2Hz5sBxX7yelEsVwJIrSZKUcoYPh3rpmzmK0XDiiXHHqRGuyZUkSUoxI0bA4c1m0KjtvtCxY9xxaoQzuZIkSSmkpATGjYs4bs3QpF2qAJZcSZKklPLBB1BeHjiu7J2kXaoALleQJElKKSNGQGZ6GUdEn8GAAXHHqTGWXEmSpBQyfDgc2WgSDXv1hmbN4o5TY1yuIEmSlCJWroQJEyrX4ybxUgWw5EqSJKWMDz6AKAoMYHjSl1yXK0iSJKWI4cOhQfomjmg4DQ47LO44NcqZXEmSpBQxfDgclTmW+scdBfXqxR2nRllyJUmSUkBxMUyeHDFg/ZtJv1QBLLmSJEkpYeTIivW4xzE8qW8C8SXX5EqSJKWA4cOhYfpGDmu7BHJy4o5T45zJlSRJSgEjRkQcHcaQ+d0BEELccWqcJVeSJCnJFRTAlCmVt/JNgaUKYMmVJElKeiNHVvyZKutxwTW5kiRJSW/ECGicvp5+B2yCtm3jjlMrnMmVJElKcsPfL+eY8g+od+KAuKPUGkuuJElSEluxAqbPSGNA9H5K7I/7JUuuJElSEhv14nIABjACbrwR8vPjDVRLLLmSJElJbMSdw2jMWg5hHMyaBbm5cUeqFV54JkmSlMRGruzF0XxEPcqgHMjLiztSrXAmV5IkKUkVFMBUetGfyj3E0tKge/d4Q9USZ3IlSZKS1KhRFX8OqDcaytMrCu7QofGGqiXO5EqSJCWpkW+tpxFf0O+mAVBWBtOmQZcucceqFc7kSpIkJakRb63nKMaReVZqXGy2NWdyJUmSklBREUxZnMWAZhOgb9+449Q6S64kSVISGvXuRgAGnJgBIcScpvZZciVJkpLQiGcX05B1HHpln7ijxMKSK0mSlIRGjq7HURmfkjnwO3FHiYUlV5IkKckUryhj8sqO9O9ZBPXqxR0nFpZcSZKkJDPqiZlEpDHg7NZxR4mNJVeSJCnJjHy5kAas57Cf9Is7SmwsuZIkSckkihgxLZsjW82mflbjuNPExpIrSZKURFaOmsKk0h4MOKYs7iixsuRKkiQlkQ/+PJWINPpfnhq3790eS64kSVISGTFsM/XDRg4/qUXcUWJlyZUkSUoW+fmMLDyAI7usoEGDuMPEy5IrSZKUJFa98BYTOJj+pzWJO0rsLLmSJElJ4sPnFlbsj3tGy7ijxM6SK0mSlAwKChgxLZvM9DIOPzzuMPGz5EqSJCWD119nBP05ovd6GjaMO0z8LLmSJElJoOTfbzOBgxngelzAkitJklS35edDjx58+NYaykmnf7dlcSdKCJZcSZKkuiw3F/LyGMEAMtnIEb/LjTtRQrDkSpIk1WV5eRBFjKQ/h/EpjWZPijtRQrDkSpIk1WX7789qmjKOQxjASOjePe5ECcGSK0mSVJddeikfcTTlpDNgn3wYOjTuRAmhSiU3hNAihPBSCGFmCGFGCOHImg4mSZKkKnj7bQY3vYbGjSOOnPEUdOkSd6KEkFHF4x4G3oqi6OwQQibQqAYzSZIkqSoWLmTByHm8kDaIn10faGRD22KnJTeE0Aw4FrgMIIqiTcCmmo0lSZKknRo8mAe5kZAWuOGGuMMklqosV+gCFAB/DyFMCCE8GUJo/M2DQgg/CiGMDSGMLSgoqPagkiRJ2koUUfT31/hr2o85//zAPvvEHSixVKXkZgB9gT9HUXQw8AVw2zcPiqLoiSiK+kVR1C87O7uaY0qSJOlrJkzgz7MHsq68ITffHHeYxFOVkrsYWBxF0SeVH79ERemVJElSTNb/4188wnWceuImevWKO03i2WnJjaLoc2BRCOHLTdcGAtNrNJUkSZK2r6yMp/8ZsYK23HJHZtxpElJVd1e4DhhcubNCPnB5zUWSJEnSjmx+Zxj3r76Kw7oWc+yxWXHHSUhVKrlRFE0E+tVsFEmSJFXFf+6dyRxO4qW7NhFC3GkSk3c8kyRJqkOiNWu598Oj2b/Zcr7/A5cqbI8lV5IkqQ4Zdc9oPivvx01XrSY9Pe40icuSK0mSVIf84YkWtEkr4JK79os7SkKz5EqSJNURU94v4I3Cw/jZdybSsLE1bkequruCJEmSYnb/rQU0piHX3O3tzXbGHwEkSZLqgEWL4LmxXbmqzWtkHdl9509IcZZcSZKkOuChOwuJCNx4zYa4o9QJllxJkqQEt2byPJ74ZwPO53n2ef5eyM+PO1LCs+RKkiQluCGnPMFamvATHoM5cyA3N+5ICc8LzyRJkhLcs0uPpwtzOYKPoRzIy4s7UsJzJleSJCmBLV0Kwziei3iWAJCWBt298GxnLLmSJEkJ7IU/F1NOOhe2fAPS0yEnB4YOjTtWwnO5giRJUgIb/OR6DuVTuk34N+y7b9xx6gxnciVJkhLU9EmbGP95By7qOdGCu4ssuZIkSQlq8G9mk04Z596xf9xR6hxLriRJUgIqL4fn3szixIYf0fbcAXHHqXMsuZIkSQlo9LP5zN+4FxedvqZiRwXtEkdMkiQpAT177xIa8QWn331E3FHqJEuuJElSgtlUtIZ/Tz+QMzpPpEmn1nHHqZMsuZIkSQnmzf/5iJW05KLrsuKOUmdZciVJkhJJFPHs82m0ySjihGtz4k5TZ1lyJUmSEkjJO58wtORYzjt2GRn1Qtxx6ixLriRJUgJ5+TeT2UgDLvr1fnFHqdMsuZIkSYlixQqe/aQrXVusoN93Gsadpk6z5EqSJCWIxQ++yIioPxddFAiuVNgjllxJkqREsHkzzz++mog0LrwhO+40dZ4lV5IkKW75+dC5M8+uPJUjMsexX8iPO1GdZ8mVJEmKW24uUxY1ZzK9uaj075CbG3eiOs+SK0mSFLeZM3mWi8iglB9E/4K8vLgT1XmWXEmSpJitadyOv3IV3+M1stOKoXv3uCPVeZZcSZKkOM2dy1/XnMdKsrgl7Y+QkwNDh8adqs7LiDuAJElSKtv0h4d4gFs57uiNHP7h6LjjJA1LriRJUlyWL2fwUxtZQkee+lXcYZKLyxUkSZJiUv7Qn7i37OccfMAGTjwx7jTJxZlcSZKkOKxezat/WkAeObzwK7zDWTVzJleSJCkG0V8e555117Ffxw2cdVbcaZKPJVeSJKm2bdjAyHvG8CmHc/P/NCDD361XO0uuJElSbXvmGe5Z+SPattzIpZfGHSY5WXIlSZJq0+bNTPi/13mbk7nh5kwaNIg7UHKy5EqSJNWmIUO4d9H5NGtYyjU/8WqzmmLJlSRJqi1RxJzfPMuLnMM116bTvHncgZKXJVeSJKm2vP8+9087mXoZEdffaA2rSV7LJ0mSVNPy8yE3l8+nF/EP5nPZmWvYa6+WcadKav4IIUmSVNNyc2HGDB7mekqpx03jL4g7UdKz5EqSJNW0vDxKoqY8xk84m5fYf967cSdKepZcSZKkmtaxI/dyK6tpzm3hD9C9e9yJkp4lV5IkqSZFEfMbHcAD/JyLwrMc3GMDDB0ad6qk54VnkiRJNWnoUG6bcQlp9TK4O/8i6HhR3IlSgjO5kiRJNaW8nNE3/pt/cR433xLo2DHuQKnDkitJklRDygc/z43519G+5Tpu+WV63HFSiiVXkiSpJmzaxPO/GMunHM7v/9iAxo3jDpRaLLmSJEk1YN1fnua2gp9zyP4lXHyplau2eeGZJElSdVu3jj/eUcxi9ua5v0Wk2XFrnUMuSZJUzZb+/h/cs/annNW/kO8cG+KOk5KcyZUkSapOJSXcfl8LykImf3jKhbhxcSZXkiSpGo37xXP8c9MFXH9xMV26xJ0mdVlyJUmSqkN+PlG37tz4t55khwLu+Pm6uBOlNEuuJElSdcjN5ZXZvfiAY/lffk3zCwbFnSiluSZXkiSpGhTNLODnvEFPpnJl9ATkxZ0otVlyJUmS9tDmsojzeZ7Pace/+QEZaRF0z4k7Vkqz5EqSJO2hO86bw7vlA3my2Y0c/sW4ioI7dGjcsVKaJVeSJGkP/PvZTdz7cleubvECP1zxB6j3YNyRhBeeSZIk7bYpU+DyK+AoPuLh57KhXr24I6mSJVeSJGk3rFwJZ+SW0bysiJdOeYrMUwbGHUlbcbmCJEnSLtq8GS68EBYughH1LmCvx/4edyR9gyVXkiRpF/361/Dmm/AXfspRv+wPnTrFHUnfYMmVJEnaBUOGwO9+B1dmvcyPGr8Ft8yIO5K2wZIrSZJURTNnwiWXwOGdPufR+RcQnngOGjWKO5a2wQvPJEmSqiCK4OqroX5mOS+tPIH6xx8DZ54ZdyxthzO5kiRJVfDq458zcmQ7/sxP6cg0uPl+CCHuWNoOZ3IlSZJ2YtMmuPn6jfRkKlfy14py+4tfxB1LO+BMriRJ0k48+ijM2bQvb3ESGWyGCMjLizuWdsCZXEmSpB0oLIT//V84JfM9TuKdik+mpUH37vEG0w5ZciVJknbgrrtg7dqI+6NfQNOmkJ4OOTkwdGjc0bQDLleQJEnajhkz4M9/jvhxmyEcsGERTJ8F7drFHUtVYMmVJEnajptugiaZm/jNsh/DPx604NYhllxJkqRteOcdeOMNuC/zLrJP7ldxFwjVGZZcSZKkbygrg1/8ImK/hku5Lu1xeHyCe+LWMZZcSZKkb/jb32Dq1MDLXEf9x34L++wTdyTtoirvrhBCSA8hTAghvF6TgSRJkuJUUgJ33rGZY9M/4ozvFMGPfxx3JO2GXdlC7HpgRk0FkSRJSgR3/z6isCjwQPrNhL89WbEnruqcKv1XCyF0BE4DnqzZOJIkSfGZN3IhD/5hE5fwNIe0zK/YE1d1UlV/NHkIuAUor7kokiRJ8bpl0HQyKON33AEFBZCbG3ck7aadltwQwiBgRRRF43Zy3I9CCGNDCGMLCgqqLaAkSVJt+GBkOS+tPZlbuZcOLIXycsjLizuWdlOIomjHB4RwN3AxUAY0AJoBr0RRdNH2ntOvX79o7Nix1ZlTkiSpxpSXw2H7Lmf54k3kkUMj1lWsxc3JgWnT4o6n7QghjIuiqN+2HtvpTG4URb+MoqhjFEWdgPOA93dUcCVJkuqaZ367gHGL23LPgYNpdECnirW4OTkwdGjc0bSb3CdXkiSltLXLv+CX/9uQw+uN5/z3r4Ls2+KOpGqwS3tiRFE0IoqiQTUVRpIkqbb94eT3Wba5DQ8+EJGW3SruOKombvwmSZJS1sI/v859E0/gvB6TOPLaQ+KOo2pkyZUkSalp4UJ+ef16CIF7Xjsg7jSqZpZcSZKUesrK+Dj3dzxXeg6/uGY9++5fL+5EqmaWXEmSlDry86FnT6J6mdww+XLaNf2C2+5tGXcq1QBLriRJSh25uTBjBs9zHp9wBL9v+nuaNIk7lGqCJVeSJKWOvDzWRQ24jXs4mPFc+vkf4k6kGuI+uZIkKTWUlkL9+ty97jYWsQ/PhktIy+kWdyrVEGdyJUlSarjpJv6y7mJ+y51cHJ7h2B4F3tEsiVlyJUlS8nvmGZ760xqu4S+cdho8ueFimDYNunSJO5lqiCVXkiQlt/HjGfzD97mSJ/nuieW89BJkZsYdSjXNkitJkpJXYSEvnvQkl5Q+yYCjyxjynzQaNIg7lGqDJVeSJCWnsjJeHfgnLih8mCN7r+O1tzJp1CjuUKotllxJkpSU3jzvn5wz+X/o23kVb4xq6n64KcYtxCRJUvLIz4fcXN6b0YEzotc4sOUS3hrXmWbN4g6m2uZMriRJSh65uXw2ownfi/5DN2bxbusLaOlde1OSM7mSJClpFM9cwTnRf2nDCt7jBFrlF8cdSTGx5EqSpKQQlazmcp5iKe35kGNok1YE3XPijqWYuFxBkiTVfaWlPHD4C7xWnsv9LX/HYenjISfHO5qlMEuuJEmq26KI0Wf9kVvzruCsvvO4ruguKCvzjmYpzpIrSZLqtMI7H+bcoReyb4sS/vZ+Z0KIO5ESgWtyJUlSnVX+7HNc8rturEhrx5j3MmjePO5EShTO5EqSpLpp1Cj+cNk03uRUHnwQ+h7iFK6+4kyuJEmqW/Lz4bvfZdTc9tzBcM797kquuc7NcPV1zuRKkqS65eSTWTF3NefxAvsxlycWnOQ6XH2LM7mSJKnuWLGCKbPrczVDWElL3uQUms2ZFncqJSBnciVJUp0w6YPVnNVtCgcxhSn04imuoHfaVOjePe5oSkCWXEmSlNAmTIDvDyqjz7HNGFZyCL/63kQWdD+J89Nf9IYP2i6XK0iSpIQ0bhzcdVdFh22RsZ7fhAe4/p99aXFxLjAm7nhKcM7kSpKkhPPLX0K/fvDhhxH/230w88s68ut/dqksuNLOOZMrSZISyv33wz33wA+vKOeBVT+k2Sv/gMceg4svjjua6hBLriRJShhPPw033ww/aPYmjz+VSzqb4dZb4Zpr4o6mOsblCpIkKSH8979wxRUwsPEYnl59RkXBDcELy7RbLLmSJCl2Y8bAOedA794RQ744ifpsrHggiiAvL95wqpMsuZIkKVbTp8Npp0GHDhFvdrmWpqxhyy3M0tLcB1e7xZIrSZJis2gRnHQS1K8f8Xa3n9HmpcfgxhuhRw9IT3cfXO02LzyTJEmxKC6uKLirV0eM7PtzurzxKNx3H9x0U9zRlAQsuZIkqdatXw+DBkF+fsRbvW+jz4iH4JFH4Npr446mJGHJlSRJtaq8HC49ey0fj2nEi5zDgE9fgd/9zoKrauWaXEmSVKt+/Wt48Y0m/IFbOItXKi4yGzw47lhKMpZcSZJUa555Bn77W/ghT/IL/ljxSbcJUw2w5EqSpFrx4Ydw5Q/LOa7+aB7jpwS3CVMNsuRKkqQal58PZ+SWsu/mfF5qejmZ//m324SpRnnhmSRJqlGrVsGg/mvYvKqU1/e5hqz334D99oPTT487mpKYJVeSJNWYsjL4wVGLmb24Le90v4VuowZDmzZxx1IKcLmCJEmqfvn5RAf05Gf1/sy7Mzryl24PctzY+yy4qjWWXEmSVP0GDeKRGQP5M9dwM3/gh+n/gCZN4k6lFOJyBUmSVL0KC/nvjM7cyIOczn+4m1/CrBB3KqUYZ3IlSVL1mTqVSb0v4Vz+RR8mMpgLSU/DLcJU6yy5kiSpegwdytLDz2DQ50/SIiudoV1/QeP0jW4Rpli4XEGSJO2ZKIL77uOLW/+X3EafsbJBOz4clkb7PiPiTqYUZsmVJEm7Lj8fcnMrbsfbpAmbS9ZwYYePmbgsh9deC/TpE3dApTpLriRJ2nW5uTBzJpSXQ0kJtzV8lFeXHMrDD8Npp8UdTnJNriRJ2h1fFlzgCa7i/vU/5dpr4Wc/izmXVMmSK0mSqi6K4N57txTcdzmBn/AYpzQZxYMPxpxN2orLFSRJUtWsWgWXXQavvkp08im8MXUfLlh8LwfUn8u/PtqbDFuFEogvR0mStHMTJ8LZZ8OCBYz9+XPcMuE8hi8OdOsGr7/bnKb7xB1Q+jqXK0iSpG3Lz4eePSE9Hfr2JX9VFucPWMqhD5zP1KmBRx6BKVNgHwuuEpAlV5Ikbdtpp8H06RSWt+T66EFyij7ktdHZ/M//wJw5cO21kJkZd0hp21yuIEmSvm38eJbOLOFJ7uSP/IK1NOGH4Sl+M/sq2rePO5y0c5ZcSZK0xcb15bz2o9f5++B6vM0iykknl9e4J9zOAT0iaH9V3BGlKrHkSpKU4qIIxo+Hvz/6Bc89u5mVZd+jY8NCbjtvBZeNuoKu89+F7t1h6NC4o0pVZsmVJCmFvfUW3HJLxQVk9UnnjPQ3uPz6xgy8/xTSMwLwZtwRpd1iyZUkKUW9+/QyTr+sFV2iuTzGnziv+0RavvqPillbqY6z5EqSlII+/BC+f1lzcqIZDOc4ssIqSMux4CppWHIlSUox40dv4LSBER2ixbzDd8liJUTArFlxR5OqjfvkSpKUQqY/N5HvHrueFptW8F7zs2ibVljxQFqas7hKKpZcSZJSwfr15F91NydemE29qJRhf1/EPuNfhZycijua5eS4e4KSissVJElKZvn5cMIJLJm3kYF8yIbMZowcGdj/iGMqHp82Ld58Ug1xJleSpGS1Zg0ceigr5q3lBN6jiFa83eGHHHhEk7iTSTXOkitJUjJ69VWiHgcwprgb3+UdFrAv/+U0+i18Je5kUq1wuYIkSclk0SI2/PQX/GtoQx5p+Bbj6ElzVjGEM/hO2mjonhN3QqlWWHIlSarL8vMhNxfy8ljcshd/WXUeT5Q9SgFt6LFvxGMXFHDxc6fSZPaEioLrxWVKEZZcSZLqstxcxkxvzoM8xyuFZ1JOGoMGrudnt8HAgYEQsuHOz+JOKdU6S64kSXXUqNdW8b/TH2YYJ9CCldzAQ/wk7XG6vDc77mhS7LzwTJKkOiSKYPh7mxnQdTH9T2/BFHpxPzexiL25P+1WuuRkxh1RSgiWXEmS6oAognffhWN7r+L4E9OZNSeNh7o9xrz/TOYXB7xJk/QN3tBB2orLFSRJSmT5+Xxw/K+5ZcFP+Zgj6MgaHs16gB/+uR8NzrkGQoDTvaGD9E2WXEmSEtTnn8PNB0/l2dXPsDcL+Qs/5rLW/6X+wtnQsGHc8aSE5nIFSZISTFkZPPTHzXTvvJF/rz6JO/gtM8nhxzxB/ZWfW3ClKnAmV5KkBDJqZMS1l65myoLmnMS7PNLgFrpunFqxKDctDbp3jzuiVCfsdCY3hLB3CGF4CGFGCGFaCOH62ggmSVLKyM9nWbf+XByepf+AQMmCVbzS/lreHLKRrlOHQI8ekJ7uhWXSLqjKTG4Z8IsoisaHEJoC40II70ZRNL2Gs0mSlBLeOua3nL/sP6yjEXfwW25v93cazZ8J9epVHDDNC8ukXbXTkhtF0TJgWeX7a0IIM4AOgCVXkqQ9EM1fwP0/+ITblv2VA5nKi5xDN2ZDQfpXBVfSbtmlC89CCJ2Ag4FPtvHYj0IIY0MIYwsKCqopniRJSWjZMtZffSMXdRnNLZ/9gLMyhzI6HFNRcF13K1WLKpfcEEIT4GXghiiKVn/z8SiKnoiiqF8URf2ys7OrM6MkSXVbfj707AkZGdC6NYs6fYfvPH4Rz0fn8rtbVvGv6QfRuMc+rruVqlGVdlcIIdSjouAOjqLolZqNJElSkjn1VJg1C6KID4tyOCu8wvrGrXn1+TRyc1sALVx3K1WzquyuEIC/ATOiKHqg5iNJkpQkVqyAW2+FvDyIIv7KlRzP+zSPVvHJZ2nk5sYdUEpeVZnJPRq4GJgSQphY+bnboyh6o8ZSSZJUB61eDQsXwpLJRSz5xzssHj6HJWX7sSTjTRaWtWcKB3ESb/N899/QsseYuONKSa0quyt8CIRayCJJUp20YQP8+vpV/PGJJmwmA2gFnA9AdlYZHdptZp+FY7jwi19yU87rpL/+aqx5pVTgHc8kSdoDn34Kl52/kRn5LbicpziJt+nAEjrsW4/2ecOpXz+Div/dDqh8uzvOuFLKsORKkrQbNm6E3/xwEX8Y3J4OLOdtruS7vPvVAYvToX58+aRUt0v75EqSlLK22gbss7aD6NtiLvcM3pvL6z/PlNue47vdF1bscQvudSslAEuuJElVMWgQG6fP5Y7Nd3Hkiv9QsqE+b1wzlCeLz6T53bfBG29U7HHrXrdSQnC5giRJ27BpE8yYAZNHrmTy81OZPONBxtOXQrK5nKd4IO1mWjxW9NUTunRxr1spgVhyJUkC1qyBV1+Ft96CyZNhxvRyyjanAS2pz6H0TJvBoPLXOY8XOCntvYrZWkkJy5IrSUpZmzZVlNrnnoPXXi1n/YY02rGMQ9ImMqh8Ar0bzOKgH3Sn6+0/IKNec8i9v+LGDt1djiAlOkuuJCmlbN4MH3xQUWxfeglWroTWWZu5nH9yAX/jSMaQVh5B27YVt+Jt1uyrJ7scQaozLLmSpJTx3//CtdfC/PnQuHHEGUct54L1T3HCmP+j3uYNXz+4sPDrBVdSneLuCpKkpFdUBBdfDIMGQeOls3k+XMCKTS155t29OGXa/dS74aew//5uASYlEWdyJUlJ7eWX4Sc/iSguivhV/fu4feOvqM8mKAXat4c5c6BhQ/jJTyA3t3LNbXfX3Ep1nCVXkpSUli+Hn161kZeH1qdvg+m8s/kCem+e/O2DGjaseN8twKSk4nIFSVLyyM8nOqAnz6ZdwgF7FfP60Ii7uY1Pel1F779eV7Htl0sSpJRgyZUkJYf8fEb2vZHjZzzKxdHTdI9mMrH5AG6bdAEZn46GK6+suPLMu5JJKcHlCpKkumvNGqJ/v8iIP03mrsnfZySv0o5lPMpPuZq/kL42wEEHfXW8SxKklGHJlSTVHfn5FVsk5OURNWnK++uP5K7SX/IBV7BX0zU83OD/uKrwHhpG61yOIKU4lytIkhJfFMH48XDooZTPmMnb5SfwndWvc0Lpm+RnH84jf4rIX9GUn318IQ17dHI5giRnciVJCSY//6utvDp3htxcil8fzTuzO/EGD/IWJ1NAGzqyiP8XruWKhY/SoEHlc12OIKmSJVeSlFhOOYXyWXOYSG/emHMqbz54Ch9zH+Wk0yp9JSdvfoPTeJ0zw3+o36MLNNj5KSWlHkuuJCl2a+cV8NnDoxnzWgFj5t3PxxxBIdkA9OMz/udX6ZxyChzaaiXp3/+9N2yQtFOWXElS7crPZ+nJV/DunM6MaXgcH3MkU9Z1oZzTAeiePptBm/9Lf0ZwSnibtj2y4K4vlyC4HEFS1VhyJUm1YtEiePnvq3npt4WMLn2fiDSarSvh8PAp/3PMIo44vzOHn9uJrJJ0yL3P2VpJe8SSK0mqfpUXj82buZGXs67kpcwL+WTp3kAzDqI+d/FrTudVDmQqaWlp8EHZV89t5WytpD1nyZUkVZsoggnjI4aeOJTXVj7NeA6BQjiEsdzd7hXO+kE6Xf/7EMybB+Xl7mUrqcZYciVJe2TdOhj2zmZef2oFr7/fiKVfNCdwHUfwMfdxE2fxMp3TF8Gyytna60/9aoswlyNIqiGWXEnSLls/fR4vnvA4Ly47hvcYyAYa0pTGnBTeZVCv+Zy6/O9kF87YarY256snu5etpFpgyZUkVdn0txby+N3FPD1qX1ZxD52Yx494gtzG73PsXy8mc9BJ0LQp5J/hbK2kWFlyJUnfttVdxzZ02I+XOt/M45/24cP1/cikLWfxMj/mcY5lFAFgQzqc/+pXz3e2VlLMLLmSpK9ZM2Mxc469mtmFPRnDlTy98BKKF7Zi/4ZLuO97H3Dp7R3IvuJ3MHMmlOPFY5ISkiVXklJUaV4+w757LxMXtWJ2o97MbnIws4ta8nlZR+AdADIo5QyG8OO0Jzlu7TukpXWoePLQoS5HkJTQLLmSlEKiCMaPi3j6oSKef64ZBdHjALT7Yhldv5jNKe3z6dqnMV3H/YuuKz5i/2gWjdM2QE4OpG11IpcjSEpwllxJSlZbratdtM/RDO7xW57+qAszSjqQSVO+x2tczDMcx3CashbS02FJ5TZf+R0gdzjkbazYGcGZWkl1jCVXkpJQWeEqJh71Mz5c/l2G8ieGzzuOaF4aR2d8zOP93uCcixvQ8tH/g7lzt31TBmdqJdVxllxJqssqZ2tXz1zKx61O46O9z+PD2W35eM0BrON1APZnNr/mLi4Kz7HfxjxIO6LiuYOOdl2tpKRlyZWkOqi0ZB1jnprBm3d8yDvrn2YifSgvSCetYDO9G8/likOncsy8Zzi66DU6RosqZmpzcir+/JKztZKSmCVXkhLZVutql7Tty1vdfsabkzvwXvHBlHAIGRzEUYzmf/gtx/Ahh6eNpdnalZXPza5cV5vuTK2klGPJlaQEVLa8iBn/msyEO15kwtoreZ/jmby0NyyF9vVWcPZBszllUDonvPhjms8d/9W62hxvnytJYMmVpHjl51M66AzG5zVmQrNjGd/qu0xY0oYpG/ZnI8cBx9GA9RzOJ9zDrZyS9g69Nk4ghDYVz//hv1xXK0nbYMmVpJq21ZIDuncnevkVZk4p5d0Xinj3tXWMKPuQtTSFVdByVTEHZy3k2n4TOXhACw5+7ma6zXubjKj0q5nasNW5na2VpG2y5EpSNSkvh7Kyr95KSyvfP/EyNuRv5BN+wLvTT+TdHo1ZQkegYueDi3mGgQyjH2PZJ20Joajsq5Ne/rAztZK0Gyy5krQHNmyAfz+6gkfvXM5nG3pt56hRW97LooiBDOPEH3XhxEvb0+mq78PMmdteUwvO1ErSbrLkStJuWLQI/vxYxF+fKKewuA05FHEHv6UR68igrOKtaUMyOu9DxrzZZKxZSS8m0zdMJL1HN3i8srgOHepMrSTVAEuuJFVFfj7RKacycnZ7Hqn3c/6z6RQAchnKdTzC8bz/1VLZtLSKO4ntuy+E8K01uV8rss7USlKNsORKEnytiK7tejALfvsM82dtYsFnK5g/Yx0L8jYyKfoPeeSQtamIm9Me4OozltPp+C5w/zxYkPb1JQedOn11bousJNU6S66klLd++WqGHP4w/yh8kPH0pWhmazj7q8cz2cg+LKQL+dzKvZzHCzQMpfBS5QViJ5/skgNJSjCWXEmpYeslAx07Ep17HuMnBJ769ECeKzmVVTxMJ+ZxNi/Rifl0YgH7/vFndPpuN9oe0Iq0Xt//+gVi3b3pgiQlMkuupKSx+KMFvHb202xeXkhW+/q0uuESWlFE1sKJtPr7/TRfu5hishi84Hs89YfzmEQfGqRt5KyDZnHF5z9hQMGLpEWbv1py8PMjvzq5F4hJUp0Soiiq9pP269cvGjt2bLWfV5K+acUKeOmFUl74x0Y+mNBkh8emU7G8YDMZ9OMzrnjsUM4/H1q0YNsXh3XpUvPfgCRpt4UQxkVR1G9bjzmTK6numDsXTj2VVbMLGNLkYl5ocBnDCnqxmXocwCz+j+f5Af+mJSspJouikE3xH56kqFlnin71MMXLS4GIc8LL9O6xCa7ZaomBSw4kKalYciUllvx8Sk49n8Wz17Mkuw9LjjqHJUsDSxdtZsnSwBIGM4VebFpTny5r5nLrfi9z3nHL6XVCW7jjBZg3D8rLyU4rhpwMuKly7ewJZ7vcQJJSiCVXUu2qXBZQOnMu8zofz8wbHydvThp5E9aTNzedvCVNKYg+qTh2OTCk4t2sjBI6sIAOLOE6HuEH/JtD08YT5mx1C9xDD3U/WkkSYMmVVN22s7a1YNEGRr34OSN/NYJRXwxmGj0pm1sPrq14WhuW051ZnM5MujGLvVlEB5bQIe1z9lqdR8PGzaHn+d4CV5JUJZZcSdXrtNMgL4/lUTYjpx/IyJxhjAxlTNvUDehEI87laD7iVN6gO3nkkEe3v91Ky8O7QdfD4eCrv11kG6dVnNsdDiRJVeTuCpJ2XeUFYOtmL2FOq8OZfeAZzF7SkFlLmzL7i72YTVeW0w6AJqzhmLZz6N+zkP4DAv3+/lPqLZjz9RK79eyruxxIkqrI3RUk7ZLi8fOZcuavWbAwUJjVjaJ+J1G0Mo2i5WUUrkyjaHU9ChnOMtpDITCi4nntGqyia70ZnFb6Bj2YzrHhQ/rmrCNj+uSvTn7hmzuejXXJgSSpGlhypVRUOVu6aWY+MzudzOQL7mHK7PpMnpLGlIXNWLK2E/DPimOLIP3tMrIopnXGKlo13kAX8jmUz+jMPLoym65p+ey/8jOaNWsB+W0h90pLrCQpVpZcKVnl51M+6Ht8nlfCvH36M/+y3zBvUQbz8jYx7+PlzC8bykL2YXN+BvwWMtlID2ZwXPpHHMR4ejGF/ZlDawppFtaStqoYmnWrOPe2LgBrVvl1LbGSpATgmlypjlo1cT4Tvn8XsxY2YEXrAyg45WIKSjJZsWgjBSsiCpaUUhhlUUa9rz2vHcvoxHw6M48u5HMgU+kVptHtrT9R74Cu0L499Or17RLrullJUoLZ0ZpcS66UyCrLZNHMAsZ3/B7jz7mbcdMbMn5iGnOXNfraoS1YSTYF33rryGI6M49OYSGdhjxIwwP3+6qgWmIlSXWYJVdKVN8ok9FrQ1m0LouJ7xYy4ZONTPjPQiZsOoCF7LvlKZ3Jpy/jOYRx9GU8PZlGG1aQGcrg0Udh//0r3gYN2n6RtcRKkpKAJVeKy7bKZLt2RPMXsHz8EuZc/whzirOYxgFM4GAm0ociWgMQKKcbsziYCfRlPH0Zz8FMJOup+2G//eDKKyu28nI2VpKUoiy5Uk36Rpnc8OxLFK+pR9H05RTf8UeKi2E5bcinC3Poyly6MJf9+IImW05Rnw0cyFQOZiIHn9mZgw/NoFf/LJpc8QOYNcvZWEmStsGSK+2Jby4pGPwcnxdnMuXDEqZOLGPKGwuZtnF/lrEXRbRiPY22eZr6VGy9tV/3euzfuYz9etRn/3//nv2XjmLfaB710sqdjZUkaRd4MwhpJzbPzmf5aVewam4Razp0Z/Wl17Hmi3RWL17NmtdHsnr991nGXkydfiBTDu5IMa22PLcdyziQqRzIVFpRRBYrybr8dFp1bkbW43eTtXQqraMVtA+fk9aj+9dL7LW3Q+4YyIuge457ykqSVE2cyVVKiObms/SUHzJvbjnzWx3C/CPPY35BY+YvzWR+UVMWrs2ilMwdnqMpqyu222IqB17Wj16HNuDAY7Nofe7A7W+35UysJEk1xuUKSg47Koxz5rD51FwK5pQwO/soZh95CbOXN2P20kbMLmrFnC/asY7GXztdO5bRKSygU6MCOn0xlX1YSBbFNGUNzcJamr74FM26tqXpOSfTdPZ46kcbvMBLkqQEYslV3bCNsrh53y4sW7CJRROLWHT171hcUJ8VZFNINoUZbSlstA+F6xpTWNacYrKISNtyugxK6ZKxkG5NP6fryk/pyiy6kE9n5rFPWEyDZfMgO7uiuPbs6WysJEl1jCVXiSM/n02DzqQor5DCTv0ouuluCtc2oHDBFxT+878Ura3P57RlEXuziH1YQns2f2PpeCYbaU1hxVuLMlo3K6X1wvG0rrz5wf7MoWtaPvuun0lGZmXp3VGJrcxlkZUkqW6x5Kr67KQMls3KZ+Fp17Awv4zFrXqz6NgLWFzcmMXL0llc1IDFhQ1YEbXZ7umbspo2rGDvypq7d5/W7N0xYu/OGez94gPsvXwsLaJiwjeLqiVWkqSUY8nVrtmqEJZ3y6H05dfY1LELm9ZsZOMxA9k0bwmraUo++zG38UHM3fc45hY0Y+7qbBZsbEcZ9b52upYU05HFdMxYTseyeXRgCW1YUTETG4pp/cdf0rprS7JuuIT683azqFpiJUlKOZZcfV1lIVwzcwnz9+3PvGv+wPzCJsybXcq8BWnMn7yaBWUdWEuTbxXWbWkZVrJfg6Xs17yA/T4fzX7MYV8WsDeL6BCW0XjuZNhrL2jQwLWvkiSp2lhyk9F2CmEUQfEns1l07k0sWVTO5y178PmhuSxf04jPC9NZvqo+nxekszxqQwktvnbKRnxBZ+bRmXnsywKaU0Imm8iklMyj+pHZsjGZHwwjc3UhjVlLlzCf/bqm0TLv469O4rIBSZJUSyy5ddE3y+DgwZTVb8zyvFUsnbWWJfc8w9KVDVlCexaxN4vT9mVReicWl7ZhAw2/dbrmrKJdWE7bzFW02ziftiynPUsrSm1YQOcn76B1tyxCu7YVX3d3byVriZUkSbXEkhuX7RW+KKK0aDXFIyZTfP1dFC/byMoWnVl5xCmsWl+flcXlrJz+OSvLmlBEK5bSnqW0Zzltv7ZFFkA6ZXRgCR1ZzN4doWOr9ew96XU6sogOLGEvltE2rZAGqz6HJk0gBGdbJUlSUrDk7okqXOxUPnMWKzv3peCW+yhY34QVCzdQsLSUFa99TMG6xhTSiiJaURxaUZSWTdHmFqyh2Q6/bBPW0JKVtKKosuIuo/33D6P9Phm079KA9g/cxF6LPqVttIz0NNxpQJIkpRxLbhWvyi/vlsMXz7xCSVpLVi8qoWTJWlb/6n5KCkspIosC2lCQ2YEVrXIo+KIRBasbUFC5Y+s393L9UgtW0npLzS0ia78sWrUsp1VryHr7eVpFhbRkZcVb2mpaLp9J8+ZQr88eFFVLrCRJSgHJX3K3Veo6d6Zs1VrWL13JhtPO4osFhXxOW5bQkSWNu7EkZyBLihqwZOFmlpS3YzltWUPTby0H+KYWrCS7/mqyG6ylTckssitvQJBNAW1CIdn33kSbLk3I3q8Zrc8/kcxZU7dfVN1pQJIkabclf8nt2ZOLp9/GB3yH9TTc8raj7a8y2Uj7jBV0KFtAB5bQluW0YBXNWEPz80+lWXYmzds1otmjv6fZsjxaRQW0DsXU67F/9S0LsMhKkiTtth2V3G3/jr2uycsjh5mkUV5ZbzfQ8IjeNGxWj4bNM2n47ms0XLWMtnxOh7CMDvs3pHXeR4SwN/Q8+dtF9bn7vjr3ufdWFtGCr4rol4YO/XZJ3VqXLl8vvd+0s8clSZK0W6o0kxtCOBl4GEgHnoyi6J4dHR/HTK7rVyVJklLLHs3khhDSgf8HnAgsBj4LIbwWRdH06o25B/ZkRtXZVEmSpKRTleUKhwFzoijKBwghvACcDiROybWoSpIkaSs73kqgQgdg0VYfL678nCRJkpSQqlJywzY+962FvCGEH4UQxoYQxhYUFOx5MkmSJGk3VaXkLgb23urjjsDSbx4URdETURT1i6KoX3Z2dnXlkyRJknZZVUruZ0DXEELnEEImcB7wWs3GkiRJknbfTi88i6KoLIRwLfA2FVuIPRVFkVd5SZIkKWFV6WYQURS9AbxRw1kkSZKkalGV5QqSJElSnWLJlSRJUtKx5EqSJCnpWHIlSZKUdCy5kiRJSjqWXEmSJCUdS64kSZKSjiVXkiRJSceSK0mSpKRjyZUkSVLSCVEUVf9JQygAFlT7iXeuNVAYw9etqxyvXeN47RrHa9c4XrvG8dp1jtmucbx2TVzjtW8URdnbeqBGSm5cQghjoyjqF3eOusLx2jWO165xvHaN47VrHK9d55jtGsdr1yTieLlcQZIkSUnHkitJkqSkk2wl94m4A9Qxjteucbx2jeO1axyvXeN47TrHbNc4Xrsm4cYrqdbkSpIkSZB8M7mSJElS3Sy5IYSTQwh5IYQ5IYTbtvF4TghhTAhhYwjhpjgyJpIqjNeFIYTJlW+jQwi948iZKKowXqdXjtXEEMLYEMIxceRMFDsbr62OOzSEsDmEcHZt5ks0VXh9DQghlFS+viaGEH4VR85EUZXXV+WYTQwhTAshjKztjImkCq+vm7d6bU2t/DuZFUfWRFCF8WoeQhgaQphU+fq6PI6ciaIK49UyhDCk8v+Rn4YQDowj5xZRFNWpNyAdmAt0ATKBScAB3zimDXAo8Dvgprgz14HxOgpoWfn+KcAncedO8PFqwldLfQ4CZsadO5HHa6vj3gfeAM6OO3cijxcwAHg97qyJ8FbF8WoBTAf2qfy4Tdy5E3m8vnF8LvB+3LkTebyA24F7K9/PBoqBzLizJ/B43Qf8uvL9HGBYnJnr4kzuYcCcKIryoyjaBLwAnL71AVEUrYii6DOgNI6ACaYq4zU6iqKVlR9+DHSs5YyJpCrjtTaq/BsMNAZSeWH7Tser0nXAy8CK2gyXgKo6XqpQlfG6AHgliqKFUPHvfy1nTCS7+vo6H3i+VpIlpqqMVwQ0DSEEKiY4ioGy2o2ZMKoyXgcAwwCiKJoJdAohtK3dmF+piyW3A7Boq48XV35O27ar4/VD4M0aTZTYqjReIYQzQggzgf8CV9RStkS00/EKIXQAzgD+Uou5ElVV/z4eWfnr0TdDCD1rJ1pCqsp4dQNahhBGhBDGhRAuqbV0iafK/96HEBoBJ1Pxw2eqqsp4PQr0AJYCU4Droygqr514Cacq4zUJOBMghHAYsC8xTpzVxZIbtvG5VJ5J25kqj1cI4TgqSu6tNZoosVVpvKIoGhJFUQ7wfeD/ajpUAqvKeD0E3BpF0eaaj5PwqjJe46m4TWVv4BHgPzUdKoFVZbwygEOA04CTgDtDCN1qOliC2pX/P+YCH0VRVFyDeRJdVcbrJGAi0B7oAzwaQmhWs7ESVlXG6x4qfuicSMVv8CYQ48x3RlxfeA8sBvbe6uOOVPyEpW2r0niFEA4CngROiaKoqJayJaJden1FUTQqhLBfCKF1FEWpeI/zqoxXP+CFit/20Ro4NYRQFkXRf2olYWLZ6XhFUbR6q/ffCCE85utri229vhYDhVEUfQF8EUIYBfQGZtVOxISyK/9+nUdqL1WAqo3X5cA9lUvU5oQQ5lGx1vTT2omYUKr679flAJVLPOZVvsWiLs7kfgZ0DSF0DiFkUvEX9bWYMyWynY5XCGEf4BXg4iiKUvF/DFurynjtX/mXlxBCXyoW4KfqDwY7Ha8oijpHUdQpiqJOwEvAT1K04ELVXl/ttnp9HUbFv9O+vrb/7/2rwHdCCBmVv4I/HJhRyzkTRZX+/xhCaA70p2LsUllVxmshMBCgcm1pdyC/VlMmjqr8+9Wi8jGAK4FRW//gXtvq3ExuFEVlIYRrgbepuNLvqSiKpoUQrq58/C8hhHbAWKAZUB5CuIGKKwBjG+i4VGW8gF8BrYDHKv/fWhZFUb+4MsepiuN1FnBJCKEUWA+cu9WFaCmliuOlSlUcr7OBa0IIZVS8vs7z9bX98YqiaEYI4S1gMlAOPBlF0dT4UsdnF/4+ngG8Uzn7nbKqOF7/B/wjhDCFil/X35qiv1Wp6nj1AJ4OIWymYteTH8YWGO94JkmSpCRUF5crSJIkSTtkyZUkSVLSseRKkiQp6VhyJUmSlHQsuZIkSUo6llxJkiQlHUuuJEmSko4lV5IkSUnn/wMagoa2iztXLgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 8))\n",
    "plt.plot(np.linspace(0.1, 0.9, 81), theoretical, color=\"red\", marker=\"o\", label=\"theoretical\", markersize=4)\n",
    "plt.plot(np.linspace(0.1, 0.9, 81), actual, \"b\", label=\"actual\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f7b8f71",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}