notebooks/BANDIT_General_linear_improved.ipynb

{
    "cells": [
        {
            "attachments": {},
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "今回は線形関数近似の探索でよく出てくる[Improved Algorithms for Linear Stochastic Bandits](https://papers.nips.cc/paper/4417-improved-algorithms-for-linear-stochastic-bandits)を学びます。\n",
                "\n",
                "前回も線形バンディットをやりましたが([BANDIT_contextual.ipynb](BANDIT_contextual.ipynb))、詳しい証明やよりよいバンディットについては触れてませんでした。\n",
                "線形バンディットやLilnear MDPでタイトな性能を達成するためには、今回紹介する\"self-normalized bound for vector-balued martingale\"ってやつが大事になってきます。\n",
                "\n",
                "--- \n",
                "**定理：Self-Normalized Bound for Vector Valued Martingales**\n",
                "\n",
                "$\\left\\{F_t\\right\\}_{t=0}^{\\infty}$をフィルトレーションとする。 \n",
                "$\\left\\{\\eta_t\\right\\}_{t=1}^{\\infty}$を$\\eta_t$が$F_t$-可測な実数値確率過程とする。また、$\\eta_t$は条件付き$R$-sub-Gaussianとする ($R \\geq 0$)。つまり、\n",
                "$$\n",
                "\\forall \\lambda \\in \\mathbb{R} \\quad \\mathbf{E}\\left[e^{\\lambda \\eta_t} \\mid F_{t-1}\\right] \\leq \\exp \\left(\\frac{\\lambda^2 R^2}{2}\\right)\n",
                "$$\n",
                "\n",
                "$\\left\\{X_t\\right\\}_{t=1}^{\\infty}$は$X_t$が$F_{t-1}$-可測である$\\mathbb{R}^d$の確率過程とする。\n",
                "$V$を$d \\times d$の正則行列とする。任意の$t \\geq 0$について、\n",
                "$$\n",
                "\\bar{V}_t=V+\\sum_{s=1}^t X_s X_s^{\\top} \\quad S_t=\\sum_{s=1}^t \\eta_s X_s .\n",
                "$$\n",
                "すると、任意の$\\delta>0$について、最低でも確率$1-\\delta$で、任意の$t \\geq 0$について、\n",
                "$$\n",
                "\\left\\|S_t\\right\\|_{\\bar{V}_t^{-1}}^2 \\leq 2 R^2 \\log \\left(\\frac{\\operatorname{det}\\left(\\bar{V}_t\\right)^{1 / 2} \\operatorname{det}(V)^{-1 / 2}}{\\delta}\\right)\n",
                "$$\n",
                "が成り立つ。\n",
                "\n",
                "---\n",
                "\n",
                "ここで、マルチンゲールのぶれ$\\left\\|S_t\\right\\|_{\\bar{V}_t^{-1}}^2$が$\\bar{V}_t^{-1}$で測られていることに注意しましょう。\n",
                "つまり、マルチンゲールのぶれ自体がマルチンゲールから出てきた量で測られているのでself-normalizedと呼ばれるんですね。\n",
                "証明はちょっと厄介なので論文を参考にしてください。(式展開を追うとなるほど！となりますが、分かりやすい説明はTODOです。)\n",
                "\n",
                "これを確認してみましょう。"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 102,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "text/plain": [
                            "<matplotlib.legend.Legend at 0x7f83be65fa30>"
                        ]
                    },
                    "execution_count": 102,
                    "metadata": {},
                    "output_type": "execute_result"
                },
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABXdElEQVR4nO3dd3wUZf4H8M/29A3pCSQkJIHQewmggKCAFfQ85VBAUU8OThErp9gVPBueovzUE7xTingIVhRCk15DlZ4QQnrbTd06vz8mWbIkIbvJbmaTfN6v175md8rud4cyn33mmWdkgiAIICIiIvJgcqkLICIiImoMAwsRERF5PAYWIiIi8ngMLEREROTxGFiIiIjI4zGwEBERkcdjYCEiIiKPx8BCREREHk8pdQGuYLVakZWVBX9/f8hkMqnLISIiIgcIgoDS0lJERUVBLr92G0qbCCxZWVmIjo6WugwiIiJqgkuXLqFTp07XXKdNBBZ/f38A4hcOCAiQuBoiIiJyhF6vR3R0tO04fi1tIrDUnAYKCAhgYCEiImplHOnOwU63RERE5PEYWIiIiMjjMbAQERGRx2NgISIiIo/HwEJEREQej4GFiIiIPB4DCxEREXk8BhYiIiLyeAwsRERE5PEYWIiIiMjjMbAQERGRx2NgISIiIo/XJm5+SERE1C4IAiBYAasFsJoBwSI+r5lX89pqFh8WE2A1VU9rvzbXml/7tVGcZzZUPzdWPzcBMhkw/g3JvjoDCxERkTOsVsBYBhhKaz301Y9SwFguLjeWA8YKcWoqB0yVVx7mSjEk2AJHTdC4KnQI1urwUSugSEWhYWAhIiJyO0EQWwyMtcND9XNjufioKADK8oHyPKA8H6jSXxVMSgFjqdTfpGEyOSBTAHIFIFcBCmX1VAXIldXTq+df9VqhBpQacVr7udJL0q/GwEJERJ7DahVDRE2AsLVU1DyvAEw1j8qGw4ep8sp6xoorr13ZQiFXApoAQONfa+onTlU+gNoPUPuKD5UPoPK+MlV6VYcFRa2AUeu5TCG+v7ye1zIFIJeLr23bKqvDisx138/DMLAQEZF7CQJQWQyUXARKMoDSXLH1wvYoqG7RKBBPq7QEuao6VNQECV9x6hsC+IYCfmHi1CuwOojUflSHE6WmTQcET8PAQkRETScI4mmSyiKgNAfQZQK6S9XTTKDkkhhUjGXOva9Mbt9CofYVQ4Xap7qVwsf+ucq7eh1vx5YrVO7ZH+Q2DCxERCQyVYnBo6LIflpZXP28uP55jp5m8QsHAmMA/0ix9cI3VGzRqGnNsLVo+ImnTNh6QbUwsBARtQdWK1CSDhRfFFs+9JerW0IuV7/OEvuANJXSSwwk2mhA2wnQdgQCOgKBncWQEhgttm4QNREDCxFRW2ExAVU6MXwUpwFFF4DC80DeH+LDkUAiUwA+QYB3B8A7qPp5EODToe487w5X1mUYITdjYCEiag6zUeyfYa66MtiW3dQgnmoxVz9MlbWeV4njcZgNV+abql/XPLcY6w4KJljEFpPar40VjQcSpRfQIU5s/dB2AgI61XreUTw9owngqRjySAwsROT5rFbxoG33sNR6LlwZ7bM5y41lYp+M2n01KoqAqpLqAcCqL7E11XpuNUu9d+ryCRaDSVAXICgOCE0CwnsBwfHiJbBErRADC1FbZjHbj2NhqrAftrv20N6159V5Xd+8a2xnMdqPg1Hz3FI9FHjN8N8WY63WgobChFXqvegYmcJ+gC2lWhwZVFnz8AZUXtXLvKqfXz3P+6p1veuO1WE3rTVImMq7usNqgDgIGFEbw7/VRK2FoQwoOA3knwHKcsQxKyqKxHErjGXi8qsH27IYpK66hcnEy2HliupBtORXRv6UObBM5XOlT0btfhzegXUvsVX7VV86W/2al8kSuRUDC5GnEQSgOB3IOgTkHAfyTgK5JwFdRtPfU66sPsB6iwNm1YyMaXtc/VreyPLaI282sE7tAbnUPmJrgbJ6qG/bEOCqWiN0XvWoEyyutawmdLDvBVFbxcBCJDWrFcg9DqT/DqTvADL2iH0n6uMbBoR2EztJ+gSLjzq//v3qTpXqFv1KRESuxsBC1NKsVrHVpCagpO8QO3XWJlcBEb2AyL5iZ8mwHkBYd/E0BRFRO8TAQuRuVguQcwzI2A1c3Amk76zbgqL2A2KSgbjrgM4jgIjeYudLIiICwMBC1Hxl+eIgXRWF1R1hC8Vb1OuzxYG7Cs7WvR29yheIGQbEjgTirgci+/HKDiKia+D/kETOqGktqWkpydwn3m22MZoAIHoo0DkZiL0OiOrPq0qIiJzAwELUGN1l4OxvwNmNYn8Tg+6qFWRiJ1jfULETrG+IOPULEwfuCk4AQrpywC4iomZgYCG6mqFUDCYXtgLnt4hjn9Sm9hdbSjqPEB/hPcSrcYiIyG0YWIisViDnCHDmN+DCFiBz/1XDrcuAToOAruOB+LFARB/2NyEiamH8X5far8oSYN9nwP7PxZFja+sQC3QZA8SPEfuc8HJiIiJJMbBQ+1N4Xgwph/575eodla8YThLGikElKE7aGomIyA4DC7UPVivwx/fAgX8DaduvzA/tDlz3JNDjdo57QkTkwRhYqO0rOAd8/3cgY1f1DBmQeBMw9BGgyw3ifXOIiMijOf0/9fbt23HbbbchKioKMpkM69ats1s+Y8YMyGQyu8eECRMafd8lS5YgNjYWXl5eGDp0KPbt2+dsaUT2LCbg9/eAT4aLYUXlC1z3FDD3KDD1GyBhHMMKEVEr4fT/1uXl5ejbty+WLFnS4DoTJkxAdna27bFy5cprvufq1asxb948vPTSSzh06BD69u2L8ePHIy8vz9nyiET5p4HPbgBSXgEsBiD+BuBvu4GxC4DAGKmrIyIiJzl9SmjixImYOHHiNdfRaDSIiIhw+D3fe+89PPzww3jggQcAAEuXLsVPP/2EL774As8995yzJVJ7d+xb4PvHAFM54BUITFgI9J0CyGRSV0ZERE3klvbwrVu3IiwsDN26dcOsWbNQWFjY4LpGoxEHDx7EuHHjrhQll2PcuHHYvXt3vdsYDAbo9Xq7BxEAYMf7wP9mimEl7npg9l6g318YVoiIWjmXB5YJEybgP//5D1JSUvDWW29h27ZtmDhxIiwWS73rFxQUwGKxIDw83G5+eHg4cnJy6t1m4cKF0Gq1tkd0dLSrvwa1NoIAbHkT2PSy+HrE48D96wB/x1v6iIjIc7n8KqF7773X9rx3797o06cP4uPjsXXrVowdO9YlnzF//nzMmzfP9lqv1zO0tGdWK7DhWWDfp+LrsS8B18279jZERNSquP0SiS5duiAkJATnzp2rd3lISAgUCgVyc3Pt5ufm5jbYD0aj0SAgIMDuQe2U1Qqsm1UdVmTAxLcZVoiI2iC3B5bMzEwUFhYiMjKy3uVqtRoDBw5ESkqKbZ7VakVKSgqSk5PdXR61dhsXAEdXAXIlcNfn4tgqRETU5jgdWMrKypCamorU1FQAQFpaGlJTU5GRkYGysjI8/fTT2LNnD9LT05GSkoI77rgDCQkJGD9+vO09xo4di48++sj2et68efjss8/w5Zdf4o8//sCsWbNQXl5uu2qIqF6nfwF2V/89mrQU6P0naeshIiK3cboPy4EDBzBmzBjb65q+JNOnT8cnn3yCo0eP4ssvv0RJSQmioqJw00034bXXXoNGc2XY8/Pnz6OgoMD2+p577kF+fj5efPFF5OTkoF+/ftiwYUOdjrhENrrLwLq/ic+HzQb63C1tPURE5FYyQRAEqYtoLr1eD61WC51Ox/4s7YHVAnx5O3BxBxDZF5i5kfcBIiJqhZw5fnNccmp9ti4Sw4raD/jTMoYVIqJ2gIGFWpfTvwDb/yk+v/V9IDhe2nqIiKhFMLBQ61F4Hlj7V/H5kEeAPn+Wth4iImoxDCzUOpRcAr66EzDogOihwE1vSF0RERG1IJePdEvkcsUXgS9vBUoygA6xwN1fAkq11FUREVELYgsLebbSXGB5dVgJ6gLM+BkIqH8QQiIiarvYwkKey2wAVk8FdAwrRETtHVtYyDMJAvDTk0DmfsArEJj6LcMKEVE7xsBCnunwV8Dh/wIyOfCnL3j5MhFRO8fAQp7n0n7g56fF52OeBxLGSlsPERFJjoGFPEv+GWDF3YC5EkgcD4ycJ3VFRETkARhYyHPos8SxViqLgY4DgbuXAXL+FSUiIgYW8hTFF8UbGuouAcEJwF/WAGpfqasiIiIPwcuaSXqXDwIr/wKU5QABnYD71gK+wVJXRUREHoSBhaRTXghsfxvY9ykgWICwHsB9/wMCoqSujIiIPAwDC7UcQQDyTgKnfgLO/ApkHQIEq7is553i3Ze9AyUtkYiIPBMDC7mW2SB2ni3LAyoKgcoiQJ8N5B4XA0pJhv36kX2Bca8A8WOkqZeIiFoFBhZyjNkAlOeLQaQ8v/p5rhhO9FmA/rI4Lc+/9vsoNGI46XazOL6KtlPL1E9ERK0aAwvVZTEBp38GMvYC2UfE1pGqEse3V3oBfuGATzDgEwT4hgFh3YHwnkD0UEDj57bSiYiobWJgoSvK8oHUr4C9nwKlWXWXy1WAbyjgFypOfcPEDrIBUUBAR3Gq7QR4dwBkspavn4iI2iwGFgLMRmDH+8Dv7wAWozjPNwzoOQmI6g9E9BaDiFcggwgREUmCgaW9u7AV+OU5IP8P8XVUf2DII0CvuwClRtLSiIiIajCwtFcZe4HNrwHpv4uvfUKAm/8pXl7MVhQiIvIwDCztTfZRYPPrwNlfxdcKNTBwBjB6vthBloiIyAMxsLQX+iwg5VXgyCoAAiBTAP3+Aox6BgiMkbo6IiKia2JgaQ/SfgdW33fl0uRedwFjngeC4yUti4iIyFEMLG1dSQawaipg0AGR/YBb3gM6DZS6KiIiIqcwsLRlFjPwv4fFsNJpMDD9R0DlJXVVRERETpNLXQC50e/vApf2AGp/4K7PGVaIiKjVYmBpqzL2ANsWic9vfR/oECtpOURERM3BwNIWVZaIp4IEK9DnXqDP3VJXRERE1CwMLG2NxQSsmQ7oMsRWlZvflroiIiKiZmNgaUsEAfjhcXG4fZUvcPeXgFeA1FURERE1GwNLW2G1AD88BqR+LQ4K9+cvgah+UldFRETkErysuS2wmIC1jwAn1gIyOXDHR0DijVJXRURE5DIMLK2dqRJYMwM4swGQq4C7PgN6Tpa6KiIiIpdiYGnNDKXAyiniHZeVXsA9X7FlhYiI2iQGltbKYgZW/UUMK2p/4C+rgNiRUldFRETkFgwsrdWml4C07YDaD5i2nvcHIiKiNo1XCbVGx74Fdn8kPp/0McMKERG1eU4Hlu3bt+O2225DVFQUZDIZ1q1bZ1tmMpnw7LPPonfv3vD19UVUVBSmTZuGrKysa77nyy+/DJlMZvdISkpy+su0C9lHge//Lj4f+QTQ4w5p6yEiImoBTgeW8vJy9O3bF0uWLKmzrKKiAocOHcKCBQtw6NAhrF27FqdPn8btt9/e6Pv27NkT2dnZtseOHTucLa3tK7oAfHUXYKoA4m8AblggdUVEREQtwuk+LBMnTsTEiRPrXabVarFx40a7eR999BGGDBmCjIwMxMTENFyIUomIiAhny2k/BEG8fLk8DwjvBfxpGSBXSF0VERFRi3B7HxadTgeZTIbAwMBrrnf27FlERUWhS5cumDp1KjIyMhpc12AwQK/X2z3avNO/ANlHxE62U78FvAOlroiIiKjFuDWwVFVV4dlnn8WUKVMQENDwPW2GDh2K5cuXY8OGDfjkk0+QlpaG6667DqWlpfWuv3DhQmi1WtsjOjraXV/BMwgCsG2R+HzIw0BApLT1EBERtTCZIAhCkzeWyfDdd99h0qRJdZaZTCbcddddyMzMxNatW68ZWK5WUlKCzp0747333sPMmTPrLDcYDDAYDLbXer0e0dHR0Ol0Tn1Oq3HmV2DFnwGVDzD3GOAbInVFREREzabX66HVah06frtlHBaTyYQ///nPuHjxIjZv3ux0iAgMDETXrl1x7ty5epdrNBpoNBpXlNo67PtUnA6eybBCRETtkstPCdWElbNnz2LTpk0IDg52+j3Kyspw/vx5REby1AfKC4DzW8TnA2ZIWgoREZFUnA4sZWVlSE1NRWpqKgAgLS0NqampyMjIgMlkwp/+9CccOHAAX3/9NSwWC3JycpCTkwOj0Wh7j7Fjx+Kjjz6yvX7qqaewbds2pKenY9euXZg8eTIUCgWmTJnS/G/Y2p1cDwgWILIvEJIgdTVERESScPqU0IEDBzBmzBjb63nz5gEApk+fjpdffhnff/89AKBfv352223ZsgWjR48GAJw/fx4FBQW2ZZmZmZgyZQoKCwsRGhqKkSNHYs+ePQgNDXW2vLbn+P/Eaa8/SVsHERGRhJrV6dZTONNpp1XRXQbe7wlAAOYeBwLb+NVQRETUrjhz/Oa9hDzZie8ACEBMMsMKERG1awwsnkoQgKOrxee97pK2FiIiIokxsHiqtG1AzlFA6QX0vFPqaoiIiCTFwOKpdv5LnA6YDvg6f2k4ERFRW8LA4onKC4AL1WOvDP2rtLUQERF5AAYWT/TH94BgBaL6A8HxUldDREQkOQYWT3TsW3HaY5KkZRAREXkKBhZPU3QBuLgTkMmB3ndLXQ0REZFHYGDxNKkrxWmXMYC2o7S1EBEReQgGFk9itQJHqgNL/6nS1kJERORBGFg8Sdo2QHcJ8NIC3W6RuhoiIiKPwcDiSVJXiNNefwJUXtLWQkRE5EEYWDxFlU68nBng6SAiIqKrMLB4il0fAuYqIDQJiBogdTVEREQeRSl1Ae2a1QLs+xT440fg4g5x3nVPAjKZtHURERF5GAYWKf00Dzi4/MrrMc9z7BUiIqJ6MLBIJfOAGFZkcuCGBUD324CQRKmrIiIi8kgMLFI59B9x2uce4Lp50tZCRETk4djpVgqCAJz+RXze5x5payEiImoFGFikUHQBKM8DFBqg83CpqyEiIvJ4DCxSSK++IiiqP6DUSFsLERFRK8DAIoWzv4nT+BukrYOIiKiVYGBpaYIAZOwRn8ePkbYWIiKiVoKBpaXpLwMVBYBMAUT0lroaIiKiVoGBpaXlHBenod0Albe0tRAREbUSDCwtrfCcOA3pKm0dRERErQgDS0uzBRaOaktEROQoBpaWVnBGnAYnSFsHERFRK8LA0pIEAcg9IT4P6yFtLURERK0IA0tL0mcBVSXiFUKh3aSuhoiIqNVgYGlJeSfFaUgiR7glIiJyAgNLS8qtvqQ5vKe0dRAREbUyDCwtqeCsOA3tLm0dRERErQwDS0vSXxangdHS1kFERNTKMLC0pNIcceofIW0dRERErQwDS0sqzRan/pHS1kFERNTKMLC0FGMFUKUTnzOwEBEROYWBpaXUtK6ofAGNv7S1EBERtTIMLC2ldv8VmUzaWoiIiFoZBpaWUtPCEhAlbR1EREStkNOBZfv27bjtttsQFRUFmUyGdevW2S0XBAEvvvgiIiMj4e3tjXHjxuHs2bONvu+SJUsQGxsLLy8vDB06FPv27XO2NM9m63DLK4SIiIic5XRgKS8vR9++fbFkyZJ6l//zn//Ev/71LyxduhR79+6Fr68vxo8fj6qqqgbfc/Xq1Zg3bx5eeuklHDp0CH379sX48eORl5fnbHmei5c0ExERNZnTgWXixIl4/fXXMXny5DrLBEHA4sWL8cILL+COO+5Anz598J///AdZWVl1WmJqe++99/Dwww/jgQceQI8ePbB06VL4+Pjgiy++cLY8z6XPEqe8QoiIiMhpLu3DkpaWhpycHIwbN842T6vVYujQodi9e3e92xiNRhw8eNBuG7lcjnHjxjW4jcFggF6vt3t4PFsLCwMLERGRs1waWHJyxINyeHi43fzw8HDbsqsVFBTAYrE4tc3ChQuh1Wptj+joVjDUPQeNIyIiarJWeZXQ/PnzodPpbI9Lly5JXdK1CQI73RIRETWDSwNLRIR4MM7NzbWbn5uba1t2tZCQECgUCqe20Wg0CAgIsHt4tKoSwFzd6ZgtLERERE5zaWCJi4tDREQEUlJSbPP0ej327t2L5OTkerdRq9UYOHCg3TZWqxUpKSkNbtPq1PRf8e4AqLykrYWIiKgVUjq7QVlZGc6dO2d7nZaWhtTUVAQFBSEmJgZz587F66+/jsTERMTFxWHBggWIiorCpEmTbNuMHTsWkydPxpw5cwAA8+bNw/Tp0zFo0CAMGTIEixcvRnl5OR544IHmf0NPwCuEiIiImsXpwHLgwAGMGTPG9nrevHkAgOnTp2P58uV45plnUF5ejkceeQQlJSUYOXIkNmzYAC+vKy0L58+fR0FBge31Pffcg/z8fLz44ovIyclBv379sGHDhjodcVstXiFERETULDJBEASpi2guvV4PrVYLnU7nmf1Ztr8NbH4d6HcfMKn+AfeIiIjaG2eO3063sFATcJRbImqDrFYrjEaj1GWQh1OpVFAoFM1+HwaWlsDAQkRtjNFoRFpaGqxWq9SlUCsQGBiIiIgIyGSyJr8HA0tLqOl0yzs1E1EbIAgCsrOzoVAoEB0dDbm8VQ7pRS1AEARUVFTY7g0YGdn0vpwMLC2BLSxE1IaYzWZUVFQgKioKPj4+UpdDHs7b2xsAkJeXh7CwsCafHmIsdjerBSirHhSPVwkRURtgsVgAiONoETmiJtiaTKYmvwcDi7uV5wOCBZDJAd8wqashInKZ5vRHoPbFFX9XGFjcreYeQr5hgIJn4IiIiJqCgcXd2H+FiIio2RhY3K2iUJz6hkhbBxERUSvGwOJulSXi1LuDpGUQEZF7XLp0CaNHj0aPHj3Qp08frFmzRuqS2iR2qnC3ymJxysBCRNQmKZVKLF68GP369UNOTg4GDhyIm2++Gb6+vlKX1qYwsLhbTWDxCpS0DCIico/IyEjbgGgREREICQlBUVERA4uL8ZSQu1WViFO2sBARSW7r1q2IjY1123scPHgQFosF0dHRzfoMqouBxd14SoiIyKNlZmbi0UcfRUJCAry8vBAeHo6bbroJx44dc+p9ioqKMG3aNHz66aduqrR9Y2BxNwYWIiKPlZ6ejv79+6OwsBD//e9/cerUKXz77bfo0aMHNBqNw+9jMBgwadIkPPfccxg+fLgbK26/2IfF3WyBJVDSMoiIqK4PP/wQvr6+WL16te0mjrGxsbjuuuscfg9BEDBjxgzccMMNuP/++91VarvHFhZ342XNREQeq7i4GFVVVcjIyGjye+zcuROrV6/GunXr0K9fP/Tr18/p00nUOLawuJPVAlTpxOcMLETURgmCgEqTRZLP9lYpmnWfmjlz5iAlJQVdunTBwIEDMXbsWEybNg09evSwrXPhwgUcP34ct99+e73vMXLkSFit1ibXQI5hYHGnKh0AQXzOy5qJqI2qNFnQ48VfJfnsk6+Oh4+66YeyAQMG4MKFC9ixYwd+++03rFmzBu+++y7Wrl2L2267DQDwyy+/oLS0tMHAsnfvXqSlpdW77N57760z7/Dhw1i0aBESExPx+uuvN7n29oaBxZ1qLmlW+QJK3oadiMgTKRQKjBo1CqNGjcLLL7+M2NhYrFixArfddhu2bduGBQsWIDg4GKtXr8aOHTvqbK/VavHBBx9g9+7dtnnZ2dkYMmQIbrnlFvj7+9ut379/f7z11lv4/PPP3f7d2hIGFneqYIdbImr7vFUKnHx1vGSf7UpWqxUGgwGhoaEAgFGjRqFPnz5Yvnx5g2OvxMfH12lhWbRoEf7+978jPT0d8+fPt81XKBRYv369S2tuLxhY3El3SZwGRElbBxGRG8lksmadlpHK/fffjx49euCGG25AeHg4Lly4gDfffBMAMHfuXNt6GRkZ1xxsTqVSITAwEIWFhQgODkZWVhbWr1+PkydPwsfHBz/++KObv0n7wKuE3Km4OnF3iJO2DiIiqmPAgAH48ccfceutt6J79+549NFHkZSUhNTUVHTp0gWAOKhcVFTjPzqTkpJw+vRpAGLryhNPPAEfH596101LS8Pzzz+PX375BV999ZXrvlAb1/oicWthrABOfi8+D+oibS1ERFTHE088gSeeeOKa61y8eNF2n6BrSUpKwpkzZxAbG4uff/4Zx48fb3DduLg4fP31107X296xhcVdtr0FZB0Sn0cPkbYWIiJqkl69euHChQvo3bs3Tp482eB6NS0sCxcuxFNPPQUvLy/bsuLi4pYotc1jYHGXC1vEadKtQMJYaWshIqIm0Wq1OHjwII4dO2Y3NsvVkpKSsHXrVmzatAkzZ860W/b000+7u8x2gYHFHcoLgZzq5sCJ/5S2FiIisomNjbXrUOuq90hKSsKePXvwzDPPQKVS2eb/8MMPOHfuHJYuXdqszyT2YXG9Uz8BG54DBAsQ2Q/QdpS6IiIiquauwBIYGAhBEOqsGxwcjAcffBDTpk1r1mcSW1hcSxCAn54ESqrvSTF6/rXXJyKiNu3YsWPo3bu31GW0CQwsrlScBpRmi8/nHAS6TZC2HiIiklRoaCg+/PBDXLp0SepSWj2eEnKlnOq7c0YNAEISpK2FiIgkd+edd+LOO++Uuow2gS0srqTLFKcdYiUtg4iIqK1hYHGlmsDCjrZEREQuxcDiKmYDkHlAfK6NlrYWIiKiNoZ9WFxBdxlYNuHK1UFRA6Sth4iIqI1hC4srHPjiSli5dTEQPVjScoiIiNoaBhZXKM0Rpze8AAx6QNpaiIiI2iAGFleoKhGn3h0kLYOIiKitYmBxhcoSceoVKGUVREREbRYDiyvYWlgCpayCiIiozWJgcQVbCwtPCRERUfszY8YMTJo0ya2f4fLAEhsbC5lMVucxe/bsetdfvnx5nXW9vLxcXZZ7sYWFiMjjjR49ut47NS9fvhyBgYEtXo+j0tPT7Y6RarUaCQkJeP311+u9Q3Rb5fJxWPbv3w+LxWJ7ffz4cdx44424++67G9wmICAAp0+ftr2WyWSuLst9LCbAWCY+Zx8WIiJqIovFAplMBrm8/raETZs2oWfPnjAYDNixYwceeughREZGYubMmS1cqTRc3sISGhqKiIgI2+PHH39EfHw8Ro0a1eA2MpnMbpvw8HBXl+U+Vborz7200tVBREQuUXN645VXXkFoaCgCAgLw6KOPwmg02tYZPXo05syZgzlz5kCr1SIkJAQLFiywa/EwGAx46qmn0LFjR/j6+mLo0KHYunWrbXlNy87333+PHj16QKPRICMjo8G6goODERERgc6dO2Pq1KkYMWIEDh06ZFtutVrx6quvolOnTtBoNOjXrx82bNhgW75161bIZDKUlJTY5qWmpkImkyE9Pd2upl9//RXdu3eHn58fJkyYgOzsbNs2FosF8+bNQ2BgIIKDg/HMM8+0SEuPW/uwGI1GfPXVV3jwwQev2WpSVlaGzp07Izo6GnfccQdOnDhxzfc1GAzQ6/V2D8nU9F9R+wMKDhxMRO2QIADGcmkebjpQpqSk4I8//sDWrVuxcuVKrF27Fq+88ordOl9++SWUSiX27duHDz74AO+99x4+//xz2/I5c+Zg9+7dWLVqFY4ePYq7774bEyZMwNmzZ23rVFRU4K233sLnn3+OEydOICwszKH6Dhw4gIMHD2Lo0KG2eR988AHeffddvPPOOzh69CjGjx+P22+/3e7zHFFRUYF33nkH//3vf7F9+3ZkZGTgqaeesi1/9913sXz5cnzxxRfYsWMHioqK8N133zn1GU3h1iPsunXrUFJSghkzZjS4Trdu3fDFF1+gT58+0Ol0eOeddzB8+HCcOHECnTp1qnebhQsX1vmLIxn2XyGi9s5UAbwZJc1n/yMLUPu6/G3VajW++OIL+Pj4oGfPnnj11Vfx9NNP47XXXrOdsomOjsb7778PmUyGbt264dixY3j//ffx8MMPIyMjA8uWLUNGRgaiosR989RTT2HDhg1YtmwZ3nzzTQCAyWTCxx9/jL59+zZa0/DhwyGXy2E0GmEymfDII49g2rRptuXvvPMOnn32Wdx7770AgLfeegtbtmzB4sWLsWTJEoe/u8lkwtKlSxEfHw9ADF6vvvqqbfnixYsxf/583HnnnQCApUuX4tdff3X4/ZvKrS0s//73vzFx4kTbH1Z9kpOTMW3aNPTr1w+jRo3C2rVrERoaiv/7v/9rcJv58+dDp9PZHpcuXXJH+Y7hGCxERG1O37594ePjY3udnJyMsrIyu+PNsGHD7M4eJCcn4+zZs7BYLDh27BgsFgu6du0KPz8/22Pbtm04f/68bRu1Wo0+ffo4VNPq1auRmpqKI0eO4JtvvsH69evx3HPPAQD0ej2ysrIwYsQIu21GjBiBP/74w6nv7uPjYwsrABAZGYm8vDwAgE6nQ3Z2tl3LjlKpxKBBg5z6jKZwWwvLxYsXsWnTJqxdu9ap7VQqFfr3749z5841uI5Go4FGo2luia7BFhYiau9UPmJLh1Sf7aCAgADodLo680tKSqDVurYPYllZGRQKBQ4ePAiFQmG3zM/Pz/bc29vb4QtNoqOjkZCQAADo3r07zp8/jwULFuDll192aPualqHa/U1MJlOd9VQqld1rmUzmEVcjua2FZdmyZQgLC8Mtt9zi1HY1yTQyMtJNlblYZbE4ZYdbImqvZDLxtIwUDyeuKu3WrZtdJ9Uahw4dQteuXe3mHTlyBJWVlbbXe/bsgZ+fH6Kjo23z9u7da7fNnj17kJiYCIVCgf79+8NisSAvLw8JCQl2j4iICIdrvhaFQgGz2Qyj0YiAgABERUVh586dduvs3LkTPXr0ACBeFAPArgNtamqqU5+p1WoRGRlp993NZjMOHjzYxG/hOLe0sFitVixbtgzTp0+HUmn/EdOmTUPHjh2xcOFCAMCrr76KYcOGISEhASUlJXj77bdx8eJFPPTQQ+4ozfUKzojTgI7S1kFERNc0a9YsfPTRR3jsscfw0EMPQaPR4KeffsLKlSvxww8/2K1rNBoxc+ZMvPDCC0hPT8dLL72EOXPm2F1ynJGRgXnz5uGvf/0rDh06hA8//BDvvvsuAKBr166YOnUqpk2bhnfffRf9+/dHfn4+UlJS0KdPH6d/zANAYWEhcnJyYDabcezYMXzwwQcYM2YMAgICAABPP/00XnrpJcTHx6Nfv35YtmwZUlNT8fXXXwMAEhISEB0djZdffhlvvPEGzpw5Y6vXGY8//jgWLVqExMREJCUl4b333rO78shd3BJYNm3ahIyMDDz44IN1lmVkZNj9gRcXF+Phhx9GTk4OOnTogIEDB2LXrl22ROjRrFbg/Gbxeefh0tZCRETX1KVLF2zfvh3PP/88xo0bB6PRiKSkJKxZswYTJkywW3fs2LFITEzE9ddfD4PBgClTptQ59TJt2jRUVlZiyJAhUCgUePzxx/HII4/Yli9btgyvv/46nnzySVy+fBkhISEYNmwYbr311ibVP27cOABiy0pkZCRuvvlmvPHGG7bljz32GHQ6HZ588knk5eWhR48e+P7775GYmAhAPNWzcuVKzJo1C3369MHgwYPx+uuvX3OctPo8+eSTyM7OxvTp0yGXy/Hggw9i8uTJ9Z5ucyWZ4AknpppJr9dDq9VCp9PZkmaLOLsR+PpP4iXNTxxnPxYiaheqqqqQlpaGuLi41jcyuQNmzJiBkpISrFu3rsF1Ro8ejX79+mHx4sUtVldr1tDfGWeO37yXUHNkp4rTpJsZVoiIiNyIgaU5Tv0kTkO6Xns9IiIiahYOzdpU6TuArMPicwYWIqI2Y/ny5Y2uU3uIfWoZbGFpqt3VowYGxQMJ46SthYiIqI1jYGkKqxW4WH2t+12fA2rHBy4iIiIi5zGwNEXRBfEuzUpvIKK31NUQEUmiDVxkSi3EFX9XGFiaouSiOA2KAxSqa69LRNTG1Aw1bzQaJa6EWouKigoAdYf9dwY73TaFLlOccnRbImqHlEolfHx8kJ+fD5VKZTcYKFFtgiCgoqICeXl5CAwMrHNfJWcwsDTF/s/FqbaTtHUQEUlAJpMhMjISaWlpuHjxotTlUCsQGBjY7HsoMbA4q7IEyDkqPg+Ov+aqRERtlVqtRmJiIk8LUaNUKlWzWlZqMLA4qzjtyvMB06Wrg4hIYnK5vE0OzU+eiScenVVc3fzZaTDg1YL3LSIiImrHGFicVXRenAZ2lrYOIiKidoSBxVlHvxGnHQdKWwcREVE7wsDiDEEA8k+Lz7vfJm0tRERE7QgDizNMlQCqR+vz7iBpKURERO0JA4szjGVXnqt4/yAiIqKWwsDijJrAovYDOLIjERFRi+FR1xnGcnGq9pW2DiIionaGgcUZhlotLERERNRiGFicwRYWIiIiSTCwOMNYKk41/tLWQURE1M4wsDjKUAqsmSE+ZwsLERFRi2JgcdT5LVeeB3WRrg4iIqJ2iIHFUVUl4lSuAsa9ImkpRERE7Q0Di6Oq9OK052RAxdupExERtSQGFkcZ2OGWiIhIKgwsjjJUt7B4BUhbBxERUTvEwOKomlNCbGEhIiJqcQwsjqppYdGwhYWIiKilMbA4ioGFiIhIMgwsjjJWiFO1j7R1EBERtUMMLI4SLOJUrpS2DiIionaIgcVR1urAIlNIWwcREVE7xMDiqJoWFhl3GRERUUvj0ddRgiBO5dxlRERELY1HX0fxlBAREZFkGFgcZet0y8BCRETU0hhYHMUWFiIiIskwsDiKLSxERESScXlgefnllyGTyeweSUlJ19xmzZo1SEpKgpeXF3r37o2ff/7Z1WU1n9UqTtnCQkRE1OLc0sLSs2dPZGdn2x47duxocN1du3ZhypQpmDlzJg4fPoxJkyZh0qRJOH78uDtKazpbCwsbpYiIiFqaW46+SqUSERERtkdISEiD637wwQeYMGECnn76aXTv3h2vvfYaBgwYgI8++sgdpTUd+7AQERFJxi2B5ezZs4iKikKXLl0wdepUZGRkNLju7t27MW7cOLt548ePx+7duxvcxmAwQK/X2z3cjn1YiIiIJOPywDJ06FAsX74cGzZswCeffIK0tDRcd911KC0trXf9nJwchIeH280LDw9HTk5Og5+xcOFCaLVa2yM6Otql36FebGEhIiKSjMsDy8SJE3H33XejT58+GD9+PH7++WeUlJTgm2++cdlnzJ8/Hzqdzva4dOmSy967QWxhISIikozbbz0cGBiIrl274ty5c/Uuj4iIQG5urt283NxcRERENPieGo0GGo3GpXU2qmZoft5LiIiIqMW5/ehbVlaG8+fPIzIyst7lycnJSElJsZu3ceNGJCcnu7s051h580MiIiKpuPzo+9RTT2Hbtm1IT0/Hrl27MHnyZCgUCkyZMgUAMG3aNMyfP9+2/uOPP44NGzbg3XffxalTp/Dyyy/jwIEDmDNnjqtLax6eEiIiIpKMy08JZWZmYsqUKSgsLERoaChGjhyJPXv2IDQ0FACQkZEBea2xTIYPH44VK1bghRdewD/+8Q8kJiZi3bp16NWrl6tLax52uiUiIpKMTBBqOme0Xnq9HlqtFjqdDgEBAe75kFc6AIIVePI04N9w/xoiIiJyjDPHb3bIcIQgiGEFYAsLERGRBBhYHFETVgD2YSEiIpIAA4sjavqvALxKiIiISAI8+jpCqBVY2MJCRETU4hhYHGHXwsLAQkRE1NIYWBzBFhYiIiJJMbA4gn1YiIiIJMWjryNqD1XDU0JEREQtjoHFEXanhLjLiIiIWhqPvo7gsPxERESSYmBxBG98SEREJCkGFkewhYWIiEhSDCyOYAsLERGRpBhYHGHljQ+JiIikxMDiCFsLC3cXERGRFHgEdgT7sBAREUmKgcUR7MNCREQkKQYWR7CFhYiISFIMLI4QajrdcncRERFJgUdgR1zcKU7Z6ZaIiEgSPAI7IvOAOFX5SFsHERFRO8XA4girWZwOelDaOoiIiNopBhZHCII4VWqkrYOIiKidYmBxRE2nW8gkLYOIiKi9YmBxSHULC68SIiIikgSPwI7gZc1ERESS4hHYEQwsREREkuIR2BEMLERERJLiEdgRDCxERESS4hHYEbbAwquEiIiIpMDA4oiacVgYWIiIiCTBwOIIgZc1ExERSYlHYEewDwsREZGkeAR2BAMLERGRpHgEdgQDCxERkaR4BHYEAwsREZGkeAR2BAMLERGRpHgEdgTv1kxERCQpBhZHsIWFiIhIUjwCO4QDxxEREUnJ5YFl4cKFGDx4MPz9/REWFoZJkybh9OnT19xm+fLlkMlkdg8vLy9Xl9Z0HDiOiIhIUi4/Am/btg2zZ8/Gnj17sHHjRphMJtx0000oLy+/5nYBAQHIzs62PS5evOjq0pqOp4SIiIgkpXT1G27YsMHu9fLlyxEWFoaDBw/i+uuvb3A7mUyGiIgIV5fjGgwsREREknL7EVin0wEAgoKCrrleWVkZOnfujOjoaNxxxx04ceJEg+saDAbo9Xq7h1sxsBAREUnKrUdgq9WKuXPnYsSIEejVq1eD63Xr1g1ffPEF1q9fj6+++gpWqxXDhw9HZmZmvesvXLgQWq3W9oiOjnbXVxAxsBAREUlKJgg1PUpdb9asWfjll1+wY8cOdOrUyeHtTCYTunfvjilTpuC1116rs9xgMMBgMNhe6/V6REdHQ6fTISAgwCW12/lXf6DoAvDgr0DMMNe/PxERUTuk1+uh1WodOn67vA9LjTlz5uDHH3/E9u3bnQorAKBSqdC/f3+cO3eu3uUajQYajcYVZTqGVwkRERFJyuVHYEEQMGfOHHz33XfYvHkz4uLinH4Pi8WCY8eOITIy0tXlNc01TgntvVCIlfsyWrggIiKi9sXlLSyzZ8/GihUrsH79evj7+yMnJwcAoNVq4e3tDQCYNm0aOnbsiIULFwIAXn31VQwbNgwJCQkoKSnB22+/jYsXL+Khhx5ydXlNIzQ8cNw9n+4BAHQN98PAztfuWExERERN4/IWlk8++QQ6nQ6jR49GZGSk7bF69WrbOhkZGcjOzra9Li4uxsMPP4zu3bvj5ptvhl6vx65du9CjRw9Xl9c0DnS6zSyutHttMFtwNLMEbuwiRERE1G64vIXFkQP01q1b7V6///77eP/9911dius0EFis1ivfVaWwX7Z8ZzoW/nIKM0fGYcGtHhK8iIiIWin2InVEA4GlymyxPVfI7U8XfbRF7DD87x1p7q2NiIioHWBgcURNYIF9KCk31Aos1f1bBEHAznMF6Bzsc2VznhYiIiJqFrdd1tymNNDCUmE0256bLOI6vxzPwd++PmS3ntFihUapcG+NREREbRhbWBxS/zgstVtYjNWBZX3q5TpbV9Raj4iIiJzHFhZHONDCsudCITb9kYf96cV1Ni8zmNHBV+3WEomIiNoyBhZH1BNYrFYBb204ZXu9ct+lBjevMLKFhYiIqDl4SsgR9Qwc9/Px7HpbU+pTXqslhoiIiJzHwOKIelpYskuqHN6cfViIiIiah4HFEbbAcqWFpbjC6PDmbGEhIiJqHgYWR9TTwnK5pLKBleuqYGAhIiJqFgYWRwh1L2suKhdbWEL9NY1uXs5TQkRERM3CwOKIelpYyg1iq0kHH1Wjm9cMKkdERERNw8DiiHoCS82lyoHejY+vwsBCRETUPAwsjqgnsFSaxMCiraeFxV+jxIu17tBssvBeQkRERM3BwNIYQUB9Q/PXtLAE1zOC7awx8XhwZBzuGxYDADCa2cJCRETUHBzptjG177Rcu4WlOrAE1Qosz0zohlA/De4a0AkAoFKI6/OUEBERUfMwsDRGqBs2BEGwXaoc7HflKqEuIX6Y0CvC9lrNwEJEROQSPCXUqLotLAazFdbq2bVPCWmU9rvzSgsL+7AQERE1BwNLY2q3sFQHli2n8myzap8SUjcQWIxsYSEiImoWBpbGXBVYBEHArK8P2Wb5e105q1anhUUpDuVvYqdbIiKiZmFgacxVgaW4wmS32Ed9JbBc3cLCPixERESuwcDSmKsCS1G5wfZy9SPD4K1S2F43dEqIfViIiIiah4GlMXUCi9jCEhfii6FdguGtvhJYagLK1a/Zh4WIiKh5GFgaYxdYZLabHtbcQ6h2YBGuakhRKar7sDCwEBERNQsDS2OuGjiuuEIMLDVXB9U+JeSluqoPi5J9WNzl+GUdjlwqkboMIiJqIRw4rjFXBZYrLSxiYFHIZXjtjp7QV5nRqYOP3aa2Pixm9mFxJbPFils/3AEAOLzgRnSo5/YIRETUtjCwNOaqU0LF5fYtLABwf3JsvZuyD4t71L5S62S2HiMSQiSshoiIWgJPCTXmqjs1F1WfEnLkVz37sLhHTSsXAPyRrZewEiIiaikMLI25KrCUVP+6D/JpPLBwHBb3KKx1afmFgnIJKyEiopbCwNKYq1tYqn/dB1ZfJXQtqupOt0aOdOtSxeVXTgmt2JuBc3llElZDREQtgX1YGmPrwyKe3rn6KqFrqbmCqNJkcUtp7VXtwfsAYNx725C+6BaJqiEi8gxVJgsOXiyGn0YJqyAgLsQXgVedDdBVmuCvUcJsFXA6pxT/2nwW/l5KhPprcOKyHnEhvogO8oYMMsSG+GJkQgi8VHLIZDKJvtUVDCyNEarDhlwMH7arhBwJLNVjtFQYGVhcqaDMWGeeyWKtM3AfEVFrUGWyYNMfuTh2WYfTOaWoNFpQabIgyFeNED8NAODGHuEY1TUUGqUYHorLjTiXX4b8UgP+u/sicvVVKDWYkV9q/4NudLdQTBkSg2FdgvFhyln8e2caVAo5zBYrrPVcwLrjXIHda1+1Al4qBe7o1xHPTUyqM6J7S2JgaYyxuo+E2hcmixWlVWYAjvVh8akOLJUMLC5TabTgg5SzdeZnlVSic7CvBBURETmnoMwArbcKMgAr91/C2xtOQV99bGnItwczbc8Vchks9aWNemw9nY+tp/Pt5l3dTSHASwl9lRkRAV4oKjfaXdlabrSg3GjBFzvTkFtahY+m9JestYWBpTGGUnGq9rOdDpLLgADvxvuw+KjE3Wu2CjCarZIm07bi3zsu1Dv/YmEFAwtRO1ZaZcK61CyczyvDdYkh6Bruj04dvFv04HqpqAL/t/089qUV4cERcbh3SEyd5R9vPYdvDmTWCRxKuQzJ8cEYENMB8WF+8FLKUVJhQn6ZAekF5dhwPAelBjHUXL3tX4bGIMxfgxA/DW7tEwmttwpZuiq8veEULpdUYn96sW3dp8d3g7+XEpeLKzE4NgjXdw2FWilHaZUJPmoljGYrCssN8FErYTBb8H1qFrafzcfOc4XYdDIXFwrKER/q56Y9eG0MLI0xVHfo1PjbOnsG+qihkDf+j6D2sP2VRgsDSzMJgoB3fjtT77JsXSXO5pbi52M5mJbcmYPJUatWZjBjf3oR8ksNKCwzorTKBH2VCRVGC5Ii/PGngdFQKWRYuS8DZQYLZo2Kt/v/pi0wmC0oqzIjuPqUCADklVZhyeZz6B4ZgJGJIfDXqKD1UWF96mW8/P0J2xhNy3elAwDGdAtFv+gO8NUoML5nBKKDfFBYZsArP5zE90eyMGt0PIbGBWFEQki9p5SrTBb871AmxveMsJ2aEQTBLgT9fjYfM5cfqDPe1nNrj8FkFXD/sM4AgA3Hc/DoVwfr/a639InE8zd3R1Sgd4P74807e2PVvgysPXwZJ7L06BEZgEV39UZUoDcCvOr+gO4Y6I3F9/YHAPx2Igdv/3oat/eNwqxR8ZDXc/zy97pyu5lO6iuDoP51VDz+Oioeaw5cQp9OgZKFFQCQCcLVd8BpffR6PbRaLXQ6HQICAlz75ie+A9bMADqPwO7r/4spn+1BfKgvUp4c7dDmic//DJNFwO75NyBS2/BfRmrc90ey8NjKw7bXE3pGYMOJHADA30bH4+Ot5wGIvyBmj0mQpEaipjJbrFi1/xKW70p3+sq3KK0X5t7YFX8eFO2m6lqG2WLFP389jc9+v2AbZPy2vlGYd2NXmCxW3PnxLpQZ7E+dBPmq7e7xNqprKL4/klWnf4ZcJv7YrD2OU20f3NsPd/TraHtttQr461cHsfFkLgCgbyctqkxi68Obk3ujtMqMP7L1+M+ei3anWJIi/BHkq8au84UNfs9x3cOhUcqRX2bAsxOSMLBzB4f3ESCGOY2ybQRUZ47fbGFpTPUpIYPCB1M+2wMAqDI5fpmyt0oBk8WMckP9/ViOX9bhmwOX8NjYRFuCp/qdv+o/8SVTB+CDTWfwr83n7M7vnsjStXRp5OEMZgt0FSYE+2lQWmWCVQCUClm9v0ydUVRuxLKdadifXoRygwWju4WipMKEkkoTpiV3xuDYILv1a36d55caMPPL/SgsM2J0t1B4qRRYtjPN7iAb7KtGbIgvOgf5QOujgr+XChqlHOtTL+NM7pV/Cwq5DFm6Kjzz7VH8ciwbn08f7FALcG1mixXKq1oYzuaW4lJxBVbvv4RKkxWzRsUjwFsJhVyGF9efQKi/BqO6hiIxzA8ymQy9O2qd/tzaLuSXYdZXh3A6t9Ru/g9HsvDDkawGt6sJIEkR/vj3jMHoGOiNZycm4dPtF7BsZ7ptPauABsMKADy95ihignywZMt5eKnk6NTBxxZWAOBI5pX/Vx75b92WErkMePfPfTG5fycIgoC/fLYXuy/Yh5bkLsH4+9gEDIsLrreVw1FtJaw4i4GlMdWnhPIMV04xdI90vBXHRy12ZqrvH4rZYsW0L/ahqNyI/FIDPrlvYPPrbSNS/shFkK8a/WOu/PKo6UMEAD/MGQmFXGZrQs2r1TP+52M5OHix2OlfLeTZisqNKKkwIkdfBX2lGQazBWaLgIIyA8qNFkRpvTAoNggmixWVJgv2XijC2bxSnMzS40xuaZ1f3Aq5DMG+avSMCsCzE5OQFGH/71oQBJzPL0N8qB/ySw0I8lVDqZCjuNyIF9YfR36pAfvTi+xuN3bs8pWD2rUOsrV9vTfD7vVfr++C2/tFISHMr94D08yRcVi67TzSC8oxe0wCgnzVeOS/B3HwYjG2nM7HlM/24MsHhjh0ikhXacIDy/bhUEYJhsQFwV+jRMqpPIT4qetcjbf9TH6d7X86mm17PqZbKD65byC8VNf+3Av5Zfjb14cQFuCFT6YOwP70IpzLK8O7v52xGwLivmEx6N1Ri+fWHrPtY4Vchq1PjUZplRnfHc6ERqlAhdGCxHA/3DMo2hYCIrXeeOm2nrh3cAzUSjlignzw1Z6L+OFIForKjfjqoaEID/DC/w5lorjciJ+OZeNopg6TP95Vp967BnSCAAGVRgsGxQbhhyNZSK1149X5E5Nwz+Bou8uHZTIZlkwdgBV7L+LYZR2iO/hgWJdgjOoWyqsZm4GnhBqz7Z/AljdwquNdmHD+LgDA93NGoE+nQIc2j33uJwDiQHOpL95kt+xElg63/GuH7XXqizfWuWa+NcvWVSIiwMvpTm9nc0tx4/vbAYj/CQ7tEoxHR8Vj6ud7sPNcIaYMicHCO3sDEM8f3//vfXXe47a+UfhwSv/mfwlyiwPVBymzVUB6QTlO55Yiv9QAo8UKi1WAxSqgtMoMQRAQ4qeB0WJFZnGl2+pRyGV47899bacEjmXqsGxnGtYevtzotsG+agyLD0ZplRkhvmqEBmiw50JRo3cTjw32wYCYDrbPuCEpDO/c3dehMZ6uJggCvtyVjrd/PY1yowV9OmkR3cEHz05IQkywT73b/HQ0G8+tPWq78vFa5DLUewlsfeJDfaH1VqFvdCAW3NLDFiI2n8pFucGCpdvO40RW/bfU6BzsgwdHxGH68FjbvNIqExasO44fjmbj46kDML5nhGOFOCGvtApD3kixm6dWyHHXwE54fVKvOi1H5/JKUWWyoldHrctraW94SshVzEZgyxsAAL3VCwAwY3isw2GltpIKE7JKKu06VZ3Ktm/6XPTLKSy6q0/T63UjXaUJu88XINTfCz5qBSICvOw6tpYbzPjmwCUs25mOjKIK23xvlQJ+Xko8eWNX3DM42qHwcrLW/YG2nM7HltP5GB4fjNSMEgDAtOTOtuVdw/1tz2UywE+tRKnBjB+OZGHKkGh8mHIOgT4qvHJ7T4QFeDVnF7R6giAgS1eFM7mlgADEBPsgNtgXCrkMVSYL8ksN6BjoDXn16x1nCyCTiWMOXcgvh5dKDoPJirxSA6zVv3OySiqRWVyJAG8VvFVydAn1w8XCcgT6qPHIdV3QwVeNwxnFyNZVQeutwrP/O+p08Kh9uWfN370gXzU0KjlkkCHYTw1vlQJ/ZOvtmu0BYMqQaPTqqMWAmA7oGu6PzOIKFFeYEB/qixxdFXaeK8BX1aMlP74qFR9tPoezDvYf8VUrcPegaMwaHY/wq/5uWawCvt57Eav3X8KJLD2u7xqKA+lFtjGZZo+Jx9xxXaFSyHHP4GgE+2kQH+rb5CtaZDIZZoyIQ5dQP0z7Yh+OZupwNFOHSpMF/54+CDKZDOUGMx5fdRgdfNTI1lXZjbfRNdwP3SMDsPdCEXL0VQCAkQkhmDkyDkqFDINjg/Dh5rMoN1jw8PVdEOSjxqXiCvhplNBXmfDZ9jT875B4WvZ8vjgUxKGMEizbmQ61Qu7QDWA7Bnrj17nX12mh8fdS4b0/98Mrd/SC1oGrM5sizN8Lb0zuhQXrjmP2mAQ8Ma4rrIJQ5zRZjYQw/3rnk3u5rYVlyZIlePvtt5GTk4O+ffviww8/xJAhQxpcf82aNViwYAHS09ORmJiIt956CzfffLNDn+W2FhazEXg9FAAw1/g3rLOOxJM3dsXfxyY6/Bav/3gSn+9IAwA8OyEJf72+i+0Xx8vfn8DyXem2a+C13iocXnCj3bnNgjIDDGYrOl6j97g7ncsrw392p+OrPRftfmFplHLc3jcKZQYzCsoMdpfNNeSxGxKgVMgxrns4lu1Mw6DYDrhncEyd9d7beAb/umqslbsHdsKag5nw91Ii9cWbbL94BEFA3PyfAYhhct5NXdHn5d/qvGdShD82zL3ema/eahSVG3E2txRD4oLsDnjpBeVILyxHucGCH45kYdMfuTDX8zPZX6NEudEMqyCGPkAc19nRX9RNNTw+GFZBQHiAF4bHByNS6w2NUg6lQgaFXA4ftQJymThQoK7ShP7RgdcMnYIg4GxeGWKCfFBptDh8pZjZYsXT3x7Fd/W0piR3CYaXSo4+nQJxIksPfZVYx33DOiM6qP6WC6mt2JuB1386aQtHY7qFIkLrhZX7LtVZV6OU46fHrkNCmP2VH1ar4HQfi0tFFfjT0l3I1RuuuV4HHxVWPDwM7/52GkXlRnw8dSB+P5uPvtGBdj9ApFBaZYKfRukRo7q2F84cv90SWFavXo1p06Zh6dKlGDp0KBYvXow1a9bg9OnTCAsLq7P+rl27cP3112PhwoW49dZbsWLFCrz11ls4dOgQevXq1ejnuS2wWK1AzhGM+egg0qyRAIA3J/fGX4bWPchey32f77X9mrlvWAwm9opEcpdgJC9KQa7egE+mDsC8b46g0mTBb09cb/tHa7UKuO6fW5Crr8KyBwbjusRQp79CucGM07mlGBDjXH+OI5dKsPZQJr7cfdHpz5wxPBZ+GiUKygyQy2VYcdU5+tr+O3NIne9158c7cai6NeVqo7qG4ssH7YPvsp1p2HOhEO/c3Rf+Xirc/++9+P1sQZ1t9z0/FmH+Vw54Z3NLEeSrtrtssqlMFisOZ5QgvaAcYQEaRGq90bGDN/w07mnEPJ1Tit/P5uNIps7WV6JnVABGJIQg2Ff89fv13gzU96872FeNIF81Mosrr3nbiDB/DSpNFpQbzEgM84dCLoNMBoQHeEFfaYJSIUOvKC06dfBGaZUZp3JLseNsAQbEBOJiUQUu5F+5MWWInxqF5UYIAnBn/46YMjSmTodUKZktVqw9fBn/O5iJ8/nluLl3BGaPSajTctKafL33Ip7/7niDy+8bFoPX7ujl0oOzyWJFlckCP40SX+5Kx2e/p8FPo0RxhRFD4oJw18BO6Bbuf83Ld6l9kTywDB06FIMHD8ZHH30EALBarYiOjsbf//53PPfcc3XWv+eee1BeXo4ff/zRNm/YsGHo168fli5d2ujnubMPi9UqoMs/fra9/s+DQ3B9V+eCw3P/O4pV++1/3UwdGoOv92bA30uJAy+MwwPL9tsug/thzkj0jApAamYJ7qzVCezVO3rilt6RkMtk6OCrhiAIMJit0FWaIIM4mF3t5tQ9Fwrx8JcHUGow4//uH4gRCSGoMJhRWG5ElNYbOfoqlFaZMCg2CJdLKvHub6cxplsYMosr8daGU3W+x6t39MRfhsTAKgC/nsjByWw95DJg86l8ZOsq8cx4sfPZ1ed7rVYBD365v85oi4DYFP2/WcPh76WCvsqEpVvP2y5P7hLqi5ggH7vtHLlkuaTCiF+O5+DjredwqejK6Qe1Qo7fnrgesSG+2HwqFw8uPwCttwrrZo9AXEjTBp0zmq347PcL+M/u9Dq/LEP8NHj+liRM6texyQeFmsGc5DLxNENhuRFv/XIK36VerjeM1CfAS4k/D4rGxN6R6BEZYOuMabJYcSJLjwqDGfFhfujgo0auvgqa6lM/HQO9YRUEmK1Cox0pr1ZlsuCz7RdQajDjtj5R6N1Ji0qjBRqlvFlXR5BzFm86g8WbzkKlkMFkEZAY5ocP/9IfscG+Tv+ZErmDpIHFaDTCx8cH3377LSZNmmSbP336dJSUlGD9+vV1tomJicG8efMwd+5c27yXXnoJ69atw5EjR+qsbzAYYDBcOTjo9XpER0e7PLCYLVZMX7YPO8+JQeK1O3pi6tDOTv+H+/7GM/UOJw8AfxrYCe/c3bfe0yDXkhThj1M59n1g1Eo5xnQLha9aic2n81BSYWpga3vjuofhaKbO7kobAIgO8kaInwaf3j8IGpW8WZeAVhot+OV4NsL8vbDpj1x0DffHqz+esF0ifmOPcGw+lWcbwfHO/h3x3j39YLZYMeKtzcjVGxAeoMFvc0dB68CdsmvbdDIXD/3ngO11/5hAnM8rs/WNiA32wYiEEFgFoLD6FJxVEBDgpUJJpREapQLBvmqolXJbKAnxU6PSZMHRTB3SCq60JCSE+SGzuMLu0vfYYB9EB/kgSusNpUKGSqMFZQYz/DRKeKkVyNFViadG/L1wPEuHSqMF/aIDYTBb8euJHMhlsnr7AKgUMkQH+WD+xO4YEBOIH45kYV1qFgRBQFSgN3p11OL+5M5QyeVtblAxctzZ3FJ0DvaF2WqFQi5rt5fEkmeStNNtQUEBLBYLwsPD7eaHh4fj1Km6v9oBICcnp971c3Jy6l1/4cKFeOWVV1xT8DVYBdjCyj2DonF/cmyT3mdIXMNN3zWnl27vG+VUYLk6rADir/1fT+TazQvx06Cg7EoQkckAP43S7sqATX/k1Xmvacmd8eodjZ+Oc5S3WoE7B3QCAIxMDAEARAV6Yc6KwygzmO3GO+jTSWvrfKxUyLFx3iiYLQICvVVN+nU+rkc4Xp/UCy+sE5vHD1efbvJSySEIQHphBdILGz5t1Rh/jRKPjo7HjOGx8NWIQ1sbLVZ8uSsdH24+W/3+FY2/US0XaoUgwP43RXSQNxZO7oMRCcEQBNj2yYwRcZgxIq7J34PapsTqU8xq8HJaat1a5VVC8+fPx7x582yva1pYXE0uE3vze6sUuH9YbJPfZ0RCCJY/MBglFSZofVTQequw/Uw+BnbuYOtbkhDmh03zRmHnuQL4apQ4n18Gq1XAvUNiIAOQVlCOkkpxXIQKozj+RK6+CnEhvkgM94fWW4WMogoculgsXgrqr0F0Bx8MTwjG+bxypF4qwXWJIQgP8IJaKccf2XpolHIUlhuRmlECuVyGyf07wmIVx54Y1AJjmIzuFobNT43CngviJa6xwT6ICPBCv5hAu9sYNHdwLwC4b1hnjEwIwZHMEuToqhAT5IOx3cNRZjBj3eHLV25q6SeeagOA0iozFHIZfNQKlBnMMJis8NUo4K1WoqjMCB+1AqH+GozuFmp3ObpaKYdaKcfsMQm4d3A0TmbrkV1ShcsllZDJAJVCjgBvFSoMZpQbLQjz10AhlyFPb0CEVrwfyIks8XRbj6gAaL3VKCwzINRfA5VCjp5RAbZTTOwbSETthcsDS0hICBQKBXJz7X/p5+bmIiKi/uvnIyIinFpfo9FAo3H/qLBKhRxPj09yyXuN7mbf2bi+TrAJYX51euvXiHWgj0VciC9G1dO/pkdUAHpE2Te11Qx+1yUUdTo/hvq33Ii7Yf5euL1vVIt8VmyIb539GKRU48GR7muVCPbTNKmz9Nju4Y2vRETUjri8jVCtVmPgwIFISbkyCI/VakVKSgqSk5Pr3SY5OdlufQDYuHFjg+sTERFR++KWU0Lz5s3D9OnTMWjQIAwZMgSLFy9GeXk5HnjgAQDAtGnT0LFjRyxcuBAA8Pjjj2PUqFF49913ccstt2DVqlU4cOAAPv30U3eUR0RERK2MWwLLPffcg/z8fLz44ovIyclBv379sGHDBlvH2oyMDMjlVxp3hg8fjhUrVuCFF17AP/7xDyQmJmLdunUOjcFCREREbR/vJURERESScOb4zevciIiIyOMxsBAREZHHY2AhIiIij8fAQkRERB6PgYWIiIg8HgMLEREReTwGFiIiIvJ4DCxERETk8RhYiIiIyOO5ZWj+llYzWK9er5e4EiIiInJUzXHbkUH320RgKS0tBQBER0dLXAkRERE5q7S0FFqt9prrtIl7CVmtVmRlZcHf3x8ymcyl763X6xEdHY1Lly7xPkVuxP3ccrivWwb3c8vgfm4Z7trPgiCgtLQUUVFRdjdFrk+baGGRy+Xo1KmTWz8jICCA/xhaAPdzy+G+bhnczy2D+7lluGM/N9ayUoOdbomIiMjjMbAQERGRx2NgaYRGo8FLL70EjUYjdSltGvdzy+G+bhnczy2D+7lleMJ+bhOdbomIiKhtYwsLEREReTwGFiIiIvJ4DCxERETk8RhYiIiIyOMxsDRiyZIliI2NhZeXF4YOHYp9+/ZJXVKrsXDhQgwePBj+/v4ICwvDpEmTcPr0abt1qqqqMHv2bAQHB8PPzw933XUXcnNz7dbJyMjALbfcAh8fH4SFheHpp5+G2Wxuya/SqixatAgymQxz5861zeN+dp3Lly/jvvvuQ3BwMLy9vdG7d28cOHDAtlwQBLz44ouIjIyEt7c3xo0bh7Nnz9q9R1FREaZOnYqAgAAEBgZi5syZKCsra+mv4rEsFgsWLFiAuLg4eHt7Iz4+Hq+99prd/Wa4n523fft23HbbbYiKioJMJsO6devslrtqnx49ehTXXXcdvLy8EB0djX/+85+u+QICNWjVqlWCWq0WvvjiC+HEiRPCww8/LAQGBgq5ublSl9YqjB8/Xli2bJlw/PhxITU1Vbj55puFmJgYoayszLbOo48+KkRHRwspKSnCgQMHhGHDhgnDhw+3LTebzUKvXr2EcePGCYcPHxZ+/vlnISQkRJg/f74UX8nj7du3T4iNjRX69OkjPP7447b53M+uUVRUJHTu3FmYMWOGsHfvXuHChQvCr7/+Kpw7d862zqJFiwStViusW7dOOHLkiHD77bcLcXFxQmVlpW2dCRMmCH379hX27Nkj/P7770JCQoIwZcoUKb6SR3rjjTeE4OBg4ccffxTS0tKENWvWCH5+fsIHH3xgW4f72Xk///yz8Pzzzwtr164VAAjfffed3XJX7FOdTieEh4cLU6dOFY4fPy6sXLlS8Pb2Fv7v//6v2fUzsFzDkCFDhNmzZ9teWywWISoqSli4cKGEVbVeeXl5AgBh27ZtgiAIQklJiaBSqYQ1a9bY1vnjjz8EAMLu3bsFQRD/gcnlciEnJ8e2zieffCIEBAQIBoOhZb+AhystLRUSExOFjRs3CqNGjbIFFu5n13n22WeFkSNHNrjcarUKERERwttvv22bV1JSImg0GmHlypWCIAjCyZMnBQDC/v37bev88ssvgkwmEy5fvuy+4luRW265RXjwwQft5t15553C1KlTBUHgfnaFqwOLq/bpxx9/LHTo0MHu/41nn31W6NatW7Nr5imhBhiNRhw8eBDjxo2zzZPL5Rg3bhx2794tYWWtl06nAwAEBQUBAA4ePAiTyWS3j5OSkhATE2Pbx7t370bv3r0RHh5uW2f8+PHQ6/U4ceJEC1bv+WbPno1bbrnFbn8C3M+u9P3332PQoEG4++67ERYWhv79++Ozzz6zLU9LS0NOTo7dvtZqtRg6dKjdvg4MDMSgQYNs64wbNw5yuRx79+5tuS/jwYYPH46UlBScOXMGAHDkyBHs2LEDEydOBMD97A6u2qe7d+/G9ddfD7VabVtn/PjxOH36NIqLi5tVY5u4+aE7FBQUwGKx2P0HDgDh4eE4deqURFW1XlarFXPnzsWIESPQq1cvAEBOTg7UajUCAwPt1g0PD0dOTo5tnfr+DGqWkWjVqlU4dOgQ9u/fX2cZ97PrXLhwAZ988gnmzZuHf/zjH9i/fz8ee+wxqNVqTJ8+3bav6tuXtfd1WFiY3XKlUomgoCDu62rPPfcc9Ho9kpKSoFAoYLFY8MYbb2Dq1KkAwP3sBq7apzk5OYiLi6vzHjXLOnTo0OQaGVioRcyePRvHjx/Hjh07pC6lzbl06RIef/xxbNy4EV5eXlKX06ZZrVYMGjQIb775JgCgf//+OH78OJYuXYrp06dLXF3b8c033+Drr7/GihUr0LNnT6SmpmLu3LmIiorifm7HeEqoASEhIVAoFHWupMjNzUVERIREVbVOc+bMwY8//ogtW7agU6dOtvkREREwGo0oKSmxW7/2Po6IiKj3z6BmGYmnfPLy8jBgwAAolUoolUps27YN//rXv6BUKhEeHs797CKRkZHo0aOH3bzu3bsjIyMDwJV9da3/NyIiIpCXl2e33Gw2o6ioiPu62tNPP43nnnsO9957L3r37o37778fTzzxBBYuXAiA+9kdXLVP3fl/CQNLA9RqNQYOHIiUlBTbPKvVipSUFCQnJ0tYWeshCALmzJmD7777Dps3b67TTDhw4ECoVCq7fXz69GlkZGTY9nFycjKOHTtm949k48aNCAgIqHPgaK/Gjh2LY8eOITU11fYYNGgQpk6danvO/ewaI0aMqHNp/pkzZ9C5c2cAQFxcHCIiIuz2tV6vx969e+32dUlJCQ4ePGhbZ/PmzbBarRg6dGgLfAvPV1FRAbnc/vCkUChgtVoBcD+7g6v2aXJyMrZv3w6TyWRbZ+PGjejWrVuzTgcB4GXN17Jq1SpBo9EIy5cvF06ePCk88sgjQmBgoN2VFNSwWbNmCVqtVti6dauQnZ1te1RUVNjWefTRR4WYmBhh8+bNwoEDB4Tk5GQhOTnZtrzmctubbrpJSE1NFTZs2CCEhobycttG1L5KSBC4n11l3759glKpFN544w3h7Nmzwtdffy34+PgIX331lW2dRYsWCYGBgcL69euFo0ePCnfccUe9l4b2799f2Lt3r7Bjxw4hMTGxXV9ue7Xp06cLHTt2tF3WvHbtWiEkJER45plnbOtwPzuvtLRUOHz4sHD48GEBgPDee+8Jhw8fFi5evCgIgmv2aUlJiRAeHi7cf//9wvHjx4VVq1YJPj4+vKy5JXz44YdCTEyMoFarhSFDhgh79uyRuqRWA0C9j2XLltnWqaysFP72t78JHTp0EHx8fITJkycL2dnZdu+Tnp4uTJw4UfD29hZCQkKEJ598UjCZTC38bVqXqwML97Pr/PDDD0KvXr0EjUYjJCUlCZ9++qndcqvVKixYsEAIDw8XNBqNMHbsWOH06dN26xQWFgpTpkwR/Pz8hICAAOGBBx4QSktLW/JreDS9Xi88/vjjQkxMjODl5SV06dJFeP755+0uleV+dt6WLVvq/T95+vTpgiC4bp8eOXJEGDlypKDRaISOHTsKixYtckn9MkGoNXQgERERkQdiHxYiIiLyeAwsRERE5PEYWIiIiMjjMbAQERGRx2NgISIiIo/HwEJEREQej4GFiIiIPB4DCxEREXk8BhYiIiLyeAwsRERE5PEYWIiIiMjjMbAQERGRx/t/OAvrcNrjyyYAAAAASUVORK5CYII=",
                        "text/plain": [
                            "<Figure size 640x480 with 1 Axes>"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                }
            ],
            "source": [
                "import numpy as np\n",
                "import matplotlib.pyplot as plt\n",
                "\n",
                "d = 10\n",
                "length = 1000\n",
                "R = 1\n",
                "lam = 0.1  # V用のパラメータです\n",
                "delta = 0.9  # 信頼パラメータです\n",
                "\n",
                "# etaはR-sub gaussianにします\n",
                "etas = (np.random.rand(length) - 0.5) * 2 * R\n",
                "# X用に適当な確率過程を作っておきます\n",
                "F = np.zeros((length, d))\n",
                "for i in range(d):\n",
                "    F[:, i] = etas * i\n",
                "\n",
                "# Vは適当な正則行列です。今回はlambda * eyeとします。\n",
                "V = 1 * np.eye(d)\n",
                "S_t_norms = []\n",
                "upper_bound = []\n",
                "Xs = np.zeros((length, d))\n",
                "for t in range(length):\n",
                "    # Xは適当に時刻t-1から予想できるF_{t-1}-measurableな値にします\n",
                "    if t == 0:\n",
                "        pass\n",
                "    else:\n",
                "        Xs[t] = F[:t].sum(axis=0)\n",
                "\n",
                "    V_t = V.copy()\n",
                "    S_t = np.zeros(d)\n",
                "    for i in range(t):\n",
                "        X = Xs[i].reshape(d, 1)\n",
                "        V_t += X @ X.T\n",
                "        S_t += etas[i] * Xs[i]\n",
                "    assert S_t.shape == (d,)\n",
                "\n",
                "    # ノルムを計算します\n",
                "    S_t_norms.append(S_t @ np.linalg.inv(V_t) @ S_t)\n",
                "\n",
                "    # 上界を計算します\n",
                "    y = np.log(np.sqrt(np.linalg.det(V_t)) / np.sqrt(np.linalg.det(V)) / delta)\n",
                "    upper_bound.append(2 * R ** 2 * y)\n",
                "\n",
                "plt.plot(S_t_norms, label=r\"$\\|S_t\\|^2_{\\widebar{V}_t^{-1}}$\")\n",
                "plt.plot(upper_bound, label=r\"Upper Bound\")\n",
                "plt.legend()"
            ]
        },
        {
            "attachments": {},
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "さて、このバウンドを使った線形バンディットについて考えてみましょう。\n",
                "\n",
                "行動の集合として$\\mathcal{D}_t \\subset \\mathbb{R}^d$ を考えて、その中から一つ行動$X_t \\in \\mathcal{A}_t$ を選択することを考えます。そして、報酬関数が次の形をしている状況を考えます:\n",
                "\n",
                "$$\n",
                "Y_t = \\langle X_t, \\theta_*\\rangle + \\eta_t\n",
                "$$\n",
                "\n",
                "ここで$\\eta_t$は$R$-subgaussianです。以下、次の擬リグレットを考えます：\n",
                "\n",
                "$$\n",
                "R_n =\\sum_{t=1}^n \\left\\langle x_t^* - X_t, \\theta_* \\right\\rangle\n",
                "$$\n",
                "\n",
                "ここで$x_t^* = \\mathrm{argmax}_{x_\\in D_t} \\langle x_t, \\theta_* \\rangle$です。\n",
                "リグレットを定義したので、以下では[BANDIT_contextual.ipynb](BANDIT_contextual.ipynb)と同様にしてパラメータの信頼区間を求めて、その一番大きい値をスコアにします。\n",
                "\n",
                "$$\\mathrm{UCB}_t(x)=\\max _{\\theta \\in \\mathcal{C}_t}\\langle x, \\theta\\rangle$$\n",
                "\n",
                "そして$X_t=\\operatorname{argmax}_{x \\in \\mathcal{X}_t} \\operatorname{UCB}_t(x)$でアームを選択します。\n",
                "$C_t$はどうやって作ればよいでしょうか？ここでself-normalized boundが出てくるわけですね。\n",
                "\n",
                "### $\\theta_*$の推定と信頼集合\n",
                "\n",
                "$\\theta_*$の推定値を$\\widehat{\\theta}_t=\\left(\\mathbf{X}_{1: t}^{\\top} \\mathbf{X}_{1: t}+\\lambda I\\right)^{-1} \\mathbf{X}_{1: t}^{\\top} \\mathbf{Y}_{1: t}$とします。ここで、$(\\mathbf{X}_{1: t})$ はその行が $(X_{1}^{\\top}, X_{2}^{\\top}, \\ldots, X_{t}^{\\top})$で定義された行列です。また、$(\\mathbf{Y}_{1: t}=\\left(Y_{1}, \\ldots, Y_{t}\\right)^{\\top})$とします。\n",
                "このような推定値に対して、次の定理が成立します：\n",
                "\n",
                "---\n",
                "**定理：信頼楕円**\n",
                "\n",
                "$V=I\\lambda$, $\\lambda>0$とし、$(Y_{t}=\\left\\langle X_{t}, \\theta_{*}\\right\\rangle+\\eta_{t})$とする。また、$(\\left\\|\\theta_{*}\\right\\|_{2} \\leq S)$とする。このとき、確率$1-\\delta$以上で$\\theta_*$は\n",
                "\n",
                "$$\n",
                "C_{t}=\\left\\{\\theta \\in \\mathbb{R}^{d}:\\left\\|\\widehat{\\theta}_{t}-\\theta\\right\\|_{\\bar{V}_{t}} \\leq R \\sqrt{2 \\log \\left(\\frac{\\operatorname{det}\\left(\\bar{V}_{t}\\right)^{1 / 2} \\operatorname{det}(\\lambda I)^{-1 / 2}}{\\delta}\\right)}+\\lambda^{1 / 2} S\\right\\} .\n",
                "$$\n",
                "\n",
                "また、$\\bar{V}_t$の固有値を$\\alpha_1, \\alpha_2, ..., \\alpha_d$とすると、$\\operatorname{det}\\left(\\bar{V}_t\\right)=\\prod_{s=1}^t \\alpha_s$および$\\operatorname{trace}\\left(\\bar{V}_t\\right)=\\sum_{s=1}^t \\alpha_s$であり、算術平均と幾何平均の関係を使うと、$\\operatorname{det}\\left(\\bar{V}_n\\right) \\leq\\left(\\operatorname{trace}\\left(\\bar{V}_n\\right) / d\\right)^d$です。\n",
                "これを使うと、$\\|X_t\\|_2\\leq L$のとき、以下の信頼区間に変形できます：\n",
                "\n",
                "$$\n",
                "C_{t}^{\\prime}=\\left\\{\\theta \\in \\mathbb{R}^{d}:\\left\\|\\widehat{\\theta}_{t}-\\theta\\right\\|_{\\bar{V}_{t}} \\leq R \\sqrt{d \\log \\left(\\frac{1+t L^{2} / \\lambda}{\\delta}\\right)}+\\lambda^{1 / 2} S\\right\\}\n",
                "$$\n",
                "\n",
                "一個目の信頼区間は\n",
                "$$\n",
                "\\begin{aligned}\n",
                "\\widehat{\\theta}_t & =\\left(\\mathbf{X}^{\\top} \\mathbf{X}+\\lambda I\\right)^{-1} \\mathbf{X}^{\\top}\\left(\\mathbf{X} \\theta_*+\\eta\\right) \\\\\n",
                "& =\\left(\\mathbf{X}^{\\top} \\mathbf{X}+\\lambda I\\right)^{-1} \\mathbf{X}^{\\top} \\eta+\\left(\\mathbf{X}^{\\top} \\mathbf{X}+\\lambda I\\right)^{-1}\\left(\\mathbf{X}^{\\top} \\mathbf{X}+\\lambda I\\right) \\theta_*-\\lambda\\left(\\mathbf{X}^{\\top} \\mathbf{X}+\\lambda I\\right)^{-1} \\theta_* \\\\\n",
                "& =\\left(\\mathbf{X}^{\\top} \\mathbf{X}+\\lambda I\\right)^{-1} \\mathbf{X}^{\\top} \\eta+\\theta_*-\\lambda\\left(\\mathbf{X}^{\\top} \\mathbf{X}+\\lambda I\\right)^{-1} \\theta_*,\n",
                "\\end{aligned}\n",
                "$$\n",
                "を使って推定値を真の値とノイズに分解してself-normalized boundを使えば簡単に出てきます。\n",
                "\n",
                "---\n",
                "\n",
                "このバウンドを実験してみましょう。論文のFigure2を再現します。\n",
                "今回は簡単のためにアームの数が有限であり、\n",
                "\n",
                "$$\n",
                "X_t = \\text{argmax}_{x \\in \\mathcal{X}_t} \\langle \\hat \\theta_{t-1}, x \\rangle + \\sqrt{\\beta_t} \\|x\\|_{V^{-1}_{t-1}}\n",
                "$$\n",
                "をアーム全て調べて選択します($\\beta_t$は信頼集合の右辺の値です)。\n",
                "今回の$\\beta_t$は前回[BANDIT_contextual.ipynb](BANDIT_contextual.ipynb)で使った$\\beta_t$ とほぼ同じなので、ほとんど同じ実装を使うことにします。"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 107,
            "metadata": {},
            "outputs": [
                {
                    "name": "stderr",
                    "output_type": "stream",
                    "text": [
                        "100%|██████████| 10/10 [01:42<00:00, 10.27s/it]\n"
                    ]
                },
                {
                    "data": {
                        "text/plain": [
                            "<matplotlib.legend.Legend at 0x7f83be6b5370>"
                        ]
                    },
                    "execution_count": 107,
                    "metadata": {},
                    "output_type": "execute_result"
                },
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABWz0lEQVR4nO3deVxU5f4H8M/srDMICgMqSm5I4pIr2m2l0LRbad2rPzMs0zK0XLKym1a2YJYtWuatW2I3zZstVpaaaVkBKuKSK24obgMqMgPi7Of3xzAjo6gM25nl8369zothzpkz38NR5sNznuc5EkEQBBARERF5EanYBRARERFdigGFiIiIvA4DChEREXkdBhQiIiLyOgwoRERE5HUYUIiIiMjrMKAQERGR12FAISIiIq8jF7uAurDb7Th58iTCw8MhkUjELoeIiIhqQRAElJeXIy4uDlLp1dtIfDKgnDx5Eq1btxa7DCIiIqqDY8eOoVWrVlfdxicDSnh4OADHAarVapGrISIiotowGAxo3bq163P8anwyoDgv66jVagYUIiIiH1Ob7hnsJEtERERehwGFiIiIvA4DChEREXkdn+yDUhuCIMBqtcJms4ldSsCSyWSQy+UcCk5ERB7zy4BiNptx6tQpVFZWil1KwAsJCUFsbCyUSqXYpRARkQ/xu4Bit9tRWFgImUyGuLg4KJVK/gUvAkEQYDabcfr0aRQWFqJDhw7XnJSHiIjIye8Citlsht1uR+vWrRESEiJ2OQEtODgYCoUCR48ehdlsRlBQkNglERGRj/DbP2n517p34HkgIqK64KcHEREReR0GFCIiIvI6DCgBavTo0bj33nvFLoOIiKhGDCheZPTo0ZBIJJBIJFAoFIiJicEdd9yBTz/9FHa7vUHf67333kNWVlaD7pOIiKihMKB4mYEDB+LUqVM4cuQIVq1ahVtvvRVPPfUUhgwZAqvV2mDvo9FoEBER0WD7IyIi/6DTGzHio404WFIhah0BEVAEQUCl2drkiyAIHteqUqmg1WrRsmVL3HDDDXj++efx3XffYdWqVa4Wj7fffhvJyckIDQ1F69at8cQTT6Ci4uI/pKysLERERGDNmjXo3LkzwsLCXMHHiZd4iIjoUgW6cty3IBu5h8/ima921OlzrKH43TwoNblgsSFp5pomf989s9IQoqz/j/i2225Dt27d8M033+DRRx+FVCrFvHnzkJCQgMOHD+OJJ57AM888gwULFrheU1lZibfeegv//e9/IZVK8eCDD+Lpp5/GkiVL6l0PERH5n5xDZ/DYf/NRbrTiuhaheG94D1EnOg2IgOIPEhMT8ddffwEAJk2a5Hq+bdu2ePXVV/H444+7BRSLxYKFCxeiXbt2AIAJEyZg1qxZTVozERH5hhXbTmDaVztgsQno3bYZPn6oFyJCxL1FSUAElGCFDHtmpYnyvg1FEARXkv3ll1+QmZmJffv2wWAwwGq1wmg0orKy0jV7bkhIiCucAEBsbCxKSkoarB4iIvJ9giBgwW+H8OaaAgDA4K6xmPtANwQ14OdXXQVEQJFIJA1yqUVMe/fuRUJCAo4cOYIhQ4Zg/PjxeO211xAZGYk///wTY8aMgdlsdgUUhULh9nqJRCLqtUQiIvIuVpsdM7/fjaWbigAAY/+WgOmDOkMq9Y771/n2p3aAWL9+PXbu3InJkycjPz8fdrsdc+fOdU0j/+WXX4pcIRER+ZJKsxUTl27Dun0lkEiAF4ckYfSABLHLcsOA4mVMJhN0Oh1sNhuKi4uxevVqZGZmYsiQIXjooYewa9cuWCwWzJ8/H3fffTeys7OxcOFCscsmIiIfcbrchDGL8/DXcT1UcineG94DA7toxS7rMgExzNiXrF69GrGxsWjbti0GDhyIX3/9FfPmzcN3330HmUyGbt264e2338Ybb7yBLl26YMmSJcjMzBS7bCIi8gGHTldg6IfZ+Ou4Hs1CFFg6tp9XhhMAkAg+2DHBYDBAo9FAr9dDrVa7rTMajSgsLERCQgKCgoJEqpCceD6IiLzDliOlePSzLSirtKBNVAiyHu6DhOahTVrD1T6/L8VLPERERH7up52nMOl/22G22tG9dQT+k94LzcNUYpd1VQwoREREfuw/fxzGaz/thSAAdyTFYN7wHghWij+M+Fo86oNis9kwY8YMJCQkIDg4GO3atcMrr7ziNnxVEATMnDkTsbGxCA4ORmpqKg4cOOC2n9LSUowcORJqtRoREREYM2aM21TtREREVD82u4CXf9iNV390hJOHUtpg4YM9fSKcAB4GlDfeeAMffvgh3n//fezduxdvvPEG5syZg/nz57u2mTNnDubNm4eFCxdi06ZNCA0NRVpaGoxGo2ubkSNHYvfu3Vi7di1WrlyJ33//HePGjWu4oyIiIgpgRosNGUu2YlH2EQDA83cl4uW/Xw+Zl8xxUhseXeLJycnBPffcg8GDBwNwTLP+xRdfYPPmzQAcrSfvvvsuXnjhBdxzzz0AgM8++wwxMTFYsWIFhg8fjr1792L16tXIy8tDr169AADz58/HXXfdhbfeegtxcXENeXxEREQBpfS8GY8uzsPWojIoZVLM/Uc33N3N9z5bPWpB6d+/P9atW4f9+/cDAHbs2IE///wTgwYNAgAUFhZCp9MhNTXV9RqNRoO+ffsiNzcXAJCbm4uIiAhXOAGA1NRUSKVSbNq0qcb3NZlMMBgMbgsRERG5O3r2PIZ9mIOtRWVQB8nx3zF9fDKcAB62oDz33HMwGAxITEyETCaDzWbDa6+9hpEjRwIAdDodACAmJsbtdTExMa51Op0O0dHR7kXI5YiMjHRtc6nMzEy8/PLLnpRKREQUULYfK8OYrDycPW9Gy4hgLH6kN9pHh4tdVp151ILy5ZdfYsmSJVi6dCm2bt2KxYsX46233sLixYsbqz4AwPTp06HX613LsWPHGvX9iIiIfMnaPcUY/lEuzp43o0tLNb59or9PhxPAw4Aybdo0PPfccxg+fDiSk5MxatQoTJ482TWTqVbrmI2uuLjY7XXFxcWudVqt9rK76lqtVpSWlrq2uZRKpYJarXZb6OpGjx6Ne++91/X9LbfcgkmTJjX4+7Rt2xbvvvtug++XiIhq57+5R/DYf7fAaLHj5o4t8L9xKYhW+/7EmB4FlMrKStcN6pxkMhnsdjsAICEhAVqtFuvWrXOtNxgM2LRpE1JSUgAAKSkpKCsrQ35+vmub9evXw263o2/fvnU+EF8mCAJSU1ORlpZ22boFCxYgIiICx48fF6Gya8vLy+MILCIiEdjtAjJX7cWM73bDLgDDe7fGf9J7IVTlH1OceXQUd999N1577TXEx8fj+uuvx7Zt2/D222/jkUceAQBIJBJMmjQJr776Kjp06ICEhATMmDEDcXFxrr/mO3fujIEDB2Ls2LFYuHAhLBYLJkyYgOHDhwfsCB6JRIJFixYhOTkZ//73v/HYY48BcHQ6fuaZZ/Dhhx+iVatWIlfpzmw2Q6lUokWLFmKXQkQUcExWG55e/hd+2HESADD1jo6YcFt7SCS+M4z4WjxqQZk/fz7uv/9+PPHEE+jcuTOefvppPPbYY3jllVdc2zzzzDOYOHEixo0bh969e6OiogKrV692uw/LkiVLkJiYiNtvvx133XUXbrzxRnz00UcNd1Q+qHXr1njvvffw9NNPo7CwEIIgYMyYMbjjjjvw22+/uSbH69SpE9577z2319psNkyZMgURERGIiorCM888g2vdYunHH3+ERqPBkiVLAAA7d+7EbbfdhuDgYERFRWHcuHFuk+c5Lxm99tpriIuLQ6dOnQDwEg8RUVPTV1rw0Ceb8cOOk5BLJZj7QDdMvL2DX4UTwMMWlPDwcLz77rtX/UCSSCSYNWsWZs2adcVtIiMjsXTpUk/eun4EAbBUNt37OSlCAA/+waSnp+Pbb7/FI488gqFDh2LXrl3YvXs33n//fSxfvhxRUVHIycnBuHHjEBsbi3/84x8AgLlz5yIrKwuffvopOnfujLlz5+Lbb7/FbbfdVuP7LF26FI8//jiWLl2KIUOG4Pz580hLS0NKSgry8vJQUlKCRx99FBMmTEBWVpbrdevWrYNarcbatWvr9WMhIqK6OX6uEg8vysOBkgqEqeRY+GBP3NihudhlNQr/uFB1LZZK4HURLh89fxJQenanyI8++gjXX389fv/9d3z99ddo0aKF2xDrhIQE5Obm4ssvv3QFlHfffRfTp0/H0KFDAQALFy7EmjVratz/Bx98gH/961/44YcfcPPNNwNwBBaj0YjPPvsMoaGOet9//33cfffdeOONN1zDxkNDQ/Gf//wHSqXSs58DERHV264TejyclYfT5SZo1UFY9HBvdI7130EjgRFQfEh0dDQee+wxrFixwtVv54MPPsCnn36KoqIiXLhwAWazGd27dwcA6PV6nDp1yq2DsVwuR69evS67zPPVV1+hpKQE2dnZ6N27t+v5vXv3olu3bq5wAgADBgyA3W5HQUGBK6AkJycznBARiWDD/tN44vN8nDfb0CkmHFmP9EasJljsshpVYAQURYijNUOM960DuVwOudxxapYtW4ann34ac+fORUpKCsLDw/Hmm29ecdbdq+nRowe2bt2KTz/9FL169fL4emX1AENERE3jy7xjmP7tTtjsAga0j8KHD/aEOkghdlmNLjACikTi8aUWb5GdnY3+/fvjiSeecD136NAh12ONRoPY2Fhs2rQJN910EwDHvDL5+fm44YYb3PbVrl07zJ07F7fccgtkMhnef/99AI6RVVlZWTh//rwrhGRnZ0Mqlbo6wxIRUdMSBAHv/HIA89YdAAAM7dESs4d1hVLu0fgWnxUYR+nDOnTogC1btmDNmjXYv38/ZsyYgby8PLdtnnrqKcyePRsrVqzAvn378MQTT6CsrKzG/XXs2BG//vorvv76a9fEbSNHjkRQUBDS09Oxa9cu/Prrr5g4cSJGjRp12W0LiIio8Vlsdkz76i9XOJlwa3vM/Ue3gAknAAOK13vssccwdOhQ/POf/0Tfvn1x9uxZt9YUAJg6dSpGjRqF9PR012Wg++6774r77NSpE9avX48vvvgCU6dORUhICNasWYPS0lL07t0b999/P26//XZXCwsRETWdcqMFj2Tl4av845BJJXj9vmQ8ndbJ74YRX4tEuNaEGV7IYDBAo9FAr9dfNu290WhEYWEhEhIS3OZeIXHwfBAR1Z5Ob8TDWXnYe8qAEKUMH/zfDbg1MfraL/QRV/v8vlRg9EEhIiLycgW6coxetBmn9EY0D1Nh0ejeSG6lEbss0TCgEBERiSzn0Bk89t98lButaNciFFkP90HryLqNBPUXDChEREQiWrHtBKZ9tQMWm4A+bSPx0UM9ERHCOacYUIiIiEQgCAIW/HYIb64pAAAM7hqLuQ90Q5BCJnJl3oEBhYiIqIlZbXbM/H43lm4qAgCMu+k6PDcwEVJpYI3UuRq/DSg+ODjJL/E8EBG5O2+yYuIX27B+XwkkEuClu69Hev+2YpfldfwuoCgUjul/KysrERzs3/cp8AWVlY67SDvPCxFRIDtdbsIjWXnYeUIPlVyK94b3wMAuWrHL8kp+F1BkMhkiIiJQUlICAAgJCQm4yW28gSAIqKysRElJCSIiIiCT8ZoqEQW2Q6crMHrRZhwrvYDIUCX+k94LN8Q3E7ssr+V3AQUAtFpHGnWGFBJPRESE63wQEQWqvCOlGPvZFpRVWtAmKgSLH+6Dts198x5xTcUvA4pEIkFsbCyio6NhsVjELidgKRQKtpwQUcD78a9TmPzldpitdnRvHYFP0nshKkwldllezy8DipNMJuMHJBERieY/fxzGaz/thSAAdyTFYN7wHghW8nOpNvw6oBBRw9FfsODJL7bBYLQgc2gyErVXv48GUSCz2QW8+uMeLMo+AgBIT2mDmXdfDxmHEdca72ZMRNekv2DBQ59uxob9p7GtqAz3vJ+N/248ymHkRDUwWmx4Ykm+K5w8f1ciXvo7w4mnGFCI6KoMRkc42XGsDBEhCvytQ3OYrHbMWLEL4z/fCn0l+3kROZWeN+P/Pt6INbuLoZRJMX9ED4y7qR1Hk9YBL/EQ0RUZjBaM+uRiOFn6aD90jg3HJ38W4o3V+7B6tw47T+jx3vDu6NU2UuxyyQcZLTZk/rQX32w7AbtdgFQigUQCSKUSSCXOBa6vEokEUikgq1onca27+FgmrbZttfVS6eXbutZL3betvq9Lt3XuVya9fF+/7C1B4ZnzUAfJ8fFDvdD3uiixf8Q+SyL4YButwWCARqOBXq+HWs3r4ESNoaZwkhR38f/bzuN6TPxiK46crYRMKsHk1A4Yf0t7NmNTrRXoyjHxi63YX1whdikNqmVEMBY/0hvto8PFLsXrePL5zYBCRJcxGC146JPN2F4VTpY82hfXx2ku267CZMWMFbvw7bYTAICU66Lw7vDuiFEHNXXJ5EMEQcDnm4rw6so9MFntaBGuQuZ9yegYEw67IMAmCBAEAXYBsAsC7Paqr9Wec623V9vu0vU1vM4uoOq1Amz2S/ZVbb3N7r7txf2ial0N6+0CQlVyDO/TGtHh/D9QEwYUIqqz2oaT6r7OP44Z3+1CpdmGyFAl5j7QDbcmRjdRxeRLyirNeOarv/DznmIAwC2dWuCtB7qhOecFCQgMKERUJ3UJJ06HTldg4tJt2HPKAAB49MYEPDMwEUo5++KTw6bDZzHpf9txSm+EQibBc4M64+H+bXkH3wDCgEJEHjMYLUj/dDO2FTnCyedj+qJLy9qFEyeT1YbMn/YhK+cIACC5pQbzR/TglN4BzmqzY976g3h//QHYBeC65qGYN6KHx/++yPcxoBCRR8qrhhJvKyqDJtjRclKfD4+1e4ox7asdKKu0IFQpw2v3JePeHi0bsGLyFcfPVWLSsu3YcvQcAOCBnq3w0t+vR6iKg0gDEQMKEdVaQ4cTp1P6C3hq2XZsLiwFAAy7oRVm3cMPpkDy085TeO7rv2AwWhGukuPV+7rgnu4MqoGMAYWIaqWxwomTzS5g/voDmLfuYtP+/P/rUet+LeSbLphtmLVyD77YXAQA6N46AvOG90B8VIjIlZHYGFCI6JrKq/qcbG2kcFLdpsNn8dSy7dAZjFDKpHj+rkSk92/L2TX90N5TBkz8YhsOllRAIgHG39wOk+/oCIWMnaWJAYWIrqEpw4nTufNmTPvqL/yy1zG8NLVzDN68vyuahSob9X2paQiCgP9uPIpXf9wLs9WO6HAV3vlndwxo31zs0siLePL57VGkbdvW8RfPpUtGRgYAwGg0IiMjA1FRUQgLC8OwYcNQXFzsto+ioiIMHjwYISEhiI6OxrRp02C1Wj08RCKqKzHCCQA0C1Xi44d64uW/Xw+lTIpf9hZj0Ht/YOPhs43+3tS4Ss+bMfazfMz8bjfMVjtuS4zGqqf+xnBC9eJRQMnLy8OpU6dcy9q1awEADzzwAABg8uTJ+OGHH7B8+XJs2LABJ0+exNChQ12vt9lsGDx4MMxmM3JycrB48WJkZWVh5syZDXhIRHQl5UYLRi/Ka/Jw4iSRSJDevy2+zeiP61qEQmcw4v8+3oh31u6Hze5zjbkEIOfQGQx673f8stdxc7yX7k7CJ+m9EMWJ16ie6nWJZ9KkSVi5ciUOHDgAg8GAFi1aYOnSpbj//vsBAPv27UPnzp2Rm5uLfv36YdWqVRgyZAhOnjyJmJgYAMDChQvx7LPP4vTp01Aqa9fUy0s8RJ6rMFmR/ulm5B89B3WQHEvH9hN1HorzJite/H43vso/DgDokxCJ94Z3R6wmWLSaqPYsNjve++UAPvjtIAQBaNfCMbcJO0DT1TTaJZ7qzGYzPv/8czzyyCOQSCTIz8+HxWJBamqqa5vExETEx8cjNzcXAJCbm4vk5GRXOAGAtLQ0GAwG7N69+4rvZTKZYDAY3BYiqr1Lw8mSR8UNJwAQqpLjrQe64d1/dkeoUobNhaUY9N4fWLun+NovJlEdK63EP/+di/d/dYST4b1b44eJNzKcUIOqc0BZsWIFysrKMHr0aACATqeDUqlERESE23YxMTHQ6XSubaqHE+d657oryczMhEajcS2tW7eua9lEAaemcJLcyns+SO7t0RI/Pvk3JLfUoKzSgrGfbcFL3++G0WITuzSqwcq/TuKueX9ga1EZwoPkeP//emD2sK4IUXJ+G2pYdQ4on3zyCQYNGoS4uLiGrKdG06dPh16vdy3Hjh1r9Pck8gcVJitGVwsnnz/a16vCiVPb5qH4enx/jP1bAgAgK+cIhi7IwaHTFSJXRk6VZiue/eovTFi6DeVGK26Ij8BPT/4NQ7o2/mcABaY6BZSjR4/il19+waOPPup6TqvVwmw2o6yszG3b4uJiaLVa1zaXjupxfu/cpiYqlQpqtdptIaKrc4aTLdXCSddWEWKXdUVKuRT/GpyERaN7IzJUiT2nDLh7/p9YvuUYfHA2BL+y+6QeQ+b/if9tOQaJBJh4W3t8+VgKWkdy4jVqPHUKKIsWLUJ0dDQGDx7seq5nz55QKBRYt26d67mCggIUFRUhJSUFAJCSkoKdO3eipKTEtc3atWuhVquRlJRU12MgoktUmKx4eJHvhJPqbq0aotq/XRQqzTZM++ovTP7fdlSYOB1BUxMEAYuyC3HfBzk4fPo8YtQqLH20H6be2QlyTrxGjczjUTx2ux0JCQkYMWIEZs+e7bZu/Pjx+Omnn5CVlQW1Wo2JEycCAHJycgA4hhl3794dcXFxmDNnDnQ6HUaNGoVHH30Ur7/+eq1r4CgeoitzhpO8I+cQHiTHEh8KJ9XZ7AIWbjiEt6uGILeJCsH8ET188lh80dkKE6Z99RfW73P8QZnaOQZz7u+KSE6sR/XQqDPJ/vzzz0hLS0NBQQE6duzots5oNGLq1Kn44osvYDKZkJaWhgULFrhdvjl69CjGjx+P3377DaGhoUhPT8fs2bMhl9e+gxUDClHN/CWcVJd/tBRPfrEdJ8ouQCGT4NmBiXhkQAKkUk6T31iyD57B5P9tR0m5CUq5FC8M7oxR/drw1gRUb5zqnigAnTdZMbpaOPl8TF90ax0hdlkNQl9pwbNf/4XVux2j/W7p1AJzH+jGycAamMVmxztr9+PDDYcgCED76DDMH9EDnWP5e5YaBgMKUYDx53DiJAgClmwqwqyVe1z3enn3n93Rn9OpN4iis5V4ctk2bD9WBgAY0SceM4ckIVgpE7cw8isMKEQB5LzJiocX5WHzkVK/DSfV7dMZMGHpxbvlPnFLO0xO7chOm/Xw3fYT+Ne3u1BhskIdJMfsYV1xV3Ks2GWRH2JAIQoQl4aT/47pi+5+HE6cLphtmLVyN77Y7JgTqWebZnhveHe0asZhr5649HYDvds2w7vDe6BlBG83QI2DAYUoAJw3WfFwVh42F5YiXCXHfx8NjHBS3Q87TuL5b3aivOov/zeGdcUg/uVfK7tO6DHxi20oPHMeUgkw8bYOmHhbe7ZEUaNiQCHycwwnFx0rrcSEL7ZhR1XfiZF94zFjSBKCFOw7URNBEPDJn4V4Y/U+WGwCYjVBePef3dH3uiixS6MAwIBC5McqzVaMXsRwUp3FZsdbPxfg3xsOAwA6xYRj/v/1QMeYcJEr8y5nKkx4evkO/FZwGgCQdn0M3hjWFREhnNuEmgYDCpGfujScfDamD3rENxO7LK/x+/7TmPLldpypMCNIIcVLd1+Pf/Zuzfk7APxx4DQm/28HzlSYoJJLMWNIEkb2jefPhpoUAwqRH2I4qZ2SciOmfrkDfxw4AwAY3DUWmUOToQ5SiFyZOMxWO+auvdi61DEmDPNH3IBOWrYuUdNjQCHyM5Vmx2idTQwntWK3C/joj8N4a00BrHYBrZoFY96IHrghwH5mR8+ex5NfbMOO43oAwIP94vHCYPbPIfEwoBD5kUvDyeIxfQLug7authWdw8QvtuH4uQuQSyWYemcnPHbTdQExTf6KbSfwwgrH3CaaYAXeGNYVA7tc+a7xRE2BAYXIT1SarXgkKw8bD5cirKrlhOHEMwajBdO/2Ykf/zoFAPhbh+aY+49uiA4PErmyxlFhsmLmd7vwzdYTAIA+CZF495/dEce5TcgLMKAQ+QGGk4YjCAL+l3cML/2wG0aLHc3DlJgxJAlto0IRFiRHuEqOsCA5ghUyn+40+tfxMjz5xTYcOVsJqQSYlNoRGbe2hywAWozINzCgEPm4SrMVY7K2IPfwWYaTBnSguBwTlm5DQXF5jeulEiBUdTGwhKnkCAtSIFwlR6hKhjCVwi3QhKmc21x8HB4kR6hKDkUTTnhmtzvmNpmzxjG3ScuIYLw7vDt6t41sshqIaoMBhciHXRpOFj/SBz3bMJw0FKPFhrk/F+C3gtM4b7Ki3GRFhcmKhv5NqJJLER50aYBRIEwlq/pecXG9yhFqqm9f21adS0ctDeqixeyhXaEJCcxRS+TdGFCIfNQFsw2PZOUxnDQxQRBQabZdDCxGR2gpN1pxvirAOL+vMFlw3mRzPa6otn2FyQqjxd6gtV2rVWf9vhLXvC8v3n09hnPeF/Jinnx+y5uoJiK6BoYT8UgkEoRWtWJE13NfFpvdEXSqQsuload6mHEFIZMVFUZH8HEEIUfwsQuAXQDKjY79QV/zeyZqwzF/RA904My55EcYUIi8wOXhpDfDiY9SyKSICFHWe/p4QRBwwWJzDzVGR5g5X61FRx0kxwO9WnNuE/I7DChEIrtgtmHMYkc4CVXKqsIJOzcGOolEghClHCFKOaLZMEIBiPfVJhKRM5zkHHKEk8/G9GE4ISICAwqRaC6YbXj0M4YTIqKaMKAQicAZTrIPOi/rMJwQEVXHgELUxC6YbRj72Ra3cNKLE2oREblhQCFqQs5w8ufBMwwnRERXwYBC1ESMFvdwksVwQkR0RQwoRE3AaLHh0cWOcBJSFU54nxQioivjPChEjaxAV44Z3+3C5sJShFRd1mE4ISK6OgYUokZSeOY83v1lP77fcRKCAIYTIiIPMKAQNbDj5yoxf91BfLX1OGx2x70470rWYsodndA+Okzk6oiIfAMDClEDKTYY8cGvB/HF5iJYbI5gcntiNCbf0RFdWmpEro6IyLcwoBDV09kKExZuOITPco/CZLUDAAa0j8KUOzrxhn9ERHXEgEJUR/pKCz7+4zA+zS5EpdkGAOjZphmm3tkR/ds1F7k6IiLfxoBC5KEKkxVZ2YX46PfDMBitAIDklhpMvbMjbu7YAhKJROQKiYh8n8fzoJw4cQIPPvggoqKiEBwcjOTkZGzZssW1XhAEzJw5E7GxsQgODkZqaioOHDjgto/S0lKMHDkSarUaERERGDNmDCoqKup/NESNyGix4ePfD+OmOb/irZ/3w2C0omNMGBY+2BPfTxiAWzpFM5wQETUQj1pQzp07hwEDBuDWW2/FqlWr0KJFCxw4cADNml28zj5nzhzMmzcPixcvRkJCAmbMmIG0tDTs2bMHQUFBAICRI0fi1KlTWLt2LSwWCx5++GGMGzcOS5cubdijI2oAJqsN/8s7hvfXH0RJuQkAkNA8FJNSO2BI1zjIpAwlREQNTSIIglDbjZ977jlkZ2fjjz/+qHG9IAiIi4vD1KlT8fTTTwMA9Ho9YmJikJWVheHDh2Pv3r1ISkpCXl4eevXqBQBYvXo17rrrLhw/fhxxcXHXrMNgMECj0UCv10OtVte2fCKPWG12fL31OOatO4gTZRcAAC0jgvFUagcM7dESchknYiYi8oQnn98e/Yb9/vvv0atXLzzwwAOIjo5Gjx498PHHH7vWFxYWQqfTITU11fWcRqNB3759kZubCwDIzc1FRESEK5wAQGpqKqRSKTZt2lTj+5pMJhgMBreFqLHY7AK+234CqW9vwLNf78SJsguIUavwyr1d8OvTt+AfvVoznBARNTKPLvEcPnwYH374IaZMmYLnn38eeXl5ePLJJ6FUKpGeng6dTgcAiImJcXtdTEyMa51Op0N0dLR7EXI5IiMjXdtcKjMzEy+//LInpRJ5TBAErNmtw9tr92N/saNPVGSoEk/c0g4P9muDIIVM5AqJiAKHRwHFbrejV69eeP311wEAPXr0wK5du7Bw4UKkp6c3SoEAMH36dEyZMsX1vcFgQOvWrRvt/SiwCIKA3wpOY+7aAuw64WidUwfJ8djN7TC6f1uEqjjYjYioqXn0mzc2NhZJSUluz3Xu3Blff/01AECr1QIAiouLERsb69qmuLgY3bt3d21TUlLitg+r1YrS0lLX6y+lUqmgUqk8KZWoVnIOnsFbPxdga1EZACBUKcOYGxMw5m/XQROsELc4IqIA5lFAGTBgAAoKCtye279/P9q0aQMASEhIgFarxbp161yBxGAwYNOmTRg/fjwAICUlBWVlZcjPz0fPnj0BAOvXr4fdbkffvn3rezxEtZJ/tBRzf96PnENnAQBBCinSU9risZvbITJUKXJ1RETkUUCZPHky+vfvj9dffx3/+Mc/sHnzZnz00Uf46KOPAAASiQSTJk3Cq6++ig4dOriGGcfFxeHee+8F4GhxGThwIMaOHYuFCxfCYrFgwoQJGD58eK1G8BDVx64Tesz9uQC/FpwGAChkEvxfn3hk3Noe0eogkasjIiInj4YZA8DKlSsxffp0HDhwAAkJCZgyZQrGjh3rWi8IAl588UV89NFHKCsrw4033ogFCxagY8eOrm1KS0sxYcIE/PDDD5BKpRg2bBjmzZuHsLDa3emVw4zJUwW6cryzdj9W73Z0xJZJJXigZytMuK09WjULEbk6IqLA4Mnnt8cBxRswoFBtFZ45j3d/2Y/vd5yEIAASCXBv95Z46vYOaNs8VOzyiIgCiief3xyeQH7p+LlKzF93EF9tPQ6b3ZHB70rWYlJqR3SMCRe5OiIiuhYGFPIrxQYjPvj1IL7YXASLzRFMbkuMxpQ7OqJLS43I1RERUW0xoJBfOFthwsINh/BZ7lGYrHYAwID2UZhyRyf0bNPsGq8mIiJvw4BCPk1facHHfxzGp9mFqDTbAAA92zTD1Ds7on+75iJXR0REdcWA4qd+KyhBudEKTbACESEKRAQroQlRIFwlh9QP7r5bYbIiK7sQH/1+GAajFQCQ3FKDqXd2xM0dW0Ai8f1jJCIKZAwofujPA2cwelFejeukEkAdrEBEsAKaEKUjwLhCTA3PhSiqtldCKRf/BnlGiw3/zT2KDzccQul5MwCgY0wYptzRCWnXxzCYEBH5CQYUP7R+n+NWArGaIGiCFSirtEB/wYILFhvsAlBWaUFZpQU4W+nRfkOUMkQEVwWWqlaZiBAFNCGKqlCjdAUd1zYhSoQqZfUODiarDf/LO4b31x9ESbkJAJDQPBSTUjtgSNc4yPygVYiIiC5iQPFDOYfOAACev6sz7u52cXZeo8UGwwULyi44AosjqJihr/79Bcdzzu3KKi0wGC0QBKDSbEOl2YaTeqNH9cilEmiCHUHG0TKjvCToOJ67dBt1kOOf59dbj2PeuoM4UXYBANAyIhhPpXbA0B4tIZeJ36pDREQNjwHFz5ypMGGfrhwA0L9dlNu6IIUMQQqZx1O62+0Cyo1WlF0wu0KM/oIF+kr37x0tNWZXi01ZpQVmmx1Wu4Cz5804W3VJxhNKuRTmqlE50eEqTLytPf7RuzVUcpnH+yIiIt/BgOJncqtufpeoDUdUWMPcAVoqlTgu44Qo0Cbq2ts7CYIAo8WOsgvmai02jhBTvcVG7ww0znBTaUG5ydHx1Wy1IzJUiSduaYcH+7VBkILBhIgoEDCg+Bnn5R1vGGIrkUgQrJQhWBmMWE2wR6+12uwwGK3QX7AgVhPEYEJEFGAYUPxM9kFHC8qA9h40dXghuUyKyFAlIkOVYpdCREQiYA9DP3KstBJFpZWQSSXokxApdjlERER1xoDiR5z9T7q10iA8SCFyNURERHXHgOJHsr2o/wkREVF9MKD4CUEQkFPVgtLfx/ufEBERMaD4iQMlFThdboJKLsUN8bx7LxER+TYGFD+RfdBxead320gOySUiIp/HgOIneHmHiIj8CQOKH7Da7Nh4uCqgsIMsERH5AQYUP7DrpAHlRivCg+RIbqkRuxwiIqJ6Y0DxA87+J/2ui4JMKhG5GiIiovpjQPEDzgnaBrRj/xMiIvIPDCg+zmixIe9IKQCgf3v2PyEiIv/AgOLjthadg8lqR4twFTpEh4ldDhERUYNgQPFxOQedo3eiIJGw/wkREfkHBhQf57z/zgAOLyYiIj/CgOLDyo0W/HVcD4ATtBERkX9hQPFhmwtLYbMLiI8MQatmIWKXQ0RE1GAYUHxYdlX/kwFsPSEiIj/DgOLDcqr6n3B6eyIi8jceBZSXXnoJEonEbUlMTHStNxqNyMjIQFRUFMLCwjBs2DAUFxe77aOoqAiDBw9GSEgIoqOjMW3aNFit1oY5mgBypsKEfbpyAI4RPERERP5E7ukLrr/+evzyyy8XdyC/uIvJkyfjxx9/xPLly6HRaDBhwgQMHToU2dnZAACbzYbBgwdDq9UiJycHp06dwkMPPQSFQoHXX3+9AQ4ncDhnj03UhiMqTCVyNURERA3L44Ail8uh1Wove16v1+OTTz7B0qVLcdtttwEAFi1ahM6dO2Pjxo3o168ffv75Z+zZswe//PILYmJi0L17d7zyyit49tln8dJLL0GpVNb/iAIEL+8QEZE/87gPyoEDBxAXF4frrrsOI0eORFFREQAgPz8fFosFqamprm0TExMRHx+P3NxcAEBubi6Sk5MRExPj2iYtLQ0GgwG7d++u77EEFHaQJSIif+ZRC0rfvn2RlZWFTp064dSpU3j55Zfxt7/9Dbt27YJOp4NSqURERITba2JiYqDT6QAAOp3OLZw41zvXXYnJZILJZHJ9bzAYPCnb7xwrrURRaSVkUgn6JESKXQ4REVGD8yigDBo0yPW4a9eu6Nu3L9q0aYMvv/wSwcHBDV6cU2ZmJl5++eVG27+vcfY/6dZKg/AghcjVEBERNbx6DTOOiIhAx44dcfDgQWi1WpjNZpSVlbltU1xc7OqzotVqLxvV4/y+pn4tTtOnT4der3ctx44dq0/ZPi+b/U+IiMjP1SugVFRU4NChQ4iNjUXPnj2hUCiwbt061/qCggIUFRUhJSUFAJCSkoKdO3eipKTEtc3atWuhVquRlJR0xfdRqVRQq9VuS6ASBAE5VS0onN6eiIj8lUeXeJ5++mncfffdaNOmDU6ePIkXX3wRMpkMI0aMgEajwZgxYzBlyhRERkZCrVZj4sSJSElJQb9+/QAAd955J5KSkjBq1CjMmTMHOp0OL7zwAjIyMqBScahsbRwoqcDpchNUciluiG8mdjlERESNwqOAcvz4cYwYMQJnz55FixYtcOONN2Ljxo1o0aIFAOCdd96BVCrFsGHDYDKZkJaWhgULFrheL5PJsHLlSowfPx4pKSkIDQ1Feno6Zs2a1bBH5cdyDjou7/RuG4kghUzkaoiIiBqHRBAEQewiPGUwGKDRaKDX6wPucs/Yz7Zg7Z5iPDOwE564pb3Y5RAREdWaJ5/fvBePD7Ha7Nh4uKr/CTvIEhGRH2NA8SG7ThpQbrQiPEiO5JYascshIiJqNAwoPiS7qv9Jv+uiIJNKRK6GiIio8TCg+BDnBG0DePdiIiLycwwoPsJosSHvSCkAoH979j8hIiL/xoDiI7YWnYPJakeLcBU6RIeJXQ4REVGjYkDxETkHnaN3oiCRsP8JERH5NwYUH5FTdf+dARxeTEREAYABxQeUGy3YcVwPgPffISKiwMCA4gM2F5bCZhcQHxmCVs1CxC6HiIio0TGg+IDsqv4nA9h6QkREAYIBxQc4+59wensiIgoUDChe7kyFCft05QAcI3iIiIgCAQOKl3POHpuoDUdUmErkaoiIiJoGA4qX4+UdIiIKRAwoXo4dZImIKBAxoHixY6WVKCqthEwqQZ+ESLHLISIiajIMKF7M2f+kWysNwoMUIldDRETUdBhQvFg2+58QEVGAYkDxUoIgIKeqBYXT2xMRUaBhQPFSB0oqcLrcBJVcihvim4ldDhERUZNiQPFSOQcdl3d6t41EkEImcjVERERNiwHFS2Xz8g4REQUwBhQvZLXZsfFwVUBhB1kiIgpADCheaNdJA8qNVoQHyZHcUiN2OURERE2OAcULOae373ddFGRSicjVEBERNT0GFC+U45zenncvJiKiAMWA4mWMFhvyjpQCAPq3Z/8TIiIKTAwoXmZr0TmYrHa0CFehQ3SY2OUQERGJggHFyzgv7/RvFwWJhP1PiIgoMDGgeBlnB9kBHF5MREQBjAHFi5QbLdhxXA+AE7QREVFgq1dAmT17NiQSCSZNmuR6zmg0IiMjA1FRUQgLC8OwYcNQXFzs9rqioiIMHjwYISEhiI6OxrRp02C1WutTil/YXFgKm11AfGQIWjULEbscIiIi0dQ5oOTl5eHf//43unbt6vb85MmT8cMPP2D58uXYsGEDTp48iaFDh7rW22w2DB48GGazGTk5OVi8eDGysrIwc+bMuh+Fn8h2Di9m6wkREQW4OgWUiooKjBw5Eh9//DGaNbt4p129Xo9PPvkEb7/9Nm677Tb07NkTixYtQk5ODjZu3AgA+Pnnn7Fnzx58/vnn6N69OwYNGoRXXnkFH3zwAcxmc8MclY9y9j/h9PZERBTo6hRQMjIyMHjwYKSmpro9n5+fD4vF4vZ8YmIi4uPjkZubCwDIzc1FcnIyYmJiXNukpaXBYDBg9+7dNb6fyWSCwWBwW/zNmQoT9unKAThG8BAREQUyuacvWLZsGbZu3Yq8vLzL1ul0OiiVSkRERLg9HxMTA51O59qmejhxrneuq0lmZiZefvllT0v1KblVdy9O1IYjKkwlcjVERETi8qgF5dixY3jqqaewZMkSBAUFNVZNl5k+fTr0er1rOXbsWJO9d1Ph5R0iIqKLPAoo+fn5KCkpwQ033AC5XA65XI4NGzZg3rx5kMvliImJgdlsRllZmdvriouLodVqAQBarfayUT3O753bXEqlUkGtVrst/ibnEDvIEhEROXkUUG6//Xbs3LkT27dvdy29evXCyJEjXY8VCgXWrVvnek1BQQGKioqQkpICAEhJScHOnTtRUlLi2mbt2rVQq9VISkpqoMPyLcfPVeLo2UrIpBL0SYgUuxwiIiLRedQHJTw8HF26dHF7LjQ0FFFRUa7nx4wZgylTpiAyMhJqtRoTJ05ESkoK+vXrBwC48847kZSUhFGjRmHOnDnQ6XR44YUXkJGRAZUqMPteOKe379ZKg/AghcjVEBERic/jTrLX8s4770AqlWLYsGEwmUxIS0vDggULXOtlMhlWrlyJ8ePHIyUlBaGhoUhPT8esWbMauhSfkc3+J0RERG4kgiAIYhfhKYPBAI1GA71e7/P9UQRBQJ/X1+F0uQlLx/ZlSCEiIr/lyec378UjsoMlFThdboJKLsUN8c2u/QIiIqIAwIAisuyDjss7vdtGIkghE7kaIiIi78CAIrLsquHFvHsxERHRRQwoIrLa7Nh4uCqgsO8JERGRCwOKiHafNKDcaEV4kBzJLTVil0NEROQ1GFBE5Bxe3O+6KMikEpGrISIi8h4MKCJyTtA2gHcvJiIicsOAIhKjxYa8I6UAgP7t2f+EiIioOgYUkWwtOgeT1Y4W4Sp0iA4TuxwiIiKvwoAiklzn8OJ2UZBI2P+EiIioOgYUkTgnaBvA4cVERESXYUARQbnRgh3H9QA4QRsREVFNGFBEsLmwFDa7gPjIELRqFiJ2OURERF6HAUUEOVX9Twaw9YSIiKhGDCgicPY/4fT2RERENWNAaWJnKkzYpysH4BjBQ0RERJdjQGlizuHFidpwRIWpRK6GiIjIOzGgNLGcQ7y8Q0REdC0MKE2MHWSJiIiujQGlCR0/V4mjZyshk0rQJyFS7HKIiIi8FgNKE3LevbhbKw3CgxQiV0NEROS9GFCaUDb7nxAREdUKA0oTEQTB1f+E09sTERFdHQNKEzlYUoHT5Sao5FLcEN9M7HKIiIi8GgNKE3HOHtu7bSSCFDKRqyEiIvJuDChNJJuXd4iIiGqNAaUJWG12bDxcFVDYQZaIiOiaGFCawO6TBpQbrQgPkiO5pUbscoiIiLweA0oTcA4v7nddFGRSicjVEBEReT8GlCbgnKBtAO9eTEREVCsMKI3MaLEh70gpAGBAe/Y/ISIiqg0GlEa2ragMJqsdLcJVaB8dJnY5REREPsGjgPLhhx+ia9euUKvVUKvVSElJwapVq1zrjUYjMjIyEBUVhbCwMAwbNgzFxcVu+ygqKsLgwYMREhKC6OhoTJs2DVartWGOxgvluKa3j4JEwv4nREREteFRQGnVqhVmz56N/Px8bNmyBbfddhvuuece7N69GwAwefJk/PDDD1i+fDk2bNiAkydPYujQoa7X22w2DB48GGazGTk5OVi8eDGysrIwc+bMhj0qL+KcoG0AhxcTERHVmkQQBKE+O4iMjMSbb76J+++/Hy1atMDSpUtx//33AwD27duHzp07Izc3F/369cOqVaswZMgQnDx5EjExMQCAhQsX4tlnn8Xp06ehVCpr9Z4GgwEajQZ6vR5qtbo+5TeqcqMF3Wethc0u4M9nb0WrZiFil0RERCQaTz6/69wHxWazYdmyZTh//jxSUlKQn58Pi8WC1NRU1zaJiYmIj49Hbm4uACA3NxfJycmucAIAaWlpMBgMrlaYmphMJhgMBrfFF2wuLIXNLiA+MoThhIiIyAMeB5SdO3ciLCwMKpUKjz/+OL799lskJSVBp9NBqVQiIiLCbfuYmBjodDoAgE6ncwsnzvXOdVeSmZkJjUbjWlq3bu1p2aJw3r14AKe3JyIi8ojHAaVTp07Yvn07Nm3ahPHjxyM9PR179uxpjNpcpk+fDr1e71qOHTvWqO/XUJz9Tzi9PRERkWfknr5AqVSiffv2AICePXsiLy8P7733Hv75z3/CbDajrKzMrRWluLgYWq0WAKDVarF582a3/TlH+Ti3qYlKpYJKpfK0VFGdqTBhn64cgGMEDxEREdVevedBsdvtMJlM6NmzJxQKBdatW+daV1BQgKKiIqSkpAAAUlJSsHPnTpSUlLi2Wbt2LdRqNZKSkupbilfJrbq8k6gNR1SYb4UrIiIisXnUgjJ9+nQMGjQI8fHxKC8vx9KlS/Hbb79hzZo10Gg0GDNmDKZMmYLIyEio1WpMnDgRKSkp6NevHwDgzjvvRFJSEkaNGoU5c+ZAp9PhhRdeQEZGhs+1kFyLs/8JL+8QERF5zqOAUlJSgoceeginTp2CRqNB165dsWbNGtxxxx0AgHfeeQdSqRTDhg2DyWRCWloaFixY4Hq9TCbDypUrMX78eKSkpCA0NBTp6emYNWtWwx6VF3BO0MYOskRERJ6r9zwoYvD2eVCOn6vEjW/8CplUgu0z70B4kELskoiIiETXJPOg0JU5717crZWG4YSIiKgOGFAaQfYhDi8mIiKqDwaUBiYIwsUOsux/QkREVCcMKA3sYEkFTpeboJJLcUN8M7HLISIi8kkMKA3MOXts77aRCFLIRK6GiIjINzGgNLBsXt4hIiKqNwaUBmSzC9h4mBO0ERER1RcDSgPadUKPcqMV4UFyJLfUiF0OERGRz2JAaUDO4cX9rouCTCoRuRoiIiLfxYDSgJwTtA3g3YuJiIjqhQGlgRgtNuQdKQUADGjP/idERET1wYDSQLYVlcFktaNFuArto8PELoeIiMinMaA0kBzX9PZRkEjY/4SIiKg+GFAaiHOCtgEcXkxERFRvDCgNoNxowY7jegCcoI2IiKghMKA0gLwjpbDZBcRHhqBVsxCxyyEiIvJ5DCgNINs5vJitJ0RERA2CAaUBOPufcHp7uiKrCbCaxa6CiBqbIIhdgd+Qi12ArztTYcI+XTkAxwgeCjCWC0C5zrFU6Ko9LgbKTwHlVV+NZY7tgzRAaIuqpXnV1+hqj6utC24GcEQYkfcRBKCiBDizv2o5AJwpcHzVHwOkCkARAiiCL1mqnpMHVVsfAiiCLll/yfaX7sP5epnCr39HMKDUU27V3YsTteGIClOJXA01GPP5moNG9e8rdIBR79l+jXrHcvbgtbeVyi8JMtcINoqguh0rEdXMZgXOHakWRKotV/u/b7cAJr1jaUwS2TUCTvUAFFKLkBTiHoKC1IAqvHGP4SoYUOop5xDvXuxTTBWXt3a4gke11g+Tofb7lAcD4TFAeCwQVvU1PAYI0wLhVUtYjGPb82eA86erLWeA8yXVHlc9b9QDdmtVGDpVuzqU4RcDS9ilrTLVH0c7WmekvMJLBAAwlVe1gux3bxU5e8gRNmoikQIRbYDmHYHmHRxfW3RyPCfYHa2rlkrHV+sF9+/dlkrAarzyOrfXVz0n2B01CDbAXO5YGkOPUcA97zfOvmuBAaWenBO0sYOsiATB8QvGFTxqaO1wfm+uqP1+FSFV4aJa0HB9Xy2QBGlq38waEgm06Hjt7aymaoGlFqHGZr74i+pc4bX3L5ECIc0vDy9hLWpurVGG1u74iLyVIDh+RzgvxZzZD5yuelx+8sqvU4QAUe0vBhBnGIlsJ06rpSAANku1QOMMODUFoGrf1yUkKYKb/viqYUCph+PnKnH0bCVkUgn6JESKXY5/MlUAhhM1XGa5pN+HpbL2+1SGubd0XKnlQxUu3vVduQrQtHQs1yIIjhaf82cc18XdgkwNoebCOcdfYOdLHEttKEIuhpaQSEClrmr+rf41oobn1I6WHbbWUFOxWYDSwxdbQk5XaxG5WktDaLR7AGneAWjeCVC39K5/vxIJIFc6luCIxn0vu71x938NDCj14Lx7cbdWGoQHKUSuxscJgiN4nPoL0P0F6HY6vpYerv0+VOqqoFHtskp47OUtHyJeU20UEomjFSdIA0S1u/b2NgtQedY9vFTUcJnJGWqcf1mVHXUsnhfo+Jmr1FV11hBiXF81NW+jDPOuDwkSn1F/eUvImf2OFkS7tebXSGRAZIJ7AGneEWje3nHZk9yJ/H+OAaUeLt5/h/1PPGK3OYLHqR1VQaQqjJw/XfP2Kk1VyIi5wuWWqoWXIWpHprj4M7sWQXBcFqseZC6cc7TYGA3VvlZ1/nV7zuC49ISqFh6TATAcr2PRkouBJUhzlYCjvhjWLmvJCfPrEQ9+SRAcLajVA4hzqSi+8uuUYdVaQqotkQmO1knyCQwodSQIArKdHWTZ/+TKLBeAkj2OEHKqqmWkeDdgOX/5thIpENUB0CYDsV0dX7VdHf0gSBwSZ+tHOBB5neevtxivHWJcX/UXv6/+2G6BI+RUjYrQH6vjsUirjkVzcXSCMtRx+cr1NQRQhFZ9DXF80F32XLVtlWGOwEees1kdrXNWk+OrsawqhDiDSAFw5mDNvyucwmMvCSEdHJdpwmMZRv0AA0odHSypwOlyE1RyKW6IZ9MgAKCy9OLlGWcYObPf0dP8UvJgIOb6amGkKxCd5PilT/5DEeRYwqLr9npBqPrwurS1plqouVrQcX5vtzr63TiHeTfk6E+p4iohJtQ7w4/dDthMVeHA5B4UnF9tV1lnNV/heaOj1azGdZe8X02/F2r8+codHVKrj5Rp3sHxx0yQunF+PuQVGFDqyDl7bO+2kQhSyESupokJAlBWdHkYuVLzfUiUI4Bok4HYbo6vUe0BaYD93MhzEsnF+RzCY+q2D0FwtORdFmLKHX1rzOervlY6/lo3O5+renzpNs51zn4OdsvF4NPQahN+ZIqrBAZTzUHD5mWzGksVjmNxjZap1irSrC1bqQIUA0odBczlHZvFcf330jBypQmImiVcvDTjvEzD5lYSk0RS1SIRUrt+N7VlNV8eYlxB5tKvlY6+PJcGocteV9l04cdF4piUS666yterrQsCZMorrK/++iu9ViV6Z0zyTgwodWCzC9h42A8naDMaHP1DdNVG0pTsrfmvLakCiO7sHkRirnd0TiQKBK6hno1wibd6+Llia05VsLFZrh0g5MpLAkO1r1I5/4Agr8SAUge7TuhRbrQiPEiO5JY++IHsnLBItxPQVRtJc6UhvSpNVatItc6rzTs5fukRUcNrzPBD5CMYUOogu2p4cb/roiCTevlfHnabY7rm6nOL6HZeeUivumW1/iJVXyPa8C8sIiJqUh4FlMzMTHzzzTfYt28fgoOD0b9/f7zxxhvo1KmTaxuj0YipU6di2bJlMJlMSEtLw4IFCxATc7GDW1FREcaPH49ff/0VYWFhSE9PR2ZmJuRy38hLzgnaBnjr3YstF4Df3wQKf68a0lvDLKsSqaMDWvUwEpMMhHrpMRERUUDxKBFs2LABGRkZ6N27N6xWK55//nnceeed2LNnD0JDHZNkTZ48GT/++COWL18OjUaDCRMmYOjQocjOzgYA2Gw2DB48GFqtFjk5OTh16hQeeughKBQKvP766w1/hA3MZLUh70gpAGBAey/tf/LT08C2zy9+rwgBYrq4X6aJThL9PgtERERXIhEEQajri0+fPo3o6Ghs2LABN910E/R6PVq0aIGlS5fi/vvvBwDs27cPnTt3Rm5uLvr164dVq1ZhyJAhOHnypKtVZeHChXj22Wdx+vRpKJXX7tdgMBig0Wig1+uhVjftOPjcQ2cx4uONaBGuwubnb4fE2y597Pgf8O04RwvJXW8CCTc7JtjikF4iIhKZJ5/f9Rrbpdc7hr5FRjpulJefnw+LxYLU1FTXNomJiYiPj0dubi4AIDc3F8nJyW6XfNLS0mAwGLB79+4a38dkMsFgMLgtYrk4vX2U94WT0/uBlZMdj29+Fuj9qGNCI4YTIiLyMXUOKHa7HZMmTcKAAQPQpUsXAIBOp4NSqURERITbtjExMdDpdK5tqocT53rnuppkZmZCo9G4ltatW9e17HpzTtA2wNuGF1suAMtHO4YjJtwE3DRN7IqIiIjqrM4BJSMjA7t27cKyZcsasp4aTZ8+HXq93rUcO1bHe3HUU7nRgh3HHa1GXjdB26pngZLdjluGD/0PW02IiMin1WnYzIQJE7By5Ur8/vvvaNWqlet5rVYLs9mMsrIyt1aU4uJiaLVa1zabN292219xcbFrXU1UKhVUKvHvQJl3pBQ2u4D4yBC0auZF94z5azmwdTEACTDs47pPCU5EROQlPGpBEQQBEyZMwLfffov169cjISHBbX3Pnj2hUCiwbt0613MFBQUoKipCSkoKACAlJQU7d+5ESUmJa5u1a9dCrVYjKSmpPsfS6LKdw4u9qfXkzEFg5STH45umAdfdImY1REREDcKjFpSMjAwsXboU3333HcLDw119RjQaDYKDg6HRaDBmzBhMmTIFkZGRUKvVmDhxIlJSUtCvXz8AwJ133omkpCSMGjUKc+bMgU6nwwsvvICMjAyvaCW5Gmf/E6+Z3t5idPQ7MVcAbf8G3PKc2BURERE1CI8CyocffggAuOWWW9yeX7RoEUaPHg0AeOeddyCVSjFs2DC3idqcZDIZVq5cifHjxyMlJQWhoaFIT0/HrFmz6nckjexMhQn7dOUAHCN4vMKa6UDxTiCkOTD0Y/Y7ISIiv1GveVDEIsY8KCv/OokJS7chURuO1ZNuapL3vKpdXwNfPQJAAjz4NdD+drErIiIiuqommwclkDj7n3jF5Z2zh4Dvn3I8/ttUhhMiIvI7DCi15JygTfQOsq5+J+VAfH/gluni1kNERNQIGFBq4fi5Shw9WwmZVII+CZHiFvPzC447EodEAfd/Ash84waLREREnmBAqQXn3Yu7tdIgPEghXiG7VwB5Hzse3/cRoI4TrxYiIqJGxIBSCxfvvyNi/5PSw8D3Ex2Pb5wMdEi9+vZEREQ+jAHlGgRBQPahqg6yYvU/sZqA5Q8DJgPQuh9w6wvi1EFERNREGFCu4WBJBU6Xm6CSS3FDfDNxivh5BnBqOxDcjP1OiIgoIDCgXINz9tjebSMRpBBhIrQ93wOb/+14fN+/AU2rq29PRETkBxhQriFHzMs7544A301wPO7/JNAxrelrICIiEgEDylXY7AI2HhZpgjaruarfiR5o1Qe4fWbTvj8REZGIGFCuYtcJPQxGK8KD5EhuqWnaN//lReDkViAoArj/U0Am4vBmIiKiJsaAchXZVcOL+10XBZlU0nRvvO9HYGPVDRbvWwhEtG669yYiIvICDChX4ZygbUBT3r343FFgxXjH45QJQKdBTffeREREXoIB5QpMVhvyjpQCAAa0b6L+J1az4w7FRj3Qsidw+4tN875ERERehgHlCrYeLYPJakeLcBXaR4c1zZuuexk4sQUI0gD3LwLkyqZ5XyIiIi/DgHIFF6e3j4JE0gT9TwpWAbnvOx7fswBo1qbx35OIiMhLMaBcgXOCtgFNMby47Bjw7eOOx33HA52HNP57EhEReTEGlBpUmKzYcVwPoAkmaLNZqvqdlAFxPYA7ZjXu+xEREfkABpQabC48C5tdQHxkCFo1C2ncN1v/CnB8M6BSs98JERFRFQaUGmQ7hxc3duvJ/jVA9nuOx/e8D0QmNO77ERER+QgGlBo4+5806vT2+hMX+530GQck3dN470VERORjGFAucabChH26cgCOETyNwmYFvh4DXCgFYrsBd77aOO9DRETkoxhQLuG8OWCiNhxRYarGeZNfXwOKch39Th7IAuSN9D5EREQ+igHlEs7+J412eefAL8Cfbzse/30eEHld47wPERGRD2NAuYRzgrZG6SBrOAl8O87xuNcY4Pr7Gv49iIiI/AADSjXHz1Xi6NlKyKQS9EmIbNid26zA148ClWcBbTKQ9nrD7p+IiMiPMKBUk3PIcXmnWysNwoMUDbvz3zKBo9mAMgx4YDGgCGrY/RMREfkRudgFeBN1kAJ92kY2/N2LD60H/pjreHz3e0BUu4bdPxERkZ9hQKlmYBctBnbRNuxODaeAr8cCEICeo4Hk+xt2/0RERH6Il3gak90GfDMWqDwDxHQBBs4WuyIiIiKfwIDSmDa8ARz5A1CEOuY7UQSLXREREZFPYEBpLId/AzbMcTy++12geQcxqyEiIvIpHgeU33//HXfffTfi4uIgkUiwYsUKt/WCIGDmzJmIjY1FcHAwUlNTceDAAbdtSktLMXLkSKjVakRERGDMmDGoqKio14F4lfLii/1OeowCuv5D7IqIiIh8iscB5fz58+jWrRs++OCDGtfPmTMH8+bNw8KFC7Fp0yaEhoYiLS0NRqPRtc3IkSOxe/durF27FitXrsTvv/+OcePG1f0ovIndBnzzKHC+BIhOAgbNEbsiIiIinyMRBEGo84slEnz77be49957AThaT+Li4jB16lQ8/fTTAAC9Xo+YmBhkZWVh+PDh2Lt3L5KSkpCXl4devXoBAFavXo277roLx48fR1xc3DXf12AwQKPRQK/XQ61W17X8xvHbG8BvrwOKEGDcb0CLTmJXRERE5BU8+fxu0D4ohYWF0Ol0SE1NdT2n0WjQt29f5ObmAgByc3MRERHhCicAkJqaCqlUik2bNtW4X5PJBIPB4LZ4pcLfgQ1VI3UGv81wQkREVEcNGlB0Oh0AICYmxu35mJgY1zqdTofo6Gi39XK5HJGRka5tLpWZmQmNRuNaWrdu3ZBlN4yKEsdU9oId6D4S6D5C7IqIiIh8lk+M4pk+fTr0er1rOXbsmNglubPbgW/GARXFQItE4K43xa6IiIjIpzVoQNFqHbOwFhcXuz1fXFzsWqfValFSUuK23mq1orS01LXNpVQqFdRqtdviVf6cCxz+FZAHO+Y7UYaKXREREZFPa9CAkpCQAK1Wi3Xr1rmeMxgM2LRpE1JSUgAAKSkpKCsrQ35+vmub9evXw263o2/fvg1ZTtM48ifwa9WdiQe/BUR3FrceIiIiP+DxvXgqKipw8OBB1/eFhYXYvn07IiMjER8fj0mTJuHVV19Fhw4dkJCQgBkzZiAuLs410qdz584YOHAgxo4di4ULF8JisWDChAkYPnx4rUbweJXzZy72O+k2wtH3hIiIiOrN44CyZcsW3Hrrra7vp0yZAgBIT09HVlYWnnnmGZw/fx7jxo1DWVkZbrzxRqxevRpBQUGu1yxZsgQTJkzA7bffDqlUimHDhmHevHkNcDhNyNnvpPwU0LwjcNdbgEQidlVERER+oV7zoIjFK+ZB+eNtYN3LgDwIGLseiLlenDqIiIh8hGjzoASMo7nA+lcdjwfNYTghIiJqYAwonjp/FvjqEUCwAcn/AG54SOyKiIiI/A4DiifsdmDF40D5SSCqPTDkbfY7ISIiagQMKJ7ImQcc+BmQqRzznajCxa6IiIjILzGg1FbRJmDdLMfjQbMBbbK49RAREfkxBpTaqCy92O+kyzCg58NiV0REROTXGFCuRRCAFeMBw3Eg8jpgyLvsd0JERNTIGFCuJfd9YP9qQKZ09DsJ8rL7ABEREfkhBpSrOZYH/PKS4/HATCC2m6jlEBERBQoGlCtx9juxW4Gke4FeY8SuiIiIKGAwoNREEIDvMgB9EdCsLfD3eex3QkRE1IQYUGqy8UOg4Kdq/U40YldEREQUUBhQLnU8H1g70/H4zteAuB7i1kNERBSAGFCqu1AGfDUasFuAzncDfcaKXREREVFAYkCpbv2rQFkRENEG+Pv77HdCREQkErnYBXiV214ALpQCKRlAcITY1RAREQUsBpTqgiOA+z8VuwoiIqKAx0s8RERE5HUYUIiIiMjrMKAQERGR12FAISIiIq/DgEJERERehwGFiIiIvA4DChEREXkdBhQiIiLyOgwoRERE5HUYUIiIiMjrMKAQERGR12FAISIiIq/DgEJERERexyfvZiwIAgDAYDCIXAkRERHVlvNz2/k5fjU+GVDKy8sBAK1btxa5EiIiIvJUeXk5NBrNVbeRCLWJMV7Gbrfj5MmTCA8Ph0QiabD9GgwGtG7dGseOHYNarW6w/VLd8Hx4F54P78Lz4V14PmpHEASUl5cjLi4OUunVe5n4ZAuKVCpFq1atGm3/arWa/8C8CM+Hd+H58C48H96F5+PartVy4sROskREROR1GFCIiIjI6zCgVKNSqfDiiy9CpVKJXQqB58Pb8Hx4F54P78Lz0fB8spMsERER+Te2oBAREZHXYUAhIiIir8OAQkRERF6HAYWIiIi8DgNKlQ8++ABt27ZFUFAQ+vbti82bN4tdkl/KzMxE7969ER4ejujoaNx7770oKChw28ZoNCIjIwNRUVEICwvDsGHDUFxc7LZNUVERBg8ejJCQEERHR2PatGmwWq1NeSh+afbs2ZBIJJg0aZLrOZ6PpnXixAk8+OCDiIqKQnBwMJKTk7FlyxbXekEQMHPmTMTGxiI4OBipqak4cOCA2z5KS0sxcuRIqNVqREREYMyYMaioqGjqQ/F5NpsNM2bMQEJCAoKDg9GuXTu88sorbveR4floRAIJy5YtE5RKpfDpp58Ku3fvFsaOHStEREQIxcXFYpfmd9LS0oRFixYJu3btErZv3y7cddddQnx8vFBRUeHa5vHHHxdat24trFu3TtiyZYvQr18/oX///q71VqtV6NKli5Camips27ZN+Omnn4TmzZsL06dPF+OQ/MbmzZuFtm3bCl27dhWeeuop1/M8H02ntLRUaNOmjTB69Ghh06ZNwuHDh4U1a9YIBw8edG0ze/ZsQaPRCCtWrBB27Ngh/P3vfxcSEhKECxcuuLYZOHCg0K1bN2Hjxo3CH3/8IbRv314YMWKEGIfk01577TUhKipKWLlypVBYWCgsX75cCAsLE9577z3XNjwfjYcBRRCEPn36CBkZGa7vbTabEBcXJ2RmZopYVWAoKSkRAAgbNmwQBEEQysrKBIVCISxfvty1zd69ewUAQm5uriAIgvDTTz8JUqlU0Ol0rm0+/PBDQa1WCyaTqWkPwE+Ul5cLHTp0ENauXSvcfPPNroDC89G0nn32WeHGG2+84nq73S5otVrhzTffdD1XVlYmqFQq4YsvvhAEQRD27NkjABDy8vJc26xatUqQSCTCiRMnGq94PzR48GDhkUcecXtu6NChwsiRIwVB4PlobAF/icdsNiM/Px+pqamu56RSKVJTU5GbmytiZYFBr9cDACIjIwEA+fn5sFgsbucjMTER8fHxrvORm5uL5ORkxMTEuLZJS0uDwWDA7t27m7B6/5GRkYHBgwe7/dwBno+m9v3336NXr1544IEHEB0djR49euDjjz92rS8sLIROp3M7HxqNBn379nU7HxEREejVq5drm9TUVEilUmzatKnpDsYP9O/fH+vWrcP+/fsBADt27MCff/6JQYMGAeD5aGw+ebPAhnTmzBnYbDa3X64AEBMTg3379olUVWCw2+2YNGkSBgwYgC5dugAAdDodlEolIiIi3LaNiYmBTqdzbVPT+XKuI88sW7YMW7duRV5e3mXreD6a1uHDh/Hhhx9iypQpeP7555GXl4cnn3wSSqUS6enprp9nTT/v6ucjOjrabb1cLkdkZCTPh4eee+45GAwGJCYmQiaTwWaz4bXXXsPIkSMBgOejkQV8QCHxZGRkYNeuXfjzzz/FLiVgHTt2DE899RTWrl2LoKAgscsJeHa7Hb169cLrr78OAOjRowd27dqFhQsXIj09XeTqAs+XX36JJUuWYOnSpbj++uuxfft2TJo0CXFxcTwfTSDgL/E0b94cMpnsslEJxcXF0Gq1IlXl/yZMmICVK1fi119/RatWrVzPa7VamM1mlJWVuW1f/Xxotdoaz5dzHdVefn4+SkpKcMMNN0Aul0Mul2PDhg2YN28e5HI5YmJieD6aUGxsLJKSktye69y5M4qKigBc/Hle7feVVqtFSUmJ23qr1YrS0lKeDw9NmzYNzz33HIYPH47k5GSMGjUKkydPRmZmJgCej8YW8AFFqVSiZ8+eWLdunes5u92OdevWISUlRcTK/JMgCJgwYQK+/fZbrF+/HgkJCW7re/bsCYVC4XY+CgoKUFRU5DofKSkp2Llzp9t/+rVr10KtVl/2y52u7vbbb8fOnTuxfft219KrVy+MHDnS9Zjno+kMGDDgsmH3+/fvR5s2bQAACQkJ0Gq1bufDYDBg06ZNbuejrKwM+fn5rm3Wr18Pu92Ovn37NsFR+I/KykpIpe4fkzKZDHa7HQDPR6MTu5euN1i2bJmgUqmErKwsYc+ePcK4ceOEiIgIt1EJ1DDGjx8vaDQa4bfffhNOnTrlWiorK13bPP7440J8fLywfv16YcuWLUJKSoqQkpLiWu8c1nrnnXcK27dvF1avXi20aNGCw1obSPVRPILA89GUNm/eLMjlcuG1114TDhw4ICxZskQICQkRPv/8c9c2s2fPFiIiIoTvvvtO+Ouvv4R77rmnxmGtPXr0EDZt2iT8+eefQocOHTistQ7S09OFli1buoYZf/PNN0Lz5s2FZ555xrUNz0fjYUCpMn/+fCE+Pl5QKpVCnz59hI0bN4pdkl8CUOOyaNEi1zYXLlwQnnjiCaFZs2ZCSEiIcN999wmnTp1y28+RI0eEQYMGCcHBwULz5s2FqVOnChaLpYmPxj9dGlB4PprWDz/8IHTp0kVQqVRCYmKi8NFHH7mtt9vtwowZM4SYmBhBpVIJt99+u1BQUOC2zdmzZ4URI0YIYWFhglqtFh5++GGhvLy8KQ/DLxgMBuGpp54S4uPjhaCgIOG6664T/vWvf7kNn+f5aDwSQag2JR4RERGRFwj4PihERETkfRhQiIiIyOswoBAREZHXYUAhIiIir8OAQkRERF6HAYWIiIi8DgMKEREReR0GFCIiIvI6DChERETkdRhQiIiIyOswoBAREZHXYUAhIiIir/P/5jdPhRueQ4wAAAAASUVORK5CYII=",
                        "text/plain": [
                            "<Figure size 640x480 with 1 Axes>"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                }
            ],
            "source": [
                "import numpy as np\n",
                "import matplotlib.pyplot as plt\n",
                "import seaborn as sns\n",
                "import pandas as pd\n",
                "from tqdm import tqdm\n",
                "Array = np.ndarray\n",
                "\n",
                "\n",
                "def lin_ucb_yadkori(bandit, n_step, lam=0.1, delta=0.0001, seeds=10):\n",
                "    # Yadkori 2011の信頼集合を使った手法です\n",
                "    regrets = []\n",
                "    k = len(bandit.arms)\n",
                "    d = bandit.feature_dim\n",
                "    for _ in range(seeds):\n",
                "        bandit.reset()\n",
                "        arms = bandit.arms  # k x d\n",
                "        Vt = lam * np.eye(d)\n",
                "        theta_hat = np.zeros(d)\n",
                "        AX = np.zeros((d, 1))\n",
                "        pseudo_regret = 0\n",
                "        for i in range(n_step):\n",
                "            # UCBスコアの計算。\n",
                "            sqrt_beta_t = np.sqrt(lam) + np.sqrt(2 * np.log(1 / delta) + d * np.log((d * lam + (i+1)) / (d * lam)))\n",
                "            beta_t = sqrt_beta_t ** 2\n",
                "            theta_a = arms @ theta_hat  # k\n",
                "            Vt_inv = np.linalg.inv(Vt)\n",
                "            beta_a_norm = beta_t * np.array([arm @ Vt_inv @ arm for arm in arms])\n",
                "            ucb = theta_a + beta_a_norm\n",
                "\n",
                "            # UCBスコアが最大のアームを選択\n",
                "            arm = np.argmax(ucb)\n",
                "            mean, rew = bandit.pull(arm)\n",
                "\n",
                "            # 結果の記録\n",
                "            pseudo_regret += bandit.optimal_mean - mean\n",
                "            arm = arms[arm].reshape(d, 1)\n",
                "            Vt += arm @ arm.T\n",
                "            AX += arm * rew\n",
                "            theta_hat = (np.linalg.inv(Vt) @ AX).reshape(-1)\n",
                "            assert Vt.shape == (d, d)\n",
                "            assert AX.shape == (d, 1)\n",
                "    \n",
                "        regrets.append(pseudo_regret)\n",
                "    return np.mean(regrets)\n",
                "\n",
                "\n",
                "def lin_ucb_dani(bandit, n_step, lam=0.1, delta=0.0001, seeds=10):\n",
                "    # Dani 2008の信頼集合を使った手法です\n",
                "    regrets = []\n",
                "    k = len(bandit.arms)\n",
                "    d = bandit.feature_dim\n",
                "    for _ in range(seeds):\n",
                "        bandit.reset()\n",
                "        arms = bandit.arms  # k x d\n",
                "        Vt = lam * np.eye(d)\n",
                "        theta_hat = np.zeros(d)\n",
                "        AX = np.zeros((d, 1))\n",
                "        pseudo_regret = 0\n",
                "        for i in range(n_step):\n",
                "            # UCBスコアの計算。\n",
                "            sqrt_beta_t = max(np.log((i+1)**2 / delta), np.sqrt(d * np.log(i+1) * np.log((i+1)**2 / (delta))))\n",
                "            beta_t = sqrt_beta_t ** 2\n",
                "            theta_a = arms @ theta_hat  # k\n",
                "            Vt_inv = np.linalg.inv(Vt)\n",
                "            beta_a_norm = beta_t * np.array([arm @ Vt_inv @ arm for arm in arms])\n",
                "            ucb = theta_a + beta_a_norm\n",
                "\n",
                "            # UCBスコアが最大のアームを選択\n",
                "            arm = np.argmax(ucb)\n",
                "            mean, rew = bandit.pull(arm)\n",
                "\n",
                "            # 結果の記録\n",
                "            pseudo_regret += bandit.optimal_mean - mean\n",
                "            arm = arms[arm].reshape(d, 1)\n",
                "            Vt += arm @ arm.T\n",
                "            AX += arm * rew\n",
                "            theta_hat = (np.linalg.inv(Vt) @ AX).reshape(-1)\n",
                "            assert Vt.shape == (d, d)\n",
                "            assert AX.shape == (d, 1)\n",
                "    \n",
                "        regrets.append(pseudo_regret)\n",
                "    return np.mean(regrets)\n",
                "\n",
                "\n",
                "class StochasticLinearBandit:\n",
                "    def __init__(self, k: int, d: int, noise_scale: float):\n",
                "        self.k = k\n",
                "        self.d = d\n",
                "        self.noise_scale = noise_scale\n",
                "        self.theta_star = np.random.uniform(0, 1, size=(d,))\n",
                "        self.arms = np.random.uniform(0, 1, size=(k, d))\n",
                "        self.optimal_mean = np.max(self.arms @ self.theta_star)\n",
                "        self.feature_dim = self.arms.shape[1]\n",
                "\n",
                "    def reset(self):\n",
                "        self.theta_star = np.random.uniform(0, 1, size=(self.d,))\n",
                "        self.arms = np.random.uniform(0, 1, size=(self.k, self.d))\n",
                "        self.optimal_mean = np.max(self.arms @ self.theta_star)\n",
                "\n",
                "    def pull(self, arm):\n",
                "        A = self.arms[arm]\n",
                "        assert A.shape == (d,)\n",
                "        assert self.theta_star.shape == (d,)\n",
                "        mean = self.theta_star @ A\n",
                "        noise = np.random.randn() * self.noise_scale\n",
                "        reward = mean + noise\n",
                "        return mean, reward\n",
                "\n",
                "\n",
                "n_step = 1000\n",
                "d = 5\n",
                "dani_regrets = []\n",
                "yadkori_regrets = []\n",
                "num_arms = []\n",
                "for k in tqdm(range(2, 1000, 100)):\n",
                "    num_arms.append(k)\n",
                "    bandit = StochasticLinearBandit(k, d, 1.0)\n",
                "    dani_regrets.append(lin_ucb_dani(bandit, n_step, seeds=10))\n",
                "    yadkori_regrets.append(lin_ucb_yadkori(bandit, n_step, seeds=10))\n",
                "\n",
                "plt.plot(num_arms, dani_regrets, label=\"Dani\")\n",
                "plt.plot(num_arms, yadkori_regrets, label=\"Yadkori\")\n",
                "plt.legend()\n",
                "\n"
            ]
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "Python 3.9.12 ('shumi-VTLwuKSy-py3.9')",
            "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.16"
        },
        "orig_nbformat": 4,
        "vscode": {
            "interpreter": {
                "hash": "d6b7cac5e0d2ff733f340da4d53ae5ecfef7f7ad39623f5982b029a09306b36b"
            }
        }
    },
    "nbformat": 4,
    "nbformat_minor": 2
}