Skip to content

Instantly share code, notes, and snippets.

@vincentarelbundock
Created July 20, 2020 17:30
Show Gist options
  • Save vincentarelbundock/0c31d8abc6f285f79b5be272152be013 to your computer and use it in GitHub Desktop.
Save vincentarelbundock/0c31d8abc6f285f79b5be272152be013 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:numexpr.utils:NumExpr defaulting to 4 threads.\n",
"INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_9ba014013ce5b118c52a635c64466e4d NOW.\n",
"WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
"To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
"WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
"To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n",
"WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.\n",
"To run all diagnostics call pystan.check_hmc_diagnostics(fit)\n"
]
},
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fc318f73850>]"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9d5gd5Xn3/3lm5vRdafuqrdAKSaiABJKQaDHVFCGawSBcQhx+JjZgO34Tx3H8viHlda44xUleQxKTGHdJYAOiG1ODsQGxqqBitUValS1n+54+M8/vj5nTds9KW7VFz+e6uHZ3zpxznrMrvnPPfd/P9xZSShQKhUIxedHGegEKhUKhGF2U0CsUCsUkRwm9QqFQTHKU0CsUCsUkRwm9QqFQTHKMsV4AQEVFhZwzZ85YL0OhUCgmFFu2bAlLKStPdd64EPo5c+ZQV1c31stQKBSKCYUQ4vBAzlOpG4VCoZjkKKFXKBSKSY4SeoVCoZjkKKFXKBSKSY4SeoVCoZjkKKFXKBSKSY4SeoVCoZjkKKFXKBSKSY4SeoVCoRgDUpZNc1f8tLzXuNgZq1AoFGcKSdOmpSdBeySJrgmqpvhH/T2V0CsUCsVpIGFaNHcl6IylON2D/ZTQKxQKxSgST1m0dCfoiKbGbA2nzNELIWqEEG8IIfYIIXYJIb7iHi8TQrwihNjvfi3Nec43hBAHhBC/E0JcN5ofQKFQKMYjCdOioS3K/qaeMRV5GFgx1gT+REq5CLgIeEAIsRj4c+A1KeV84DX3Z9zH1gFLgOuBfxdC6KOxeIVCoRhvpAV+X+PYC3yaUwq9lPKElHKr+303sAeYCdwC/Mg97UfAre73twAbpZQJKWU9cABYNdILVygUivHEeIrgezOoHL0QYg5wAfAeUC2lPAHOxUAIUeWeNhN4N+dpR91jvV/rPuA+gNmzZw923QqFQjEuSJjZHPzpLrIOlAELvRCiCHgS+GMpZZcQot9TCxzr8/GllI8CjwKsXLlynP56FAqFoi/xlEU0aRFNmuNa4NMMSOiFEB4ckf+ZlPIp93CTEGK6G81PB5rd40eBmpynzwKOj9SCFQqFYizoiqfojpt0x1OkzHGu7L0YSNeNAL4P7JFSfifnoWeBe9zv7wGeyTm+TgjhE0LUAvOBzSO3ZIVCoTh99CRM9jd1czgcpa0nOeFEHgYW0V8KfBb4QAix3T32F8DfA08IIe4FjgCfBJBS7hJCPAHsxunYeUBKaY34yhUKhWIU6Y6naO5OEE1MfPk6pdBLKd+mcN4d4Op+nvMt4FvDWJdCoVCMCV3xFM1dCWLJiS/wadTOWIVCccbTGU3RnUgRSVgkTXuslzPiKKFXKBRnJNJtlemKmRxpi47xakYXJfQKheKM5Gh7jO64iWVPvOLqYFFCr1AozgiklDR3JxBAezQ1KVM0/aGEXqFQTGqiSZMTnXF0IeiOm2O9nDFBCb1CoZi0tHQnaOqKj/udq6ONEnqFQjGp6IymiKUsehLmpGqRHA5K6BUKxYRHSknKkhxo7kHTmJC7V0cTJfQKhWJCIqVECIFtS9qiSTqiKSxbYp05NdYBo4ReoVBMODqjKY51xCjyGXTGxpf3+3hkIBOmFAqFYtxgWjYnumJYtpzwIh/uSXCopWfU30dF9AqFYtySsmwMTRBLWaQsSXskOSlaJA+3Rnhq2zHe2tfC5Qsq+f4fXDiq76eEXqFQjEtSlk1LdwIhoLUnOeFbJKWUfHisk6e2HaPucDs+Q2PNedP52nXnjPp7K6FXKBTjBtuWJC0bXRMcaolMit2rli1551ArT209yv7mHqYGPHx69WzWnDudsiIvNWXBUV+DEnqFQjHmpA3Gwj0JokmLpGVPeJGPpyxe29vMpm3HaOyKM32qn/uvOJurFlbhM/TTuhYl9AqFYkxpjyRJmDbt0SSWLSd8iqYzluKFncd54YMTdMVNzqku5nOXzmF1bTm61u+s7VFFCb1CoTgtSOmIuKY5ve/HOmJI6QijEEx4gT/RGWPT9uO8uqeJpGmzak4Zn1g+k8XTp+BMZB07Tin0QojHgLVAs5TyXPfY+cB/An6ccYH3Syk3u499A7gXsIAvSylfHqW1KxSKCYJlSxKmRSxpEUtZdMZS2DmZmYks8vuaunlq61HeOdSKJgRXnlPFbRfMPC2594EykIj+h8DDwI9zjv0D8NdSypeEEGvcn68QQiwG1gFLgBnAq0KIBWpmrEJxZtPcHSdlTvy+9zStPQnqDrfz5u+a+fB4FyGvzicumMVNy2ZQFvKO9fL6MJCZsW8JIeb0PgxMcb+fChx3v78F2CilTAD1QogDwCrgnRFZrUKhmDDEUxamLWloi2Jacszy0wBIGyPWihmogCGmUaSU7D7RxXM7T/DOwTC2hIoiH/deWsu1S6oJesdvJnyoK/tj4GUhxD/h7K69xD0+E3g357yj7rE+CCHuA+4DmD179hCXoVAoxhNSOu2RkYTFcTcHn2bMJjlJm9oX7iLUtIVI9Qrqb3wcxOBMAXaf6OI//+cg9eEIIZ/Ozctm8vHF1dSUBsY8/z4Qhir0XwS+KqV8UghxJ/B94Bqg0Ccu+NeVUj4KPAqwcuXKCZyhUygUpmXTEUvh0TXaIkkiCXPc5N2NWCuhpi0IaRJq2uJE9sHKAT03mjT55YeN/OTdw1QU+XjwynlcvqASv+f0tkcOl6EK/T3AV9zvfw78t/v9UaAm57xZZNM6CoVikpGybKJJi5buBPGUhRDkFVnHA2aggkj1ikxEbwYqTvmc4x0xnt95nFf3NBNLWSyfXcKfXnsOxX7PaVjxyDNUoT8OXA68CVwF7HePPwusF0J8B6cYOx/YPMw1KhSKcULaGjidhmmLJGntSWZ+Hi9RfB5CUH/j4yfN0Vu2JNyT4HBrlJd3NfL+R23omuCy+RXcvHQG86uLx2DhI8dA2is3AFcAFUKIo8BDwOeBfxNCGEAcN9cupdwlhHgC2I3TdvmA6rhRKCYP8ZRNRyxJV8wxFktZ9vgU994Ird90zZNbj/LTdw9juherqQEPd15Y41gUjMMOmqEwkK6bu/t5aEU/538L+NZwFqVQKMYX6T74w61OB81EJ5o0+atnd7GnsRuAlWeVctHccqZP9bNw2hS8xuRycB+//UAKhWJMSUfrpm3zUTiK36NNeJE3LZv36tv477frCfckAPjkilncvWo2Hn1yiXsuSugVCkUG25aEIwlsG2wpSZo2kaSJbUMkMcpZ2BHode+PjmiSn757mJd3NwFQWezj9uWz+PTqyS3waZTQKxSKDO3RJE2dTqR7Wv1nRqDXvRAnOmP8x5sH2dbQgSZgxVmlfGx+JZcvqDy9G7hG8SI2EJTQKxRnMKZlk7RsPLrGiY44CTMbtZ/OIutwet17E0tavHMozOt7m9l5tBO/R+eulTWsqi1jQW73zOkS31G6iA0GJfQKxRlIyrLRhOBYR4xIwkLTIGWOXf59KL3uvWnujrOzoZNHf32IWMpi2hQ/d6+azbWLqykv8uWfnBHfOiLVK/sX35NdDAZ4oRjJi9hQUUKvUJxBxJIWlpQ0d8WRQNTNu1tD3eRUSOyGEikPoNe9v/dusafwT6/sY+exTgACHp2/umkJy2eX9GtPYMRaCDW+hwBCje9hxFowg9V9Xr/fSHwQUfpIXMSGixJ6heIMwLYl3QmTo+1R9+cReNFCYgdDT1OcpNe9N7Zt8U8/WM+vY7V4MUm6Uva5S+Zw1cIqSoKn6H+XvaPzvheEk0Xig4rSh3IRG2GU0CsUkxjLlvQkTDqi2U1OI0UhsQNGNU0hpeTD4138YvNBtsZqAbhBe5ebbvs01dNr+nlS3zsMM1BOtOICguGdTpRdYI15kXjVckA6hQshBh+lD+IiNhoooVcoJhlSSmIpi564iSUlHdHUqPS/9yd2gxLAAaZ5DrdG+MtndhHyGzS0RZnq1/lzYz2f0l+jWMTYO+UPMQu9Xu+7jjUbMGKt1LzxIMHwTqKVy9y7jgLvnY7Eoy3UvPEgCzdclHeXMtZR+mBQQq9QTBIs2+l7b48maYskAUbPZMwV1Po1GzHibXliN2ABHEBBtLUnwabtx9i03fFGTFo2X7qkmsvPncPCl/+ZUFOSSPUqJ1oukEoyYi2EmuoQ0iLUWMfc5+4g2LIDsBBAsGUHvrbfYQUrMAOVfdcrNBBa4buUMY7SB4MSeoViApOybLrjJinLRggyPfBpRqVFsmAhMkcg0wIobYxouF/B768gmjRtfnswzGt7m9ne0AHApWeXc22tj9v2f53ibVuIHF9B/Y0bMeLtmdc3ouF8QY62UPP6AyAtp/BcuYxgeCcC52eJhmUEmL/pegAi01YXvNiMRjHVYwjKQl7KTlVLGCGU0CsUE5RIwiSSMGnqSpz65BGkT24+2uKIe6ACkE7E6yth7vOfzOTACxZlexVAj3emeOytvfzmQBiA6ik+1i6dzrWLqrjqvT8k9HYd6Ug81LTFEfmciLq3ICMEoeatCECic+Tq71Hz5oPO45UX0LTiT5j70qczQzRCTXX5nyVzcZI0XPmwm5svEPUPgqBPpyLkY0rAOK0DS5TQKxQTCCkl8ZRNVzzl5N7HwPy9d5Gy5o0HswVLIQg11mEZfnQzkhXlAkVZM1hJpHoVNH7I80W38+dPfpR57KG1i1l+VimaEBjRFufCko7EhV44shaC+jUb8LUfIFG6AITIj8RDVXk597kvfRrLCKCbPQBEqlZkP8vJuogKzlfqHyGgJOihosg3ZgNLlNArFBOEhOkUWI93xNG00zDgo1Ch1DbxtR+g/oYNGIl2kBYLN17s5MCb6kAIBFZG5CUQrVhaON0hBN+Z+S88dvgwuDclt+q/4W9n1fHRWY9niqlImXdhabjq4cItm7bp5OBz7iLSwp5+LSPe5lyM3AuHbsbYf+svsYLO+hZuuGjEuogCXp2ykJeSgAetP7sF24ZoGELDu1M4FUroFYpxjmnZtEWTdEZTxFOOup8Oke8TyUqbRT9Zhp7qxvIUs+cz26h98VOZHDjY2HoRmhnBNkJoZoxoxVKOXPM9wCkWCwG7j3fxjac/yLxVVbGPmcUaP2v7FEFiyGbDvcCUZ9dQtZy9d7+D6a+g9sV1fQu40nYLrVvz7yIC5U6U3liH5QmimzEiVcvzIv1E+aKMsc9wu4g0DUqCTu494D1F9G7b8KO10PAe1KyGe553XmAUUEKvUIwz0umZWMoimjSxbDm8Hvgh7FTt3a1iRMPo8Tb0VDcC0FPdBE/U5eTAnYSGZvZgG0E0M0q0YhlS97Bww0UcLr2Ya5sfINmrzfOPPjaX6xZXseCldQSIISEjqEYsTKixDoFFqGkLoGHEwwULuEaslWDLjj53EUas2fkcWJm1h5q2sHfdO6DpmP4yjFi2YFyoY2ggXUQBr055yMvUk0XvvYmGHZG3TedrNAxFVQN77iCZ/P6cCsUEojOa4nBrlCNtUY61x2iPpAYm8tJ2UhS922zcyHzhhtXUvnCnkwoZwGvVvPZATqRuUfPGAyRK5mEZRaTfoXrLPxOtOA+JjqUHc8Q+isAmGN5JqnEPr1rLuOLE/RmR9+iC73/2Al66q5y1504jkGzPXDBAp+GK7zri6yvF8jiva3mCmP6yfne0moEKItNWItGJVi7n0E1PgbSY/av7Mp8ju3aLmje/hOkvpfbFdSxcv4q5z94GtpXtGCrURdRL5IWAsiIv86uLmFdVRGnIO3CRByddU7MaNMP5Ghq9Vs2BjBJ8DFgLNEspz805/iXgQZyRgS9IKf/MPf4N4F7AAr4spXx5NBauUEwWTMvG0DWOtkfpiKYG3xIpbWqfvzObZln7RCaHbcRas1Fx4/tOu2Po5FGjEWvtE6mHGuswYmEOrX2C+c/c5ET6Le8DGtGKZRy79G+Z/8zazGvskzX8BQ9Sl3B2q/pIcom2i3/1/ifH173K/CdXOSkgI8Sez+xw0yN1TkH0zS8RatpCtGIpeiri3kFE3X79cixPcSZ9ZAbKnTfsHY0jmfvsHQTD2zJdN4fW/pz5z6x16wlb8LUfyPxugi1bmfvc7Ry6+alTWjYIAcV+g8piH0HvMJIiQjjpmnGSo/8h8DDw4/QBIcSVwC3AUillQghR5R5fDKwDluAMB39VCLFAzY1VnPH0U3RLj+crCXhoj6T6Pm8AaRenK2Wzm5bYjBFtwQw5Bl2mvyyns8Sm5vX78y4Ehd4j01XTWJfzXIv5T16Nnoo6ue5UBLDdyH0b859ZS0rqPGNdytfML+St79vGo6zR36VIxIlUr0JPdGRTQGaEs5+9FdsTBCkRqSjB9j0IaRFs2YZtBNDMaCaiN+Jt6GYsI/6+9v0kys4BaTndNiXzMGJhkBbBlu15qZxE6TlEK5ZlirWJ0gVOb72b1w+Gd5602FrkNygNeij2e9CRzt/TM0yB1rRRS9fkMpCZsW8JIeb0OvxF4O+llAn3nGb3+C3ARvd4vRDiALAKeGfEVqxQTDQKFN26kha6EBxpc2awNqUK9MIP1CGxz27O7M9GPIxuRrO94s29OkesJGc/exuBtj1OR8vVj2AGsm2Is1+5LxMVZ/PzEWKliwi070ICMenl+9Ya/tm8M28ZW3x/RLnoztwVSKDhqocxA1VYRijTmRNo2+0sG5tg24fYRhGaGctr0dTNaGYHbuYi5Akyf9MaIpUX4G/fi57qRgodISWWEQTsTJpJah5qX7qbYMsOx/ZgzUbQNA7d9CRzn7s963nTT7E14NWZUx50et9PYxF1pBjqfccC4PeEEN8C4sCfSinfB2YC7+acd9Q91gchxH3AfQCzZ88e4jIUiglANIxseA9hm8iG9zhxvIFWStA1gWX3n6cZqEOiGagkMm11phPFDLjnpHPtOYLnPO6KmW06XTSZfvfNrp9LtpvF2UnqiLRtFKFZcaLl5xEMb+OgPYM/S93HVrkgbz2/9n6ZGi3s7j7tjQYC4qWLCLXU5XTrBNGsaCbHf+j69cz95aezEXn50syr1d/4OL62fczftMb53TRvwbm7AKTlXhh68lNPzVtA4KRpwjudi0awEjSdQzc/1e9dU8CrMzXgyd/gdBqLqCPFUIXeAEqBi4ALgSeEEHMpvJOg4L9kKeWjwKMAK1eunNgThxWKfjjWEcOnTyFQtYJg0xYiVStolVNB0L/Ip1Mp/vKBtfb10y3SO9e+/5aXSFQszjzuaz+Q1+8O5HXZmMFKItNWEmp8n2j5Eg6tfQpfVz0HogHuf6Yhbwlni2P81Pt3TBPt2FogU/PN7i2FaMUFmMFKp5smvD1PiDUriq0H0KwYYFO97TvuZ9/srLXzQJ6pWKLsnOzvpldEj8RNL0XzWirT/fN9fpe9PGuEgKkBD+VF3sI5+HQRNR3Rj2IRdaQYqtAfBZ6SUkpgsxDCBirc47leobOA48NbokIxMbFsSWc05Qj6oI2+8t0W+x3qkfNzRqxyLxSVFxBqfh+AuS/eyZ7P7gDh/G+fKF2AZRShmz1Ywk+84lxCLY7NQM2rX6BpxVdp+Ni/MPuNBwm27KDhR5/j/vj9ecv9qvFzbtLeobq0gsaLH2bv1LOZ/doXc4qgDtHy85xOmHQNoGo5oaZ0DSDiin0ChO5cbJq3sv+W55n/zI0ImdMa2ViHr20fibJz8i9uuTn6RHsmn5/5mmvP0M/fwOfRKAl4KA15Tz4w/DQWUUeKoQr9JuAq4E0hxALAC4SBZ4H1Qojv4BRj5wObR2KhCsVEIJ6yaI0kMS2bnoSZ3dg0QKfD/C6Zuoylbp7FQGaH6HepeePL/W/Zr1pO0/I/zqRA9FS3I4ZlCzAiTeixMImSeQTD29FlHH/7HmepQKiljrm//DQAL1kX8k/mtzkkZwBwvbaZz95xB4vfup+SDje109nI3Jd/P1PszKRcSpeA4SUY/oDal9Zl1ogQIATx0kUE2vegmREsI0i8fDGh5q1OsbRsIZHqldnNTm6EPn/Tmsznzf5ONaxAudMb7x7r/RVEn79Bkd8gkjCprQgR8g1CDk9TEXWkGEh75QbgCqBCCHEUeAh4DHhMCPEhkATucaP7XUKIJ4DdOG2XD6iOG8WZQMK0EAg+ao0Ma/aq6S9zRa3biXZj+Y6MCDe90rSZhRsuJp2bzm7Zt7MbnZo2M/eXn3HTGRaWp5hEyTxqn/tkJsqHbHolN41zVFbwb+bt/MK6PHPeddpmvmQ8zRLtMHuLPkPLFd+mZNMN2deRVqbYmf565JrvuRYJhWwFLEIt2zKfQTejNFz5XfREp+NVo+V4vvvL8LXvz+blc2sW6bugxjqilcs4dNOToJ18V6qmwYypAXweDSl9gxP5CchAum7u7uehz/Rz/reAbw1nUQrFRKEnYRLuTtCTMNHEyYurfSjQOum0D0YzBcX5m9Y4dgJWjEjVioxpmOPiaOebfPnLqH3+zpyNTiDcqUiHrt9AsnQeRqKNUPOWvGJa+lzLU8yBRAnfNT/BC/bqzOM3+Hfx5/Z/c5bW5Nr9Lsf0l4MtMz3tzsVEEi1fzKE1j+Prqu9rLOZOaTL9FQXz65YRpOaNL2Ui+nRBOB2F5+Xlc/LsuXdBp+qH1zVBRZFjDVwaOj0WweOByX0ZUyhGiXjKoqU7QUc02/tuDWanUz+tk07fuxPRgyPUmtlDtPy8jO+7EQ1T88YDTpG0bDFHrn0MM1SNEW0m1Px+Nm1SeQHB8AdEqpZTvf1fM8JqaX50O5pZSrT8fPZe+T2+/nIjB7ojmeNrtXf4q/LX8AdChJqbMgXVI1f/h+s3s4VIxfk0rfgqyeKzWPD09QRbP2DJjxcicNwgG65+xBlOEgs7U5rWr3ai7hsfx9dZj+UvdSJ+QDdjzkXI3dDUp8uon6KzGag4ZT/8lIBBeZEPQxNj5iA5lgg5KpMJBsfKlStlXV3dWC9DcSYySPfAnoRJU1ecaKJXRnKQfjJGtIWFG1YjpIkUBnvvfs/pSIm2sHD9qowlL6R70PV8f5ZIIwue/Dia2ePsLv3sToxEOwvXr8wI/d67N4PQAZlxZZTogE2P9NMip7LBuor/stbmrW295/9yib7bfX+NdGpFohGtPD9vQlNv9cjtpJEAQidSvZKGK/8fCzdcnLkLSbdq9q475HbG1N/4xMALnbaV1w+ffq6uCUqCHmaUBAb2OhMMIcQWKeXKU52nInrFmcsgNr7EkhbHOmLEkgVKTv1tbDqJ+Pc3tSjTkeLm0G09hGZFyNuZqvkyEXnu7tKDtzyb308frHbe17aIViwl2LKdaNki3k6dwx+13J63nmKi/IfnX7lE350RY+f1s+mhaPlSguEdOROayKwB92eZ8z1k2zVnv/IF0v38afMzp799q3MBc+9mjFh4aAM+evXDh/wGPo9OScAz6fPvA0H9BhRnLgU2vpiBCgxdI56yaOyMU17kpStu0taTzH9upoXRKRL22diUa7FbaFdrP2kIkAi3f8FpOYxlWg4zO1PtaCZiTgtnoG03Rrwtr3hpRJsdY6/X7ifYso0uGeIHzQv4VzMr8rdqb/Ntz6N4hdknGs9N/TRc9XCe0ZmlB9GteN5FIVqyhGDnXoS0stF/eKdzkQnv6HMHkKktBKsAOawBH0Gfjkfz0K1XMqcsSLHfM+DnngkooVecuRTY+FLf3INH1+hJmEgJ3fECzpE5XR59NubkWuymxT+9Aam3mViB2ap97XaXIQ1vznt157+Ee16sbHHWq8ZfytxnP0EwvI0OGeIfzdv5sfW/sHPMaitpZ5Pv/zBDtOW9Ti5Hrv5eJlXka9+f4zAJuhVzLgQt27E1P5odJdixyynOur+HzOBwXymLfross/a0k+T+T76euevoM+91EAM+SoIeasqCgLN3QR+Mg+QZghJ6xZmLu/HF7mlBK66iI5YiYdqZ4R4ZeqVgMtYEOR7neirK/tteIlE637HY9Ze7m4I2k7b5LTw3tfcGqY3ubtR0q+BTIMjMYfW17mPec2vBjfrTQ64P3rTJ2TTUuo+Zb/8FsmU3P7Cu52/M38+81WzRxD94vsdF2t5e+X83/16xDKl5CDWnR+9Vk4m08wzOIDJtFfVrNjotj0/fkL0ApCLsv+2XjtGYEJm6g56K5t0x6FbcqR/kFFQHO4C7yG9QFvIyxZ+VMSXyhVFCrzijsREcThZhhyN9C6xQMP/e21grHcEmSs7OKwg2XPVdFm68pHAXibQxIk142/fnDPh4HyPSXHAQdTYVVEek/HzazrmDmt98E4ENCIx4mNpfXEs0ZbEs8e+YOf9re0nytPchlmiH89I96e8j1atouPrfMf3lfXLkmUi719i99OPZlkdnX2Rk2oUZkU+T9ooPNb6PFMLp6TcCjr98mn5TWWS99oXAmFLNrLIQzV1RKrVOgv7pE2Jn6lijhF5xRtIRTZIwbcI9CWf3aj+F0/6MxfJy4e5W+3S6JL2BCZyOk7wo1RWtmle/6NoNOKTNvc554hKnNbFyOccv/VsSpQtB05yJT+nJSuEthMJbMs+zPEHCrS3c0PMfeZ/xXuOX3FDexIrOV7CMANKkT24/6yhZUTBH3jvSzozd67XytC99/ZqNBd00M0ZkbvSf8ZfPa58ssHs447W/2Vnx7IsR9zzHWc/dNaHcI8caJfSKMwbbnVna0BajM5bj/X6SwR1ZoasrnFJwxcmINOX5u0QrzsMMlOdH55k0yPtk3BbTS4A898VQ8/vMf/r6zGxWPRLOf1ucYVI75Nn8Ino5P32uK/PYElHPL7x/jV8koTPdmRMjVn4egdZdfd4bCb62fYVz5KeItH1t+9zcvU2w9YO+4p3ze+pvw1N/CAHBVHvOBi8JR9+D8L4J5x451iihV0x6kqZNU1cc05ZIKYn0StGcbHAH7s5SJwyWmZ/zttyv/Tl6LF+Ij1yVs6mo8gKOX/q3WP5yNw2S45Oe85x0rhycGkHan+bsZ24l0PZh5pyk1FlvXcNfm/dknnvtwgpubvkvbuj5BZqQmedn0jPTLnRy6m37mL/pusxrRcvPdyY69UpDnczd0Xlyr4J0oef1uUuSrscN2Z8Rfc4PlE5jVlkQQ9dATkHOWg0Nv3XOrbkIKhdOOPfIsUYJvWLSYtmSlu4E4Z7EycfznWxwR8bq1y9fTDUAACAASURBVHFUzHi15Gy5X/yj89Dcvvb028x+9fMEW3chsLPRuVHkFmi3ZAqb0aK5BCOH3ZZEQaxyGcGWbXnLCbTtQgANdgW75By+kPpfeY+/57ufnkXfY+7hpxFCIoH41Hn4Ow9kNls1XPkIaDqJ8oVEqlcRaqojWrGMIx9/1N2ZamULyr1y7BlyhDu/IF3geQVqG9k0WOGaxdwX7yLYuAV71ir0z72Q+VuIz70AkWZAOJH7BHSPHGuU0CsmD+4uVxmsoDWaorkrke8905/Fr7+i8OAOenWDVF7gOj7mj6DT7OwEp3SMGmz9ANDyult0s4fjq/8SPRUhWrGUec9/gmD77xzLAzPm9pt/kB+Jly/noFnG15o/zj6ZdQA/WxzjV94/QxMSyyiiedoqt0D8PrYewN95wPV3T7i96u5nEoL6tU/kzFYlc/GJVC3HCpT1+7vq0x2Um7vvdXEoWNs4SWfNVNlJ0D1fP7Y5Px2jaVA8Lf9vPcHcI8caJfSKiUuOfUFrT5zghlvxn6gjNm0FJ9b0amUs4PNe++Ld2Z9v+Bm+jkMZI670c4xYq9MPHm1m/pPXMP/p65y8+ae3MveFu9zo25FzCdhaAM2OuWJt54m2pYeY+9Ldeb3w6QsAQkfqnoxYSyHYZ03nR80XsiGVdZC8Jrifh8pfZ1bL21hGCKwo8fLFoAnq12xg7qabCbZ9mN1slf0FUHgDksykUfzte53hHn2sCB7vK9y5m7MGuvO3QL4/6NOpnuJHJ4RQ6ZhRQwm9YmLi2hfIhvdITFtJ8xXfZeGJOoQ0CTT23XDTW6h87QfyNjTNff7OHJ+UQr7uX83pme/G11nPoZufwtf6u0zOGyBedg7B8Ha3G8bdPFR5Accv+Rv0eAdzf3l3Xz8YyKQz3rrmWf6+s5M97dmPOqtY57aVc1hTa2AGLqWDe4i17mX+puuzdYVYC6ARdNM8kJX1vDmxvS54DVc9nKkbZIZ75NkhnyQaF3393TMUEHUhAKGhTakmqAuKfQaVxb7siD6Vjhk1lNArJiRtLccodeew+k7UgRQn7ejo3T2TKF2QOd9JmezsxzPddHzdX3LcutO+L4mpczFirU4niZv2yR+6obP/9tcyO0trX1znRup63man3Ih/g7yGv3i+M2/dr/v+hMhnfgO6QXaPrsDq9fn0aBuJsnOIVpzfZ7pT7pzY3he8zO8tt6haaOzeybpv+iOniOv3aFQV+zFtm9KgN/N0kfs6Kh0zaiihV0woEqbF8Y44PclivFU5wp7b215QiHp1zwhnyLQRbQFwpjg11hEtPw+k7QzcdgUw16kx3QI594WcO4A1G5zdn/F2ZrzzUP7OUiHcrh4nakYKDl2/kWj1Cs5+9iZ2tsK65F/mrfSvjR+xzngDH0nH2TLR3idyNoOVmaKq5Qkx/5kbiVSv4NBNv2Duc3cQDO8gWr6UI9f+l+Ml098O1EJ7AvobuzfAKVlphACvoeE3dKaX+E8+nk8xqiibYsW4JpIwCXp1hBDEkhaHwj3ZDU7ubknTX54VqH4izTxbYHT23/oCVqDMGXThFlqFlSTYusN532mr3XmtbRnvd0fos/YIAhwhvus3zH/yGnd4RhH773gFNE92Z6uU1D7/yZzdo6tZv+jf+buXfpe3xnX66/yp8TglHpk3Uq9fu163jz09dSljdxwozym2FhDsdKHVV4Kv42BmmtNIUuQ38Hs0Kop8aK5dsGLkGTGbYiHEY8BaoFlKeW6vx/4U+EegUkoZdo99A7gXsIAvSylfHsL6FQq64yka2mIEvTqmLYmnLLeV3c7Ln/cuHBaaLJRnW2AEmL/p+sxjzgalLeRuJAo11WHE2zFDVTm7Oq/Pt+R13Rf1REc2f2/2cNYr9xFo25O3noarH2Hu+kv536k/YPvheez/KCvyP/D/K1e4o5UFIE2dhqsecSLo/i5erljnpqD65M77s08WGqa/lEU/cYzGLE+xMzRcG/4NfmnIQ1nIS8Cj56dlFGPKQP6yPwQeBn6ce1AIUQN8HDiSc2wxsA5YgjMc/FUhxAI1N1YxGJKmTTRpcrwjjmXLPg6S+XnmOqfQ199UojQFtuFDVrAdr/VtmWPRimXZPH9mV+fKbFSe9ocJVILMjtSzjBCBtj2ZIq+vbR/x0gW8cVTwf+T/pc6ahYHF8poSVswp4+ZlM8C+iP1t+5nxzl9m0z456ZY+FOggMuLtA7ZvAPC1H8gvLrcfIFG+cEh/LyEg4NUpD3kp9ntU9D4OGcjM2LeEEHMKPPQvwJ8Bz+QcuwXYKKVMAPVCiAPAKuCd4S9VMZlJmBbtPQm6204Q95SfVOSQMscZ0sYyigvvzOxN3jb8XoLtr6D2RceeIFp5vusambOGdP95Ol2UOxhDCPZ8docjliXzqX1pHaHGOl7iEtY//gZv286uWZ8xm3VFW/hj8wcU6WdTv9Tp7kEzSFQson7twIqdfVsd++bw4eSOkInSBdmLk6fYSd8MAiFgasBD0KujCUGR3xhcDn6Qk70Uw2NI92pCiJuBY1LKHb1uz2YC7+b8fNQ9plD0QUrpDPWIJOmJJal94S6qT5aCyYlko+XnkR3IcYodnbmcRLALFnPzNg5pOdYIvdCMTET81iU/4L9f+5CtTc6dyMe0Hcw870rWLStn6c8/g8BENnUUmIk6sGLngC19T9Ypo2nZi9MgcvRFfoOKIi+GpuEzNLShRO+DmOylGBkGLfRCiCDwTeDaQg8XOFaw2iuEuA+4D2D27NmDXYZighNPWXzUGiFlOv88+k0zFNp2L01nclHlskzny4BEPk1/gt1baPvLcffDzqMdPL3tGHWH20nr3y7f5wiKBHuXf9LtlDmFQA9k9uxgWh1PdvHIuTidDEMXWLZEE4LpU/3DH65dYLKXaqscXYYS0Z8N1ALpaH4WsFUIsQongq/JOXcWcLzQi0gpHwUeBafrZgjrUExA4imLpGXT0BZ1umdc+kSp/jKMnhPMfu2Ljpi74+xyz6m/YX22a2QUbv9PluNOc6w9xvrNh3lrf9bU7Lol07h56XQu/+09BJtMItWrTu0ECYO7sAyy1XGoFPkNpk/1o2uO2A9b5KHgZC/F6DJooZdSfgBkLr9CiI+AlVLKsBDiWWC9EOI7OMXY+eC2EyjOWEzLpj2aQgho6U5gWgWu67ki6C9zBS/biRJq2szCjRcTqVrB3rvfcXPq6wYcbRfkFNHzyVIk4Z4EGzYf4dU9TaTtdD65YhY3LZ1BacgL0nY7Z3rn8/sX6IFcWE4XXkMj4NEpCXky4j4SGg8oU7IxYCDtlRuAK4AKIcRR4CEp5fcLnSul3CWEeALYDZjAA6rj5swmlrQ43hkrPL2pN2lv98wGI4c8m4DmrYCGEW8bnigOJHouEIF3xlL8YksDL3xwAilh7dIZfHLFLEqC3pO/9gAGXQ9lnN5IomsCTQNNCOZWhByb4NFC7YI9rQyk6+buUzw+p9fP3wK+NbxlKSY63fEU4Z4kPYWGa58CM1DhdNU0vw+ApQeJV5yb2UCUcV0chigOOHoWGslABduPdPDImwdo7k6gCbhqYRV3Xzibqin+ob92n/cags3ACOA1NLyGRnmRF4+mYUk5uiKvOO0oCwTFiNGTMImnLDqiSWJJ+9RP6A8haLj6YRZuuNgx27LimUlNuS2NwxHFgUTP8ZTFO4da+c4r+zLHzq8p4b7fm0tNWXBYr90vpyn3rmkQ8Oj4PTpVxT50TagNTpMYZYGgGBFMy2ZfU4/j/z6QzpGTkTcn1CFSvWpAO2AH+z6F1tnYFee5Hcd5ZXcTsZSTcrp9+Uw+vfqsgfeKD/d3MIpoGlRP8VMS8KjIfYIzYhYICkUhbFsS7kkQ8hm09iTpiqf62hPkCnKuN03g5AW47FSnrJFYQevc4Ua+vaLn9+pb+eWHjWw90o4QgsvmVXDFOZVcUFM6+N2epykyHwjpX7Xf4+xe9Xt0At6RqqwqJgJK6BUDJmXZRJMWsaRFezTpdM/IGEasFelGroUnC5XnReiRaatPGpHn+dL0ts5trCNasRTTXz5in2nX8S62N3Tw5NajeHTBJy6Yxdql0ykv8o3Ie4wVuiYwdMEUvwevoVEW8p76SYpJiRJ6xYDovcEpHaHXvPFgXvReKD9txMJ5XTShprqTR+S9Wi0zzpTSci14d1L74l3DSt9YtuQn7x7mya1HM8cunlvO1647Z1LY6U4NeJhe4idl2fgMHWU/c2ajhF7RL/GUE7nbEtojyeyA7XR6JserPTed0rtImhV/J6KPViw7dUSek/pIfzVi7X0GhAw2PWJLyW8OhPnhbz+iuTtBTWmAKxdWcfHccmaWBCZ8QVLXBLPLgwQ9OpomJsVFSzF8lNArCtIZTdHQHiWvVu8WGEG6EbqVZ9frRN2FipDS2TwkLWa/dr8Tkb9wJw1XP4IZOIlLYy+G082yt7GLlz5o5HdN3RzriFFTFuTPrjuHy+ZVTHhxB/B5NGaUBEiZNkU+g/HQZKEYP6iuG0Ue6Si+tSfZR+QLesBXLc+x682mVgrNXs0d2ScBhE6keuXgUjCD6GaRUrLtSAc/39LAh8e7AJhZEuDuVbO5bF7FhLfTDXh1/B6NkqAXn6Gp6P0MRHXdKAaMbUs6YilCPp3mrgSdsVSfc/KKrM1b2bvunfzBGNJ2RL5la14qB7KzV4MtOxwjspYdgDW0DppTdbNIGxEJ8+sTgl9sPcrBlgjlIS/3XlbL5QsqmTIJ/NKFgCKfwbSRMBhTnBEooT8D6UmYRJMmSdMmkrBIWbYzRtWZeteXtAe82wkTrVjqRPA51rJGrJVgy45MS6RzTgXYJrGyRdmJS2s2YsRanfF8I7jVX0rJbw+GOfHbDfxP9ww+ktOYWeLny1fN44pzqiZ8tDs14Ml8P6PEj1Dj+RSDQAn9GUY0aXK4NeLMye4l6v2JfCZlU3E+sfLFBbtezEAFkWkrnQtB5TJncIe0WPTT8zOTl+pvWA+anhnPNxIbilKWzc6jnWzafoztDR34WcJicZiveX/OWbd+FznB/VSK/c7/orPLg9i2dDcHK4FXDA4l9GcAti3RNEHCtDjSyx74VOSlbFqcOkomNRNtyUvf9BZvX+u+nFmqEcdSOO1/PswNRUfbo7x7qI3H644QT9kEvTpf+Nhc7j/8FYrdcXz1E9D+Nn1XpWtZ7/f05qYhDflQKFBCP6mxbUlXPEV33ERKCubeT4XpL3OKqOm8Om5qpnQhNa8/4M44zRZUc8W74Li6YVoDxFMWGzYfYdP2Y9jS8Wv5y7WLWTarBK+hcfi8028KNpJUTfGhCYEmhGN3rFCMAEroJylSSnaf6ELaBYR1oGIrbWpfXOcUUSuWInUPoaY6bCNIsG03YDvRfeN7GLEWzGCvqU29x9UJBjWxKZeeuMkTWxp4bsdxTFty3ZJprJpTxvSp/nyDsXFkPTBYppf4mRrwTPh6gmL8oYR+kpCybDy6hm1LWnoS9CRMpF3IF72A2EJW+JEZTxqkTaipDoFFMLyTveveQU90MH/TGgR2/oxI2c8FI2dcXcZnfoAbnqSU7G3s5pcfNvL2gTBJy+aqhVVct2Qai6dPGYHf2tjj92iUF/nojqcoD3lV/l0xKiihn+BYtkTg7FwN+gxOdMSIp5wkfCHfGSD/WK6NQdVy93FnB6tlFIG0XEG3qHnzS9Sv2Zj1oTECjg9N9YoRHWotpaTucDubth9j59FOAh6dqxdVccO506mtCA33VzbmBLw6PkPD79GpKHLEvTToUSKvGDWU0E9gpJQ0tEUxbYllS5q7E3mdM/0Ja+ZY1XL0WNiJ2qVFqOl9QGQ8aXSzp4+DpBFvK+xDg8SIhoc11FpKyZYj7fzs3SMcaOkB4A8umcOac6dPaLdFjyFImRJdExT7jYJe9krkFaPJKXfGCiEeA9YCzVLKc91j/wjcBCSBg8DnpJQd7mPfAO4FLODLUsqXT7UItTN24KQsm7ZIkq5YCtOWheev5lIoH582JHv9gTzP9944hdRo1kGyegX1Nz4xvKHWhd7Hdnrgn952jP3NPVQV+7j7wtl8bEElXmNi56s9hqCmNEhbJMnMkoDqnFGMKCO5M/aHwMPAj3OOvQJ8Q0ppCiG+DXwD+LoQYjGwDliCMxz8VSHEAjU3dvh0xVNEExbhngQeXSNpDrBHslBxUmggtD5zWbPfC/bf8iJWoBQ90UmiZD5Gor3faH2oo/PiKYtXdjfxzI5jNHUlmDHVzwNXzOPqRRN/g9Os0gCmLSkJOsVVv2syplCMBQOZGfuWEGJOr2O/yvnxXeAO9/tbgI1SygRQL4Q4AKwC3hmR1Z5hSOkUVlOWpK0nmTk+YJE/Cb0dJS2jCN100iWR6guZ8e5DWf/46lU0XP3IAF5rYDtd2yNJnv/gBC9+cIKehMmi6VO497K5rK4tQ5vgKQwhwND7tkaqXayKsWQkcvR/CDzufj8TR/jTHHWPKQaJbUuOtEXpHsJw7X7plcapX/uE02GDpOb1Bwk11RGtWEbDVel5rQ6hps0s3HBR/wZkA5jfWh/u4TcHW3l7f5gTnTGkhIvmlvOJC2aycBJ00AS8GtVT/CRMm5BXlb4U44th/YsUQnwTMIGfpQ8VOK1gElkIcR9wH8Ds2bOHs4xJQ9K0aeyMM73Ez+HWyMAHbA+kL76fPLoZqnbaHpu3IrAJhneCFHnRPgxghF/vFJG7ppi3jMfrjvLzLQ1ICUtmTOGSs8u5ZlE1M0oCg/jtjE9mlPgp8hvoQmDoGkVSqsKqYtwxZKEXQtyDU6S9WmYrukeBmpzTZgHHCz1fSvko8Cg4xdihrmOykLJsDjQ7w7UjSfPURdY0AyyEniyPbgYqiFQtd4XdoubNB6m/cWOmHTN3itSADMikjXz6C/yg6SyesC4H4OqFVfzhpbVMyTHnmuhUFvv6jBtUIq8YjwxJ6IUQ1wNfBy6XUkZzHnoWWC+E+A5OMXY+0H9bxxlM0rQxNEEsZdHSnchL0QxY5Bl4IbRPHt1fhhFtydwFNFz9CAs3XJSN3OPtmCFnp+tgDMiOtEV5fWc9zx7/FF5S1IoTfOaqFaxaNHfAn2k8Ul7kpS2SpKYsSLFPpWYUE4tT/osVQmwArgAqhBBHgYdwumx8wCtuBPOulPILUspdQogngN04KZ0HVMdNPpGEScK0iaUsp0VyEKJeiAEXQoWgfs0Gx46gZB61L67rNeu1ikj1ysKvcwpbgfQO1me2H+M3B527gFuL9vFN8xEC0xZQv/COfp873vEaGpqA6il+ykJe5f+umJCoCVOnkXjK4kBzT2E74OGQm6NHFo6+c1I8eZOehMHeu99zhHyQ05s+PN7FwZYeXt3dxOG2KAGPzi3nz+D6JdMoD3kmtLkYQJHfoKY0QMqSE3rDlmLyoiZMjRO64im8ukY06YzoG5Xrajri7jXur+HqRzD9FRixMHq8LWfS03aiFcuckX9Vy7PRu9AwA+UYsZPvcG2LJNn4/hFe+rARgHlVRTx45Tw+Nr8yTxAnqrlYwKszfaqfkJuiMZTGKyY4SuhHiZ6ESVcsRWtO//uokBOF5+frN7Nw/WosI5jpj8/2ytv4Og84PVIix+TgJIXdhGmx/r0jfHCsk8NtUVKmTW1FiAevnMeC6uLR/YynCa+hURL0UFXsU0VVxaRCCf0IkXaPlFLSFklyojM+OtF7Lr2F+Yb1jtFYqhsAgZ3xqwFn+AdoznF3IEhuAbdQYfe4WczbB8L89L3DpNx6whULKvnU6tlMnzrx2yPTBLw6M0r8BFUPvGISov5VjwC2LTnQ3INH10hZ9rALrHmcJG9uxFpyDMm24Os4iJ6KZozInP80wOnHj1StxN++Fz3VjRROPiK38Jpb2H2u+JP84JUmtjX8LvN+f33zEpbPLh25zzaGCOGM59OFIGXZTA0o90jF5EUJ/RCwbEk8ZZE0bVKWTafbPWNaI9xgdLIeeWlT89oDGRvhSOUFWP5Sd27rZhzHGgCbQ9f/jGTpAhDCaaEEkLD/tpdIlJ2TvYAIwQvL/4uXdxzh9YNdlAQjXLOoilW15SyePiVvQPVERNcEfo/jOzMl4KFItUkqzhDUv/Qh0DDS1gT9kJdKaaxzbIBDVfmP4UTugbbdLNx4MZHKC4iWLSbYtivzOnNf/n3HviDtJe9eONIib9mSLYfbeK++jV/tbiLg0bnxvOnce1nthDcXA/B5NKqKfSpqV5yxKKEfIKZlc7wjPrhdq8N9z947Vt94IGtdEKggWrmMYMtWBKCZESfn3ryF9Ii/9Cqzm6Da8jY/HeuI8159K6/saeJoewyAy+ZV8JWr50+KfnGv4Qh8iRrqoTjDUUI/QE50xgc3XHuYQ7CBwjtWoy0ZoT9005PMfe52gi07Mp7x0fLzCIa3ZV4iWnkBwfAH2Vy8EHzY6eOnv/qQD493YdmSmrIgX7lqPufOmsq0Kf6hrXUcoQReochHCX0/SClJWZKOaJKOWIpEahDWwMMcxJFL3o7VquV5vjP1Nz7OoZufyp/25C+n9oU7CTXVZdwm01Og9jX38G+v7edIm1OwXTmnlAeumNfHr2WiIgRUFfuoVO2RCkUeSuhdbFuiaYKEadETN2nqSmDZQ0vRDHUQR0FyLIBButF9/utmzMncTVPuEwHoSVg8sa2Hpq4wvz3Yit+jsWTGFB5au2RS7PYsCXrweTQMTWOK38CYBDUFhWKkUUKP00VTH45g2o7RWDxlD6sHfrCDOE5JZuerdF+3rt/XNWKthJq3kpSC1457+eZP6uhx7Y7XXVjDbRfMnBS94romqJrioyTgUeKuUJyCif9//DCQUhJLWUQSFrGk0xqZKmyfPzgGMIhj0LhzXpHSqbLaNkasGTNQlff6TVYxr/g/zy86FnBEVlNZ5GHFnCnce2ltn6lHE41iv4FpOxfhmrLgpCgYKxSngzNS6CMJk56ESdK0SVp2RuRHlFM4Pvahv+JtepD3Gw8SaqwDLLe75v3M1Kd912/gpV3N/PZgmP3NPSTN36NmqpeHfm8ey88qnfDj+QCml/ipmCS1BIXidHPGCL2Ukp6ESUc0RWcsNfr2BIOhv+Jt+niOwOftdrUtfnNc8tDPttDY7XQEXXp2OZfNr2R1bdmE7oH3ezQkoAkoC/kom+B3IwrFWHJGCL1lSw63RogkxoE1foHIvb/ibXZTlJUReMsIIlIxnhQfZ2PqMrba85ld5OF/3ziPGVMD1JQFx/bzjQBBn86c8pCzP0AN1VYohs2kFXopJV1xk4TpDPgY8PzVUV1UP3Nb+yne5h2vWs6b5zzE66+9wFPWZQAUewVfuKiW68+djj5JBLHYbzC7LKgEXqEYQSad0MdTFm2RJB3R1JDbI0eL3pYGvrZ9GRuCgsVbITi4ZiPv7z3MEx92se9XnYAj8l+Y+h6X3/Fliia4/4wQEPTqFPkMivwGAY+ueuAVihFm0gi9adkcbY+dFg+aoZKJ0BvrsDxB5j99A9HKZRy66UnQ9OyUp2iYlL+c/9kf5j/+5yDRpMUUn87nL6vl6oWVTLU7MQOXTtjJTeAsPT1/VUXvCsXocspqnRDiMSFEsxDiw5xjZUKIV4QQ+92vpTmPfUMIcUAI8TshxHWjtfBcIgmT/c0941rkgUzkvv+2l1w7YYtgy1bmPvsJjGgT2BahTZ8j8dO7+JvvP8E/v7IPn6Hxv0p/zVbxGb5y9I8J+QzngjABRV6I7LJryoJMDXiUyCsUp4GBRPQ/BB4Gfpxz7M+B16SUfy+E+HP3568LIRYD64AlwAzgVSHEgtEcEN4ZS9HQFh1fXTQnQ2gkys7JMyQLhrcxdf0N/IP2eTZFPw9ACd186eIKrjunhCWPP4rAyt9lOxJeOqcRjyGYUx7Co2v0JMwJb3msUEwkThnRSynfAtp6Hb4F+JH7/Y+AW3OOb5RSJqSU9cABYNUIrbUgSXN4u1jHBCE4dNOTRCqWs9VewEOpe7gk/v/YFD0PgBXid7zu+xNumm0x540Hs57z6UKtW9RduGE1tS/cmWN7MP4IeHVKgh7OKgvh9+jomlAir1CcZoaao6+WUp4AkFKeEEJUucdnAu/mnHfUPdYHIcR9wH0As2fPHuIyJiaWLXl+5wnein+TfckEPmGyQBzl3rKd3J58Bj3VjeUppilYQah5q9s/r9Nw5SMgBEY0PHJeOqNEwKszqzSgdq8qFOOAkS7GFsohFIy3pZSPAo8CrFy5cqLF5EPCsiVf/OkWTnTFM8e+UvIbLrnjK0yxuzEDd7FH/g2+9gMk3IlQeW2XafOyQu2Y4yCVo2uCYr9BwKtTHvKq7hmFYpwwVKFvEkJMd6P56UCze/woUJNz3izg+HAWOBmwbMkLH5zgv359KHPs41odXzCeY3l8P3vtP8AMVjsPCINE+ULne2nTcOXDIARmoDKv7TKvHRM5YrbIg0EI0ITAawjKQ8r/XaEYrwxV6J8F7gH+3v36TM7x9UKI7+AUY+cDm4e7yImILSWv7G7i3UOt7GvqpsvtCPoj/Xm+bmxACJmdAlWoyFBoc1XuDVOOl85YpXKKfAY1ZcFJs1lLoZisnFLohRAbgCuACiHEUeAhHIF/QghxL3AE+CSAlHKXEOIJYDdgAg+MZsfNeOXlXY08/MYBAGaWBFh5Vhkr55Ryydxy5j//T4gW2SvH1Tf6Hoyn/YjbIp8CjyGYNsXPFL9qj1QoJgKnFHop5d39PHR1P+d/C/jWcBY10bBsyYHmHsI9CX65q5HtDR0A3L58FvdcfFZeOuPQTU8y99nbCIa3AxCpXlVQwAcl3qNhi9z3LZga8OAzNMqLfCqKVygmEJNmZ+xY0dgZ519e3cfuE12AM8puzXnT+cQFM6kuNH9VCKTuBQTRsiXUr9kISIxIs5OL95dnRv8NSrwHa4s8QCqKvbRFkiyaNkVF7wrFBEUJ/RAI9yR4e3+Yt/a3sL+5B4Dza0q4dnE1GDz6BwAAETNJREFUl5xdkR/t2ma2i0bTMGIthJrqEEiCbR8y97nbkbqHUJNTyrA8xehmLJOXH6u2yYBXo6LIR0nQS2nQq0ReoZjAKKEfBFJKthxp5+9e3EPKkpxdGeJzl8zhsnkVVBWK3m2TRT9ZlumL3/OZbdS89gBgI8HdFbsTyObs9VS3M1hkDPrje9sTpFG98ArFxEYJ/QB56cMTPFHXQLgnyfSpfh5au4SZpYH+n2CbFH/0Ska49VQ3xYdfc6N5d4CI0IlUrQAoGNGPdlE1jaELfIbGjBLn8yhhVygmF0roT8Gu4508+tYhDoUjLJkxhc+sPovVc8sp8vXzq5M2RqSJ+U9eg57qdg7hiPpZr38BkJkhIr+767eYoWmAdObB9srRj+bGJ10T+D0aAa/udM8IoQReoZikKKEvQGtPgu+9dYjWSIJ9TU4O/uOLqnngynkn7zaxTeY+dwfBlm2k0zESaLjs29T85psIV+SdV7DRE52YRdMBgRmqzrzMaKZrNA3mlIdcoVfCrlCcCSihz+HlXY1sfL+BcE8ic+zey2q5fsm0U4uitF2R35rdCOUSmXklkarlmfSMxEnRJEoXjPhnOBlCQG1FiKBX/dkVijMJ9X88UHe4jf/8n4M0dTkCf/OyGZw7cyrLZ5fgMwYW9RqxVoItO/JEXgKWHsIMVdJw5b+xcOPFmcf33/ayE16PMh5DEPIaTAl41JAPheIM5YwVeiklDz27i23u5iZNwPyqIm49fyYfWzDI1Im0QUoi01YSaqwjWr6YYOtuBBa6FceIt6HHu/KectZrf8TBW58bFU+atAeNpsHZlUV49NG/oCgUivHLGSX0tpRsOdzO1sPtbDnSzolOx0XypqXT+f2L5wwtZ53rSVO1nL2fehfTX0nti3fl7Wo1/eVYRgjdjCCAQNueUWmfLA15CHoNiv0GpiWVyCsUijND6JOmzSt7mnh2+zGOd8bxGRpLZkzl1vNncsU5lcPKWWc2QEmLUPNWQANBX9dJIdjz2Z2c/eytBNr2jFj7pKELLFsytzKELcnrBlK1VoVCAZNc6JOmzSu7G/n5lqO0RpIsqC7ia6vP4ZKzyzFGItKVtrMBKj0Bqmo5pr+M2ufvJNS8hUjVSurXuq6T0sZIdHDwlmdHtH1yZmmAeMpSBVaFQtEvk04dpJQcbInw+t4m3tofpjOWYvH0KXz14wtYOnPq0P3SCwz2yDhM4hRYhW1iRJsINW12d7e+hxFpxAxNK+AXP3SRD/p0pJTUlAXxGTrF/fX0KxQKBZNM6Bs74/z7mwfY1tCBoQlW15Zxw7nTWTprGAIP2f748M68wR5moCJ/yHfLDrzt9XlPPeuVz3P4useG7RcvhDOeL+Q1sKWkvMib6QhSwz4UCsXJmBRCf7wjxtPbjvHa3iYMTePeS2u5ZlE1Rf4R+Hi9+uPzhNod8j332U8QDG8DLKq3fQdLD6Jb0UzRFSmG5Rc/NeDB59Eybpi2LVWbpEKhGDATWugPt0b42+d38/reZjQhuHphFXevmk15kW/E3sOIteRtgopWnJcv1JrOkWsfZeGGi5yCbMs29t75Nme9+nkCrbuJViwdvOUwzilVxY57pNfQkDlTqJTIKxSKwTChe+8Egvc/aucTF8zisXsu5MGr5o+oyAMg80X1yNWP9hFq019BtGIZUhhOxB6q5vDH/9tJ64R3UvviXc55wcpTinzIp+PzaJxdWUTVFD9ew/kTqfSMQqEYKsOK6IUQXwX+P5xg9wPgc0AQeBz4/9u7/9i46/uO48/X+WzHd05Cfjgm5AekrJRm2hAhC6HdKrRUW5uiZZOg4o+1ZkNi1dR12z8jW6X1X7pNkzap3YQ2uqyiVJS1I6XpBkNFm7SSEULCjwXKj9AkzS8HRgt2SGLfe398v47t89m+8/nufF+/HpL1/Z7v67v3vW2//fH7Pt/P9xrgTeDTEfF/dUU5jY2rCuz7wq/w9tDFRjw8kBTnof5t422X4prJB0SJTfvupDB4mOG+Gzj6yW+wad+dFE8fAEarXnK4K5/jikIn/cuWEBEu7GY2b+Zc6CWtA74AbI6I8+m1Yu8ENgNPRsR9knYDu4F75yXaCsZGvA0jcfS2hye3XSbMwBmfeTNK4dzzdL/z+uXbl5cinqEv35ETXXlxbV/v5eLuIm9m86neHn0e6JF0iWQkfxL4U5KLiQPsAZ6igYW+KSZepm/imbD9N3F05zcnvdF6YcV147fXbOH4jq+OnzRFsrzNks4OchIrCp0s7+l0YTezhppzoY+In0j6K+AYcB54PCIel9QfEafSY05JWlPp6yXdA9wDsHHjxrmGMf8qzJef6PIIfmyq5PtvT3mjdbo3Xnu6OlhZ7GJlsYuR0dL8nLRlZjaLOVcaSSuAXcAm4CqgKOm3q/36iLg/IrZGxNa+vtZcF3WKdLR+/UM3s+l7n04WKysz0rOaof6bxt947Vk9PuIfK+rlt0muwXr1qgIri10ALvJm1jT1tG4+DhyNiEEASd8GPgKckbQ2Hc2vBc7OQ5xNMWndmuneQJ1hxF52GEs6O1hR6KS7s4NCZ4enRZpZS9RT6I8B2yUVSFo3O4ADwBAwANyXbh+tN8imqLBuDQREcuG/SYV9Ys++zJLOHCuKXRS6Orz+jJktCPX06PdLegQ4CIwAzwH3A73Aw5LuJvljcMd8BNpo+fNvUTw7dmJUDpVGuP6h7UnBl8rWqZncdsnlYNmSTi6Olri2r7c1L8DMbBp1DTkj4kvAl8o+fYFkdN9WxnrvxTPPMrzqFygMHkKUKJ45ANJ4O2d48PI6N2Ntm/UrCvR257k0OrWnb2bWau4tjBnrvQ8PsuHJ3wdKyaUA8wXeX7WZ4tmDDK3ZwoYffJ7imWcZWfdLXPrMXro7Oy9fMLwj5wXgzWzhcaGfSDlQ7nILB6Bj9DzHf/Urabsm0jVtRug8+Qydl96B7oqzR83MFgzP8StzefokyboOQ/1bGSmsIb+8n+t/7lrYcDPk8sm2uECmhZqZzcAj+vITpMaWPBgevHwpwJ7uDjasLKBcDu56DIbPJUXeZ7SaWRtY3IV+bDmD0wcY7ruBN277FvkL77B09VpGC+tY3tNJb3d+8slNuRz0ul1jZu1jURf6/PlBiqefQZQoDB5k89d/kVzpAtpwMww8lhR1M7M2t3gr2dgJUunsGgG5kSFUGoHj+5P2jJlZBizaQp9//y2Kg+Oza4Dx/atu9ButZpYZi6rQJ+vP5FjWk+fqjVej9dtAHWj9drjqpgkH5tOlD8zM2l/me/S5HPQt7aa3O8+S/ISFxUZH4PQLEKNw5kXou278i47/EN49A8vXtiZoM7N5lOkRfVc+vfbq0iUUuvKTV4889yO4+G6yf+k9OPnc5C/+1meg5CUNzKz9ZW5EL0Ghq4NVvd0s7+mc/sC+66F7GVz4Gcnfu7KifvK55A1ZT6U0szaXmUK/rCfPymIXxfKR+3RyOfiTo3Dsafj6rnT0Llh7Q9LK8ZmvZpYRbV/oe7vzXHdlL935KhcUK5XGz2ztyMM1H02K+vH9yXbguzD8ls98NbPMaPtC39NVw4qRpRLsuW28qH92L5x/e2pxd7vGzDKk7Qt9TYbPJUW+NAI//iE88OtJL37DzXDX9zyCN7NMyvSsmymKfbB+W3qjBD85kEyvPPbfMNQ2l7Y1M6tJXYVe0hWSHpH0sqQjkm6RtFLSE5JeTbcr5ivYmpVK8N7Z5OSnUgmGBuH2B2DKBUIEeDRvZtlUb+vmb4B/i4jbJXUBBeDPgCcj4j5Ju4HdwL11Pk/tJvbjx0bxJ/4n2V+fvvnaVYRLQ7Bhu/vyZpZZcy70kpYBHwPuAoiIi8BFSbuAW9PD9gBP0YpCP7Eff+xpyAlKo0mx/+OXkitGFVZ5ho2ZZV49rZsPAIPA1yQ9J+kfJBWB/og4BZBuKw6VJd0j6YCkA4ODg3WEMY3yfnxncfzKUL39yQg+15FsXeTNLMPqKfR5YAvwdxFxIzBE0qapSkTcHxFbI2JrX18DTkyS4PavjffjLw7B7/2XZ9eY2aJTT6E/AZyIiP3p7UdICv8ZSWsB0m3rprMs7U/677k8bNwOaz7sIm9mi86ce/QRcVrScUkfiohXgB3A/6YfA8B96fbReYm0GqVSOk1S4y2ZAV/j1cwWt3pn3fwB8GA64+YN4HdI/kt4WNLdwDHgjjqfozqlEvzTp5I58Qg23pK0aXyNVzNb5Ooq9BFxCNha4a4d9TzunIzNsgEg4MR+rz5pZkYWzowdOymqsDqZUQOAkt68V580M2vztW7KFykb+G56UW952qSZWaq9C/3Ek6KO709Wolx6ZaujMjNbUNq7dVPsS0byYydCuVVjZjZFe4/oPX3SzGxW7V3owdMnzcxm0d6tGzMzm5ULvZlZxrnQm5llnAu9mVnGudCbmWWcC72ZWcYpIlodA5IGgR83+GlWA+ca/BxztVBjc1y1cVy1cVy1K4/t6oiY9UzRBVHom0HSgYiotNJmyy3U2BxXbRxXbRxX7eYam1s3ZmYZ50JvZpZxi6nQ39/qAGawUGNzXLVxXLVxXLWbU2yLpkdvZrZYLaYRvZnZouRCb2aWcZkt9JL+UtLLkp6X9B1JV0xz3CckvSLpNUm7mxTbHZJeklSSNO1UKUlvSnpB0iFJBxZQXE3NmaSVkp6Q9Gq6XTHNcU3J12yvX4m/Te9/XtKWRsVSY1y3Svppmp9Dkv68CTE9IOmspBenub8luaoytlbka4OkH0g6kv4u/mGFY2rPWURk8gP4NSCf7n8Z+HKFYzqA14EPAF3AYWBzE2L7MPAh4Clg6wzHvQmsbmLOZo2rFTkD/gLYne7vrvS9bFa+qnn9wE7g+4CA7cD+JnzvqonrVuCxZv08pc/5MWAL8OI09zc9VzXE1op8rQW2pPtLgR/Nx89XZkf0EfF4RIykN58G1lc4bBvwWkS8EREXgW8Cu5oQ25GIeKXRz1OrKuNqRc52AXvS/T3Abzb4+WZSzevfBfxzJJ4GrpC0dgHE1XQR8Z/A2zMc0opcVRtb00XEqYg4mO6/CxwB1pUdVnPOMlvoy/wuyV/AcuuA4xNun2BqUlspgMclPSvpnlYHk2pFzvoj4hQkvwjAdJcUa0a+qnn9rchRtc95i6TDkr4v6ecbHFM1FvrvYMvyJeka4EZgf9ldNeesrS8lKOk/gCsr3PXFiHg0PeaLwAjwYKWHqPC5eZlvWk1sVfhoRJyUtAZ4QtLL6SiklXE1JGczxVXDw8x7viqo5vU37OdqBtU850GStVHek7QT+Ffggw2OazatyFW1WpYvSb3AvwB/FBE/K7+7wpfMmLO2LvQR8fGZ7pc0ANwG7Ii0uVXmBLBhwu31wMlmxFblY5xMt2clfYfk3/O6Ctc8xNWQnM0Ul6QzktZGxKn0X9Sz0zzGvOergmpef8N+ruqJa2LBiIh9kr4qaXVEtHIBr1bkqiqtypekTpIi/2BEfLvCITXnLLOtG0mfAO4FfiMihqc57Bngg5I2SeoC7gT2NivGmUgqSlo6tk/y5nLF2QFN1oqc7QUG0v0BYMp/Hk3MVzWvfy/w2XR2xHbgp2OtpwaaNS5JV0pSur+N5Pf/rQbHNZtW5KoqrchX+nz/CByJiL+e5rDac9bMd5Sb+QG8RtLHOpR+/H36+auAfROO20nyzvbrJO2LZsT2WyR/lS8AZ4B/L4+NZPbE4fTjpWbEVk1crcgZsAp4Eng13a5sZb4qvX7gc8Dn0n0BX0nvf4EZZlY1Oa7Pp7k5TDJB4SNNiOkh4BRwKf3Zunsh5KrK2FqRr18macM8P6F27aw3Z14Cwcws4zLbujEzs4QLvZlZxrnQm5llnAu9mVnGudCbmWWcC72ZWca50JuZZdz/A8bJrf96EiNkAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import pystan\n",
"import numpy as np\n",
"import pandas as pd\n",
"import patsy as pa\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# data \n",
"d = pd.read_csv('../rethinking/data/Howell1.csv', sep = ';')\n",
"d['weight_s'] = (d.weight - d.weight.mean()) / d.weight.std()\n",
"d['weight_s2'] = d['weight_s']**2\n",
"d['weight_s3'] = d['weight_s']**3\n",
"\n",
"def prepare_data(formula, data):\n",
" y, X = pa.dmatrices(formula, data)\n",
" out = {'y': np.asarray(y.ravel()),\n",
" 'X': np.asarray(X),\n",
" 'N': X.shape[0],\n",
" 'K': X.shape[1]}\n",
" return out\n",
"\n",
"data_lin = prepare_data('height ~ weight_s', d)\n",
"data_qua = prepare_data('height ~ weight_s + weight_s2', d)\n",
"data_cub = prepare_data('height ~ weight_s + weight_s2 + weight_s3', d)\n",
"\n",
"# Stan model\n",
"model = '''\n",
"data {\n",
" int<lower=1> N; // number of observations\n",
" int<lower=1> K; // number of coefficients (including intercept)\n",
" vector[N] y; // outcome\n",
" matrix[N, K] X; // regressors\n",
"}\n",
"parameters {\n",
" real<lower=0,upper=50> sigma;\n",
" vector[K] b;\n",
"}\n",
"transformed parameters {\n",
" vector[N] mu;\n",
" mu = X * b;\n",
"}\n",
"model {\n",
" // probability model\n",
" y ~ normal(mu, sigma);\n",
" // priors\n",
" b[1] ~ normal(178, 20);\n",
" b[2] ~ lognormal(0, 1);\n",
" if (K > 2) {\n",
" for (i in 3:K) {\n",
" b[i] ~ normal(0, 10);\n",
" }\n",
" }\n",
" sigma ~ uniform(0, 50);\n",
"}\n",
"generated quantities {\n",
" // predicted outcome\n",
" vector[N] muhat;\n",
" for (i in 1:N) {\n",
" muhat[i] = normal_rng(mu[i], sigma);\n",
" }\n",
"}\n",
"'''\n",
"\n",
"# compile model\n",
"sm = pystan.StanModel(model_code = model)\n",
"\n",
"# fit models\n",
"fit_lin = sm.sampling(data = data_lin)\n",
"fit_qua = sm.sampling(data = data_qua)\n",
"fit_cub = sm.sampling(data = data_cub)\n",
"\n",
"# extract posterior \n",
"pos = fit_lin.extract('muhat')['muhat']\n",
"d['ub_lin'] = np.quantile(pos, .975, axis = 0)\n",
"d['lb_lin'] = np.quantile(pos, .025, axis = 0)\n",
"d['muhat_lin'] = np.mean(pos, axis = 0)\n",
"\n",
"pos = fit_qua.extract('muhat')['muhat']\n",
"d['ub_qua'] = np.quantile(pos, .975, axis = 0)\n",
"d['lb_qua'] = np.quantile(pos, .025, axis = 0)\n",
"d['muhat_qua'] = np.mean(pos, axis = 0)\n",
"\n",
"pos = fit_cub.extract('muhat')['muhat']\n",
"d['ub_cub'] = np.quantile(pos, .975, axis = 0)\n",
"d['lb_cub'] = np.quantile(pos, .025, axis = 0)\n",
"d['muhat_cub'] = np.mean(pos, axis = 0)\n",
"\n",
"# sort data for plot\n",
"d = d.sort_values(by = 'weight_s')\n",
"\n",
"# plot\n",
"plt.fill_between(d.weight_s, d.ub_lin, d.lb_lin, alpha = .2)\n",
"plt.scatter(d.weight_s, d.height, 5)\n",
"plt.plot(d.weight_s, d.muhat_lin)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fc33a5d46d0>]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydeXwc9Xn/39+Z2VuydawkX5IvbIzBNj6wIW2aBMJlbghgAgnNUZoE0qT9JWnTNCFpoGmbJk1ztqS5AzaEQDDmyMFRSAK2ZYOBgPElbPmQtCvJlrT3zHx/f8zsJe1KsizZlvx9v16w0uzs7ncl+ZlnP9/n+TxCSolCoVAoJi7aiV6AQqFQKMYWFegVCoVigqMCvUKhUExwVKBXKBSKCY4K9AqFQjHBMU70AgDC4bCcNWvWiV6GQqFQjCu2bNkSlVLWDXXeSRHoZ82aRXNz84lehkKhUIwrhBB7h3Oekm4UCoVigqMCvUKhUExwVKBXKBSKCY4K9AqFQjHBGTLQCyF+KIToEEK8VnDsbCHEi0KIl4UQzUKIlQX3fVYIsUsI8aYQ4uKxWrhCoVAohsdwMvofA5f0O/bvwJeklGcDX3C/RwixEFgDnOk+5rtCCH3UVqtQKBSKo2bIQC+lfA7o6n8YmOR+PRk46H59FbBOSpmSUrYAu4CVKBQKxYnEtqGvA05Rt96R1tF/Evi1EOI/cC4Wb3OPTwdeLDhvv3tsAEKI24DbAJqamka4DIVCoRgC24afXA6tG6FxFdy6AbRTa3typO/2o8DfSikbgb8FfuAeFyXOLXkJlVLeI6VcIaVcUVc3ZGOXQqFQjIx41AnytuncxqMnekXHnZEG+luBh9yvf0FentkPNBacN4O8rKNQKBTHn1Cdk8lrhnMbOvUSy5FKNweBdwDPAucDO93j64H7hBBfB6YB84BNx7hGhUKhGDlCOHJNPOoEeVFKeJjYDBnohRBrgXcCYSHEfuBO4K+A/xJCGEASV2uXUv5JCPEA8DpgArdLKa0xWrtCoVAMD02DivoTvYoTxpCBXkp5U5m7lpc5/27g7mNZlEKhUChGj1Nr61mhUChOQVSgVygUigmOCvQKhUIxwVGBXqFQnBxM5O7VE/zeVKBXKBQnnmz36tfPgB9f5nx/DMiCgGrZEimd/3qTGfflJMmMhWUfh8A7yu9tJJwUowQVCsWpjYxFEP27VwvKIW1bkracAGloAl0T9KZMpIRYykQCAY9OyrTImJK+lInX0PAZGomMRchnIKUklrKYXi3Y3x3Hq4HVF8GobMBjaPg9On6PjpQS05KYtsS0bUxL4vNoVAW86JpAEyCEIJmx6E2aGJrA0AWGprm3AlFYq1+qM/c4l3qqQK9QKI4bpmVj6FrR9wB9WhWhaedgHNxMZto57O3zY/b2EK7w0dGbLEqCA14dW0pSGRshyqshibRFIu208aQyacDplWqJxEDazHjsRkLtzcQaVtBy2f0gSggc0sZIdHI4EKb9SCp3WAiQtnOfGQgPaMLSNYFHFxi6hiFC1E89B++hzVjTziFlVGOYFh5NQ9OOT/OWCvQKhWLk2PbAjtMSx5IZC6+u8VZnnAqfgRDQm8yQzNhoQmDZEnHJWvREJ6Y/DKYTvduOJAe8ZDZ4I230eOlAW47sRcFIRAi1bUQAobaNGIkIZrCh38k2sx+7kVD7FmINy4suBtIufx84cpFlS8g4V6jDl6zNXxSi8dx5XkPj9CmVw1r7saACvUKhGBn9XCF71zyMZUuCa6/G42bm7dc8SMKUpE2boFcvyrKzWG70lWiYgWH60OSC8BAZednHi8G/B4xEJ6H2LQhpEmrf4gTqYN2Q95VEaCXvt4/T5qzajFUoFMOmcJOzK3IA6WrPct9GDh7YT0fbATwHNiNsE8+BzfR1tpHK2K6WPnpuKLmMXFq5jLz8om2MeKRI4zEDtcTDS5HoxBpWlgzCZiBMrGE5UhjE6pcBMvccRfc1LHcy9ZMYldErFIqSSOlsSHp0DSklPUmTWMrEsiUZyyZlTcLXsJxgmyNfpH21AMQaluckjTELgOUyctlPN+8vv6x2JJTGZ+4gGH2FeN0S99NACelHCFouux8jHqHxmTtYsPbcIpmm5bL7y2r0Jxsq0CsUigHE0ybd8Qw9CaccUdcEGcseUBm4Z/XAYHdUAbB/YB4mZrCOWMPK/AUlWFdSUy+SWNqamfPoewhGtgEWAghGX8FIdpWXXYQGQist05SRY05GVKBXKE5x4mmTRNrCtCW6JuhJZIinraJqFtPqpyUXBOgBwW64AXAkOnvB67Zc/kDRRcKIRwcGZH8N8fBiJ3t3bwUWEpBCdyQZaWHEO5z9gRIXm6xMM+afUsYQFegVilMI25akTNupSZfQHU/TmzSP7kkGqUbpf95g2fqwKl8Kn8uVUApft/CCMiAg+2uY/fgagpFtxOuWsOfyB5n9xE3O/fXLaD3/2zQ+fQcL1jpzk2JTVpV+L1kJZ5zINKVQgV6hmOAkM84m6OF4hoxlcySROaZO/AEVJ/GIk8UHwoB0AqKvijkbricYfaX8xWAYlS/OcffC0tZMVnIpWenSLyAbCTfDx3IkmlR3/n5/Db7une79DqH25uL3Isro/uMQFegViglGynRkl/aeJH6PTqQ3NaoWK0WZc/2yfJZdvwyEINTWjGX40c1Y+aBMCZ09UIsRjwwIqLkLS6HkMgwJpaTkIgRmoDZ34bCMALrZB0CsfvmATwzA8D69nOSoQK9QTACklAgh6IqlifalsGynjb8ncZSyTNGTlslkCzJnkCxYe24uu0eAwMoFeQnEw4tLB2Uh8jq7K7MMCKjSBimLLiyt53+7rFTUPyiXklwKLxy6mWDn1U9iBZ31Fb4X5/1xdPXyJykq0CsU45iuWJrOvhS2dCpj+jcjjZhyOrxt4uveRap6vhPwbMvZ5IxsIx5ejNS9hNqbsTxB9Ezc0caveKi85OFu3BrxyMCAms283QC//aYXMP3h0hcEXM2/vdmprS8IymagFiMRxfTXOBU2/tqiTD9VewZZL4VSm65juhErXVfLMZ5lqwK9QjFOkFISS1v0JU1sKYmnLZIZa0ycb0t2fvqrOeNnS9AzvVieSt645SVmP3EzwcjL2EaQYHQb8fAStq/5Q64s0fTXYiSH1rcLq2OyATWnsUuTUMdWQMNIdpXOsKVN41O3g3TknVxQLtD3LU8Q3Uzk6+mT3cXrKrPpOmYbsdKm6dEboX0LNK5yBphrYyMLqUCvUIwD9nXGSVt2LmMfzMxrACPYTMwFXjdTN/21+Lp3oGd6EYCe6SV4qNmVQGw0s8+tS3+J+Q++G81KFmn28fBi9l34P05VTf81SLuoOqZl9TpXSw8Tq1+Wy+izgdsyAs7Fxghg+msA98LUsdWVi3Ra3/Udt+Syw8nysXJrD7U1YyS6MEP1bjVPNP+zKVUaOkb18kaik2DbFpBj72o5/nYVFIpThCPxDHs7Y+zq6ONIIlMkywwI8iXa/LPHZz92IwvWrmL2Yzc4UsFQ5ALvy9hGgGB0G7Mfv5FU1WlYRgXZV2jY8jXi4UVIdCw9iAQEOEHfzbhDbZvdqpeXWLB2ZX4Ntomvczu4DpBF1TGJqPtebDf44l4cJEayCz0Tcy82MYxkF9DPkmDKCldWMmn6zW25LD+/dovGZ24H23R+NvetZM76a8AePYuG4WAGwsSnLAfNcDL60Nhp/yqjVyhOErKe67GUSSJj0ZcyyZjDSNulzewNN+R168sfKNCtiztDjXjUyWQHIR9485l6qH0LRqqbndc/xYJ1b3N08MhmQCMeXsKBP/sy8x65PL8koROrW4awkgQ7Xy0uYYy1M++X73az8hBv3LKtZBVP9hOFwHKzcGfT1vKEXPkolMvoB8guSKcLNvpSLsvfc/kvmPfI5TkN39e9i1Cbk+0HI1uZ8+h17LnyoeNXVSME+654gDMmpU+8Ri+E+CFwOdAhpTyr4PjHgTsAE3hMSvkZ9/hngQ8BFvA3Uspfj8XCFYqJgmVLYmmT1q44hqaRsex8Yj4M2cXZyNzkBuRNToliyGk8MgNhYnVLCXVsJpvJDlUimCtLbCvYVA0vxvRVYyS7iDWsKKhptwlGX3KDvEC6OXO8agEIQbDztaLjsfrl6KnDeQnIjDH30WvYfeXD+A7vwfJXs2DdeQhp5rR/zewjY4Q4kA7SYHdCJklG6uxMhTnY8gbxypkcjqVpj3bSI330pVpYWG3z2qEq+ljNARlGBmp5+fc2deJOZqR346mopuUljVbra1xiP08bNQQPpkn+4U38wUoCXp2g1yDo1Ql6daqDXqqCHgIevXioyLEitOMyhETIIYQ+IcRfAH3AT7OBXgjxLuBzwGVSypQQol5K2SGEWAisBVYC04DfAfOllIN+JlqxYoVsbm4+9nejUJyslPFo70lm6EmYpTdVh9mBasQ7WHDfilw54/b3NmMG63PPMWf9tQQjrn4tDLbftLFYcy51Mcke805i7qPXEuh8Pb+ZWb+M1nd+g/kPXohmxQByr50Ngdm3UnhcorH9pk2YgTBn/PTMXAnmEVnBW9Xn0d7ZTU/laUz3xujsirLLv5idMR+/tZYTJMVhKjA0QUAm6ZW+Yf/oBTY6Ep/PS5VP43DSIpYulrCq6MMUHmLSx2AR0Wdo1IScoB9PpOiIWVQFPVQHvVS7t1UhLzXZr93jVUFnOtWA350uOGPqpGG/lwHvTYgtUsoVQ503ZEYvpXxOCDGr3+GPAv8qpUy553S4x68C1rnHW4QQu3CC/gtHsXaFYmLRz7fdfv+jtPel6exLD7qhOlzPczNQR2zKqpxnTKGnu5GIEIy8XL6m3Uozd/01BLrecAL4Bd/BDNS73aG1zFl/LQFXesltZrZvoel3t6FZsaLAnn0NCr52vteQ2MSNyTRHdfZ0HuRHfd+nXuthHnv5g70I2twHdhW8MWcoFCvEm8wKpak++0q64hl2tPfQ4Elx0YFv0yP9zBCd1IluZop2JtPHDtmEQYZJIkGIBBUkkejE6xZT0fkKsWnL2b16HUJodMczaNjU0oMVrMPGuQAn0hZx97YvaXI4keFwPE13PE1XLE2i9SXC6cP8eciktf5CDscz7O2K8/L+wyXtmAUwdbKf0+ormFtXkbudHPSU/wMYRUaq0c8H3i6EuBtIAp+SUm4GpgMvFpy33z02ACHEbcBtAE1NTSNchkIxDohHka0bEbaJbN3I7r1vkfQOXY89bDOtcl4s2ZJDbCfI1y0rrmm3TadcMtfBusm14l2Rez5HI3cCtm1UoFlJ4rWLCrTv3CKISw9dTOI5axEd1LBDNvEbaymVJOjGnaL06Bu5R3TYk6jTJnGd9hy1eowe20eAFOdqb2Kffyd//n83Ukmf80nguhdANzAD03IlR7MfixNqf95dz87ces4Ub+UuMJYnhDB1EuHFVERfyV00va5jZU3IC4CF330ErmRjUFvm92LEIyxY+zmE10SaBtvfvqboApw2bfeikKHbvTh0xtLs7Yzx+qFentsZzZ07bbKf61c08rcXzh/ir+HYGGmgN4Bq4FzgHOABIcQc8p/cCimZs0gp7wHuAUe6GeE6FIqTlt5khkTGoifuZ2r9coLtW4jVLyfpKRdCXAodGodbw12iBLC45FCw7/z/LqrT9nXvKupgBdxBHu6mbbCO2JQVhNo2E689kz2XP4SvpwXLV8W8tau433oXnzM/zJQgWJqXSF+65NIWaXuw0PB6Ayx/+2pW1Fmc/tBFVHOk6JOA7Q2gWQkAYtu/BA0Lke2bAJt5D707XwPf3w/eX8vsx28k1LbZ2ajNyksXfNet4+/KnzMKjU9DXYC9hkb9JD/1k/wlH384nmZ3JMauSB8tkT6M4zA3dkiNHsCVbjYUaPRP4kg3z7rf78YJ+h8GkFJ+xT3+a+CLUspBpRul0SsmCrYtiWcsDsfTHElk8v7tw61lL6XLQ1E1Sf/BGmX1dX8tszdc727E4jQ5vW+bU87nLJYzfnoWutmHJfwkw2cRijj/DmMNK2k9/1sgJU1PfYyejn18zvwQT1tLCIo0CelBlsjrlk4LsKT3OT6U/hmTiZHEQ53oIV67iD1Xrgddd7LxDdc7HbRGoOBio4EQCGkhhcHOqzYw75HLnO+h/B5D0Xt2u19L/ZxH05xslJ7rpNHoy/Ar4HzgWSHEfMALRIH1wH1CiK/jbMbOAzaN8DUUinGFlJJDPUm6SmW2w2y6KeUMOcA0rGBa0uzHbxrUhOvg277EvF+tzmnsvu5dpGoXuGuCZM0ZhDo2o8sk/u43MKVGhCoeODCD53/8OHXiCIfk1WyTp+XWGJdeljYYVEU2s0Ds4zr9eeZqB0HoxPUlBDPbEMIJznr1HOKGl2D0VWY/sSa3Rqc+XpCsPoNA9xtoZgzLCJKsXUioY6tjTVCzIFfhU7gRnBvpV+S/k//5DjZEZNQan8bR0BEYXnnlWuCdQFgIsR+4E/gh8EMhxGs4Wya3SuejwZ+EEA8Ar+OUXd4+VMWNQjHe6Yo5Omwifex2BP27QbOBvcg0zJVXgoc2D9ishQITrrZmkBqWpzJnW5Cqnu9ko7F2vId3E2zfwqtyNjvkDN7IzOQH1uqi9dTKI8wWh/iwvoHr9N8zQ3QQ8Hp449ptzFn/LwQ78xu9Qlq57tbs7b53/0+uXHLgGi1CkZfA/Wygmwlaz/9OkU1woa2wkYiWHOkHTAgr4bFkOFU3N5W565Yy598N3H0si1IoxgNSSg4dSdJZRpse4bPmu0EBpF3UTJSzAfYEmfPEe4v8W7IWAdlGI8sTZN76y4nVLeXgn32ZVPUCbCT+hz/MXx26nH2yAfh57pU1HJ3pWv9W5pvbuVhrZpbWRqJ6IVLTCXbuc4K6mcFIdbPvontyQzsgax+8gpbV6/LyCaVNwYreExTdn+12BYoz53Ij/Qr8bBwTtV+Cpo/i72T8ozpjFYqjJJmxiPalMC159NOZCimRhealGytXBeMYhf3Rnb4k8XXtcOUYCz0TZ+c1T5CqOR2QzH7sRoKRl4nXLCTY/QZCWsQ69vBaT5ANm16jeX8fPakPFS3jev1ZjjRdzM1/cRb14jBNv76LYPfruUx9/9u/yrz1VxSUaC5xM+xOx0++w9lkbn3XN9FTR4ptGEpk5UWbzK49cfbCVs6ZEspvghqJzhPb4ToOUIFeoRgGR+IZIn1JbAmpzDD8YoaiTDOU6a/JmXYBuc7Tpt99JBe8UjWnu52rm4nXLiRVNc8x8Ip1EGrfhCU13or2sLPyOn7afSZb7Xnw2H4AFnrbOUSQeWI//+i5lwViH+aUZbSs/hwgaXrkI0VBPl51BtNe/BK5Es3wUvZd8L18QK5bys6rHsXyVTPvlxeiZ3qRQkfgdMFm6/ILh304Y/1+AYCRLLywFUhTw5gglb04moGwKxdtzQ/8Hqe+8WPFsKpuxhpVdaM4WelLmbQdSZBIDxHcj1IjdmqxVyGkWVRJYsQjLLhvZW6aEhRUm6x5IXcxMGJtzP/lhWhmn+MX875X6OmO8IdffJ3vWlcPeL0AKX7u/ReWij0IYSPcZ5fAzqufJFVzOr7uncx7+BKEK+Fka9HBzlXFxOvOJhjZRnakX//oUdwJC7hyTuv532bBfefm3le2Jr9og7n/ZvNlDwxfb7ct5jx6XcHowqN47AnkZK+6USgmBiWsCcAZx7e/O0G8RJfjAIY7LLuAcjKEGQi7tevNxGsXgSYIRl4hXruIxqfvKBjqESMjBb+1V7A1M497vrfRfeZ8kP/kBfM4Z1YNk/wGsx+7gVB7C7H65YDTHJVl3q8ucTds41h6AN21NQDnE0V2fF+8djHB6Lb8SL/cOe6PoeBYUV1++xawJfHahTmDs5xZWsfW/AWsVPnocNF09lz5kNqQLYMK9IpTl37WBOb71hOJZdCEOKo5q2WtCgbL8st1swpBy6X3Mnf9NQS7XiNWt5R43WKCkVcAizZZw0/iF/Hf1pVFT7dKe4ON9hn8/cXzecc0sIIFFy5p0/qub4O00FNHSE2eS7BtM3OevDmXvedMxqwYoBGvXYI0PE6pozu+r/HpO3KWv/Hw2QSjr+QuBGSPdb7q1r272X/0FceR8tmPuwZnjh2CFM5maaxhuePLU64n4GgYZyWPxxMV6BWnLvGoE+Sz1gRv7SXtH2bHZEGDTskRdMPJ8ksFJmkzZ8MNOX+ZUMdWknj4jbWUu8xbOCRrMPv9s93k+yjBKafTcvmHQQis7PpiHSAtmp76GMHIy8hsMxIaAgZk7/kM3SbY+RLbb9oIwiie9pQ9T0Jsyjm5Tx4H/vxfsAJhGp/9eEGd/zq3gsadK4tEZtuspM3Oa37tbCK7QX4iDOE+WVGBXnHKkvHXYsxYBa0bidcvJ+0bwpogS6nxdNmZpgEnOzXiBWPwhvKBL8hkC/1lDsoa/lF8kmeT+Walc8Imd/Z8Fj8Z5oiDaMLV2c/7cpGHzZz11xGMvpR7nBNcszKLXZC9Zw0S8v/Po+c+mSBl8YZn12uO5AI0Pn078x5x6u9jDSuLfg7ZubKFG8wSsDwVpKrnlag2Gt9DuE9WVKBXnJJYtqS1O0H60vuQvdHB5YJ+kkLhRKSco6M70xSkE9T9tW7j0yYG9YEvyGS761fy1PL/5pv6P/FE7PTcKYYm+MvlYS4+ezYV6S5Of2A/wu1DzOrnqep5ztSmrh1M//0/5kzHci9TcFvsMilzz5HX791PJoU16m6FTbxuKcHoqznJpTDTBwh1bHF+DgU/S2cqVLxoo1Y3E06tfbabdbgGbooRoQK94pTjSCJDpDfljuYTMFjmWEJSGDCYI9uw5K8pOrf1/G/lpzH1z1Ld0X96IoqvbSsvWPP5+r4Lad77GuAE+TOnTuKO809jRnWweB3hpUTP+iAzn/24G/AFRiKSm9oEDMjNC4/Jfl/HGlY6BmCBOvpvhhZ9Mom8VLxxmp3r2rA8t7nr2CQXB+n8BvNmpBAgraJ5r87iyuxZFPwehvT4UZRFBXrFKUPatNnfHS/pF16OcpJCUROQ2wVqJNxB1NlKE5zSwgFZqjv6b9PBNE/a5/CA9RMAAiS5Vn+eW4ynqKufyuErf0LJkYCRrWjbkrmN0Vj9sqKpTdkAbmkBUrULnC7ZIvOwwmCv03r+d50N0RLBc0Cmnd04zZG9nDgjBbODvYtwg7ivawfzHr7U9d2JF2X0znllNlNt0x0L+Ep5j5+CT0q6JrClzG2m65rAsiV+j4Z0/w4KN9q9hpYbti6RmJY8ZiuLkw0V6BWnBL3JDB29qfLlkmUyxLKSQn8TLWnT+LuP5YOvK330z1JTmQw73nyDT+27jjdlY+51/kLbxv94/hO/SDsBOLKL9MGNxKeuAq24kUoiCXS9lpdCbJNU1TwsLYhux53lAbqdzFXG6GaCRO0iZ1KU5kO3465WHsQM1JbfDB0i085bIdsEO18dGLyzZBu9ppS48BWQdVHOun5qwmbuY+/B5+4NhNq34Du8i1CHe9Hr2MJUTx9xTy1+j0bAq2NomhO8gaRp4TN0ehIZKv0GmnAuAl3xNLUhZ0pVqclPpmWTsSQZ2yZj2iRNm2TGImPZuQuBEM7EqVS/C8fJiAr0igmLZcvcRKBBu1kHGa6NELSsXuu4PlbPH2gFnJM4IoQ6NueCb+s7v0WhDNKTNPnFllYefXk/Fhr1VPA27TX+zniQFdqOvIYuHBtfKQRznrjRsRa+5SXHOz6Tnegki/TuYGQbRqqbvRfdw5wnbymQbSS27ni8Zz1ofF1vMu9Xl+Q0dT0Tw0hE0ZPd5TdDy2Xa7ibtsLV1IXjrsrV4u3eRrp6Pz6uTNm3qKn3E0xYVPoPakJfueBpNCEI+A/o68HRsy61XTF/KvDNXwBZnE100riLcMKOsfBP0OiGu2h0wAqAhqK8s7RWfxdA1DB0ClPbMyVg2mhDomkBKScq9ECQz7q1pDW+w+3FCBXrFhCOZseiMpemODT6qL8tgw7WRdnkr4IJ2fj0R7feskqr1t/Lrtgq+b19Ju+VMWDpL7OVq/Q/cov8Wv8j0k1E03rz+9wQiLzPzmY/mrIXnPnI1ga7X3Gft/yogBZjeSTRu+UbROdnGpHjtIkdS0XSsgkCc9a1pfOaOgfsNQ22GFn4C6FdxlJVBPLrESHSR8dcyOegl5NGovP9qxP5NyBkr4dYNCE3LD9u2bYhHqC1sXpvcAE3nwr4XYfoy+OBvEJoGt24o2eh2vPDoealICIHfo+P3FF8UTCv/SSCZsbBtckPSHZkI9OO0dhXoFROKRNpid6Tv6D5Kl9CUs5TS6IEiE62FP1mE5komtoRWWc/TD93DV3v+GoAVYjsfNppZ4tnPwno/ofYtWJ4gMmMRrz4DXL/2WP0yGp/7JKG24hEOga4/FXWf2kYIzcwP5UZaBNu3EormLYMT1QsIdG93SyHfyEkqZrDOMSJrbyYeXsK+C+9xbIT7G6SVCkD9ykBzP5eOrei6ju7VqfQbVAW8CGy8P78SsX8TzFiFvPVRtEQn7N8EtukcT3RChVty2q95jVs3ODqOEKWDuqblH3uSYugaFbpGhe/Eh9kTvwKFYhSwbImQNm0H9yL1muFlebmmp3DZ4dpFGr079ML0F5toaXYcJDxur+L2zCecB7rOxcvFm6zzfhlDOLruzlUb0DMx4nVnM/ex9xDsesMxBrt6A5a/1g26Ml8RE16eMzbLviPNShZvunoqiU9ZWVQJFOjega2HHNmmflnB3oKg5fIHCiZWkfe/b1g+aJDPZvCpqSuIXvcg9oyVaAec7Hxm00x8noJw0pcP6uzfiEh0OoG6cVU+mIcK5KCC5jVaNzrfZwP5OAjqJzsq0CvGNb3JDNG+NH2JNLMfu5FZg3VWFurqrqVv4bQmI9k9cMMxuxkZj9D49O0suG8VsYbl7Ln8F8zZcD3ByMs8Z53JrZnPFr3UBu8/Mk1EqRF9+YCsh5jzxE3uhqpAuFp7qGMz8351GbH65QXBOoBuJkA32LP6PuZscKtO6pfnjb/CZ9O+4lPEp6wCXXcqWzq3uxbGNpoVA6G77ycrEFGsuUubQv97I9HhSDBIPKlO7ICTRU81YrkNUH9bMzN8CfjAYxCPIkJ1+CxHC40AACAASURBVPpfHEoF9XLZebnzFaOGCvSKccuReIqD+1sHSgmlOiv7VZW0nv/tgZLMoBcGinR84p00X7COnz3/Jr/b7dSuX2Fs5EN/+VHmP3E9gc5WR54xDXfwx10gZcFGaHEfqpBWzuBLT3Tmzgu1b8RIdbHnyoeL58a6F545T97sbLS6ewfT/vh5KPCfKVvDP0CCyfvfZ6adA0LgObgZGlchbt0AorJ04C6XaR+t5DLYRUBxzKhArxh3dPQmSaQyhB+8jgUFviqDVX/0vxAgRZEkk5vLWmr2av0yDp73JQDaZDU/NS/k5+v20JOy0QR8zFjPX+pPUid6if/6eQKdr+c2aY3U4XxwjkVyY/0KsY0gmpXK2fRagX5WDFL0q3xxMvDchadtI0YiAmhO1l34UCAeXlxcw+++r8SU5Ry+/iHS087Be3ATSAshLbwHNzsZvm0VyyhHG4iPVnJREs2YoQK9YtyQzFhEelMcSWTQYxGChRl5smvQeu+81t48oMY9Z7pVbvZq+yaqH17D7NR9gDNyr6HCw9tnB7j07Fm888X/ItTeRzy8xK1bt5zhF6nDJW0EDr7tn4tsCjQzTnzSXABnHmrd0qI5r2b/wA9O8C9Aj3eRqjndXcNL+cEhtYvYc8VDIASaBpOtnpwEE2jfgt+IIT70GCIWhQc/4AT2GSudn19/GUUF4nGLCvSKkx4pJZHeFB0F1sElG5myJlqlnyVf0yZd1VxomIFajHg0p43HaxeBtDEDdcQalnPk0B6+YV7DY9a5uWd6xvt3hCfNILj3FWJJV9+PR9CT3Ux74U5ntF7/UXcFFSpWoI49Vz7EnEeuzPmzB3t2Q8/uvG+OKO0Jk6WwesbyhJj3yGXO3sEVD7pdpNuIh5fw1tUPo2kaM6qC6LqgwjspJ8GIxlWISrfTtbKhOGOXUskoEwg1YUpxUtOXMjl4ODGw4cn1inH8VoYORgMmOmWdF7OSTd1ShJUm2LkNgEN1b+dbDXfxq22HSGRsLtE28SnjF8zVDuSeU4DzXDf+IeczYxkV7Lz+KczglILmKsnsDdfn/WCmrHLkISmZ8/ClBN0yyOy/RMuoIFm70PGCH2xakrTz82Pd93Xgg1vJ+GsgFqWmfjqWhJBPx2cU1HiXGbaiGH+oCVOKcU3GsumOpemMpTGtfslI/3b91WsxEiUqZgoo+gRQt5Sm3/110Ug8x3XRpktW8nPr3fygdTU9rQc4b04tH3zbTGZaYeY+/M3iaUpCJ9awvMhnRjf70BPdOR93p6JF0HrBd1x5qHhzdM/Vj7PwZ4vRzD4gn8W3nv+dIvOwkggNM7yA+JTlBNu2kJq2gmnTGrEBWRNymnpKBfWsBGOZEN0BdQvy3gOKCcmQgV4I8UPgcqBDSnlWv/s+BXwVqJNSRt1jnwU+BFjA30gpfz3qq1ZMaPpSJq1d8YEB3qVICmnbXGx4VW5gRUGZpBPkt+Zr0YXOrsl/xrpIIz+zLiKNh/MDO7nqivcwp97paE3J012DsqxLY4Hjo5R5Td2oKJJvckO/A/WlDc50g9ff/wq+rp1Me+ELedlngHlY0VvBo4GMRZg6rYnQXz2J3RfBCIbRdI3cuy/XhAROkP/32ZDqAd8k+EwL6Crvm6gM5zf7Y+DbwE8LDwohGoELgX0FxxYCa4AzgWnA74QQ86WUw7cLVJyyZLX49p7UICe5/io5r3c7F7SHHFghNBBabrCHBDpqVvEp7+d4/i0no75e/z9uqduDdu33i7PcbKNRKblICN543zZ83buw/NVO05M7cMTXtSPXhJS90AwI4JpBKnwGLZcPYtMLBLwa4Qofk3w68ieXo+3fhHADuDapgQGXt8GakKI7nCAPzm10BzQsLP9zV4xrhgz0UsrnhBCzStz1n8BngEcKjl0FrJNSpoAWIcQuYCXwwrEvVTGR6UuZHOqOYfZEoJxcUSDZxGsXgdDd0Xh5GWUoj5asN7p+aBsfFf/I823zsOw+Lj1rChcvrOf0itPLyyVCy3vg9EczSNUuyJt8uR2q8351aVGd+4AyzsJPHyXMwzTN2RetrfBSX+l3nBb7OvJdp/0DeCGDNSHVLXAy+WxGX7dg0J/bqKP2CY4rI/qsJoS4Ejggpdwmin9J04EXC77f7x4r9Ry3AbcBNDU1jWQZiglA2rRp70lyOJYacmZooWQTjL7i2BC4w6dzMsoQQSOesbgr/FV+29FGZ9zEqws+f8VCls+sBsDs/4CjHXCR9V7v3J5veiqscx9iXJ7HEGRM6di8IJhbH8LQtGIr3eF2kQ7WhKRpjlxzIjT6wSQlxZhw1IFeCBEEPgdcVOruEsdKCq1SynuAe8CpujnadSjGN1JKIn0pOnqcksnhzAztX1KZHT49nCAspeSPuzv5/vN76IylaawJcsM5M7n4zCkl/cidB41wYLXQilwinecSmMHy4/Kyy59VEyDe3UaweioeQy+9tqPpIh2s9l03ToxcM5ikpBgTRpLRzwVmA9lsfgawVQixEieDbyw4dwZw8FgXqZhY9CYzHDycJG3mSybLDvjol1H3b4oazgDpHe29/OD3Lbx+qIfZ4RD/cOkCFkyZNOTjjmVgdb7OvWD+aon1ewxByGtQV+kjY5r4770Kf2GmWzJ3Ynw3Lylfm+POUQd6KeWrQO4vTAjxFrBCShkVQqwH7hNCfB1nM3YesKnkEylOSTKWzd7O+EAb4f5BEInR10bTUx8tkGe+gxmoL+nZUiqrbetJ8rMX3uK5nVGqgh5uf+dpXLiwoXwG349jGljd3yUyt3Hr6PCaBn6Pjs/QnJmwgD/VeWpkusrX5rgznPLKtcA7gbAQYj9wp5TyB6XOlVL+SQjxAPA6jtx5u6q4UQB0xdL0JjMkM4OMXctuRuYmPjk5QtZIbMHac4s2NsvJKn1Jk/ubW9nwykE0TXDjOY1cu3R6btpQjqH096EGVg9Giec2dIEmBDOqAwAEPHrxU55Kme54/kQyDhlO1c1NQ9w/q9/3dwN3H9uyFBMFKSUZS3LoSCI3B3QocpJJ9jnc28JmI6CkrPL09na+/3wLsZTJu89o4OZVTdRW+EosbJj6e7kxeoO+6YLqoCnLabvmQdIWzKuvQBMCrdwnCpXpKsYI1SGhGDMsW9LaFSdj2cMO8uBKJvXLCHVsdp5HD5IMn5WzBMgNzCiQVfalQtz/xx08vb2DM6dN4q//Yg6zwxVlX+NY9PehKKoOat/C3GDSsQTWh7GRqzJdxRigAr1iTMhYNm9FYyQHG8pdDiFoveDbLFh7HgIb3UqVtARouex+OtoP8dNtvTx371Z0TXDdshm879yZQ+rwx6S/D750ahqmY81YiX7AbWgK1WGo7FxxAlGBXjGq2LYkZdrs64oXVdUcFdKm8emPkx2gEa9b4tbIg5GIYgbCSOCBLQe4b+NePLrGtctmcNWSaVSHvMN7jWPR3wc+FZP8HiYHPQQ8Ol5Dy01fUhKM4mRABXrFqCGlpKUzRjyZcXT0EQZQI9FJqCPvRROMbmP2YzeAEPjbtrIudDMPBq/n5dYjnDenlo++Y+7wA3whI9Hf+6Frgjl1Iby6Vqy9KwlGcRKhAr3imMk2P/UkMiRSZulNzqOwFc7JKm3NgJXbhN0jp3B35hM8nVxGOBnjA2+bxVVnTx92ueRoIARU+g0mBzxkLEml38Dv0Yd+oEJxAlGBXnFMZEyTgwda6dGqQYjSm5yB2qJyyZwfe7ku05wBWJTGZ26Httf4iudjrO1dgobNNaFX+cD7b0PoxzfACgEza4NU+j3H9XUVimNFGUwoRkTGNDnSsY/0D1bT9OMVjrQi7Vw2LoWR2+QsLJd0auKbcyWSZREaZqieXy76H843fsS9vUu5+MypPDXtHr5uf5U5T6xxnCyPE5V+g9OnVKogrxiXqIxecdTYloX5w8uYdHAzucEdBSWKA2wKchUuTkYfDy/B9JeYg1rAkUSGH/6+haff7GB6VYB/vfYMllSlaVr74piURJYj6NOpCXqpCnoQalNVMU5RgV5xVJiWTWvrPmYdakZQxiK4/yZn1g4g1p6zNJj9+I1l5Zute7v5zrO76OhNccOKRm5c0YjX0DCzFsCjXBJZiqwHzbSqwHHdA1AoxgIV6BXD5kgsRduh/aS9NfmA298iuJStgHsMTScYfaXkUA6Azr4U3/99C3/YFWV6VYCvXLOIs6ZPzi9gFEsiSyEEVPgMpkz2qw1WxYRCBXrFkJiWTbQ3ScW6q5lfOKc12W9Oq20OHOtHgSdN/bJ+QzlWE2tYzq5L17Hh1Tbu3bgP07a5eVUT1y2b4cw87c8olET2J+jTqQp48Hl0Knzqn4Ri4qH+qhVDciSRobvjIA2F1TTJ7uKAK20nyPcb6wcFnjQdW9m+5gX0ZDfzfrUaIU12HDrMp+7fyp6uFMuaqvnIO+YwdXJg5Is9ykEhXkNj2uQAAa/K4BUTFxXoFWVJZiwOHE6QythYQ1gGGInOolms8fBi5xxpEw8vzmX5ZrAeM1jPofB5fOvg6ayz3kWNIfmHSxbwtrm1x7bheRSDQvwejaqgl3CFV22yKiY8KtArSpI2bVqiMUxL5rLkltXrMBLRkplydhZrqK2ZeN0S9lzxECCZ/fgagpGXidedTcvqdSAE21oP843uT9Blp7lyyTTeu6ppoIXwCBjKqEwIqAk5FTSj8XoKxXhB/bUrcvQmM+iaM7M00pfKBflclly31J3X+urAjLnERqkRbyfUthEBBCNbMZJRnnhL8s2ndzG9KsBX37OE+Q2Vo7b+wYzKAl6Nhkl+VQevOCVRgV6Ro70nOWAwiJGIEGpvdmwIXNvgnAYfjxQ7SvbfKJX5zH+vXc8nH23lzUiSxTMm8/nLFo5+ZUuJi82kgEG4wkdIbbIqTmHUX78CgI7eJIl0v05TadP41O0grfzwD1wNvnoBjU/fTqhjS37qUz89PDs39bmDgs+YH4Eek/ef28QVc734jTFqyhYaVqiOmpCjv/sMtcmqUKhAP9Gx7YF2ue4xGQyTsR3v+PYjqQEPLXaR1IjXnU0w8jK2ESTY9TpgO9l920aMRAQz2FD0+MOJDO/s/AzRTJr59RX8/cXzWfX8rYS2Db1ZOlI8huC0uorhDflQKE4R1L+GiYxtw08uh6+fAT++zPm+4Fjmh6vZ3dFDa3c8/5isy6SUmP4a4uHFSHTidWez5/IH2XnNk2hWEkH/7L94g7Z5bxcfX/sS0b40c8IhvnLtYqZ5YgM2S0cTIaCpJqiCvELRD5XRT2TiUWfQtG06t/Goc9w95jmwGfqiZAK1jq7tr2H242ucypnwYqTuyWfw0W3MfuImWlavyzc9GQF0M+GWTTra/JFEhrWb9vHYq4eYVRvkrqvPYmZtCABTH5upTuA0PU2vCqiOVoWiBCrQT2RCddC4ygns7kg7AGvGSrTWTU6w9dfkB1mHF7u18BbB6EuAo8lrZl9+AzbZld/w9NdgJLtyG5/Ne7v46q/fJJ62uHLJNG49b5YzbSnLGFgY1FZ4qQl5VYBXKAZBBfqJjBBw64acRm9J2N8dJ37JWuiLYvpr8HW9mauqCUZeIlF9OoHu7blNV4mG5QnlM/esI6WbwZvBOo4kMvzHb7bzcusR5oRDfOqi02msCZZZ07FbGAgBUsL06gA1I5kspVCcYqhAP9FxR9oVD+sW0G8YiHT/H+jejm1UoFkJYvXLHcMyf21R5l7IwcMJvrj+T3T19vEp42FuCHbQXn3vmL6lxpogk/yG6mhVKIbJkIFeCPFD4HKgQ0p5lnvsq8AVQBrYDXxASnnYve+zwIcAC/gbKeWvx2jtimESS5m0dsfJmPkC+cJhIOAE+uzXmhln59WPYwVqnANlsvCd7b3c+eifEEjWeu9mmXgTGTHoHCOfeF0ThCu8+D2aCvIKxVEwnPKEHwOX9Dv2W+AsKeViYAfwWQAhxEJgDXCm+5jvCiGUeHoC6exLDQjykO8ideQZsIyK3NexhhVMe+ELLFi7kgX35adHZclYNj/+YwufenAbfo/OV69bwvypVUVTpUaTyQEPM6oDnFZfQf0kv6qNVyiOkiEzeinlc0KIWf2O/abg2xeB97hfXwWsk1KmgBYhxC5gJfDCqKxWMWySGYueZKZkfTyQHwYSjwCSxqfvINTeTDy8hNbzv82CteflMvxQ2+acd/y+7gRf+82b7InGuGhhA+8/bxaTA54x8YnXNJg2OUC10uEVimNiNDT6DwL3u19Pxwn8Wfa7xwYghLgNuA2gqalpFJahyBJPm+yJxAbe0d/CV2iYoQaMeMRtjLIJdr4KaEWj/yxPiNMeXs3dvo/zgyPnMMlv8PnLzmDl7IJxgKPsE2/ogunVASYpbxqF4pg5pkAvhPgcYALZ3bdSqZwscQwp5T3APQArVqwoeY5imBR0v/alLaK9qSK/GmBQC9+cGZhbP28GwrlsX0920fjw1dyZ+QAPJM+hLmTwpasX01hdpqqm4PVGkuFrGtRV+giHfGhqhJ9CMSqMuIVQCHErzibtzVLmwsp+oLHgtBnAwZEvTzEkBZ2uqf+9lJaOXnqT5oDTSln45hCCltVridctyc1zBTBDDXSGTuOdmW/wgPUuPjipmR/cunJYQX72YzeyYO2qAfr+YFQFPZzeUEl9pV8FeYViFBlRoBdCXAL8PXCllLKgf571wBohhE8IMRuYB2w69mUqylLQ/eo9uLmsrUBu8zW7YeqvyVkdABjJ7vw8V/dC0NoV5zO/fIV2q5L3LKrmmvd9AqEN/Scz6EWlBLomaKoJ0qjsCxSKMWE45ZVrgXcCYSHEfuBOnCobH/Bbt8ztRSnlR6SUfxJCPAC8jiPp3C6ltMZq8QqwA2HMqefgObh58IqXwq7UrNVBgYxT6OUeqTuXh3ak+emLO/F7dO666iyWNFYNe02D+cL3p9JvML06UHo+rEKhGBWEHCDmHn9WrFghm5ubT/Qyxh2mZbM7EiOdMYenh7u6OUgWrD0XIU2kMNh+00ZnI1XatB44wF3PtHPwSJIzp03i0xedTm2F7+gXN4RGr2uCKZP9qrNVoTgGhBBbpJQrhjpPdcaOU6SUdPSmSJt2ccVLuQBbuBlbv6xkxv2nQ718+Yn9+AydL1y+kOUzq9FGWio5SBVOTYWXKZP86EqHVyiOCyrQj0O6Ymk6+1KkzIGDQgqDeesF38H0hzESUfRkV5Fuvv3GP6Cne0hVz0cC3/jdDp7e3kHDJB//dsl0auuqR60ePoumwSS/h2mT/aqzVaE4jqhAP86wbElnX8r1rKEogy/eBN3EgvtWYRlBdLPPeaxR4X5tMe+hC9HNBJG6c/ls8E6e3h7hmrOn8enOzzP10RfGZDDIjOogkwOqLl6hON6oQD+O6EuZ7C+0M+hfG3/pfY5HfKYXAIGN7loMA+hmDNCc45lettjz+bvW62iVEW5Y0chfLgowdd0LRdUyo+E0WV/po9LvIeBV1gUKxYlABfpxQmtXnMPxTNGx/mWMvsO70TPxAothx2YYdxpUrH4F/u7taOle/te6nH811zDV6OUrV5zFmdOrsKQc1cEgfo9GY01QecUrFCcYFehPcpLpDKkj7RxOVgzQzLOj/oLRV4g1LCdVPZ/YlBWE2poBy83kBTuvfhIr6ATt0+87l6+Y7+X71uUsn+LhM1ecR9DnyimjNBikKuihrtKnArxCcZKgAv1JjGVZaD+9gkkHNjO7v2YubWY/voZg5GXitWfScsm9GMlOWlavw4i1M/+X70Yz+7CMAFagBjNQx+5IH9emfkBCeripchs3XvNR9P7168fgWSME+D3OSD/V2apQnDyoQH+SkrFsWvftZfaBzY4009aMEY9ihuoBV7Zpa84ZkS38+dloVpJY/TKElcmN/9PNPhasO48Xqy7jn4+8j4T0cMvSWq4/76Now+hyHS66JphbH1IWwgrFSYhqRzzJsGxJyrTYE4kRM2qI1S9zXeEsGp+5PecbYwbCxOuW5AaGaGafq9U3E4y+lNPpAR4xz+XWQ9fi0eCba5Zy45+dMepBfk6dCvIKxcmKCvRjhW1DXwcDbSQHJ5Gx2NXR5zZCCVov+A4IPT+cOx5xPeRhzxW/JF63DImO5alECoN4eEnuuaSEr3g/ziczt7PY18HXbljK7HBoNN9lLsgrPV6hOHlR0s1YkHWUbN0IjaucAd3DyKCllER7U9gFfVBmoJ5Yw4pcE1TjM3cUedTsufKhnH+NkezC9Ncy+7Eb8La9zGeMT/OrnjO5ZP5k/vr88zBGMeOuCjobuGrTVaE4+VGBfiwocJSkdaPzfUV92dOllHTG0vQlzYEWwwWVMIUeNYV17tnN0+ztr1f8L99/djvbO03eu7KJNec0jmon6pTJfuoqR+B/o1AoTghKuhkLQnVOJq8Zzm1o8CqWg0eStB1JlvSRB3KVMGagrthquESd+1vRXv7l4Y1EOru4q/Y33HTOjFEL8l5Do2GSTwV5hWKcoTL6sUAIR65xpz4NVo+eMi26Y+mhpXzX6qDl0vvwHd5Nquo0jES0qN79qTfa+d6zu5hkZ7jXezfz4u1sT3x8VEb8CQEza1Xzk0IxHlGBfqzQtCHlGinh4OGkE+QHcZ004hFHm29rxvIE0TNx59ZMEGtYzp7V6/jRC/t4+KUDLJo+iW+Je5jZ2T4q3a0AQZ+utHiFYhyjAv0JIJmx2N8dJ2XazsZruXmu2eMFna56prfoNt32Jv/xxGs819LL6kVTue3tc+gVP2L7MXa3FlIb8qoh3QrFOEYF+uOMbUv2dydIpPOlNaVG75nBuvxxrJxvjeUJoWdiWJ4QXRkPN1tfpGVvHzesaOSWVU2uHi9GxYzMZ2iEfAZVQTUcRKEYz6jN2ONIxrI50B0jc7itqL5+wDxXV24pPr6S7Te9SLLmDEDwesV5XOX5H/bJer54xZm879yZo7bpGvTpzG+o5LT6CqZVBUblORUKxYlDZfTHkYPdMWofvI4Z/SWacmZi/Y4biSihjq1ssufxkUPXIf02d119FgumTBq1NdZUeKmv9KkZrgrFBEL9az5OdMXSxLva8xJNWzO+rh35zD5rJtY/Ky84bgbC3Be6hRvTX8Crw5evXjRqQV4IqJ/kY9pkvwryCsUEQ/2LHmOklBw6kuBAdyIvxaBjeYLMe/hS5qy/Bmyr4AFOlU1RvaW00WMd/OrlA/xT9CJmV/v41gcvZFa4YlTWGPTpNNYEaZikRvwpFBMRFejHECmdjddob9o54EoxO695wh0QYhGMbGXO+msx4u1gW8zecAML7juH2RuudwzMpE3ThjXc+9Pv8YM/vMW5c2r49xuWEfCNjupWFfQwt65CjfhTKCYwQwZ6IcQPhRAdQojXCo7VCCF+K4TY6d5WF9z3WSHELiHEm0KIi8dq4Sc78bRJa1diwFQohEaq5vQi58lg9CUWrD2XOY9cQ6h9EwKbUPsmjHg7Zm+Ez+8/h3utC/iIsYF/ekd4VOrZNc3J5NVmq0Ix8RlORv9j4JJ+x/4BeEpKOQ94yv0eIcRCYA1wpvuY7wohTrkum2wJ5ZFEpvQJQhQ5TwIIaRHs3FZ0WuxwF9f/fBfrrfP4e+N+7pixGztbNllK4jkKpk4OMCccQlcDQhSKCc+QgV5K+RzQ1e/wVcBP3K9/AlxdcHydlDIlpWwBdgErR2mt44KUafFmey+pjD34iZrOnisfYvtNG4lNWeWUUNavcOyGgRZ9Nh/7XRzTlnxMX8/7p7bSsnqds2vqNlItWLuK2Y/dkPOoHw7hSi/hSi81Ia/S4xWKU4SRCr0NUspDAFLKQ0KIbK//dODFgvP2u8cGIIS4DbgNoKmpaYTLOLlIZize6oxhWsPMsoWGGaovLq2UFnt2vsEXfh/DNjP80PsfvEtshajh2BAXNlL1a7AaCl0ThCt8GCqLVyhOKUZ7M7ZUBCkZ9aSU90gpV0gpV9TVHbvp1okm2pdiV0cfGXOIIJ+VXKwMvs7tjnd9QQnlpr09fPrpXgIeg3+/bgkrp3lAG6yRamg/GyGcTP70KZV4dE1l8grFKcZIM/p2IcRUN5ufCnS4x/cDjQXnzQAOHssCxwOR3hQdvclhOVA63jWbkUIgpIXlqeSN920DzaC1K87XfvsmTTVBvnjlmUwOeIbVSJWVc0qZotVUeKn0G8qrRqE4hRlpRr8euNX9+lbgkYLja4QQPiHEbGAesOnYlnhy0x1L03YkWTQVqhx57xobIfMmZUbnTv5vR4RPP7gNr6Hx6YtPd8odyzlaQnGDVRnNvmGyj+lVARXkFYpTnCEzeiHEWuCdQFgIsR+4E/hX4AEhxIeAfcD1AFLKPwkhHgBeB0zgdimlVfKJJwBHEhn2dyeGfb4ZCBOrX0ao3bn2ScBG599fEjyz401m1Qb5/GULqZ/kL+9oWYJSmr0VqqM2pAaEKBSKYQR6KeVNZe66oMz5dwN3H8uixgO9yQytXfGje5A77HvBfasQOFn3f2Wu4ZkdUa5fPoObV83MlTsezYZrVrPPXhRqG6aj65oqnVQoFIAyNRsRPckM+zrjIyphzw77DrZt4j/MG/iudTUXnlHPLefORCuQZ/oH70E3XAs0e62ynnmVfjQV5BUKhYsK9EdJtC/FocPJY3gGSa/t5VOZT/KkvZKLF9bz0XfMxZsdC4jM6fIlN2LLITSqGqbToIK8QqHohwr0R0FvMnOMQR7SPRHOa/0YAJ8z7uXPzvkCTU+sIdTWTDy8GKl7CHVszenyw6mPFwImBzxMUaZkCoWiBCrQD5OUabHvaDV5ANvE172LVPV89h1Ocvt9OwH4ovFjbjV+Q/zpQwQj2xyDs+hLgNOMMNxGKCFgTl2IoFf9KhUKRWlUdBgGsZTJgcOJYZVQFmGbnPGzJeiZXlr1Jj6R+DcA3q69wq3GbxxDs8g24rULCXa+isCpxJFCH3YjVGNNUAV5hUIxKCpCDEHW1mAkQb7yrd+i9+7JYQAAEzNJREFUZ3pJSB93xP8Kr27zgPFFlmk7CwL6ClouXcucDe8hGH2FWP0yWi/4LmagxBASF10TzKgJ4DM0fMYp5xmnUCiOEhXoByFj2Ucf5KWNEWtn3i/fjZ7pZbc9lb/K/D/2yil8T/wnS7NBHo03b/wjZmiK42Z55UPD3nidMtmvmqAUCsWwUYG+DBnLZm9nbGjvmkJyFgebAEmbrOHm9Odop4bveb/JxVpzzoMebPTUEcyKqc5js52uQ+AxhBoSolAojgo1YaoEpmWzq6OPRPro9BojESHU3oxAco95Geelvk03Fdz17qm8faqdc3eTgOWpJFU9/6ieXwiYWaM85BUKxdGhAn0/snLNsK2Gs0ibxqduB2nxrLWYr5g3A/CzKb9gyelzaL3gO4DmZvOCndc97Yx5GgZCwKSAwfSqAAGv0uQVCsXRoaSbAqSUtERjQw8NKYGR6CTUsZXfWsu5LfP/APiR56uc6fPQgsT01SCFAAkSgemvHuIZHYSA0+or8ChLA4VCMUJURu9iWjbb23pJm0cf5J0h3pJXqt/N32Y+xmyjkyd9n+Vd+kuEOrZiJDrxHd6dc6wU2Mx99JohJ0N5DY15DRX4PboK8gqFYsSoQI/T8bq9rRfTkkfvX+NuwLbcewfva7sBj7+CL7z3ImZMnVI0GCRVPR/LCOU2YwNdb2AkOss+bcCrMbcupMonFQrFMaOkG6Arlh7pjG2MRIRHD4T4bObDLBZ7+JvLLiBc6aP1Xd8GIfL18ELwxvteYe76qwl0vTFoQ1TIpzOzVm26KhSK0eGUD/SH42l6EubIHixtIo/dxRcyf8lfaNv4r8bn2N9w80Af+eyERd3D7qsfHbRevsJvMKs2qDxrFArFqHFKSzexlElr1/AHh/Rn9779fKJ9NTNEhG96vkPbu76JkexySiwLfOSB/KxYRH4yVD+ylTUqyCsUitHklA30XbE0LdHY8B+QDdSuxrOt9TCffXI/Qd3mR55/Y7Loo/HZOzC9k9zOV5DSxvROLjvqr5CqoIeZtSG8xin7K1EoFGPEKRlVOvtSHOhODF+Xt03mrL82F6hf2B3hi4/+ifpKP/92/dnM1Dsdx8m2ZoJtWwZU1ziNVFsGZvk4na4zqgM01gTH4q0qFArFqRfoM5ZNe09q+A+QNnMefQ/ByFaENHnhoMW/Pvkmp9VX8JVrFzG5dhqx+uVu16tFw9avDaiuQQpiDcuLqnDACfKn1VVQHfKO/htVKBQKl1Mq0Gcsm5ZoDMsefomNkYg4QR7Yap/GxzMfZ25dBf985VlU+j25ObAI3cnqIy+x87qnSNQuygf2YB0tl93P9ps20nLZAyAEhi6YXhXA0E+pX4FCoTgBnDJVN2nTCfJH3RAlnY3RXfY0PpT+NDWVQe684swiK4LsHNjcfNeKqU51TTyS33QtMC0TAmbVhpSdgUKhOC6cEoHedP1rRtL1agbr2FC5hs9FLsDQ4ItXLSrhHimL6+aRGPEIjc/cUVxmKTR8Ho2mmiB+jwryCoXi+HBMgV4I8bfAh3GKTF4FPgAE/3979x4c9XUdcPx7VqvXLgi0egPiVdu87BCweNhtM26dOBn8wG7jmEw8xY4nNDNxXzOehmnaYqfT1EnTzjjNk9T1kNSlpo5jMIYEm+A6jccEjAEbBAaDjUCAhMTDaBeJ3T394/cTSGIl7UPa1f50PjM7u6vfXe3Zo9XR1d177w94DpgKfAB8TlXPZhRlBuLx9PevAdjVdI7HW++htNDHvzwwj9rxpb0bdG9N3F3Ql6xl2qbPEzy1E4j1Oi1goKLOhmuMMVmXdsURkYnAnwMNqnojUAAsA1YCW1X1emCrez9nTpyLcCnNIv/bo238w8b9TKkI8MMHG64W+R5TLf2Rtl4zaorPHnbuE3OauqcFlLFV1JcHrMgbY7Iu06rjB0pFxI/Tk28GlgJr3ONrgHszfI60tX7Uybnw5bQe+5vDZ/jG5gNMrQzyjftuujozps+c+GhJqNeMms7yG3rcX8iBz2/nw7vXMSkUxGdbGhhjciDtoRtVPSEi3waOARFgi6puEZEaVT3ptjkpItWJHi8iK4AVAJMnT043jH61d3Rx6vyltB67/+QFvvXLA1xX7cyuCRZfTVPfHrz/UjtH73yu17YGPe8HSvzUlwdsIZQxJmcyGbopx+m9TwMmAEEReTDZx6vqalVtUNWGqqrBT6GXinPhLk6cTW9rg6b2i/zjxn1UjS3m7+6c3avIA0RLK6+dE989o6bPDJtgiZ8pISvyxpjcyuTD2E8CR1W1FUBEXgBuBU6LSJ3bm68DWoYgzqTE4sqpC5dov9iV1uNbLkR44rlfUxSL8Uzli8RKV1/bqE+Pvb8TeZcHC23fGmPMiJBJV/MYsFhEAuJUs9uBRmADsNxtsxxYn1mIyYnFlaNnLqZd5M9HLrNq/V7CMR8/KXqSWe2/6n+/+L49+D5qx5VYkTfGjBiZjNFvF5HngV1AFHgbWA2MAdaJyCM4fwzuH4pABxKPKx+2daR8Mu9u4a4oj294l9YLEX5a9G1m+Y7RUb0Q57x/zhZlg/Xgu5UHC6kaW5xWHMYYMxwymkevqquAVX2+3InTu88KVaXpbJiOzlhajz9zsZOv/NcuLl2O8ePCp1joO4jiQ+JRZq5dTEf1fBC5ZuFTXz4fjCstpG5caYJnMcaY3Mn7lbHtHemfOOR85DKPb9hHuCvGY5+6gUWH/OhpP+HKjxFo3YMQI3h6J4ggGruy8Kl7K4OeplQEGVOc9+k0xnhQ3lemFPYn6yUWV77/2mGazob5+j1zmDe5nKM3uB+yFo1j9n9+HF/0IjF/gEsVswm27HJm2ZSE8Idbew3jFPl9VuSNMSPWqKxOsbjy1Nb3eOP9NpbfMpV5k8udA+IjWlrB9A1/hC96EQEKYhGa/vB7zrGSENM2LbtmGGfC+JKcvh5jjBnIqCv0qsp3tx1i28FWvrBoMp+9eVKv4/5Imzts42zgE66cSzRQ7WwtHO59ApGiznZqJ0x2tis2xpgRatSt5Fm/u5lXG1t4YEE9yxZcuyI3WlpJR20DSgHhqvkcufuFK0M0PRdLddY1MH3qVMYFrMgbY0a2UdWjf3lvM0//5ii3/k4FX1joFnmN9546OdCCKPdYTcFHVNfVDzrV0hhjRoJRU+g37m3mR68fYdG0EI/dMcNZzNS9xfCpnYSr5nLkrv/B33mOaGllwpk1AGNKi6gMWZE3xuSPUVHotx1s4UevH2Hx9BB//emZFLpbBfsjbQRP7USIEWjdxeyfzsUXu9TvfPnqsmKqxxbbildjTF7x/Bj9rmNneWrrIT42cVyvIg8QLQkRrph95UTevujFq7tS9tn+YEyJn5qyEivyxpi84+lCf+j0R/zT5kYmhwL8zZJZvYo8GmfapmUE2vYRLwjSPR1f4equlK7xgULqxtkUSmNMfvLs0E3zuQhPbNxPWUkhj98955rthq/sK08cXywCUoBoDMVH023fvTIGX1NWTHWZFXljTP7yZI/+bLiLx1/aR1yVJ+6ZQ6j77FA9XJ0qWUBHTQMdVfPdXn2c+m2PIsSpD5VakTfG5D3PFfpwV5Svv7Sf9o4uVt01h0nlgX5aujtTuheJO6ccdE7mvZ0KzjM+cO0fCGOMyTeeGrq5HIvz5OYDHDlzkb+9czYzasf229YfaSV4+rdOYW/ZAQg9P2YN2UIoY4xHeKZHH1flO1sP8XbTOR79g+tYMDU08AO09+yZSGgOzs7zjuL1X4J4evvbG2PMSOKZQr/mjQ947b1WHlw8hU/Nrh20fTRQRUfNQhQfscKxlLY3Eqm40flQFqBpO4TPDHvcxhgz3DxR6F/cfYIX3j7Bkpvq+FyfTcr6JcLRu9Zx6L5fUBCNIMQobW9EJs4Hnx/qF0FwaE9abowxuZD3Y/Rb9p3i6f87yi3TK1jx+9MHX9DUa28bH52hGYRrbyZw+i2kfhEsfwnCbU6Rt8VRxhgPyOtC/9aH7azasI85E8p47I4ZFPgGL/LTXn6g137yY0qLCHxpM9KzuI+pzs4LMMaYLMjrQj+ztow/nj+Jez8+kSL/4KNQVxZJudscBC63URMMIBKw4m6M8ay8HqMPFvt57NMzGFOS3N+rK4ukKCBceRPTfvUVAv92IzyzxGbYGGM8K68LfcpEOLpkLeGquQRa9+A7/iZoDI69AR0tuY7OGGOGRUaFXkTGi8jzInJARBpF5BYRCYnIKyJyyL0uH6pgM6Jx/OFW/JfaCbTtRejZgxf3Yowx3pNpj/4p4BeqOhOYCzQCK4Gtqno9sNW9nxtucSceY9rLDzBz7SKu+99HkUkLQQqguAx8BTDlVhujN8Z4VtofxopIGfAJ4CEAVe0CukRkKXCb22wN8Brw1UyCTEvPs0dVzCbQvh/RGIUnd8Jf7XNOKhKosKmUxhjPy6RHPx1oBZ4RkbdF5N9FJAjUqOpJAPc6J13lXmePansHLSy9uhBqTI3Tg/cVONdW5I0xHpZJofcD84EfqOo8oIMUhmlEZIWI7BSRna2trRmEkVi0tJJw1dyrZ4+6HIE//TU89LIVdmPMqJJJoT8OHFfV7e7953EK/2kRqQNwrxNOZ1HV1araoKoNVVXDsNWACEfu/hldtTc7PfnJi6F6lhV5Y8yok/YYvaqeEpEmEZmhqgeB24H97mU58KR7vX5IIk1DqKyU4hWvOpuT2Ti8MWaUynRl7J8Bz4pIEXAEeBjnv4R1IvIIcAy4P8PnSI27l00gVMuEcSW2pYExZtTLqNCr6m6gIcGh2zP5vmmLR5n+0mcJnNkL9YuQhzZaL94YM+p5Z2Wsxp0i37oL0Shy3PaTN8YY8FCh90faCLTuubq+dcI820/eGGPwQqF3V79GSyqI1bsrXictgC9usWEbY4whz7cpJh5n3Lr7qDyxg8sTF+B/eCNE2m2GjTHG9JDfhT58hsLmHYhGKWre4RR5m2FjjDG95PfQTbCK2MQFqJ3j1Rhj+pXfPXoR/A9vsgVRxhgzgPwu9AA+nw3XGGPMAPJ76MYYY8ygrNAbY4zHWaE3xhiPs0JvjDEeZ4XeGGM8zgq9McZ4nBV6Y4zxOCv0xhjjcaKquY4BEWkFPhzmp6kERuoG9SM1NosrNRZXaiyu1PWNbYqqDrr3y4go9NkgIjtVNdHZsHJupMZmcaXG4kqNxZW6dGOzoRtjjPE4K/TGGONxo6nQr851AAMYqbFZXKmxuFJjcaUurdhGzRi9McaMVqOpR2+MMaOSFXpjjPE4zxZ6EflnETkgIntF5OciMr6fdp8RkYMiclhEVmYptvtFZJ+IxEWk36lSIvKBiLwjIrtFZOcIiiurORORkIi8IiKH3OvyftplJV+DvX5xfMc9vldE5g9XLCnGdZuInHfzs1tE/j4LMf2HiLSIyLv9HM9JrpKMLRf5qheRbSLS6P4u/kWCNqnnTFU9eQHuAPzu7W8C30zQpgB4H5gOFAF7gNlZiG0WMAN4DWgYoN0HQGUWczZoXLnIGfAtYKV7e2Win2W28pXM6weWAJsBARYD27Pws0smrtuAjdl6P7nP+QlgPvBuP8eznqsUYstFvuqA+e7tscB7Q/H+8myPXlW3qGrUvfsmMClBs4XAYVU9oqpdwH8DS7MQW6OqHhzu50lVknHlImdLgTXu7TXAvcP8fANJ5vUvBX6ijjeB8SJSNwLiyjpVfR1oH6BJLnKVbGxZp6onVXWXe/sjoBGY2KdZyjnzbKHv44s4fwH7mgg09bh/nGuTmksKbBGRt0RkRa6DceUiZzWqehKcXwSgv5MEZyNfybz+XOQo2ee8RUT2iMhmEZkzzDElY6T/DuYsXyIyFZgHbO9zKOWc5fXJwUXkVaA2waGvqep6t83XgCjwbKJvkeBrQzLfNJnYkvC7qtosItXAKyJywO2F5DKuYcnZQHGl8G2GPF8JJPP6h+19NYBknnMXzt4oF0VkCfAicP0wxzWYXOQqWTnLl4iMAX4G/KWqXuh7OMFDBsxZXhd6Vf3kQMdFZDlwF3C7uoNbfRwH6nvcnwQ0ZyO2JL9Hs3vdIiI/x/n3PKPCNQRxDUvOBopLRE6LSJ2qnnT/RW3p53sMeb4SSOb1D9v7KpO4ehYMVd0kIt8XkUpVzeUGXrnIVVJylS8RKcQp8s+q6gsJmqScM88O3YjIZ4CvAveoarifZjuA60VkmogUAcuADdmKcSAiEhSRsd23cT5cTjg7IMtykbMNwHL39nLgmv88spivZF7/BuBP3NkRi4Hz3UNPw2jQuESkVkTEvb0Q5/e/bZjjGkwucpWUXOTLfb6ngUZV/dd+mqWes2x+opzNC3AYZxxrt3v5ofv1CcCmHu2W4Hyy/T7O8EU2YrsP569yJ3Aa+GXf2HBmT+xxL/uyEVsyceUiZ0AFsBU45F6HcpmvRK8f+DLwZfe2AN9zj7/DADOrshzXo25u9uBMULg1CzGtBU4Cl9331iMjIVdJxpaLfP0ezjDM3h61a0mmObMtEIwxxuM8O3RjjDHGYYXeGGM8zgq9McZ4nBV6Y4zxOCv0xhjjcVbojTHG46zQG2OMx/0/vpUW0Xb34pwAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# quadratic model\n",
"\n",
"plt.fill_between(d.weight_s, d.ub_qua, d.lb_qua, alpha = .2)\n",
"plt.scatter(d.weight_s, d.height, 5)\n",
"plt.plot(d.weight_s, d.muhat_qua)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fc33a4b0950>]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydeXwc1ZW2n1tdrd4k25JakjfJtsBG2IDBK8uQAGG1sdkSMAkMJCTM5COTmcyXScJkm8yETCaZyUoyCZkQIN/YZg2YnUAgJAFsy8YLYIONZWTZ2nf1XlX3+6OqV3VLsrzKus/vl1hdXd19W8anTr33nPcIKSUKhUKhOHHRjvUCFAqFQnFkUYFeoVAoTnBUoFcoFIoTHBXoFQqF4gRHBXqFQqE4wdGP9QIAgsGgnDlz5rFehkKhUIwpNm3a1CGlrBjuvOMi0M+cOZP6+vpjvQyFQqEYUwghPhjJeUq6USgUihMcFegVCoXiBEcFeoVCoTjBUYFeoVAoTnBUoFcoFIoTHBXoFQqF4gRHBXqFQqE4wVGBXqFQKE5wVKBXKBSKY4BhWrT1R4/KZx0XnbEKhUIxXjAtScdAjPb+GC5NUFniPeKfqQK9QqFQHAUsS9IRsgO8ZR3dz1aBXqFQKI4gliXpDMVp749hWsdmdKsK9AqFQnEEOB4CfBIV6BUKheIwYlqSzlCMjv74MQ/wSYatuhFC3CuEaBNCvJVx7EwhxBtCiC1CiHohxJKM5+4UQuwWQrwrhLjsSC1coVAojidMS9LWF2VnSx+tvcc+i89kJOWV9wGX5xz7HvAtKeWZwDecxwgh5gKrgHnOa34uhHAdttUqFArFcYZhWrQmA3zf0d9oHQnDBnop5atAV+5hYILz80TggPPzVcBaKWVMStkA7AaWoFAoFMcSy4KBNpCHL8s2TIuW3ig7W/ppO04DfJLRavT/ADwvhPhP7IvFuc7xacAbGec1OccGIYS4HbgdoKamZpTLUCgUimGwLLj/Sti3HqqXwi1PgTb6XtFw3KAnnKA7HD+ug3smo/22nwW+IKWsBr4A/No5LvKcm/cSKqW8R0q5SEq5qKJi2JGHCoVCMTrCHXaQtwz7z3DHQb+FZUna+qO819rP+20hOgfGTpCH0Wf0twB/7/z8MPA/zs9NQHXGedNJyzoKhUJx9AlU2Jl8MqMPjCyxjBkmvZEE4ZhJ3LSIJcZQZM9htIH+APBh4BXgImCXc3wdsFoI8QNgKjAb2HCIa1QoFIrRI4Qt14Q77CAv8gkPNjHDpCecoC+SIHoUAns0YdIxECNY7DminzNsoBdCrAEuAIJCiCbgm8BngB8LIXQgiqO1SynfFkI8BLwDGMAdUkrzCK1doVAoRoamQXFlwadtg7EYXaH44dyvLYhpSV7c0crqDY2ce1I5d398wRH9vGEDvZTyxgJPLSxw/l3AXYeyKIVCoTjShOMGAzEDn9vFvq7IUal7l1KyqbGb+/6ylw+6wpw6pYRPnjfziH+u6oxVKBTjhphh0h81iMRt/f1oZO9J9rQP8JvX9rJlXw9TJnr5yuV1fGhOkLlTJx7xz1aBXqFQHB9Y1oh09IPFMC1MKYnGLdoHokTiR3dTtaEjxGNvNvHHd9sp9uh85vxarjhtMm6XhjiM33MoVKBXKBTHnsNc694bSeDS7CC6vztC3Di6wd2Skvq93TyxdT/bmnrx6BrXnV7KdUvnUOx1H9W1gAr0CoXieCBfrfsQm6f5iBkmRS6Nlr4oHf3x4V8gLfRIJ4YveNjuICJxk5d2trJu6wGae6MEi4u49ZwZfKbpTqbufo1Q/0Ialj8I4ugO91OBXqFQHHtGUetuWpKBmEFfJEF/1MCSEq9bG1lZpLSY9fQNBFrrCVUtKhx8h7oYZDzX1h/jqe3NvPB2C6G4ySlVJdx89gzOqS3HG+tk6pbXENIg0LrJfo3/6DaJqkCvUChGTz5dfTRa+whr3U1L0toXxet20ROKEutpywrCI9Xf9Ug7gZb1CCDQsh490o7hr8o+KXUx2ESoKicTlxYzn7qBd5t7uUe/kRcjswE47+QgK+dPpW7yhNTbGL4goaqFqfcxfMGR/U4OIyrQKxSK0ZFPV4fRa+15at1NS5IwLdr7YwBEEqbdoTrSjLwQMjc7H3xh0SOdBFo3DcrEw3GDP73VwMuNV/G2nMmEeIhrTy/jigW1+ee/CkHD8gcPu0x0MKhAr1AoRkchD5lD0NotSzIQN+gcsId2RBMmQjDIV2ZEGXmSPPKL4SsnHDwLf8c2O8vOI6VkZuIDFQt4tz3CM+/v4tVd7UQTFnPcAf6V+7l8SojWD/3v0AFcaEddrslEBXqFQjE6CunqB6G1Synpixj0ROzN0/6o4RzPPCffCwtk5LlBPVd+WbYGPdJJ9cufw9+xjXDFfOduIE+QFoJXz/0Nr+9oZOPbO9nZ2IhPJDi/bhqXzZvCnEo/7uhf0XqMsvSDQQV6hUJx8CR1+FuehHBntq4+hNYupeRAbxTLkoTiBqaRwN25m1jpnJFJPBmBPFS1JK17+yvyaup6pJ1Aaz1CmgRa6ql98qP427cCJgLwt2/F0/Uupj+I4bPXmzAtXnu/k+feauatA30I4DQR4t/0e1mpr+fAOa9g+EsAjmmWfjCoQK9QKA6OfNp8ZkB3tHZpmcR729CKK+kMJRiIJYgbMm01YBmc+tv5uBL9mO4Sdty8FbQCIUla6OF2ql/+XEYgX4se7U5l73q4I1tTD7dT/Yc7QJpIIFwxH3/HNgT2Y4mGqfuY/bg9QG9r2RXcX/llXtrZRl/UYPIEL399zgwurqtkwcs/PqabqYeKCvQKheLgyNXmQ20gNKQ/SDRhYg60YRRNwvP/VuBt20q4aiHteTZLPd27cSX6EYAr0Y+nezex8rrBn5fM1FvqSWbigdZNdpDPyKhzq1sQgkDbZgQgcdH4kV9S/Ypzoag4i5YF/0jfM9/kNetUXjAXse3ASbiaD3B2TTGXnlHDmTWlaI78s+/Cu0GIVNY/1lCBXqFQHByBCuT0pYim9cSnLEI8eCv6/o1EqhZgIQi01GPqXlxGyJZHCtSOx0rnYLpLUhl9bNLJ6OH2QZUpqeqXZCYuXPkzayFoWLbGvmCUzgEhsssaA5U0LH8Qo7eVN5+9hwefaORd+S8AnCYa+HLpH1lW/C41bX8mtH0hDTUPgpSDSyzzzlcaHT730RmprQK9QqHIj6PDJ7zlaJpGdziO1+0iEjdpv2w1RbEuEoZJ3ZqzEdLA17IJBAjMVJCXQDh4Rn65Q9PYcfNWOzBPOplZz944uGZdWiBlOmBXLmDfRXcXbG6a9Uz2eyTLGhuifn7/2l52NvexryvEQPwjzBMN/If+K2pX/BPFpZcBUON8l2Q5JZC3xPJQ8XtcVE3wUuw5OiFYBXqFQpGFZUnC8QT6AyvwNNcTr1pI8zWPEElIsAw83bsxS+cQKSoH3SQcPAN/+1bCwTOQriICrfWYbj+uRJhwxXz2rHissNyh6cTK69DD7YMDqq88nU1XLmDnja9jeIPMemZV3iamrI3X1k2IUAfPfWDxzPZm9nSEcGmCGZOKWFBTxqcH/pu/6nmS8OSFNEydhyFE9gUl447hcDY7DQrwR8jILRcV6BWKE50hgkmyISmaMAnHTSwpCcVMrL5m6po3IqSJv6WeRE8r+MrSm6d6gB03bWHWczfhb9+Cpfvxd2whXDaPnTf8GTSXXTvuLUePtAPSflxA4za8ZfYFI1nX7guiRzoItNQjsAM3aOjRrvwZtrSofukOtpk1rDPP5Y/aUvY+sAvDkpRrIb6kP8W1no1UhVsJB+az55qHeTf2NQxvGXqkIyUX5WtsOhzNTgGPi8rcDP4wG7kNhQr0CsWJQKFg7gQT6QQT8+Z1RM30CLvkKVlDN6TFrJfS1SpgUv3yHRw451vpzVMjxEmPr8DXuwuBhWYM2Hp811vUPXierY9X2rOJAq3paaKhyUvtLFxaGVo6zHpmlX1XUDGfhitW28HXU+rcGfRjuv0Y3jJ7NS4vLmMA0+WlhwlsfLeNN/c009C4ir1yMgCnBX1cNbWUCz74CctDj6EJkCZOSeVmap/6GHtWPGLfHbTUO3cej4LmGizNjLLZSQiY5HcTLPbgzafFHwYjt5GiAr1CMdaxLLhveTozvPVpLAQxwyLR20LJvvUIy0B+8Aa79jRg+IcOJnqkM6NahVSVi+kpx9QDKf3d17MTnHMszYdmRZzXWAgJgdZ6IHvrMtBajx5qZfajF6fuDHZ99A+pzVZ/+1Zqn/qY3cwUPANXIuRU5YTRo13EEgavx2bSIKfwfHwxb/xmI6YEr65xpifOp8z7Oa8qQXjlr6l98jr84c123xRgufxoZtgO9h3b8HTvTt0x+Ns3U/vkdexZ+dghO0sW6RplgSLKAkUpq+S8jHJo+WgYyczYe4ErgTYp5WkZx/8O+Bz2bNinpZRfco7fCdwGmMDnpZTPH4mFKxTjgpFouKE2ZONrdpBtfI2Olkb63UFCMRNBgJnBswi0bQRMqv9wx7C+MKkyxZYMrT14BoavnF3XvcSMFz+Dr/MdkqWOEtCsCDg/239qmLoflzFAZmNrqHIhrlhP1p3BjBduI1S5gEDbJsJBp9ZdGilJqCch+ItYwJOvd7OlsYv2xD8DUCsO8NE6L4vnzmF2VYCiyEkgLibsq0CPtOFv35LeEC47jT0rn6D26Y+l5KFY6Ry7tr59cyr4H8pma7FXp7y4iBKPPrKBIgcxtPxQGUlGfx9wN/BA8oAQ4kLgKuAMKWVMCFHpHJ8LrALmAVOBF4UQc9SAcMW4ZzSbbgU0XMuSRA0T05KE4yYDHSFqM17WMZDA8Nv/5KSUCOefX6r+PDeY5doGZGrVRRM46clr8bdv5dT/Nx+XEbE3Rm/4E3MeuQTNDKXeW2b8CRYuR85JH9fYd9HPMXzB7DuDrh2Eg2fYXgeWxftlH6K5vZPd/rN4pq+GjZZTW7+jjXlTJvAV/xrO7H2ZGZ5+XHsjhML2YO1UVc5H7qb6D38HWOkLjzsAmqDx4l9m1cPvWfGonfln7A2MCOd3ZvqDlBV7KAsU5ZdnhmOYoeWHi5EMB39VCDEz5/Bnge9KKWPOOW3O8auAtc7xBiHEbmAJ8PphW7FCMdYY7aZbhoYr962nu30/3VopRS6N/qiR1tXd5YQmL005ORq+dBDXI+3ZmW1uqaMZ56R11+Dr2uEEyZ9h+CqdjdNyatddi69ze6qpKXmxqHnxdjQzlJJlMoN85jF7hRoSC9MdwPCV22WVN23h5CdW4u1+l/bS+bzTavEr8yts3z+LPgL2G8RhumjjC/rDnFau4V/xPXs6kzwNs+s9XI8vczZl6wHLXlvbRurWnJN6nJaeHPsDJ6Cn6uE1F3tWPnZQm60uIZn1zCq8LfUwfSni1iO3iXq4GK1GPwc4XwhxFxAFviil3AhMA97IOK/JOTYIIcTtwO0ANTU1o1yGQjEGGMWmW3coTsIsYdLUxbj3byRcuZD98RIQJmFybpAL2eA6lSjJzDZcsSC71DFpQeBk1oHWDdStOTtl+6tHOvG3b00FTEsvRjOjhMtPx9/xZkb2DiCQyKxjluZHs8Ikg64rEUaGOtkX97P5mV9hDCxhm/w425tnYDihKECEO/U1lJ3/Nyx57TNU04IQgp2Xvw5WD4YMgtCIlZ2SKntMrie9EittceAO4DIiqYqevPXwI9hs9RVpVJR48bldFEU7oKXe/vtsOrKbqIeL0QZ6HSgFzgYWAw8JIWrJ3zKWd866lPIe4B6ARYsWHcVZ7ArFUWaITbekDJMcfefWBf2ROGaf3SHaetka9HB7VvDOdWdMPk4FK+cYyIxNVUHjRb/Iyjw93buzGpsAx/xrI3rIvkm3NewthMvnsefKx/D0NWB6JlG3dknqfcITTmbPynWc/PR1eLvfJVx2OkLG8HXvBAk7ZA07ZA2bPGfz7MMNdIcT2GEDzhK7uN31NGdrb1Oj91JtNaEJSWhPK0yuRrQ2A5LZj11sy0YZtfOpi5u3nFnP3ECgZWMqsNt3Jz+3yzujXelzDrIePlk5M22SL627H8VN1MOFkHk9QHNOsqWbp5KbsUKI57Clm1ecx+9jB/1PA0gp/905/jzwL1LKIaWbRYsWyfr6+lF/CYXiuCdDozcl9EUSdIfj9oapyLDizWOrm+r2rFxg+7fkey4lR5DVZISUzkYsg43DTIN5vzkZgW32LtFSP5uaH5cVtpeuF6MZYWdz1Q6iwjJysnr7AtJbfhbb4lNo6+nlDWsuz5pLCOEDwI1BdXAiZ1VP5Lwd/8ql8i/4RRyReScgXAhpIoXOrqueYvYTy+3HJKUhF7uueZZY2SnZMkvyguctswP7MKP/hpJohICAR2eSz02JV0d35ZFljlKj03AIITZJKRcNd95oM/rHgYuAV4QQc4AioANYB6wWQvwAezN2NrCh4LsoFCc4UtpujTHDomnABwMDJEyroN+6HulMNwm11NslgBkNQghSlrv+5o2Dmocgu2X/gwvvJtC2Mds4rGwOeqiVop73EWTehlupn11WOPVzskY+ubkaaNvIzhveYNpzt7K3O8YGq45vGzfbJ+/P/v6TtV5Wite5TNvAAtce9q/8IwB1O19FYDghPk2y6iZUtZBYWR2hqkVZ1T+m28/sx5cNHu2XZCgZpsBzyZucEo8d2AsG99wXHedyTSYjKa9cA1wABIUQTcA3gXuBe4UQbwFx4BZp3xq8LYR4CHgHu+zyDlVxoxhPxAwTgaAnHCccNzGcKUkjuHEG7A7RVJOQ7sP0TMzyeUGIVOCrffbj9rmOpGH4giCtlCWB6fYz4+U7kMIF0kwZh8168mOpLB/Sso2pF+MyBrLWY9ef+xBGhP0EaZHlvGKewc/u3wN8Y9D6J2hRPj/pdWr7N7GgSmPfioeZ9ex9BFrfyW8rkLrr2ESo0rYe9vS8n/KnT8szZXi6dzE7tQGb3RWbdLfMbHwaCULASRXFo6uYGUOMSLo50ijpRjFWSZgWlpT0Rw1KvDqNnWGiiZENqM4nJejhdurWLEnJFaARDs6n8ZJfOqPyJJ6u91IBL0vKQDLrqesJtNYTLpuLv3uH8z4aey7/X+KlJ4MQ1K1ekpJowNbvd139nJPptzDjuVvQehr4vbWQF82FdE75MI0H9tMo06P6Sjw6/TF7GtQP3P/NnKAPcem3cMf7iE08CU/vnvQwkXyyCjJ1zPauqbc7aTOlqdyMXUpmPX19xvMPOT707c53MtObzgUan4SACV47cxcCfEUuPPrYDfJHWrpRKMYlpiXpjyYIx01ihsWAM/oOoPlg3ijPNCTbG6bM1sIT/YBdQeLveJOaF/82FbxSFSctGwmXzyU2abYd8EJtBFo3pKwIwhVn4e/YTqhyAVVbfpTyYTc1b0p/B+iqPJsNkck0vdXCL15tYJa4nQY5Nb3W/QaniBifdj1Fv3c6iz+0jMt3fDn1fgfO/Ramp5TZj12KK9GPFC4EdnNUslwzZVCWzLqvfBgAPZocwG1mS1P56v0LVBcZvmDBxidNA13T8Be5mOB14/e4cA8ny5yAqIxeMb4Z4aZaJG7SE4nTHzWI5cvYR7jRl8TO3JfaWbnQ2Xnjegx/xaDsFNKbkDtXvW57sXjL0EMtzHn0EjRjwDYYu3kbeqybutWLUhubO2/cAMIFyJSVsMTFDmsaG6w6/sW4dcg1Xu96mS/rD1Im+lINT+GKM7NG8eVGj0GNU8JFqGoR+y78CXVrzkmVPiZLNbM2mHM3m52MfURYZlbjU8Pyhyhyu6gs8TDB5x7aimAMozJ6hWI4hmlkihkm/VGDzoE4cWMIOaZAdj5U8M+dhpTUrg1f0LEDsDV0yxVwuk9NZj/6EXtDUvOkMvKUwdi6q3n/qnXZjVP+KjpDcfZ3h3jC/VkaBnR2abPYYwwuLbzAtZ1V2otMEgOcJd6jSKS7aZPDPsLlZ+Dv2Joxio/UOSTPy/gZ0hvHNb//W5L1/IL0Bm+gbbN9AXPuZvRIR1bn6ojRXDRd+zgVWh/6hCpOcesU6eMvcy+ECvSK8UtOI5PR30aoqBywG5b6M2SZQWTozp7uXUN7qefTmws1OZFtWaCZkVTJYcofxqmISQZW20bgHayBDn53xi95cfs+3myO0Pezv2Qs+ByqRRtzrT3M1vZQKXpY6XqNReK9lOlXbjaeKf3su+huu/nK2TswXX5cZjSVoQOEJ83D37sztS8QrjgzZU5mXyCy7wBSk6L8lcDBTXJKSjJCwESnDNKlCTz6xMJ/Z+MYFegV45dABdb0JYimDcSnLGb3gBdLhod/XUaVR2bZX2b1ix7JGFTdUo8e7sAI5JTjJcv9pGU/7wsO6kYNB+cj9aKMz7K1e0sKdsgaHjY/jIGL38ultD6wK+vtL9I24yVOndbIR12vMkV05fGkGfw4SeNHfpmSijzdu1LNVwAuM2JfCNq3YGleNCuMv+dtezRgsrFp2Vp789VTavvkOGuX2BU+uz72B3uDOd9g7wLmYkW6xiS/m8oSD1KCdoJKMocbFegV45Jw3KCtL0b/JavTWbUsEDRyJJjMGabp4dZhu/qldLbtpe4ttyWY1g0k/dzz1n0PapBaS2jyooxSwcdAgNHXzpudgqZ9+7h3eyTrLTQsFs8sJxgzuHiayfmNv+Cs3udxiXToLhTU03cGdnWP1Nx2qWPVQoyAXeWTuqjpvlT5ZWjyEhqWrbVLHn93RfoCkAix65rnUg1NyX0HVyKcdZFxmVF7/yBjQzWflKW7BOWBIiRQ7NEJZAzuGIMzuo8ZKtArxg1SSvqiBp0DMUIxp71jOJ+TPPr7IBtfJ4ONlc7OGnO376KfUrf23PxVJNKyA2CkIz3+zrEe2Hfh3WxribC1Q0B9E2/t72Xb/t5BS1vhep0iGefbRQ+w94JXbI/3t5xqHVF4o1Tm/ByqWmLbBfgqSJU9Ji9qyUwbE5cRYdfVz2H6gykNPe05Y/dFhiYvHtS1aviCzsVrI1IIu6Zf96UGiSQXnJSy3BMrKfO48egaXrfLnsqU3DQvcrT746QzdaygAr3ihGcgZtAditMbSRRuXCqwcZrK3nMkhcxGnmRtuB5uyZpZCnbFSVaW6gT46pf+T2rD9X1rKm/JmWy1TuKd3z7FeuvUrKXNLPdzZd0kanf/hmXaG9SKZszSk/E6gz9CVUuyPN6TAdwUXmLBU+3mKd2X5WuTaSmw76Kf2zp5nt9BbqYdKz81J7Amf6H2HUHDsrWDA68TxD1d76Wy/+QgETNQ4dSya1SUeIgmipnoc2e/3jTgN5fDgTftTfO/XgcPrDwqI/hOFFSgV5yQxA2LnnCc7nBi6IoZsLP2pzIaca58KCWxpANdfX4zrOQdgWUw++GPpDYrQ5ULMHzl7Lvw7nQViSODuJq3sMk6mX9I/Ddx9JQtr5cYp4pG/EQJ4+Vr+m+5wbuevddvQA93MKfpWTRHOnH37ExLIWac2KTZWf40AnDJKP72N53KnAiR8tPxdb5tT4ByXmuP6CtFD7VS/fLnBm8eF9w0tn9vnq73HO3ewt+53b7o5btDcur/w1UL8bduIj51MbNrZ6HnNCsNal6yLDvINzmdvPvWQ8d7R20E34mCCvSKE4K4YRGJm0QNk1DMSEszI0APt6cajQKtG9DD7Y4+DSBtMxqJY0ojnQ7N7OYff/OGrEEbB5b+S0rG6Q+exTMnf4PXW1y8+cEqGuXfZ33+SWI/N7l+zydcL6ELmepaFYA04KQnrsbX9Vb6WM6f/o4t6JFWYhNq8GdeAMiQZyYvpmHZWvzNG6h99oYMTX3AHt2XURs/Ihvf3A3pTBuGjHP0aCfuCVUESzxoQuL36AgNPG4X5G6k5pNjwh2wf3P6nKlnQUXdmHOPPNaoQK8Y08QNi9a+KD3hxOjfJI/UkCQ9P9Uk0LY5bRyWMWt07v2nO77r6etB00s/48fdS+mQl7BzXw29+/qdd6ziU67nOEd/j7PZQrGI2h+JbUVgNyS9SSa+rrezatUz/0xS1PMB/p53C1TTuNh34c9AcxGesjTL0yYcPMv2aU/WxidLHvPZ+GZIO9kb0vZGNFWnolmSyRO8aEgmPHgN2v4N6eEc4Q5o2gCWOTgTL9TTEKiAmrOh8Q2YtgA+9YJ9/CiN4DtRUIFeMSaJGxZdoTidoRhWUpk5yO7UdC18sOCEpiyNunIBIDG82S33CTPOg+alJND5tnET5fTSGbPruafSzmWujSwS77Hc9Qb7r3oUV2I54Yqv4nrqOuh6G8vlswd6BOfj79ienYkHF6ZsEFLfKGnlS1J+KSE8ZSmhyUsItKzH0n1oRtQZhu1k2smMXNPY8dfb8HS9h+krx/BVpP1jKhew7+L/zt+slLMpvW/FgykZxpi2mNJZ8ykr9qQ7UAfaYP+G7OEcQ/m4FxrOUmiu6hhzjzzWqECvGFPEDJP93ZHB0kyh7tScc1IXgtwGnWVr0KPdgy8SSY063E71H+6gbvVSQlULeeeytfz54R9S31/KJms2UTwAeIgzW2viy9palrveICBiaXdIV4DaZ2+0/WAQKR92zQyDcCFd7oxqHh8uIwIunT3LVlP7lDMGL9P4K3gmrYu+SHjyUnC57A3Pzp224RnS7qhNlTBmDPnTdGLBuanfiT0nlozzsn9XQhMUmz3pTem2TZxSEkfc/iwi3Ik7UEFF7oUhX1Afahj2UBcBFdQPGRXoFcc9AzGDgahBOG4QjptZQzoGSQmFGm5yM9KL7h7s5T7EhUFKSUvzPl6zLuIv++bx0r2bMKzzmT7BzcT+NqbRwXe8v8X98Qepe+6/8HW+g+n2Iw3TNv4679sgJbMfv9wJt9lj94Q0U3YArkhn6rxA63r0WBd7Vv4u6yKVvPDUPveJ1Og/gKmvfZ3ModgFSzsH/d7M1CjBcOVChCbwtW7CmLoYcetTaKIEUWMHYlG9FFHiZNuFAvDBZuJDXQQUh4wK9IrjFsO0aO4toL/nazTK03CTJPdCgBRZksygihOg5qlVvN4sec59MX+2TqMj/l+ALcecNd3HwunFrDx1IqeuvdE2DJMuwi98Al/nO6lNWj3Wkw7OoXa7czTRn7U2W7qJEao4C4TA9JXnfC2hdpsAACAASURBVFeRsyFqZ+CpDeSW9Xi6dmL6ggTaNmUZBwwaCJ7ze2u88kHCkxfib6kHaSKkib/NcZG0TNwHNkK0yw7OBxuIDzYTV5n7EUMFesVxh2VJusNxmnujBeveB2Xw0a7CZYA4Az10X2qgh+ErT0kyrmgXsx9fjpAG7Qca2bxtO5s7XWxtvJ69cjLeeIwPaxs4X9/Oedp2prv6iLlOJ7B1M6HmBWnL4NLZKdMvf8c2O8jnDsYoP539f/Udpv3lqyndXTPtTld/2ybq1iwlVLEgdUEw3SUYuYEfBnXxzn58OaHKhY7Onx7xFy4/nT0rHkNoAk0IKkS/fTFwJJjaQJSiTz+LCHfAI5+0M/bpS+zfX66MogLxmEUFesVxg2FadIbidA7EMa2h7bMNb5ljlrUtncE7Lff50KNdqTZ8VyJsuyQC1S9/jgMHDvBZ8+/ZaU1lr5wCf+rHIwwWemL8X/NuLtPWUyQy5pZK0JPZdNtmdl7/J+Y8egn+7p32NKecypXM8YD+zi1Me+1r7FnxCLXrrsbfuT17lJ+0R/XheLq7jEje2nTDX0GoagmB1nrndY70c8Nr1Lz4N/g7thKtPJOW656gwluUMv3y6iUpLVxUL8U7abId1EuqsjN2KZWMcgKhAr3imBM3LDoGYnSF4iMbuSctZj2zCn/7VsIV82m4YrXtLzNEtU2qDd8ZlqE/90882DadV6zr2C2nA+Anyh2ux7nQtYXTtEYabvwTcAnxl+/A3bIRMuvbSZciuuJ9KdtdpMmey9cQL5ud9dmDBmPEethz1Tpqn7gGf+eWQevNnJ2at9RRCBqufAg93EHNy3fgb91EePJCiium03fT01hmDyVlU5lpWbginWnrAEgHdF8ptO2w69I1LTtjH0p/V4w5hh08IoS4F7gSaJNSnpbz3BeB7wMVUsoO59idwG2ACXxeSvn8cItQg0fGJ9GESXt/bGhrgjxkDe3IsMMtWG3j0NYb5s/vfMCGbW/zdmIyAHPEPj6h/4F5pTCv9w/4RCJjHN3vHNMY212y+qXPpmwLsvxhpOTU357uyELFRMvnpYzBUuvJMxgj9d6hVlzRbqa+9o2s16XG7uW5eHl1KDF70YorKA24cUe7sHxBtMzpSUP57ZsGfG8WxPrAMwG+1AAulfeNNQ7n4JH7gLuBB3I+oBq4BGjMODYXWAXMA6YCLwoh5qgB4YpMYoZJW19sdE1O0gIpM5whrVSmXMjetrk3wks72vjdlv3EDYszRIgv62u4XNtAZeVUGi/9FYY3iPnMDciWjYQrznRcI50AKzSMQCUNKx5GD7cPHowhBDtu3oqnezemt5S6teek7Ik9Xe/ZJl+aiz0rH0u9PoXQMIqnYBRPoWHFQ1l7DMnvkfRc1zSBrgkCRRqB1VcjmjICeHElgy5vhWrTwbYRiPU5fyF99uOquQf/96EYEwwb6KWUrwohZuZ56ofAl4AnMo5dBayVUsaABiHEbmAJ8PqhL1Ux1hmyigaGb3jKqBgJl5+e3TiUo4tHEyavvd/Bizva2L6/FwEsmFHKbefN5EOv/TjbBtjJcvNu5matScuwRshB04mV19kXoQxny9mPX5FV/pjXTyaJU1nj1gWBIp2JfjfRuElZoAg9M1MfaLObkIbzehmqNr2izs7kkxl9RV3+73WkUO6TR5VR3asJIVYC+6WUW0X2X9I04I2Mx03OsXzvcTtwO0BNTc1olqEYI0gpae+P0dYfG9I9criGp8xKG3/HNlv37thmd3R+5Od0MokdDV3U7+3i1V0dRBImUyZ6uensGVx0SiUVJXZTU8HqnFxPl5E0YeWSdGrs3JmuhW9Zjx5pB7SsSqGA0U3MG0QIKA948Be5MEyTYqMHrbgCoWlM8LoHf8ZQATxnLQVLIjXNlms63ktr9EeLYUY4Kg4/Bx3ohRB+4KvApfmeznMs7z9tKeU9wD1ga/QHuw7F2CAUM2jujRKJD63eFWx4ysiocy1zG5atxRXpZGtPEU+/2sLre3ZhWhKPrnHeyUEuObWKeVMnIHIzxuE86Idb03AIDTN3A1UKDH/QqVnfRGLqYmbWzERzaUgp7TVaFqwegf3uwTQXDVUS6dKPjVwzlKSkOCKMJqM/CZgFJLP56cBmIcQS7Ay+OuPc6cCBQ12kYuxhWpKeYWrhMxk0Ychbhj7QTM1Ln83I2n+WGk/XIybx8vYWXnhrP3u74wQ8LlacMYVzTwpSWxEYbHc7CgpNPRrRa1Plj5uITllE+ZTpFHvduG97BhHtoigjQKcuRAcTAMdyTftI70gUh41hq24AHI3+qdyqG+e5vcAiKWWHEGIesBpbl58KvATMHm4zVlXdnDiYlmRfV3jowdqFyBi4bUsm9tSilFWAcLGz7CL+s+RL/GlXJ3HTYrZo4uayHZx63dfwFuWROUbyeUOZoB2kUVqRrmFakuIijWKjE11zURKcihiJNCEl3Lc8HQBvffrE1a+VRn9YOGxVN0KINcAFQFAI0QR8U0r563znSinfFkI8BLwDGMAdquJm/BCOG+zrigw/6KMQjqRi+8OnW/ktCa9bc7nfvIwX9y/ErXdw8ckT+Ozez3Ma70NYZ6fxOYyig8gMR6q/DyHzJKthYoaFr8geeVfs0UFauB5YkS3BjOj7jyO/l7F8RzIGGUnVzY3DPD8z5/FdwF2HtizFWMEwLXojCXoiCcIHMexjyPf0BQlVLsDdupkXrYX82Pgo78nplNLPpyZu4ryrP0uw2MOsp4PQ+sFByypwCPq7g1sXVJZ4KQsUAaR1doCBztFr0CoAKo4AqkNCMWpihkljZ5hoYpQZfAHihsnqvtP5cexvCeHjZLGfL5xXyYdrT8LjXojh8ww94m4EHIz+7ve40DWBQKBpMMlfRKDIlbXJm7XhqzRoxXHGiDT6I43S6McWccNiX3eYSKZlcCYHOwDEYSBq8MxbzTy5pYmeqMmZYhef1Z/kIm0zscmL0j7sIy11HI4h1qlpdsmjpkEw4EHLHXs3HEqDVhwFDmdnrEKRojec4EBvBMMskCAU0r6llb+rFDvAP7F1P+u2HiAcN1lQM4nPxf6HC3oeRwh7kLXe6ljn5vNXHy05+rsQUDnBg1vTUp2oo0ZJMIrjCBXoFSMiEjdp6YsykFtNk5MV59W+feXMeur6VBVNaPJSGpY/yEDMygrw59SWs2pxNbUVxSB/wLvhf6b65TvSY/yEsDtag2dgePNY944CIexKGbdLoyxQxETfQVbuKBRjABXoFcPSHYqzvyeSLdM4GXpuS38+7VuPdGRV0dCynYdef5dH3+ohHDc596RyVi2uYVYwkH7/pL9Mpg4vTWqftEfqzXrmhlHJN5oGwWIP4bhJsUdngk+nyKUNbqpSKE4gVKBXFMSyJI35auIzBmmAOchQLHeTNBn8zZbt/MS4jsflh+nc3Mk5teXcuCQnwOeSIa/okW78HdtGXSlT4tWpLvOnB1grFOMEFegVgzAtyUDMoCsUz5ZqHJkGpJOh5zEUy7PBmTBNflL1bR5p66bXTPBhbSufmvwe5Rd8DcPnH/G6RtOpmqx1LwsUEfCo/9wV4xNVdaPIIhw3ONATIRLPKZnM3GRN6uXOz2lf9rS0kpRy6vd2c//zr7HXKGWJp5GvcC8LxHupTteUs+NIJZgRVvQU6bbmXup3Zzs/KhQnEKrqRnHQtPZFaeuL5X0ua5O1bTM7V71uyyrJgCstO8g73vAtzfv51ye2srEpxCwR5Tfu73GB2E64Yj6y3QWYo6ugGaZTdYLHRbnWR6B0iiprVCgcVKBXEDcs2gdidA3EC54zSDbxV2YFUj3Sib99K/3Sx0+Ma7jPvIKitiifOncGf9f4AyZ1vmVn78vWokc609U0o+hqzUQImOR3Ux7w4NVB3L9C2d8qFDmoQD/OSZgWjV2hwVJNJo5cYgfpjryZcsxTzgOBW/hJ52K6mMDFp1Zx89k1LHj5ZgIdmwlXnEnDsrWguQZX04wi8y7x6kz0uQl4dIp0J5gPtCn7W4UiDyrQj1MsSxI1TPa0h4a2Ec7U5ivOcoZ+bM9qhnqnuY97Xn2f9zsvZm6lj699eA4nV5Wgh1sJtKy3B2K3b0aPdmD4nQlNI/SEz0XTYHqpP3+9u7IeUCjyogL9OKWtP0bHwBATnxz0SDuB1npbT3cGY9tTk+rpaW3i3o2tvPxBjPJAEV+89BQ+NDuYrkmXOZl67uMRIgSUFxdR6i9C10ThzdXx5P6oUBwEKtCf6OR4rsQMk5aeMOGuVuRwsom0qH7pDnDmsoId5OPSxRPmedz1yNuE8fC5ohf48Me/hdfjyXp55vANW9evOCgfHLcu8LldVE3w4nWPcJCIsh5QKAahAv2JTM5sTuPmdTS0DzD9ieupKWQOlhGI9Ug7gbZNzuAPQbjiTN5ujfD1xK28K2s4TTTwU/dPmam1sGvgVmKenAHTQtBw5UPpwI4ckQd8MoMPFntwq9JIheKQUYH+RCZjNJ3ct57de/ciLfJ60aQmOz2zKuUnIzU9lc2HZBHf7r6UtfH5THH18QvtR1zi2oSGiekuIVY6J/8aMjtbwx3DesBP8rupmuBNb7AqFIpDRgX6E5lABdb0JYh9GwhVLiThscsYc2ezJrPscPAM/O1bEZj4O94EbKnmNXMuXzJuZ38syNXzSvnEeUspNhbzjmcSnp737SA/gjLGQp2tyRLJihLPYZn1qlAoslGdsScwveEETd0DaKEcTTxjNqun611mP7EcIU0kgkjpKfi6dyKAvVYlXzU+zV+s05gpWvjXqj9Ret2P8mvrI9XeM84TmqAsYEs0KoNXKA6ewzkz9l7gSqAtORxcCPF9YAUQB94HPiml7HGeuxO4DTCBz0spnx/1t1CMmoGYwb7uMFJqWLlljEIbZB0snf/3de/E0ovZFJ/OrYmvEJZFrDxjMp88czZ6yXUFg/yI5q86ny2LK6kqLqIsUKTsCRSKo8BI/pXdB1yec+z3wGlSyjOA94A7AYQQc4FVwDznNT8XQqh78aNMTzhOwzD18SlLA0jZBwvAkC7ujl7Kx+NfpYoO1hV9le/03YleUrhcMZ8HfS5C2P+bOsnL3KkTqJzgVUFeoThKDPsvTUr5KtCVc+wFKWXS1vANYLrz81XAWillTErZAOwGlhzG9SqGwDAt9naE2NcVGf5cRy+X4PxP43XzVK6M38UPjY+yXHuDx4u+yRlaA4GWjXi63qPQlSP1XkLPa2kgBMwo93PK5BLKiz1530OhUBw5Dsdm7KeAB52fp2EH/iRNzrFBCCFuB24HqKmpOQzLGN9IKdnTESI21KDuHB294cqH0MPtxPs7eeXxe/ixcS3TRAdfu3AKq/ZsJNAaQgKmO8Dsx5cVlmUKDOrWXYLygJJoFIpjzSEFeiHEVwED+N/koTyn5U0DpZT3APeAvRl7KOsY7wxE47S3NBFzlRXeCC2go3dqZXzjlQPsMT7KJdom7pryZ1pPvZqGufZFwBXtYvbjy4cf9pFjaVDqdzFVD6GVlKgOVYXiGDPqNEsIcQv2Ju0nZLp0pwmozjhtOnBg9MtTDEdPKIq470pm3r+YWU9fDzJ/Rp9PR9/a1MPfr32Tpu4Id142hx9OfZmKrk3MeuYGAIxAFbGyuiFlmUxcmiDgcVFT5mX6E9ej/Wgu3LfcbtxSKBTHjFFl9EKIy4EvAx+WUoYznloHrBZC/ACYCswGNhzyKhV5aemN0t3aRN0wTUiQXcPeU7mYH7zWxe93tlNd6uOry07l1JIogVfzjOkrIMvkIgTUVgRsqwLlIqlQHFeMpLxyDXABEBRCNAHfxK6y8QC/dwys3pBS/q2U8m0hxEPAO9iSzh1SSvNILX48s68rTE84ASMdrycEDcvW0L//Xb71RoJ3P2hnumhndfFj9FTej0Fx4fcZxmlygk+nosST9qNRLpIKxXGFapgaY1iWpD9q0NiVcSOVucmKzJ99S4v2R/6Rf2q9hIjw8X33L1muvY4UOjtvXH/QhmNgm45Vl/rzz2LNMVNTKBSHHzVK8AQkblh80BkimltZk8y4c+a67vvIzzC8QVzhdh7fcoBft17DDNHKWvddTK2sQHa4CFUuSGfvTiOVHukYMtj7ilyUeHVKvDr+ogL/CSkXSYXiuEEF+jFCNGEPCTGtnDuwLLfJzA3XDdStXkqfq5Q7I5/gGetsLnNt4fv6zykREczeXrtGStjelCBG1OFaOcFD1QTv0fraCoXiMKAC/RjAsiTNvdG8QT4rMF+xGlP34Ur0A9BgVfK52OfZKWv4ir6a211PI4RAAK5Evz1AJGPjNV9lTlKbL/bqTPDqlPqLju6XVygUh4zqYjnOsSzJ/p4IA1Fj0HO5gdnT8z6uRBgBPGsu5or4d2mSFfza/X3+Rn+KcNUiTHfA7oQVLqRwZW28FupwLfHqzAoGKC/2oGlKb1coxhoqoz+OkVKyp72PRG875NHMDW8Z4eDp+Du2EapcgOktpatyKQ8dKOeHxkepE438quiHDFzxM3aWngxCULfmbLurTcKua54lVnZK+n1zSil1XWNaqQ/fSKc7KRSK4xIV6I9jmrpCTP7dx+xBIBXz2bPiUdCcoOvINv522zfe27WD6jUXcJP1LbYYU/krbRs/df+UiVqUyudvIlS1iIZla7NKKLOCfBJnY9fvcVFT5lcTnhSKEwAV6I9T2vqjDHS2ML2l3h4E0r6Z2ievY8/Kx0BoWe6TAKG4we3x/8t2WcXX9Qe41fU8LiGRkrQWH+0atvlJCHvDtaLYkx7yrVAoxjQqXTtSWJbdIXqQfQoJ06I3nKC1N4bhCxKumJ+sicHfsQ093I4ebsfwlqfcJ/dZQa6Mf5cNso7vlD3Lp/Tn7CAPWHpxtuaeLMXME8R9RS5OriymssSrgrxCcQKhMvojQc5Qbm55akSj9pJ18jHDqZMXgj0rHqX2yetSOnz1y5/LKH9cy5ZdjfzHH1uQRYK7LpzC3NrvEHp6N4HWesJBW+7RY91D1sVrGkyZ6KMsoCpqFIoTERXojwQZQ7lH6vVimBbvtw9gmDl3AJqLPSsfc4Z5SHsz1amyeWX7Xn70l1ZmlPn50mV1TCv1AdCwfC2e7t2pWa5D2Rd43RozygNqlJ9CcQKj/nUfCZJeL5o+Iq+XaMJkT0docJBP4sgthq/CKX90cbf3M/zXn1qYN3UC/37t6akgj7SY9cyNzH58me1CWcDNEuwgX1tRrIK8QnGCozL6I4EQtlwzAq+XgZjBB52h4Z18pYUebseyJN+O38ivI+fzoVnF/MNlc3Hr6fLHoZqeMpng05k2yYdL1cUrFCc8KpU7UiS9XoYI8n3RBHs7nCDvBPJBm7fSQg+1MuvpG5i5+nz+oelD/Npcxk2uF/hN87XMeW5VVtY+3Fi/CT6dqokeasr8auqTQjFOUBn9MUJKSXNP1I7rhTxmksdb6umXRdwW/yIbZB1f0B/h71yPocHgrH0I//iy4iKmTlQVNQrFeEOldMeAhGnR0BEi7lTX5JNbMo/HpMbN8Tupl6fw7+XPc8VHbyc8eWnhqU95SignT/QybZJPBXmFYhyiMvqjTMIwaGluIsSkVCA2CgwPMXxBOiuW8H+aLmarPJl/vmAyV+59j8ATVxKqXMDOG1/H8A0tD7k0kZ78pFAoxiUq0B9FpGVi/mY50w/UU5r0i3cCdT65pSuc4Iuxr/CuNcAXLp7N+TUagfVO5t+2GdCGDPK+IhfTJvlUkFcoxjlKujmKdLTux3OgPssvvnbdNWCZGeP6JHq4ndbeCF96dCt7O8PcefkpXFqjpbthRzCo2+9xMbPcj69IBXmFYrwzkpmx9wJXAm1SytOcY2XAg8BMYC9wvZSy23nuTuA2wAQ+L6V8/oisfIzR2helzZiAr2ohgZZ6wERg2R42666l8dJ7MLxBZj19A53NH/B181sM6KV85+p5XFb/aQJ/dDZql61BjxbudE161QQDylJYoVDYjCSjvw+4POfYV4CXpJSzgZecxwgh5gKrgHnOa34uhBj3KWVbf5S2vlhKotl543rCFQsyPGzepG7N2dQ+cQ1tzY3cEP8aEVPw3UunMHdClEBrfXqjNtpd0KtGCJhR7qeyxKuCvEKhSDFsoJdSvgp05Ry+Crjf+fl+4OqM42ullDEpZQOwG1hymNY6JumL2gZlKYSGEahkz4pHnWBvXweFNGlub2NV/OvEcfO/RXdxij9E9Ut3gDSRkC3X5NTde9watRUBSrzuo/wNFQrF8c5oNfoqKWUzgPNn0shlGrAv47wm59i4JGFaHOiJ5H/S8bDZeeN6QpOX8q6cwQ3GtzCFi9VF3+YUTzemP0igbbNjRexi34U/s9N2p76+bs1SZj19PQKLaZN8hQd1KxSKcc3hjgz59IK8Bi5CiNuB2wFqamoO8zKOPXGjgElZkuRQb38FLy25l6+vewuXV+OulXPRxL3sKJ0DQmSXXTpNUbl19yf5o/g8pUfx2ykUirHEaAN9qxBiipSyWQgxBWhzjjcB1RnnTQcO5HsDKeU9wD0AixYtOjjT9uMc05I0dhUwKXMkl+o/3EGgbRMbJ13BV7v+Gq9b566rT2PqJB8x6lKn5yu7TNXdt21CVC/FVzrZtkYegbeOQqEYf4w20K8DbgG+6/z5RMbx1UKIHwBTgdnAhkNd5FhCSklDR4hIPI9LWcrSYAMgeduayWeaV+IPwF3Xns7kCd5B5+edBiUE4Rsfp9g1YAd2KUflf69QKMYHIymvXANcAASFEE3AN7ED/ENCiNuARuBjAFLKt4UQDwHvAAZwh5TSPEJrPy5p748Rief/ynqkk0BLPQLJLmsqn4j/M8VE+PYV8wjmCfL5/G+EsO0MgsUewG+fG2o/aP97hUIxfhg20Espbyzw1EcKnH8XcNehLGqs0hOO09YfK/h8cjRgvO1dvpj4LBqStUX/xkDJcxg55+bzv5HFlcwKBgY3QSX975MZ/TD+9wqFYnyhyjQOE/3RBPu6ClTYJBGCNy9ey1ceeJFWWcxP3D+ldEotPXn84nP9byiuKOxZcxD+9wqFYvyhAv1hIG5Ywwd5IBw3uOk39cAk7nH/gPOqTPYsfwzb9qDDqZGXKV0+cyO2tnwYY7Kk/71CoVDkoAL9IRIzTBo6QphWgcIhy8DTvZve4pP48mPbAfiS/iCXuuqRHS67AueVvyPQUk84eAbS5SbQtjmly8viSmaW+Qh41F+VQqEYHSp6HAKmJfmgM0zCKBzkT/3tfLR4P18wP89e42y+VvYHbgs/4TQXmNS89Lf427ciMPF3vAnYzQiB1k0EjG6mT5mpZroqFIpDQkWQUSKl5IPOELFEgWGvlkHJ3t/jSvTzC3MFTxhn89cnx7gtfG+qq0wA/vathMvnpnxvAKRwEZm8kJrqGSrIKxSKQ0Zl9KNkf0+EUCxPGaUz43X2oxfjSvSz2riI7xk3stL1Gt/c93MQVqpVWAoXoapFNFyxhtqnPoq/YxuhygW0XPoLZsyYia6Pez84hUJxGFCBfhR0DMToDiUGP5HTEPW8uZivG5/i7PII/xn6FS6sVOYugT2XryY89WwQgj0rH0OPdCJKKphVUYxbBXmFQnGYULrAQdIfTdDSG837nB5pty2Fkay36vh84g7mi/f58uXziFfNT2fygOkuITxlaboUUmiYgQqmlfrxqCCvUCgOIyrQHwT90QQfdIaTzsDZSCtlKbzDqmZV/BtMFV38dOqL6JOmsO+inwJpHX7XNb/PsinwFWnUlPuVzbBCoTjsqEA/QgzTGjLIe7reI9C2mb3WZG6L/xMAv3T/iEl6DLK2Wm1qXv4sSHsjV3cJTqooZoIK8gqF4gigNPoR0hmKFwzyti5fT0Qv4Y7o39NKGfcVfZ9TtEZk2wH0SCeZ11R7qtQ29EgnZqCC6jI/QnWzKhSKI4TK6EdAXzRBx0B+D5ukLo80+Y/I1bwjZ3Dn5XUsrhJZQ7wNfwWhysVIcKZFLcJXOpnqMj/FqhlKoVAcQVSEGYZowmRfVxgrX7l8hi7/c3Ml95uXsvKMKazaeYddH18xn4Zla1Mbrg0rHrbH/wmBKKmkttSv6uQVCsURR0WZIbAHiBQI8jgOk22becT4EN83VnFhbQl/s7DEqbwxbXkm6ozbTU2UqkSbMJnJE30qyCsUiqOCyuiHoKk7nO58zTMExPAFWVf8Mb7SvozztW38wPwjTZ4H0mWU0sIompjlLR+evBDXJ5/GW6Q2XhUKxdFBpZQF6ByI0RdxXOLNOCc9viI1jDtZLdPQGeKfulZyhtjDf7t/xKTWDfhbNiGkiQAEFic9eY2j49ve8v7WTXjj3cfuiykUinGHCvR5GIgZNCebohxjMl/ndmcISD16pJNw3OC7z+6k2KPzw2kvExBRwKRq839h6oFUQaWvawdIe8i31HSEGgyiUCiOMirQ52CYFk3d6Xp5T/duXEYoZVsQKZtLwlvOT17aRUtflC9dXkfk0u+DcNmuk+1vsuu6l4iUn47ERTh4BoYvSPPVDyP+cQfc+rQaDKJQKI4qKtBnIKVkf08ky3Y4VjoH011i2xboAd5fuY4nt7fwl/c7ueWcmcybOhHDGyQcnJ8upwxU8cEl/0O4Yj7+jm2c8vyNzAoG7MEgKsgrFIqjzCFtxgohvgB8GjvZ3Q58Enti9YPATGAvcL2UckyI0h0D8bQun0TT2HHzVjzdu4mVzmHDB93c+5cGls4q45qzptkbrc+sSpdTXrGaWc+sItBSD9hafVHzRoh2qQlQCoXimDDqjF4IMQ34PLBISnka4AJWAV8BXpJSzgZech4f98QNi85QgcHemk6svI7Gnijff+FdZgX9fPHcMgQZQ7ydckpPz/upxwBSc6mB3QqF4phyqNKNDviEEDp2Jn8AuAq433n+fuDqQ/yMI46Udr18wUlR2PNe//3ZHXh1F78s+jGLHjuXWU9fj+EtszdaHdkmWQCb7QAAD0pJREFUVjon9diqPgfxBaXLKxSKY8uopRsp5X4hxH8CjUAEeEFK+YIQokpK2eyc0yyEyKtXCCFuB24HqKmpGe0yDhnLkrT0RYnEnSEieerlpZTc/fJuDvRE+M7l1dS+8qpTgbMJPdqVNcQbIWhY/iBzJ8ZxlShNXqFQHHsORbopxc7eZwFTgYAQ4qaRvl5KeY+UcpGUclFFxbGTNdr6Y3QOxJ1F2Y1NufXyT21r5k+7Orjp7BnMq63JyuDt4K5h+CtACDQNaoLFuCZUqSCvUCiOCw5lM/ZioEFK2Q4ghHgMOBdoFUJMcbL5KUDbYVjnEcEwLfqj6UlRqcEh0rSz9Ugnb/V6+LWz+XrdgumpjD036wf7x9pgMb4iNThEoVAcPxyKRt8InC2E8AvbY/cjwA5gHXCLc84twBOHtsQjx4GeKNEMi4OkQZkEQpUL6Ikk+O5zO6koLuKL55alf1kZGXwmUyf5VJBXKBTHHaMO9FLK9cAjwGbs0koNuAf4LnCJEGIXcInz+LijrS9KbyQzm7cNyuzQrWGZJnc/+CQD4TB3B+5Nbb4m5ZxcpkzyUhYoOiprVygUioPhkOropZTfBL6ZcziGnd0ft/RHE7T2ZZdSGr4goaqFtvFY+en8a/Ni/mLN43v6PSzs+VNazgm32xl9hmxT4tUJFnuOxVdRKBSKYRl3nbFSyrSPTSaO9r5z1es8Gp7PI+aHOFd7i2u99YQqnc3XygVUv/y5rM1aTYNppb6j/0UUCoVihIw7m+K2/ljaejgXoXGg3+B7XecD8Gv3f+IyTfZd9DMQGiCpW3N2urQy0snU6hm4XePueqlQKMYQ4ypCRRMm7f0Ful+BhGnx76+2o2mCP3k+j1fECVUtwvBXYvgrMHwVWaWVpZXTmOhTvvIKheL4Ztxk9FLKLFfK9BPpBqn7X9vL7rYB/vmK+QxMfpadQmD4MqprMkorPZOqmDXRe9S/h0KhUBws4ybQt/bFiMRzJJvk5KeWep4uvoYnOq5lxamTOKc2iFGo2UloiAmV1JQHEKohSqFQjAHGhXTTHYrnlWz0SCeBlnqa5US+1nEJc7VGvt9w3ZBllELAjLIAutLlFQrFGOGEj1bRhMn+nkje5wxvGaGyufxD/A7iuPmp/iO8RFIbrbkIAbOCAdUUpVAoxhQnfKA/0PP/27v72Lrq+47j74/v9UPs2HUcP+YBsqhRCmnJSjIGY6OpGKgKUMofVdk6NVsRrCrt2mlaS8U0bX9MKkzaOiRWimAtrIg0bWhhqLRhyVKkloYESGhpQvNE2zRp7DwQkthxYvu7P+4hmOTavtfXvjc+/rwk655zz7n3fO8311+f/M7v/H5957fLw9lx5L/VPY8X4hLuzn6ThVW/O3tX7MCM1nfsLsGC1gYaaqdNa5eZpUSqC/0bvac52T+Yd1u27zB7DhzmXwY+zvVVW7g1s5G3Wtw1eIbcXCpva2+qZaaLvJlNQakt9ENDI9wYleivmcVnT99JMye4p+5RIlt3dkLv+kMvk+3rObtvTbaKNt/5amZTVGoLfc+JfgYGR55IZPXzr/Gr6ODvsmtoHjpEZrCPd/ShGXYtdn7LDPewMbMpK5WF/vTA0Kg3Rm3ae5jV247yscwGPp7dAEBfyxKCtxts5m/8DMQQ7U211Ne4ycbMpq5UFvoDx0a4AAscPtHPvz/7S97d1sAX5r5KUMVgdSMzjmynt+W9oAwCGg6+yMzBo7Q3usnGzKa21BX67jdP8WbfwIjbv/bcbs4MDvH31y9m/02Ps/OWH5AZ6MtN7n1kO72tSwll6e1cRteci9xkY2ZTXqoK/bHe84cfHu4nu3p4fs8R/laPc/WPc3Oj9Lcsfnv8ms7l7LlpLXs+sZna256hzk02ZpYCqalkQ0PB/mP5b4wCONE/wAM/2sUSvc7tmafJHKw6O7b83pWryZ46wsCMVlQlujovIpv1TVFmlg6pKfRj9bL5+o/3cuzUIF/r/AmZo6J39vuYv/5OGrpf5GTHMvbeuAYk2ht98dXM0iUVTTdj9bJ5Ye8R1v3iIB/5/bnMuOU+etuWUn9oGw3dLyAGaTj4AtneHtqbaj0doJmlTkmFXlKzpO9I2iFpu6SrJLVIelbSzuRx1kQFO5KDb54acZiD44cPcN+GnSyYXc+qy+rJnjpC/aFXEO8ctKyqKndTlAcrM7O0KbWq/Qfwg4h4D7AU2A7cBayPiEXA+mR90vSeHuCN3jPnb4ghFjz9MR781hOc6uvlKzUPcNmaq5i//tO5i69kGKxuJJSht/MPmTv/Yqqq3MPGzNJn3I3RkpqAa4C/BIiI08BpSTcDK5LdHgE2Al8sJcjRjDaWzdrfNrNxaCn/lP0GS4+uy/WP736JHbc+n5vgu66F7KkjzJzdxfw6zxRlZulUyhn9QqAH+LqklyU9JKkB6IiIAwDJY3u+F0u6Q9IWSVt6enry7VKSfWdm8uXBP2dF1VZWZdYBEMpwsmPZ2akBqcpQ3dzBnFn1E358M7MLRSmFPgtcDnw1It4PnKSIZpqIeDAilkfE8ra2thLCyPvefPVHu1G2ln/o2ARVWU52XMGOP9vE3hvWnJ0asCZbxbxZ9WTcZGNmKVZKP8J9wL6I2JSsf4dcoT8oqSsiDkjqArpLDbJYz+08xEu/foPb/2QhfZc9wo5kTljOuct13qwZ1FW7v7yZpdu4z+gj4nfAbyQtTp66FvgF8BSwKnluFfBkSREW6XhfPw89t4tF7TO54X1dubb4+rbzinx7U60nETGzaaHUSvdZ4DFJNcAe4K/I/fFYI+k24NfAR0s8RuGGBli7+mGOn1rCQy2PktH9wPnNMg21GQ9WZmbTRkmFPiK2AsvzbLq2lPcdlxjiyLc/z9qTt/LXmf9h2dEf5pps6t/Z/p+pEhe11HuwMjObNqZ+20UMke3t4fTgIP986Brmq5u/yT5Bb+tl5837CrkmG98UZWbTydQu9ENDvGvNLbT+djP31nyaPXEl36i+F9qXsOemJ85rl6/OipZ6D3FgZtPL1C70vYeo3r+Z3YPtPHxsGSsWzqRjxcPsydPDBqC9sc53v5rZtDO1C31DG/1df8CX9l5HXdUgn/zAJQzU57/IWlddxax63/1qZtPP1G6slvjvxfezeeg9rPrAe5nVkL/IZzPi4tkNvgBrZtPSlC70Pcf7uW/DLpbMaeK6Szvy7iPBgtkN1GSn9Ec1Mxu3KV39eo73095Yy50ffPeIZ+ud76pjRo3vfjWz6WtKt9FfOqeJx2+/csR5YutrM7TO9I1RZja9TekzemDUdvfOproyRmJmdmGa8oV+JDPrsh7LxsyMFBf6jiY32ZiZQUoLfWNdlvoan82bmUEKC72U62ljZmY5qSv0c5o9mYiZ2XCpKvTN9dW0NHjQMjOz4VJT6N1kY2aWX2oKfXtjLdUeZ97M7DypqIzVWfkOWDOzEZRc6CVlJL0s6elkvUXSs5J2Jo+zSg9zdJ1NHmfezGwkE3FG/zlg+7D1u4D1EbEIWJ+sT5r6mgzNnjXKzGxEJRV6SfOAG4CHhj19M/BIsvwI8JFSjjEWD3NgZja6Us/ovwJ8ARga9lxHRBwASB7b871Q0h2Stkja0tPTU2IYZmY2knEXekk3At0R8eJ4Xh8RD0bE8ohY3tbWNt4wzMxsDKW0e1wNfFjSSqAOaJL0TeCgpK6IOCCpC+ieiEDNzGx8xn1GHxFfioh5EbEAuBXYEBF/ATwFrEp2WwU8WXKUZmY2bpPRj/7LwHWSdgLXJetmZlYhE9JlJSI2AhuT5cPAtRPxvmZmVrpU3BlrZmYjc6E3M0s5F3ozs5RTRFQ6BiT1AL+a5MO0Aocm+RjjdaHG5riK47iK47iKd25sF0fEmDciXRCFvhwkbYmI5ZWOI58LNTbHVRzHVRzHVbzxxuamGzOzlHOhNzNLuelU6B+sdACjuFBjc1zFcVzFcVzFG1ds06aN3sxsuppOZ/RmZtOSC72ZWcqlttBL+ldJOyS9Ium7kppH2O9Dkl6TtEvSpE57OOyYH5X0qqQhSSN2lZL0uqSfSdoqacsFFFdZc1boPMTlytdYn1859yXbX5F0+WTFUmRcKyQdS/KzVdI/liGm/5LULennI2yvSK4KjK0S+Zov6f8kbU9+Fz+XZ5/icxYRqfwBrgeyyfI9wD159skAu4GFQA2wDbi0DLFdAiwmNxDc8lH2ex1oLWPOxoyrEjkD7gXuSpbvyvdvWa58FfL5gZXAM4CAK4FNZfi3KySuFcDT5fo+Jce8Brgc+PkI28ueqyJiq0S+uoDLk+VG4JcT8f1K7Rl9RKyLiIFk9afAvDy7XQHsiog9EXEaWE1uztvJjm17RLw22ccpVoFxVSJnZZ2HeAyFfP6bgUcj56dAczIJT6XjKruIeA44MsoulchVobGVXUQciIiXkuXjwHZg7jm7FZ2z1Bb6c3yS3F/Ac80FfjNsfR/nJ7WSAlgn6UVJd1Q6mEQlclbQPMSUJ1+FfP5K5KjQY14laZukZyQtmeSYCnGh/w5WLF+SFgDvBzads6nonE3IePSVIul/gc48m+6OiCeTfe4GBoDH8r1FnucmpL9pIbEV4OqI2C+pHXhW0o7kLKSScU1KzkaLq4i3mfB85VHI55+079UoCjnmS+TGRjmRTAH6PWDRJMc1lkrkqlAVy5ekmcBa4PMR8ea5m/O8ZNScTelCHxF/Otp2SauAG4FrI2ncOsc+YP6w9XnA/nLEVuB77E8euyV9l9x/z0sqXBMQ16TkbLS4JBU0D/Fk5CuPQj7/pH2vSolreMGIiO9L+k9JrRFRyQG8KpGrglQqX5KqyRX5xyLiiTy7FJ2z1DbdSPoQ8EXgwxHRO8Jum4FFkn5PUg25uW+fKleMo5HUIKnxrWVyF5fz9g4os0rkbMx5iMuYr0I+/1PAJ5LeEVcCx95qeppEY8YlqVOSkuUryP3+H57kuMZSiVwVpBL5So73MLA9Iv5thN2Kz1k5ryiX8wfYRa4da2vy80Dy/Bzg+8P2W0nuyvZucs0X5YjtFnJ/lfuBg8APz42NXO+JbcnPq+WIrZC4KpEzYDawHtiZPLZUMl/5Pj/wKeBTybKA+5PtP2OUnlVljuszSW62keug8EdliOlx4ABwJvlu3XYh5KrA2CqRrz8m1wzzyrDatbLUnHkIBDOzlEtt042ZmeW40JuZpZwLvZlZyrnQm5mlnAu9mVnKudCbmaWcC72ZWcr9P03PHPympWDlAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# cubic model\n",
"\n",
"plt.fill_between(d.weight_s, d.ub_cub, d.lb_cub, alpha = .2)\n",
"plt.scatter(d.weight_s, d.height, 5)\n",
"plt.plot(d.weight_s, d.muhat_cub)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.7"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment