task 4 .ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TASK 4\n",
    "* Bayesian Linear Regression from scratch with BlackJAX\n",
    "  - Implement Bayesian Linear Regression from scratch with any appropriate sampling\n",
    "    method in BlackJAX.\n",
    "  - Create your own 1d linear dataset with adde|d noise.\n",
    "  - Plot the learned predictive mean and 2 standard deviations around the mean like the\n",
    "    below plot\n",
    "    \n",
    "Here x is our feature matrix of 1D which is randomly generated using random function. Noise is variable that stores noise or disturbances. Now we have formulated line 'ys' as (feature_matrix * weight + bias + noise). Finally we have plot our line 'ys' using plot function of matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 296
    },
    "id": "L86z9jvRkZeZ",
    "outputId": "5b47e42d-8125-43dd-b881-7ad2f227649f"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7efe6cbd77d0>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2de3xU5bX3f2vvuSQkXGJQEMJFGpEmVFCjaMELqBUVsK3UniPW0/aDvPYV29OqoPIiIqcXRXtOWz1tKdpTK/bUQis3bdWCVVDQoAkmiJhSEYIipIAEwlz2ft4/9uzJzOxnz+y5Zc9k1vfzsSUzs2c/M4S1nmddfouEEGAYhmFKD8XtBTAMwzDuwA6AYRimRGEHwDAMU6KwA2AYhilR2AEwDMOUKB63F5AOAwcOFCNHjnR7GQzDMEXFtm3bDgkhTk18vKgcwMiRI9HY2Oj2MhiGYYoKItoje5xDQAzDMCUKOwCGYZgShR0AwzBMicIOgGEYpkRhB8AwDFOisANgGIbJgI7OAJr3HkFHZ8DtpWRMUZWBMgzDFAKrm9oxf9V2eBUFIV3HQ9efjRnjh7q9rLThEwDDMEwadHQGMH/VdpwM6TgWCONkSMe8VduL8iTADoBhGCYN9h3ugleJN51eRcG+w10urShz2AEwDMOkQU1VOUK6HvdYSNdRU1Xu0ooyhx0AwzBMGlRX+vHQ9WejzKugr9+DMq+Ch64/G9WVfreXljacBGYYhkmTGeOHYmLtQOw73IWaqvKiNP4AOwCGYZiMqK70F63hN3E1BEREA4hoJRHtJKJ3iegiN9fDMAxTSrh9AvgJgD8LIWYSkQ9AH5fXwzAMUzK45gCIqD+ASwB8HQCEEEEAQbfWwzAMU2q4GQI6A8BBAL8moreJaDkRVbi4HoZhmJLCTQfgAXAugJ8LIc4BcBzA3YkvIqI5RNRIRI0HDx7s6TUyDFOg9AYtHrdxMwewD8A+IcTWyM8rIXEAQohlAJYBQENDg+i55TEMU6gUmhZPR2egKEtCXXMAQoiPiWgvEZ0lhHgPwOUAdri1HoZhioNYLZ6TMDpy563ajom1A10xvoXmjNLB7U7g2wGsIKLtAMYD+IHL62EYpsApJC2eYheGc7UMVAjRBKDBzTUwDFNcFJIWj+mMzJMI0O2MiiEU5PYJgGEYJi0KSYunkJxRJrjdCMYwDJM2haLFYzqjeQk5gGLY/QPsABiGKVIKRYunUJxRJrADYBiGyZJCcUbpwjkAhmGYEoUdAMMwTInCDoBhGKZEYQfAMAxTorADYBiGSUJvFp3jKiCGYRgbilnnxwl8AmAYpiRJtbMvdp0fJ/AJgGGYkiNxZ7/w2jqMHdo/rpGr2HV+nMAOgGGYkkImJ73g2RZU+lWEdREN8xS7zo8TOATEMExJIZOTBoDOgBYX5ikk0bl8JaL5BMAwTEkh29nHEhvmken89PT0r3wmovkEwDBMSRG7s6/wqZbnE8M81ZV+jBs2ANWVfqxuasfEBzfgpuVbMfHBDVjT1J7XteY7Ec0nAIZhehVOduixO/uW/UexZN2OlHLO6Y6izMVJId+JaHYADMP0GmLDJUFNw9zJZ+LGCcOlxtJU8Bw3bACm1g9OaazTMca5CtvkOxHtegiIiFQiepuI1rm9FoZhCot0kp+J4ZJAWOCRF3fh8z9KHaqJDfPY4dQY5zJsk+9EdCGcAL4D4F0A/dxeCMMwhUO6u2jZDh0AAmE9aajGKU6nf+U6bJPPgTOuOgAiqgFwLYDvA/iem2thGCYz8lEVI4u337lyO+pO74faQX2l1ySr7slV3NyJMc5H2CZfA2fcDgH9F4B5AGxrsohoDhE1ElHjwYMHe25lDMOkJFdVMYmhHlmtfjCs45qfvmp7D3OH7veQ5blUBjidUFOqcFEh9Q+kwrUTABFNA/CJEGIbEV1m9zohxDIAywCgoaFB9NDyGIZJQbpVMXbIQj0TawdKd/NBTSS9h7lDf3rrh3h0Yxt8aupB7bmus+/oDGBEdQXWzZ2E40GtoOcEuxkCmghgBhFdA6AMQD8iekoIcZOLa2IYxiG5iHXbOZHN86fgoevPxp0rtyMYjncEqe5RXenH7Zcb1T+pQlO5cmImMmcybtiAtN+np3AtBCSEuEcIUSOEGAngXwBsYOPPMO6QidRALmLdslCPaeBnjB+KO68cbbnG6T2cVPYku3+6FKN6qNs5AIZhXCbTOH4uYt3JnEhHZwA/fmmX5ZqF19blLKSSy4RtLp1JT1EIZaAQQrwM4GWXl8EwJUe2IZDYqpgKn4rjQS0qpOaExNLKoKbhtstqAchDTBV+FWOH9s/gkzq7f6p8QTKKUT20IBwAwzDukIs4fnWlH5vaDmWcSJ0xfijqTu+HJzb/A6veascvX/k7fraxDXdeOdpiUDVd5Nyg5qrOPpfOpKdgB8AwJUwudq3ZnCI6OgNYsfVDPLrhfQQ1o8gvEDae+8HzO3FDw1Csaf7IsUHt6Aygdf9RAIT6If0c3d80/LlI1uazaSsfsANgmAIhm4aqTK91ums1398M86Q7OUu2vtVN7Zi3shmBsH1197Nv78dz377YUTnl6qZ23PmHZoQijsSjAD++YbztSSQd3aB0yFfTVj5gB8AwBUA2tejZ1rGn2rWa7y90gYAmUOY1Ep1OJ2eZ13sUQlATWDS9DlPrB2P+qu1JjT8AeFUFx4Nayt15R2cA81Zujxp/AAjrwF0rm6UnEdmp5ZEXd+HRjW1YOrN3DX5PBlcBMYzLZFM+mO/Sw9j3D0SM68mQHr1P24FjaN3/Kb7x+ZHwe8hSDRR7fWdAQzCsY8GfWrD81d3SqVyJnAxpUs3+RPYd7oKqWDuAdV1Iq3DspoKZukGFXLqZS/gEwDAuk00iNhdJ3GQnCDuBNQAQQuCq/3oF5qbbqxLmXDIqLoyy73AXPBLDvHzTP0BI3divCeDan21KuSuvqSqHplvfL6RD6kB6QjeoGOATAMO4TDaJ2GyTuKlOEMkMZSAsEBNxQUgTeHRjm2V9QU1i6HUBLeZtVQJmXTAclX6rsXayK6+u9GPu5FrL42VeI4Qke32mukG9CXYADOMy2TRUZduMlap5Kfb9fQ6shapQXMilutKPRdPrLK8LCcQ5D0UhfOmcIQhLdvGJa4ql7cAxrGzci7YDx3DjhOFSg25nzGeMH4rX7r4cd1w5Gn5P4Qu35QMSonj01RoaGkRjY6Pby2CYvJCrKiAAthU7susmPrgBJ0Pd2/Eyr4LN86fEXdPRGcDGnZ9g0ZpW6Y7axKsSttxzueV+K7bsweK1rfAohBMh+YnC51HwL+fX4Pdv7rUkh2Vruu/Zd/Dklg+jP9980XA0jDjFUtHkJKHb04Peexoi2iaEaLA8zg6AYXoPiRU7fo8CAWDR9DrMmjBCes2apnZHRlPmLBLxKMDWe68AYHVCALC2uR33r33X9voyr4J1cyfh+ZaPLWqesWtqO3AMV/znK5brX/ruJaiq8PVqY54Jdg6Ak8AM00uIjeebBCJKmgv+1IJ9/zyB2RePshhFswy0df+nAATqh1ilFswd8sJpddEB6oGwBiKK3gMAyr0erNj6If77ZSMXcDKkG6EjItz5hbPg8ySPI3kVI2afSs2zae8R6fWb2g7hnOFVbPwdwg6AYYoQWciidf+nUGCNgZv8/G+78cTmf2DpzHGWHX4yKYfEKqGF19Zh7ND+qPCpmPboprj3ORnW8NOXdiE2ghPUAUDgB8/vhE+1Xx8ABLXuBGyyhqrxNn0BS9btQB+fJye6/qUAJ4EZpsiQqXeubmrHLU824kTIPj4PGJU7d63cjld2HYxW1dhVArUdOIZXdn2CeSub455bsn4HaqrKUTuobzRBbDaHhTSBZL1dsRVBfpWgEuLq9zVdx+a2Qym/g9pBfXFDg9W4awJFI8VcCLADYJgiQmas71rZjHkrt8eFYpIRCOu49bfbos5DVl2j6QLX/PRV3PrUW5aErFmRY06+euqbF0CzKRVNhgZg/tVjEBsVCutwbLhnTRiZtEms0KWYCwEOATFMESFrzFJJsS2ftMM8Kdy1shn3T6+3JHZNSYWgZj1RBDUNLe1H8dVlrxu5AE2Hg54uC2FNYOlfdsLvUREId9/HaSNWTVU5tCRFLKVUz58pfAJgmB4kk8lbscgas8K6hqCWegcuy78GwgL3r92R1hp0ASxa0xI9hQTDOpIUBiUlpMHqfNKY+BXbA+FVCR7FmBng8yhYOC13g2N6K3wCYJgeIlvRtmglzrV1eGBdK1RSoAkdcyefiUc3vB/V6gGMnV2iTVaI4FMFEsv4nYaOTEKSzl6fqjhyQjKEMMpVnQxwTyRRyO7PLR9j8dpWeFUFS9btQF+/hxPBSXDNARDRMABPAhgE4wC5TAjxE7fWwzD5JFvN/OWv7sbjm/4Bn0fByZBRfqkqAAThlEofSKG41lodgE+l+KSrR8XFo6vx3DsHHK+7j09BWBNQFEpa/09kvR9g6PCEND3p9X18Hjw261z0L/dmVL5pVgt1dAawZP0OBDURDV1lM+C9FHDzBBAGcIcQ4i0i6gtgGxG9KIRI7zzKMEVApqJtq5va8b3fN0Vte3dMXiAU+fOSdTvwvStH4wfP7Yy7NtEYBzUdf333k7TWHQzrmHfVGMtsXq9KUAjwqWp01w4grqFs4bQ6jB0iLxeNJaTrjoa3pCIXwnilhmsOQAjxEYCPIn8+RkTvAhgKgB0A0+vIRLTN0LhvhkxLLRavYoRP7DBLNG+7rBbLXtmNQDgc/7yHcNKmdjOsAz9+aReuqhuE1c0fRR//1wuG4TuXj7b0ItjNFYgdOnMyrEEIgXKvB0FNj84ATvzs6XbzFuNMXrcpiBwAEY0EcA6ArZLn5gCYAwDDhw/v0XUxTK7IZF7svsNdUEmBUTBpT0jXMbDSZ/u8rgs89+2LUVXhw2Mvx6t1qoSUCdxgSI8z/gDw+zf34orPDkL9kP5xn8GueSsxVg8AK7Z+iMc2vo9lr+zGYy+3RXMimeZKinEmr9u4rgVERJUA/gbg+0KIPyZ7LWsBMcVOOjvbjs4APv+jvyadmuVVCY98ZRwm1g7EhT/8qzRB29fvwVOzJ2DcsAFY09SOu1Y2R0pHNQiQ9Bon9PGp0IXIqOPWToRu3dxJmPboppTidKnem7WA4rHTAnK1DJSIvABWAViRyvgzTG+gutKPccMGOJZ6/mrDsKSvUQjRJOcjXxkHv6TWMzYMcuxkGLoABAQECJREOiIVJ4Ka447bxPJXOxnqpr1HkspTOyGd77jUcbMKiAA8DuBdIcSP3VoHw7iBk11qR2cAz2zbl/R9PIqC1v2fon+5FxNrB+K1u6fg6a0f4mcb3odHVRDWdHzj8yPRuv8o3v3oGH74vJEo7t71Zx8BiDXQiZLUNVXlUp2hibUDpfH68cMGOIrjy74/3vmnj2shICKaBOBVAO+gu2T5XiHEc3bXcAiIKSbsDJLTGHfz3iO4aflWHAuELc+ZEAwd/dgaegFg3spm6DoQStEhrCJVhiE1ZV4FC6+tw5L1hkpoVygMIkKZR0VQ06HpOmJbDcyQzua2Q/FVQxGRuZb2o9H3kn0/su9PAFn1WPR2eB4Aw+SYtgPH0LT3CMYPG4DaQX3jnrMz8k4HsADO9PcTMdQ2rc1e+cCsLjIlop2us8KnYvGMekwecxoA46SQaPTNEtJE5yn7Toywl4jLlaSbN+jtFGQOgGGKlfuefQdX/OcruHPldlzxn6/gvtXvRJ9LNmd33+EuiIRdudAFWvcftUhEdM+tVdDHp0qHqydiNEHl7nMmQ9cF1s2dhGFV5VDIeS7heFDD/WtbMfHBDdjcdgg1VeVYsn5HvOLouh3SUI4sd6AqFKmW6oaF4JzBDoBh0qCjM4B1zfvjRhECwJOvf4i2A8cAJJ+zW+FT4yQbACCgCcz+zZtx8s4mIvK/YU1PW/AtV9g5Hq+q4LmWjw0Z6hRjIhNz052B7gRy6/5PHSd+ZbX+mi6gCa7/zwR2AAzjEFOH/84/NEufN6dUJWtIOh7UoqGTWIKaVcfePEkEwsLS1ZsNPjW92p8rPnsqyiTD1oOajsc2vm8pU1XJMPp9/R74PQpmTzoDdocXw/ALxw1ciQJwZV4FS2eejaUzx8U9xvX/ziiIRjCGKXRk4xYTMadUZduQFLv7TZQ2yBaPQhg/bADe+OCw42v+3CqXj5g96Qw8sekfcY95FcLjXz8f9UP6RRu9fvP6HtuwlCED0T+t7yuxqSxVFzJjDzsAhnGATGcmlotrq+MSwbLO1+a9R1BTVY6Hrj8bd61stm3wit39Ju6MsyWsi7SMvx33Xj0GU8achp//bXfc4yFdYEj/MgDAf7/cFvmMcuvv93Tv1O2Muh2yjuNkIyQZORwCYhgHyMI6sbzxwT8tzVBmQ9KmtkOY+OAG3PirLbjoh3/FsZNh/OrmBvSRTLPyqRQ1irHhjlSzdHOB30O4oaEmZbK5wqdizOn90LT3iKXxrMxrDHW3S8D28arwewh3XDkar909BRNrB0aT39zA1fPwCYBhHBAb1lFAltm7gbDA01s/xO2Xnxmt9qnwqdh/9KRlXOOCZ1tw7zVjoCeUYPs8Cp67fZLlJLHp/YN4Zls7nOJVjKEt6aQNrj9nKM4dUYUl63fAp1LShPPJkIZbnmyERyHLLAFNF6ipKsfh40FpuOyhmWfjos9Uo7rSn/V8BCZ72AEwjEPMMEXr/k8x+zdvWhKzj25swykVPjywbgeEMJKkHoJ0SPrDL+zCnVeOxsMv7oJXJWi6oamT2E/QduBYWsYfMMTdzh7aF9vbjzl6vUqGuudNT7yRNMfhUQx1UE0AWliHTPzB7Cs6HtTgVymu4smvEoad0ieq3Z/pfAQmd3AIiGHSoLrSj0tGn4or6k6zPOdRCfetaUEgrEenY9npuBGAh194D16FEArr+N4VozGiuiIujNTRGcDvG/dmtE4nxt8M9HgUwo3LLUK8UVQC7rlqtCPRiHKvJxrHp4RQEikUJxORreYPkz18AmCYNOnoDEgHqwTDOpxORTRDJ+aAlx88vxMVPhVaRF3TkHPYjlx36quRwWFetVsFNLpLt4kZaQL44V92SZ9LxExgmyEzU3lUE/GVPazdXxjwCYApGbIdyG6y73AXfKo1gVs/pF/Ka8u8CnweBZJWAByPqGvetXI77nimKXKSyN4BTDlrILyR3bj5dtK5vgrSTjYrZISGZPX3xh0octSIf19ZPT/X7vc8fAJgeiWJQmy5TDjaVQQ17T2a9DqVgGVfa8CqbXstA1ZiIQLSnNOelPNGVOHND44glERUDgCCOnDv1WchENbx07/uSjkoBjC6gdffPgnHg5pFmdNoYut+k8QYf7qln0zu4RMA02swd/grtuzBxAc3RKUVVmzZY6vNkwndGj1p7pYVQh+vktT4A/LdeaYQgKvqT3fcT/Djl3bhxgnD8fx3LoFPMlsgEZ9qlH0mlm86jfFz6ae7sANgegWmTMONv9qCBc+2xBn7xWtbLbXtChFa9yffsSdjxvih+NXNDSiXxXJsKPOoeOX9Qylfp+VQ8+fWS0ehdlBfLLy2Dk6iO6aRrh3UFw/PNEI03iQX2sXtOcZfHLADYIqe2JLC4xLNAa+qWGLpJ4JGLXus8JqT+8TmEPb+swtdkjiJXRw9qOkYV9Pf8f1ywUWfGYjVTe14YF2ro76AWCM9Y/xQbJ4/BStv/TxW/p8L8eVzhsZF8r0xTWuJcIy/OOAcAFP0pJJp0ITAoul1eGDtjriYdCAsHNeeJ+YQzAEoiXgUsk3cnlFdjvXvfJzGJ8sOjwL08SqYvXI7gknmCseycFqddMj7ii178Me3451lSBNJh9VwjL/wYQfAFD12SdkKvxptsJoxfiiGVfXBrb/dFtfFa4Y8ZMYptqM3sWlp8dpWaYw8WQftzgPHsfPA8Uw+Ytp4FODGCcNx4/KtjiuJKnwqxg6xnlA6OgNYvM7q7ABg8dodmFo/2Na4sz5PYeOqAyCiqQB+AmMy3XIhxI/cXA9TnHTXnG+HqhhdtfdNM8YLxu4864f0g57QzhTUNBztCka1aExid/yBsAYlIYfgVZW0JnX1NAoRfv/m3rTKSMM2MXrjhAUEJdd4VSOX0r/cx7v8IsTNofAqgMcAXAlgH4A3iWiNEEK+1WCYJJiDUyCMkYh9yzwYF5FnNkmUae4KhaEL4LYVb1vGNibu+BMD6CFNh+7SgBYnqArBOBQ5X+M3J54RrdKJNeQt7UdxPCh3dsGwjluebIRPVVnPpwhxMwl8AYA2IcRuIUQQwP8CuM7F9TBFSuzglBMhLRrbl5V5monNx2adC1VRonHs2NLQ1v2fQkloXPKpBDXmFBDURE5U+u3+AWar/dkV0hFw2pYcud8Tmz+wTCXr6AxIcx0+leD3KBDCmMWbi/Japudx0wEMBRArdLIv8lgcRDSHiBqJqPHgwYM9tjimeLCrOW/d/6m087e60o/+5V74VOs1K7Z+aIw4TFD7DGoip+WZAPD1zw/HDQ2FsVtWI8qeshnGid9tH5+KH3zpc5Ey2PggAuv5FBcFnwQWQiwDsAwAGhoaCvfMzbiGLAl8IhiOhCbknb+ya+xGHOYDvwqMrK6U7q4BY0fuZBVmziMdEt/7is+ehq27/xlX0WMactn3pAuByWMMMTyu9S9u3DwBtAMYFvNzTeQxhkmL2Jpzv9qteSPb0cquqfCr8HkUfHPiSIvGj9+joMxBR2wa/WAAgIAG/PC5nbbicU6DN5PPGpjejQH8x5fG4qXvXoKHZ56Nl757CR68/mxbQ56snp9r/YsfyrXaoOMbE3kA7AJwOQzD/yaAG4UQrXbXNDQ0iMbGxh5aIVMsmKGKUFjDjY+/gaBESKev34OnZk+wJIZXbNmDxWtb4VUVhHUduoiXYjAmXgnpqcDvUSwDUXoaVYFjBVIA+NalozD/6s9aHl/T1G6ZyRt7YkrUVool2XNMYUBE24QQDYmPuxYCEkKEiWgugL/AKAN9IpnxZ5hEOjoDkcHjbfCpipH0tNnQnAxrCIW16FxecyjJkvU7ENREVJZZJcOw+1QFQU3H3Mm1OKXChyXrd0AlQlDTMOuCEWgYWYUNOz/BH9/e35Mf2UI6xl8lYPbFo6TPpWraSlbPz7X+xYtrJ4BM4BNA6WG3u1zd1I55SQar2+FXCaQYEgYjqitw469et5Q4ehXgkrNOxcs7D0KNjFe85nOn47l3PoKIjFrMZPfvNK6fClPTP10UAl7490ssU8diMaqgjgIg1A/px4a9l2B3AmAHwBQsMvmFsUP7o8KnYtqjm7JqxPJ7CFPrB6dU5swWrwIoioK5k2txxsAKzP3d23m9Xyp8KuHhr4yT1uqvbmrHnX9ojobAPArw4xvGc11/L6DgQkAMkwxZM9aCZ1tQ6VcR1AQoy41LICzybvwBQICwPjLo/Zd/+3ve75eKoCbXP+roDGDeyu1x+Y+wDty1spnn9PZiWA2UKUhk9ecA0BnQEAzrccPGCxmfx9DL7+gM4JEXnY1VzDeyWv19h7viGt1MVOK6/t4MOwCm4OjoDOBoVyg6WF1GmVeBRzFi8SrgSOveDcK6QE1VOfYd7sq6uzdXyGr1a6rKpf0EmuC6/t4Mh4CYgsKM+wNImmQNaQKabgiYAchNdjUP3PmF0dESVbdLRmOHzsuqfJbOPBt3JOQAls4cx+GfXkxKB0BEtwN4SghxuAfWw/Qy0qkRj437J6IS4FEVeFRCMKzndGxivrhu/On48Yu7DEVRTYdCRkVRT9LHq0KHwMJpdRg7pH/SvwezFJSrgEoHJyeAQTCUOt8C8ASAv4hiKh1iXKG7Rv99x0qRrfuPQiF5oKSPz4NZFw7HE5v+ASqS3761zR9BF7AdVJNv/B4Fv/jaeUkNeaKDrq7045LRp/XwShm3cFQGSkQE4AsAvgGgAcAzAB4XQvRoWQOXgeafXHR1GjX62y0hjzKvgs3zp+Dw8SD+0mpMxrqqfjBqB/VNWdefrCOXicenKlAUpHS4iWW2LOXce8mqDFQIIYjoYwAfAwgDqAKwkoheFELMy+1SGbfIhUHolma27nq9ioKFz76D51oORB9b+sIuzBh3Ov7c8rF0eIlXJShEmDu5Fr/4298RgHXmL9ONTwWW/9t5qB/SP6kDl5XZOh2PyfQeUlYBEdF3iGgbgIcAbAbwOSHEtwCcB+D6PK+P6SFiDUKigFriMPRk2JVvAobaZqzxN1nT/JHt5CoSArrQ0fbJMenAd8agwqeizKvg4a+Mj4Zwkv2d2Uloc8lnaeHkBHAKgC8LIfbEPiiE0IloWn6WxfQ0ssHqpj7+f7/cBq+iIKhpmDv5TNw4YbjtLtFuPq/fQ7j+3KF4+o29kqvsMVUaeqJpqxBQFcLtkz+D//prm+U5mZSE36PgvunxCV4nJznZ3xNLOZceKU8AQohFicY/5rl3c78kxg3k+vgaHtvYFj0VBMICj7y4C5//UffEqEQSJYL9HgV3XDkar919Ob458Yye+ChFjaYLHDkRkj5nqlL7VIJXMZQ9X7t7CmZNGIFxwwZEBe7sTnKxsJQzA3AfABMhcV5uSNdx22W1WPbKbks8PxDWk8aL7ZQlqyv9uPmi4Xjy9Q975DMVK7/dIt1vwayODWoCFT4Fv37tA3z29H5xu3u7k9y+w12Wv6tUCqBM74cdABMl0SAAwKMb35e+1s6omNhJBD9w3edw84Ujo1VAXlXBIy/ugqoQgmENug54VIIujC7aUsTvURDURNLPbyqYzlu1HXWn98PxoIaaqvK0Qzss5VzasANg4og1CKub2m0bl7KJF9cO6huVJF7d1A5AIBQWCOuG8RMAFs+oAwAsXtMKInK9i7Yn0QEsvq4eD6zdASIkVT0VusA1P9sEf8zoy8STHId2GDvYATBSzFhyYsdthV+FpsvlBDK9R2xtv2nol6zbgYXX1kHTBbQS6jv0e5TupK0A7l+bfEaSIYonolPQ5q3ajs3zp2Dz/Ckc2mFSwg6gBHHS7CWLJVf4VCyeXo/JY05Ly6jY3U92DxMC8P+ebUkq8aPA+exct7Ab3uJXCToAIQTKvR5LhVVHZwAPrGu1lbzo41MR0neletAAAB18SURBVDQQlDjRPDM0ZyaFGSYZrjgAIloKYDqAIIC/A/iGEOKIG2spNZw2e8liyZoQaRv/ZPerqSqPjmJMpMvBsBcdRqNYIesCaUJevqkDuH9GPabWD5Y6x+Wv7k7a9XzDeTV4+o0PLYqpXMrJpINbctAvAhgrhDgbxmD4e1xaR0nhtEQQyE2ZYKr7bWo7lLU4WiEbfxNVIXgTtPZDmsCSdTsAwLJbX7FlD37+t9227+f3KBHjLyyPc7yfSQdXTgBCiBdiftwCYKYb6yg1nJQIxoZrsi0TtLtf6/5PAQjMW9lc0AY8mXpnzQA/9h1J3RkNGA4gLGmOU4jQuv9onPhaR2cA969tSfp+Akb1VOzpqY9XxS++dh4uGX2qozUxDFAYOYBvAvi93ZNENAfAHAAYPnx4T62pV5KqRDBZuMaUCEjHCcju1xUK45YnG6EqVPDCbnbGXyXgaJfmeDi7XQXTiaCGW55sxNKZ3TN6V2z9ECEbxQtTz3/htLro6SG6VgjUD+mXejEME0PeHAARvQRgsOSpBUKI1ZHXLIAhLrfC7n2EEMsALAMMNdA8LLVkkDV7mSEDO3GwYyfDWLJ+R0YCcYn3C2oadJF80EsyPIoRTXe7IlQTwLFAWPpcYrzfoxBUgu0Iy0C4e0YvADy20SoBAQD3Xj0GE0ZVR09iff0eLvVksiZvDkAIcUWy54no6wCmAbic5wv0HHZhHVm4RlUIi9e2IqgJqWKkk2qi2Psd7QrithVvI6TJjWcqiqExLHGFCgmAFMkz3cSKsPlUxeIgv3XpKMy59DNxj3EXL5ML3KoCmgpgHoBLhRAn3FhDKSPr/pSGhzRhiTWbxmpT2yHH0tHm/To6A1KhuAtGVuGND6wD57xKt/xBIsVQAgoAiqLgtstq8bONbdFa/USCmhYNwyV+P36PgtkXj5Jex128TLa4VQX0KIC+AF4koiYi+oVL62AiyKp+Fk2vszRhhXQdobCGu/7QHFfdc9dKazVRoox0daUfC6+ts9z77b1H4FOta0pWCdpTxr/cq+D6c4dAtZlUloqTIR1Xjx2MO68cbfuamecOAyD/O1g6k0M7TP5wqwqo1o37MsmRhRUSY803nFeDGx9/w1KCGAjreHrrh7j98jMBWBPK5kza/uVey31DmsC3Lh2FX7/2ARQQTthlQV2gK6TjrMH9oCofwaZlISWr3tqHX7/2ge3zzza1Y9Xb+6KnKA7tMD2Fo5GQhQKPhHQHM9Zf4VMx7dFNtto0fg9h/e0XY//RLtzyZKOlyqfCp6IrqEl379+6dBRmXzwKy1/dnbQG3g08ClDmVdEZiPcAyUJUsfhUBV6VUg60MUdmstFnco3dSEi3QkBMEVFd6UdNVTma9h6JVOLI0XWBq3/6Km596i1piedxG+MPAE9s/gcOHw8m3Sm7RViHxfgD8jJR2dfj9RBCCR27HgWo8MX/81NAkR4Jg3QmsTFMJrADKGGcGpjVTe2Y+OAGLFrTKjWEJiHdCOecyGB0I4Hw+8a9UJBZrN0NFIXi/gGpCnD31DGRAfbdhDSBO79wVlxsf/GMsZYeghMhoy9gTVN79Du/aflWTHzQfgAPw2QDh4BKFEuM/to6jB3a3xJ37ugMYOKDGyxhnwq/imBYh66LpM1QPhXoraN8fSpZciEelXDbpZ/BL181wlgnQzr8KoEUiuZBzO94TVM77lrZbDkt+T0EIF4Cm8NDTDZwCIiJItPoWfBsC2Yt3xK32+zoDGDjzk8sYR9TFXTLPZfj1984H+Ve+a+RAkAX8h29T83vr96/nO+sWS0bZIPsw5rATza04aq6QdAjMaKAJnAypGPJuh1xDnbG+KH41c0N6JNQAqWSAjXhO+eB7Uw+KAQpCKaHsZNhNsM7sR3AKlmTl2FdoMyronX/UQzpXw7NrkGL7Ju3bpwwDKdW+vCTv/7domiZLgoAr0eBRzFi7XdedRaG9C/HM432A23sUBUgy+UAMIbYV/iUOCchm6JWP6Q/9IRTuCZ0IMFxssonkw/YAZQgsqavWGI7gGMxwz7BsI65v3sbgCHHfM3nBmN100dxr/UolLRz939e2wMvAaEsIpAqGXH4WFnllvajWLSmJSO5CI8CLJ4xVvrZMyExtCMz4nbyHABY6oHJO+wASpBYoyPb4cs6gCt8Ku76wln4wXM74kQNQprA6qaPLJ25tqeC2GsztLEKAIoIsWmawKLVLejr96Du9H64f63c+HsUpHQK/3vLhWg4oxogYMGf4hU5kymD2qHHXOBVydaI29X+cz8Ak2/YAZQosUZn6+4OPPzCe1AVBboQuG+6VW1SEwIDK/0gkoswJD6Sz9ICPeEGYR347u+bQES2Rj6V8b/5ouE449RKNO89gqn1gwEBLF7bCq+qQBMCt11Wi0c3vC8VdfOphKsTTkGJDkMhRAXfZMhkHVjqgck37AB6ManE2qor/djUdghLX3jP0OXXDInjvn4PHrr+bNy1cruhZa/pmDLmVHzvmbcLtqLHGI1r73ZkFTt9fCpmTRiOcTUDcLQrhIkPbogLubx+z+XR7w8AHnu5zaL/7PdQVM759snH0LT3CMq8Cu75Y0ucYqhPVS3xf4ZxG3YAvRQnox87OgOYtzJ+8LsmgLtWNuO+afUABEJhgZAu8Nw7B3pk3bLxiTl5XyL4VBHnwEKajidf/wAeRYmGwWJVTzfPn4JxwwZEXx8vba1j7uTa6AxfAKgd1Be1g/pKRe84icsUIuwAeiF22v6mjLPJvsNd0rYrAmHxuh0IujCwxaNSpJuWMp4bkIhPJSydGZ9YDWo6NF1HQAMCsB5rZBU7TnV6ks1dYJhCgh1AL8TJ6EcAWLF1D05KjKym6/B7VQQzk+3PijKPih9++XP43jNNjq9JJg3tVQnPffti1A7qCwAxswlCuG3FW7aDXex27E7j8izqxhQD3AjWC0k1+hEA2g4cwzON+yzXKgDunzHWteErQU1Hv3IvvGk0in3pnCGR7lkr98+ojxp/wDDg44YNQP2QftJS2Aq/ijJvboarm/di488UKuwAeiEyXflEg9a094j02vumfxazLhwRvd6vGobVzsA6QU3j0rmTa7H3nyekypmymQEAMHnMIDzylXGWjto+PgVjh/SP/hyrfST7jr7/xbF4evaF2Dx/iuOxlwxTzHAIqJeSKgQxPia5Gcuk2lMt11f4VOw/2oVv/PrNuCIYpwlbTaRuDAMMJ3PRqFPw1WVbLM95FIJCBIIed08CcM8f34nG9GM5EdSxdXcHAKCl/ah0tjGHaZhShh1ALyZZvLp2UF/cfNFwPPn6h9HHbr5ouCVcYl5/PKihj88TFzPv41dxMqglFYMzUR04gPNHnoJ/Xb5V+n5hXUSv9yjAd68YjZ/8dReCWveAdq9qdUk/eH4nKnyqtMrHTIqz4WdKFVcdABHdAeBhAKcKIQ65uZZS5IHrPocZZw/BK+8fwiVnDjS6YG2oqSpHVyg+YRoIaXjgurG4f01LysEoiXo3Mja1dThad1g3YvV+jwfBmAHzXlWBSsKS2JaFk2RJcRMnw+4ZpjfgWg6AiIYB+AKAD1O9lskNifr/q5vacdMTb+DXmz/ATU+8IdWcN685fDwISpiLS0SYOnYwnv/OJfB5kv8qfemcoZEdem4YWFlmSeJquoBw2EVgV+XDOvxMKeFmEvg/AcxDflUDmAiJhm3Flj0WSeh5q+IHu8dec81PX7VIFJd5DEXQ40ENi6bXRROqfg9ZJmP96e123D+9Hv6Iamc2KAAu+ky1dID60pnjUOZVUGGTMa7w2Vf5yGSyE78ThulNuBICIqLrALQLIZoTd5VMatINUcgawxavbbXs2mPDIrJrEoPzXSENs3/zJvweT9xQmT+9tQ//8/qeuNeGNIFhp5TjtbunYN/hLhw6dhK3/HabVGDNrxqp3qs+OxjrWj62PH/zRSNQXelPKaK2YuueuFLXGxpqMGvCCNvvzWn/BMP0FvLmAIjoJQCDJU8tAHAvjPCPk/eZA2AOAAwfPjxn6ytUUhl3U+LBoxjaNoum12HWhBFJ31Nm2DyKtdM2Niwiu8YwzEblTVDrTsqacfgl63dg3dxJeDLB+HdDcUnXJV8ca1HdBIwBKn6V8OLOT6SVRr97cy/OHVGFGeOH2oqoAcCa5v1xj69p3o/5U8fYGnMn/RMM05vImwMQQlwhe5yIPgfgDADm7r8GwFtEdIEQwrLdE0IsA7AMMEZC5mu9hUAq/Z7YXbnJgj+1AAKYdaG9E5AZthORUYWAMW4QQFxYRJb0DWoCXomomolXUbCp7ZC0K1choH5Iv+jn2He4K6q6KUsiG6qbxv0gRNzzgbAulbaIJZPdPEs4MKVGj4eAhBDvADjN/JmIPgDQUOpVQE70e/Yd7pLGzxevbcXUsYOjr0s8RVRX+rHw2josXtsKj0I4EbGmprSxros4uYSOzgBa938aSfp2G3sB+RhEk4Cmo8wmGfx/Lh2F6kq/1Mk9/vXzcetTb0mHyZd5VNzxhTPx0F92xT2fyphnupvn3gCmlOA+gBySTfmgkx1rTVW51AB71e7XmQZWJWM84qLp9ags82DJ+h3weRQEwsbOP1bX3u/prpM3r1dAcSqhTiAhsGhtK1SF4gbCeBTC7EmjbJ3curmTELKZwxjSdUyqPRU/+vN7lseTGfNsdvPcG8CUCq47ACHESLfXkAucyC8nw8mOtbrSj0XT6yxxc00I1FSVy0NEz7bAqyY35oGwhgqfKr3eCeY5wXQqHiW+81cIgae3fohxwwZIndz+oychJH0Cfo9RrVM7qK/FmC+cVhcdkm5nrHk3zzDJcd0B9Aacyi8nw+mOddaEEZZpVebrmvcegSqpqko0/h4FEVkFw2grCmHao5tw22W10mHxyfAQkKgarVB8nkATwCMv7oLfo1jkGgynJ1DuTegy9qr4xdfOwyWjrdIULe1HsWSdVdZBBu/mGcYedgA5IFflg053rLMuHIGpYwfHTatq3nsEFT7VNpQSS1gHHr/5XHzr6bcBiOiO/9GNbUi3LUM2MsAuTxAI6/CqBL/HmJBlGu/6If0tpx8dIpo0NjG/j68uez0rZ8swjAE7gByQy/JBpztW83UrtuzB4nU74FONkMvVYwdjdfNHSa8t8yo4fCIEv6ogGFMK6lMVzLlkFB57uQ1eRcGJYNiRzk8iyUJOZR4Vj806B/3LfXFOzmm8vnX/p1ASxthwrT7DZAY7gByQj/JBJwnlFVv2YMGzRj7AHN7yfOsB+BQgmOIgMH7YAIvTCmo6xg0bgHVzJ2Hnx8fw3WeaoGXgAfQkom8hXUf9kP6Wz+Tk9LO6qR3zVm5P2r/AMIxz2AHkiFwmHJ3O8128ttVyrUchhEV8+WZiM9X5I6osidWTYQ2aruO2FW+hKxSGAOAgmiRFUcjWcSycVhf33chKVmWYeZZE4+9TwbX6DJMh7ABySC4Sjm0HjuGuPzQjqImU83y9qoKgFl87H9YFFk2vj2rfBzXdYjRfbetA24FjUafVuv8obnmyEYEYaeVs8KoKFOhxpaYAUO6NH9CSTuWULM9iwFIiDJMpPBGsgFjd1I5rfrbJkkQ1Y9yx1FSVQ5OUTi6aXodZF47A5vlT8NTsCbjn6jHSe5kTwaor/ehf7oNPtRm3lQKVYFH5tFPl1AWioZp0hddkeRbASDizYBvDZAY7gALBNIhByZB2WYw7dqRhhU+FTyV8/4tjo7pA5jzaSbUDpfeLnQhmZ1xNvCrhtstGSZ/7jy+OxbennAm/J16V8+GvjEdsU7BXJSyd2R2qMXf0cfeROLrEz+uTSEonu45hGHs4BFQg2IU4fCrZxrid5B1aP/rU0pkrm/wVmw/oCoVBRCjzqAhqGuZOPhNXjx2MX76yG7H+iQDcv3YHvCpBFwI3XTQcsyeNilPlbN3/KQBhSfxmUjk1Y/xQ1J3ezzglxSyEk8AMkxnsAAoEmUH0eRQ8d/ukOGOdDuapItb4+z0KvnP5aMtrrTOAT+L1vx/CE5s/wC9f+Tt+trENXz63Bqub2qGSgrCuQcBQFDXTBj9/eTdqBvSJCtNVV/qjjVyJZFo5VTuoLx6eyYJtDJML2AH0AE5KOu0MYjLjnyqJKguLEMG2Zr660o9NbYeiWkKmPpBp4J9p3Id7rx6DCaOqcbQriP+74i2EEpLQicJ0yci0csruOh7lyDDpwQ4gz6RT6ZKOQZTJT9y1cjsG9PGhfkg/VFf6UeFTLbo+J0O67bQsJ1pAD7/wHl6/53JHwnROyLRyKvG6bLWYGKYU4SRwHslkxKCZvE1lFGVJ1EBYx62/3RadZXs8qEU1/038KkmHpNu9ZyKxBn7R9DrL86YwXU/CoxwZJjPYAeSRdCtd0sGucudESIsawAqfCkqYH0AK2RpoY1cvdw4msQZ+1oQR+P4Xx8KnUtJZu/kmn98zw/RmSiIE5FZsOFONoHRzBgqRZZiKV1FwPKillWj9c8vHKbt/Ezt5E4Xp3Ii98yhHhsmMXu8AnIxZzJfxyqTSZcWWPRap5xnjh0bXWeEzhrfUVJVbO3ljpDlNAzhu2ABLXkH2mWN1hZIxrKqP9HO6mXTlUY4Mkxm92gGk0unvicRhOondOHG3SChm3qrtOHYyjCXrd0DoAgFNxM3wnTF+KC4ZfRqWzhxnawBjDbTsM0+sHSjVFZLxaVco4+8in/DwF4ZJn17tAJLp9APIeoiLU5zskDs6A1i8boflcQXA4nU74hqfzCqd2PU6MYB2DnHZ186T6grJuOMPTdCFKMgKG7dPIgxTbLiWBCai24loJxG1EtFD+bhHsthwrhOHHZ0BNO89knHlyb7DXVKZg6CmSx8HrOtNVUEk+8wqEXYfPB4d3xjLdeNOt+j8BMKsvcMwvQVXHAARTQZwHYBxQoh6AA/n4z6xejmmTo0ZGkkncZjKuK9uasfEBzfgpuVboyWY6VJTVS41wvOuGiN9PNl6k90j8TMfD2p4+IX3oOnGtK5YXaGf/Ou5ePzfzkcfb3zfgOl4snV6DMO4i1shoG8B+JEQIgAAQohP8nUju9CI08ShkyRyLkJJsetRFWOi1qLpdZg1YQQG9y/DvFXbpTmAdEMet11Wi0c3vg9PpEoIADoDxv/7PcDPbzo3Trenfkg/6AnKniFdR0v7UXx12evceMUwRQwJiaRw3m9K1ARgNYCpAE4CuFMI8abNa+cAmAMAw4cPP2/Pnj05XUuyKqCOzgAmPrghrjO2zKtg8/wp0dc27z2Cm5ZvjdPR7+v34KnZEzAuRnEz2/XIqoDSMf6xjiyo6bj+3KFY07w/avyTrXtNU3uco1x4bR2WrN+R9HthGKZwIKJtQoiGxMfzdgIgopcADJY8tSBy31MAXAjgfADPENEoIfFGQohlAJYBQENDQ869VbLEoZNh77muQbdbTzYJTtkpZdVb7UgcAG+37sRTlJPvhWGYwidvOQAhxBVCiLGS/1YD2Afgj8LgDQA6ALlwvYs4Me7J8gyFQEdnABt3fgJPQkewT1Uwd/KZjtcdm2DmxiuG6R24lQN4FsBkABuJaDQAH4BDLq3FFqd5gkKtQTfDPrHKniYhXceNE4bjxgnD0143N14xTO/ArRyAD8ATAMYDCMLIAWxIdV1DQ4NobGzM9/IsFKPMsCx/AQAVfhWaLnKStC3G74VhSpEezwEkQwgRBHCTG/fOBDcajLIxrmbYR6X4sE+FT8Xi6fWYPOa0nHwebrximOKmV3cCFyvZSFSY13oUa9hHEyJnxp9hmOKH5aALjGy07WOvjS3vdFOqmWGYwoVPAAVGNiWWsmsr/LkN+zAM03vgE0CBkU2JpexaTeewD8MwctgB5JFMtHKy6Sso9J4EhmEKC1fKQDPFrTLQTMh21kC2VUBcnskwjElBlYH2dnIhEJdNiSWXZzIM4wQOAeWBfYe7LNILPKScYZhCgx1AHmhpPxpXhgmwVg7DMIUHO4Ac09EZwJL11tGOC6fVcViGYZiCgh1AjpGNXazwqRg7pL9LK2IYhpHDDiDHSGvxheDwD8MwBQc7gBzDtfgMwxQLXAaaBwp1PgDDMEws7ADyBNfiMwxT6HAIKEdkIvvAMAzjJnwCyAHZyj4wDMO4gSsnACIaT0RbiKiJiBqJ6AI31pELstHvZxiGcRO3QkAPAVgshBgP4L7Iz0WJrO6fZR8YhikG3HIAAkC/yJ/7A9jv0jqyJhv9foZhGDdxywH8O4ClRLQXwMMA7rF7IRHNiYSJGg8ePNhjC3QK1/0zDFOs5G0eABG9BGCw5KkFAC4H8DchxCoiugHAHCHEFanes5DnAbAGP8MwhYrdPABXBsIQ0VEAA4QQgogIwFEhRL9U1xWyA2AYhilU7ByAWyGg/QAujfx5CoD3XVpHwcJ9BQzD5Bu3+gBuAfATIvIAOAlgjkvrKEi4r4BhmJ7AFQcghNgE4Dw37l3o5GKcJMMwjBNYCqLA4L4ChmF6CnYABQb3FTAM01OwAygwuK+AYZiegsXgChCeJ8AwTE/ADqBA4XkCDMPkGw4BMQzDlCjsABiGYUoUdgAMwzAlCjsAhmGYEoUdAMMwTIniihpophDRQQB78nybgQAO5fke+YDX3fMU69qLdd1A8a7d7XWPEEKcmvhgUTmAnoCIGmWyqYUOr7vnKda1F+u6geJde6Gum0NADMMwJQo7AIZhmBKFHYCVZW4vIEN43T1Psa69WNcNFO/aC3LdnANgGIYpUfgEwDAMU6KwA2AYhilR2AEkQERLiGg7ETUR0QtENMTtNTmFiJYS0c7I+v9ERAPcXpMTiOgrRNRKRDoRFVypXCJENJWI3iOiNiK62+31OIWIniCiT4ioxe21pAMRDSOijUS0I/J78h231+QUIiojojeIqDmy9sVurykWzgEkQET9hBCfRv78bQB1QohbXV6WI4joCwA2CCHCRPQgAAgh5ru8rJQQ0WcB6AB+CeBOIUSjy0uyhYhUALsAXAlgH4A3AfyrEGKHqwtzABFdAqATwJNCiLFur8cpRHQ6gNOFEG8RUV8A2wB8sUi+cwJQIYToJCIvgE0AviOE2OLy0gDwCcCCafwjVAAoGg8phHhBCBGO/LgFQI2b63GKEOJdIcR7bq/DIRcAaBNC7BZCBAH8L4DrXF6TI4QQrwD4p9vrSBchxEdCiLcifz4G4F0AQ91dlTOEQWfkR2/kv4KxKewAJBDR94loL4BZAO5zez0Z8k0Az7u9iF7IUAB7Y37ehyIxRr0BIhoJ4BwAW91diXOISCWiJgCfAHhRCFEway9JB0BELxFRi+S/6wBACLFACDEMwAoAc91dbTyp1h55zQIAYRjrLwicrJthkkFElQBWAfj3hJN6QSOE0IQQ42GcyC8gooIJv5XkSEghxBUOX7oCwHMAFuVxOWmRau1E9HUA0wBcLgoowZPGd17otAMYFvNzTeQxJo9E4uerAKwQQvzR7fVkghDiCBFtBDAVQEEk4kvyBJAMIjoz5sfrAOx0ay3pQkRTAcwDMEMIccLt9fRS3gRwJhGdQUQ+AP8CYI3La+rVRBKpjwN4VwjxY7fXkw5EdKpZjUdE5TCKBwrGpnAVUAJEtArAWTCqUvYAuFUIURQ7PCJqA+AH0BF5aEsxVDAR0ZcA/AzAqQCOAGgSQlzl7qrsIaJrAPwXABXAE0KI77u8JEcQ0e8AXAZDmvgAgEVCiMddXZQDiGgSgFcBvAPj3yUA3CuEeM69VTmDiM4G8BsYvysKgGeEEA+4u6pu2AEwDMOUKBwCYhiGKVHYATAMw5Qo7AAYhmFKFHYADMMwJQo7AIZhmBKFHQDDMEyJwg6AYRimRGEHwDBZQETnR+YvlBFRRUTzvWC0XhgmGdwIxjBZQkT/AaAMQDmAfUKIH7q8JIZxBDsAhsmSiCbQmwBOAvi8EEJzeUkM4wgOATFM9lQDqATQF8ZJgGGKAj4BMEyWENEaGJPBzoAxurCgZkgwjB0lOQ+AYXIFEd0MICSEeDoyL/g1IpoihNjg9toYJhV8AmAYhilROAfAMAxTorADYBiGKVHYATAMw5Qo7AAYhmFKFHYADMMwJQo7AIZhmBKFHQDDMEyJ8v8B2AEct3PfKsUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from jax import random\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "x=random.normal(random.PRNGKey(30), shape=(1000,))\n",
    "noise = jax.random.normal(random.PRNGKey(42), (1000,))\n",
    "true_w, true_b = 2, -1\n",
    "ys = x * true_w + true_b + noise\n",
    "data = pd.DataFrame(dict(x=x, y=ys))\n",
    "data.plot(x='x',y='y',kind=\"scatter\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PERFORMING BAYESIAN LINEAR REGRESSOR\n",
    "\n",
    "## setting priors \n",
    "* Now we have imported PYMC library for Bayesian statistical modeling and probabilistic machine learning which focuses\n",
    "  on advanced Markov chain Monte Carlo and variational fitting algorithms. \n",
    "\n",
    "* We have defined function xyz which takes parameter weight of feature matrix. ALso we have defined priors and created variable \n",
    "  sigma (error term) that represents half couchy distribution (halfcauchy is a cauchy distribution i.e truncated to only have     nonzero probability density for values greater than or equal to the location of the peak) and has 3 parameters that are name,   testvalue and beta value.\n",
    "* Now we have defined prior intercept with mean value equal to zero and represent it as normal distribution and standard error   or std.deviation equals to 20.\n",
    "* Again we represent x_coeff as normal distribution with mean equal to zero and std. deviation equal to 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## setting likelihood\n",
    "* we set prior term value as weight square and divided by 2\n",
    "*  We set likelihood and give label as 'y' and is normal distribution and paramters passed are mean (value is intercept +         coefficient_of_x * ( x))  and we finally return updated weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "hCV58dt7gjNV"
   },
   "outputs": [],
   "source": [
    "import pymc3             as pm\n",
    "\n",
    "def xyz(w):\n",
    "  with pm.Model() as model:  # model specifications in PyMC3 are wrapped in a with-statement\n",
    "      \n",
    "      # Define priors\n",
    "      sigma = pm.HalfCauchy(\"sigma\", beta=10, testval=1.0)\n",
    "      intercept = pm.Normal(\"Intercept\", 0, sigma=20)\n",
    "      x_coeff = pm.Normal(\"x\", 0, sigma=20)\n",
    "\n",
    "      # Define likelihood\n",
    "      prior_term = w @ w / 2\n",
    "      likelihood = pm.Normal(\"y\", mu=intercept + x_coeff * x, sigma=sigma, observed=y)\n",
    "      return -prior_term + likelihood.sum()\n",
    "\n",
    "      # Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We are downloading blackjax from package installer from python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "W9OBFMoXkxZx",
    "outputId": "250ffe3b-1f86-48cf-a8eb-2d5fd391bcd7"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: blackjax in /usr/local/lib/python3.7/dist-packages (0.4.0)\n"
     ]
    }
   ],
   "source": [
    "pip install blackjax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* After importing blackjax, jax and pymc3 ,rng_key is random seed, w0 is a initial weight, rmh is a sampling method of blackjax\n",
    "\n",
    "* We use blackjax's Random Walk RMH kernel to sample from the posterior distribution.\n",
    "* Since blackjax does not provide an inference loop we need to implement one ourselves , so we have made function \n",
    "  inference_loop which takes paramters rng_key, initial_state, num_samples. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fF45EzSRkY0O"
   },
   "outputs": [],
   "source": [
    "import blackjax\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "import pymc3             as pm\n",
    "from jax import random\n",
    "\n",
    "rng_key = random.PRNGKey(314)\n",
    "\n",
    "w0 = random.multivariate_normal(rng_key, 0.1 + jnp.zeros(1), jnp.eye(1))\n",
    "\n",
    "rmh = blackjax.rmh(xyz, sigma=jnp.ones(1) * 0.7)\n",
    "initial_state = rmh.init(w0)\n",
    "\n",
    "def inference_loop(rng_key, kernel, initial_state, num_samples):\n",
    "    @jax.jit\n",
    "    def one_step(state, rng_key):\n",
    "        state, _ = kernel(rng_key, state)\n",
    "        return state, state\n",
    "\n",
    "    keys = jax.random.split(rng_key, num_samples)\n",
    "    _, states = jax.lax.scan(one_step, initial_state, keys)\n",
    "\n",
    "    return states"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We can now run the inference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 502
    },
    "id": "IZqybNo_EXsO",
    "outputId": "5de74aaf-c23d-44ab-9947-d9fa047eb877"
   },
   "outputs": [],
   "source": [
    "_, rng_key = random.split(rng_key)\n",
    "states = inference_loop(rng_key, rmh.step, initial_state, 10)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Untitled1.ipynb",
   "provenance": []
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}