example/05-jeffreys1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit a single variable logistic regression model with Jeffreys prior\n",
    "to simulated data. Plot out both the prior and posterior distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.integrate as integrate\n",
    "\n",
    "from bbai.glm import LogisticRegressionMAP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_correlation_matrix(p, param):\n",
    "    res = np.zeros(shape=(p, p))\n",
    "    for s in range(p):\n",
    "        for t in range(0, s+1):\n",
    "            corr = param\n",
    "            if s == t:\n",
    "                corr = 1.0\n",
    "            res[s, t] = corr\n",
    "            res[t, s] = corr\n",
    "    return res\n",
    "\n",
    "def generate_design_matrix(n, K):\n",
    "    mean = np.zeros(K.shape[0])\n",
    "    return np.random.multivariate_normal(mean, K, size=n)\n",
    "\n",
    "def generate_weights(p):\n",
    "    return np.random.normal(size=p)\n",
    "\n",
    "def generate_data_set(n, K):\n",
    "    p = K.shape[0]\n",
    "    X = generate_design_matrix(n, K)\n",
    "    w = generate_weights(p)\n",
    "\n",
    "    u = np.dot(X, w)\n",
    "\n",
    "    p = 1 / (1 + np.exp(-u))\n",
    "\n",
    "    y = []\n",
    "    for i in range(n):\n",
    "        y.append(np.random.binomial(1, p[i]))\n",
    "    y = np.array(y)\n",
    "\n",
    "    return X, y, w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Jeffreys Prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_a_matrix(X, u):\n",
    "    p_vector = 1 / (1 + np.exp(u))\n",
    "    return np.diag(p_vector * (1 - p_vector))\n",
    "\n",
    "def compute_fisher_information_matrix(X, u):\n",
    "    A = compute_a_matrix(X, u)\n",
    "    return np.dot(X.T, np.dot(A, X))\n",
    "\n",
    "def compute_log_prior(X, u):\n",
    "    FIM = compute_fisher_information_matrix(X, u)\n",
    "    return 0.5 * np.linalg.slogdet(FIM)[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_log_posterior(X, y, w):\n",
    "    u = np.dot(X, w)\n",
    "    log_prior = compute_log_prior(X, u)\n",
    "    y = 2 * y - 1\n",
    "    likelihood = -np.sum(np.log(1 + np.exp(-y * u)))\n",
    "    return likelihood + log_prior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "\n",
    "n = 20\n",
    "p = 1\n",
    "\n",
    "K = generate_correlation_matrix(p, 0.5)\n",
    "X, y, w_true = generate_data_set(n, K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "wx = np.arange(-20, 20, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_prior(w):\n",
    "    u = np.dot(X, np.array([w]))\n",
    "    return np.exp(compute_log_prior(X, u))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAs5UlEQVR4nO3deXzc9X3n8ddH931atmRbPmSMjR1uYyBNSICWawkODTTQ3S1JaUm7pWke3TRLt90cdLtbmjak25LdkkCWJM1Clk2Km0AJVwLh9AHYGN8Xkq3LumXd0mf/mJ9kRYytsTU/zU/S+/l4zEO/+R0zH49n5jPf29wdERGRidJSHYCIiESTEoSIiMSlBCEiInEpQYiISFxKECIiEldGqgNIlnnz5vmyZctSHYaIyIyyZcuWY+5eEe/YrEkQy5YtY/PmzakOQ0RkRjGzwyc7piomERGJSwlCRETiUoIQEZG4lCBERCQuJQgREYlLCUJEROJSghARkbiUIERCtr+5mx9urUNT68tMM2sGyolEUW1rD5/4n6/Q3jPIoWPH+eNrVqU6JJGEqQQhEqLHNtXS2TvIh86ax8MvH6JvcDjVIYkkTAlCJCTuzhNvH+FXzprH7390Bd39Qzy3synVYYkkTAlCJCQHjx2ntrWX6z5QyWU15RTnZvLS3uZUhyWSMCUIkZC8+V47AOuWlpGeZly4pGRsn8hMoAQhEpI3a9sozM5g5fwCAC6sLmVPUxedfYMpjkwkMUoQIiF5u7aD86qLSUszAC5cUoI7bK/rSHFkIolRghAJwfCIs6exi3Mqi8b2ra4qBGBPY1eqwhI5LUoQIiE40tZL/9AIKxcUjO2rKMimODeTvU3dKYxMJHFKECIhGC0lnDW/cGyfmbFyfgH7GpUgZGZQghAJwWgpYXwJIna/kD1NXZp2Q2aEUBOEmV1nZrvNbJ+Z3RPn+BVmttXMhszslnH7LzCzV81sh5ltM7NPhhmnSLIdaO6mojCbopzMX9q/oiKf9p5B2nvUk0miL7QEYWbpwAPA9cAa4HYzWzPhtPeATwHfn7C/B/gtd18LXAd83cxKwopVJNlq23pYUpb3vv3Vwb7atp7pDknktIVZglgP7HP3A+4+ADwKbBh/grsfcvdtwMiE/XvcfW+wfRRoAipCjFUkqWpbe6kuzX3f/urSvLHjIlEXZoJYBNSOu18X7DstZrYeyAL2xzl2l5ltNrPNzc2awkCiYXB4hPqO3rHSwnjVZbGkoRKEzASRbqQ2syrgu8Cn3X1k4nF3f9Dd17n7uooKFTAkGurb+xjxE6WF8QpzMinJy6S2VQlCoi/MBHEEqB53f3GwLyFmVgT8BPgzd38tybGJhGa0dLC47P1VTBBLHLVtqmKS6AszQWwCVprZcjPLAm4DNiZyYXD+j4DvuPvjIcYoknR1QYKIV4KAWDVTnUoQMgOEliDcfQi4G3ga2An8wN13mNm9ZnYTgJldYmZ1wK3AP5rZjuDy3wCuAD5lZm8FtwvCilUkmWpbe0lPM6qKc+IeX1yaR11bLyMjGgsh0RbqkqPu/iTw5IR9Xxy3vYlY1dPE674HfC/M2ETCUtvWQ1VxDhnp8X9/VZfmMjA8QlNXP5UnSSIiURDpRmqRmai2teek1UsAizUWQmYIJQiRJKtt6x3rzhrPibEQShASbUoQIknUNzhMc1f/qUsQwQA6DZaTqFOCEEmio+2xL/1FcUZRj8rJTKeiMHvsXJGoUoIQSaKGzj6ASRufK4tyxs4ViSolCJEkaugIEkTRJAmiOGfsXJGoUoIQSaLTKUHUd6iKSaJNCUIkiRo7+ijKySAv69RDjCqLc+jsG6JnYGiaIhM5fUoQIklU39GX0OC30VHWqmaSKFOCEEmixs4+FkzS/gAn2ijUUC1RpgQhkkQNnX0nnYNpvEqVIGQGUIIQSZKh4RGau/on7cEE4xKEShASYUoQIknS3N3PiMOCBEoQeVkZFOVkqAQhkaYEIZIko1/2iVQxxc7LpV4JQiJMCUIkSUYTRCKN1BAraTSqikkiTAlCJEnGBsklmCCqinJUgpBIU4IQSZKGzj6y0tMoy89K6PwFxTkc6+5ncHgk5MhEzowShEiSNHT0saA4GzNL6Pyq4hzcoamrP+TIRM6MEoRIkjR09CVcvQQaCyHRpwQhkiSNnX1UFp98HYiJxkZTK0FIRClBiCSBu8fmYSrKTvia0e6wmtVVokoJQiQJOnoH6R8aSbiLK0BxbibZGWnq6iqRpQQhkgSjXVyrTqOKycyoKlZXV4muUBOEmV1nZrvNbJ+Z3RPn+BVmttXMhszslgnH7jCzvcHtjjDjFJmq0S/5yuLEq5hi52uwnERXaAnCzNKBB4DrgTXA7Wa2ZsJp7wGfAr4/4doy4EvApcB64EtmVhpWrCJT1Xiao6hHVRblcLRdCUKiKcwSxHpgn7sfcPcB4FFgw/gT3P2Qu28DJo4UuhZ4xt1b3b0NeAa4LsRYRaakvqMPM5hfeJoJojiXpq4+RkY8pMhEzlyYCWIRUDvufl2wL2nXmtldZrbZzDY3NzefcaAiU9XY2Ud5fjZZGaf3kaoqzmFw2GntGQgpMpEzN6Mbqd39QXdf5+7rKioqUh2OzGENnX2n3f4AJ6qkNBZCoijMBHEEqB53f3GwL+xrRaZdbBR14j2YRp0YC6EEIdETZoLYBKw0s+VmlgXcBmxM8NqngWvMrDRonL4m2CcSSWdagqjSynISYaElCHcfAu4m9sW+E/iBu+8ws3vN7CYAM7vEzOqAW4F/NLMdwbWtwF8QSzKbgHuDfSKR0zc4THvP4GnNwzSqvCCb9DSjQaOpJYIywnxwd38SeHLCvi+O295ErPoo3rUPAw+HGZ9IMjSMjYE4/Sqm9DRjQWG2qpgkkmZ0I7VIFJzuQkETVRbnqJFaIkkJQmSKRkdCn0kbROy6HLVBSCQpQYhMUf0UqpgAKotyaejow12D5SRalCBEpqiho4+C7AwKss+sSa+qOIeegWE6+4aSHJnI1ChBiExRY2cfC05jHYiJFgRdXTVpn0SNEoTIFNV39J3WNN8TabCcRJUShMgUxUoQZ9aDCcYvPaqxEBItShAiUzA0PEJjZ99YKeBMnJiPqT9ZYYkkhRKEyBQ0dfUz4lBVcuYJIisjjXkFWTR0qgQh0aIEITIF9UG10MIptEFAbCyE2iAkapQgRKZgdDW4qZQg4MRYCJEoUYIQmYLREsRUejFBbBS2RlNL1ChBiEzB0fY+8rPSKcqZ2ryXVcW5tPcM0jswnKTIRKZOCUJkCuo7ellYkouZTelxxrq6qhQhEaIEITIF9R19VJVMrXoJYo3UoKVHJVqUIESm4Gh7HwunMAZi1FiCUFdXiRAlCJEz1D80zLHu/ik3UMOJKiZ1dZUoUYIQOUONwcjnqXZxBcjPzqAwJ4NGJQiJECUIkTN0NEmD5EZVabCcRIwShMgZGhsDkYQSBMQWHFIvJokSJQiRMzQ6ijppJYiinLHHFIkCJQiRM1Tf0UtJXia5WelJebyFJbkc6+6nb1CD5SQalCBEzlB9+9QWCppoUWnssdQOIVGRUIIws4+Z2WknEzO7zsx2m9k+M7snzvFsM3ssOP66mS0L9mea2SNmtt3MdprZn57uc4uE7WhHcsZAjFoUDLg70qaxEBINiX7pfxLYa2Z/bWarE7nAzNKBB4DrgTXA7Wa2ZsJpdwJt7n4WcD9wX7D/ViDb3c8FLgY+M5o8RKKivqM3aQ3UAIuDEsSR9p6kPabIVCSUINz93wEXAvuB/21mr5rZXWZWeIrL1gP73P2Auw8AjwIbJpyzAXgk2H4cuNpik9o4kG9mGUAuMAB0JvqPEglb78Aw7T2DSa1iqizOIc2gTiUIiYiEq43cvZPYl/ijQBVwM7DVzP7wJJcsAmrH3a8L9sU9x92HgA6gPHie40A98B7wN+7eOvEJgiS12cw2Nzc3J/pPEZmysTEQSSxBZKanUVmUoyomiYxE2yA2mNmPgJ8BmcB6d78eOB/4jyHEtR4YBhYCy4H/aGY1E09y9wfdfZ27r6uoqAghDJH4Rn/lLyrJS+rjLirNpa5dCUKiIdFJ7H8duN/dXxy/0917zOzOk1xzBKged39xsC/eOXVBdVIx0AL8JvCv7j4INJnZy8A64ECC8YqEqrY11k5QXZa8KiaINVRvOtSW1McUOVOJVjE1TEwOZnYfgLs/d5JrNgErzWy5mWUBtwEbJ5yzEbgj2L4FeN7dnVi10lXB8+QDlwG7EoxVJHS1bT1kphsLCpNXxQSxEkRDZx9DwyNJfVyRM5Fogvi1OPuuP9UFQZvC3cDTwE7gB+6+w8zuNbObgtMeAsrNbB/wx8BoV9gHgAIz20Es0Xzb3bclGKtI6OraellUkkta2tQWCppoUUkewyNOY1d/Uh9X5EycsorJzH4f+A/ACjMb/wVdCLw82YO7+5PAkxP2fXHcdh+xLq0Tr+uOt18kKupae6guS277A5wYLHckSEAiqTRZG8T3gaeA/86JX/cAXfF6FYnMFbVtvVy7sDjpjzs2WK69ByhL+uOLnI7JEoS7+yEz+4OJB8ysTElC5qLj/UO0Hh9IegM1nBgsV9uqnkySeomUIG4EthAbvDa+wtWB93U9FZntatuCHkylya9iyslMZ0FRNodbNJpaUu+UCcLdbwz+Lp+ecESib/TX/eiv/WRbWpY/1o1WJJUma6S+6FTH3X1rcsMRib4TYyCSX4IAWFKex0t7NTOApN5kVUx/e4pjTjBWQWQuqWvrJTcznfL8rFAef0lZHo2dsXUhcjKTs9aEyJmYrIrpyukKRGSmqG3robosl9i8ksm3tDxWMqlt7WHlglPNhykSrsmqmK5y9+fN7NfjHXf3H4YTlkh0HW45zpKy/NAef0lQdXW4RQlCUmuyKqaPAM8DH4tzzAElCJlTRkacwy09fOTs8CaHXFoeSz6H1VAtKTZZFdOXgr+fnp5wRKKtvrOP/qERls0LrwRRmpdJQXYG77UcD+05RBKR6HTf5Wb2P8xsq5ltMbO/M7PysIMTiZpDx2Jf2svKw0sQZsaSsjzeUwlCUizRyfoeBZqBTxCbdbUZeCysoESi6uBoggixBAGxhmpVMUmqJZogqtz9L9z9YHD7r8CCMAMTiaJDx46TnZFGVVFyp/meaGl5PnWtvQyPeKjPI3IqiSaIn5rZbWaWFtx+g9g03iJzyqGWHpaW5yV9mu+JaublMzA8Ql2bShGSOpN1c+3ixBxMnwO+FxxKA7qBz4cZnEjUHGo5Tk3I1UsANRWx59jf3D3Wq0lkup2yBOHuhe5eFPxNc/eM4Jbm7kXTFaRIFAyPOO+19LB8GhLEiooCAA40qyeTpE6ia1JjZqXASmCs8nXiMqQis9nR9l4GhsPt4jqqND+L0rxM9itBSAollCDM7HeAPwIWA28RWyP6VTQXk8whB6ehi+t4NRUF7G/unpbnEokn0UbqPwIuAQ4H8zNdCLSHFZRIFO1rin1Zr1xQMC3Pt6IiX1VMklKJJoi+YP1ozCzb3XcBq8ILSyR69jZ1U5qXGdosrhPVVBRwrLufzr7BaXk+kYkSTRB1ZlYC/DPwjJk9ARwOKyiRKNrX1MVZ8wtCm8V1otHeUipFSKok1Abh7jcHm182sxeAYuBfQ4tKJGLcnT2N3dxwbtW0PWdN0JNpf1M3F1SXTNvzioxKtASBmV1kZp8FzgPq3H0ggWuuM7PdZrbPzO6JczzbzB4Ljr9uZsvGHTvPzF41sx1mtt3Mwh26KnIKx7oH6OgdZOX86Wl/gNh0G5npxj41VEuKJDpZ3xeBR4ByYB7wbTP780muSQceAK4H1gC3m9maCafdCbS5+1nA/cB9wbUZxAbl/Z67rwU+CqgiVlJmb1MXMH0N1ACZ6WnUzCtgd0PXtD2nyHiJliD+LXCJu38pmAL8MuDfT3LNemCfux8IShuPAhsmnLOBWOIBeBy42mIVvNcA29z9bQB3b3H34QRjFUm6sR5M86d3AZ9VlYVKEJIyiSaIo4wbIAdkA0cmuWYRUDvufl2wL+457j4EdBArpZwNuJk9HUwx/oV4T2Bmd5nZZjPb3NysRd4lPHsbuynMzmBBUfa0Pu+qykKOtPfSpZ5MkgKTzcX098TmYuoAdpjZM8H9XwPeCDmuDxEbe9EDPGdmW9z9ufEnufuDwIMA69at07SXEpq9TV2smMYeTKNWBUuO7mns4uKlZdP63CKT9WLaHPzdAvxo3P6fJfDYR4DqcfcX8/5Sx+g5dUG7QzHQQqy08aK7HwMwsyeBi4DnEJlm7s7O+i5uOLdy2p97VWUsQexqUIKQ6TfZkqOj7QOYWRaxqh+A3e4+WZl3E7DSzJYTSwS3Ab854ZyNwB3Epu24BXje3d3Mnga+YGZ5wACxtbHvT+yfJJJc9R19dPQOsqZq+uenXFyaS35WutohJCUSnYvpo8Qakw8Rm/q72szuONVkfe4+ZGZ3E1s3Ih142N13mNm9wGZ33wg8BHzXzPYBrcSSCO7eZmZfI5ZkHHjS3X9yZv9Ekal592gnAGsWTn+CMDPOVkO1pEiis7n+LXCNu+8GMLOzgf8DXHyqi9z9SeDJCfu+OG67D7j1JNd+jxPrT4ikzLv1nZjBqsrUzHC/urKQp95pwN2nvQ1E5rZEezFljiYHAHffA2SGE5JItLx7tJNl5fkUZCc8O35SnVNVRHvPIEc7+lLy/DJ3JZogtpjZt8zso8Htm5xowBaZ1XY2dKak/WHUBxYVA7C9riNlMcjclGiC+D3gXeCzwe1d4PfDCkokKrr6Bjnc0pOS9odRa6qKSE8z3jmiBCHTa9IyczBlxtvuvhr4WvghiUTHrqBxOJUliJzMdFbOL2C7EoRMs0lLEMEUF7vNbMk0xCMSKW/XtgOwdlFql2A/d1Ex7xzpwF3jQWX6JFrFVEpsJPVzZrZx9BZmYCJR8HZdB4tKcplfmNrJhM9dXEzL8QHq1VAt0yjRbhn/JdQoRCLqrdo2zq8uTnUYJxqqj3SwsCQ3xdHIXHHKEoSZ5ZjZ54iNVVgNvOzuPx+9TUeAIqnS0t1PbWtvJBbrGW2o3lbXnupQZA6ZrIrpEWAdsJ3Yug5/G3pEIhHxVtD+cP7ikpTGAbGG6jVVRWw93J7qUGQOmayKaY27nwtgZg8R7gyuIpHydm07aRar/4+Ci5eW8timWoaGR8hIT3gxSJEzNtm7bGxCvmC9BpE5483ads5eUEheVmpGUE908dJSegeH2VmveZlkekyWIM43s87g1gWcN7ptZp3TEaBIKgyPOG++185FS0tTHcqYi4NYthxuTXEkMlecMkG4e7q7FwW3QnfPGLed2o7hIiHaWd9Jd/8Qly6PzhoMC0tyWVicw+bDbakOReYIVWSKxPHagRYA1kcoQQBctLSUrUoQMk2UIETieONgK0vK8qgqjtaYg0uWlXG0o4/a1p5UhyJzgBKEyAQjI84bh1ojVb006vIV5QC8sv9YiiORuUAJQmSCvU3dtPcMRq56CWDl/ALmFWTzyv6WVIcic4AShMgEo7/OL6spT3Ek72dmfHBFOa/sb9HEfRI6JQiRCV7c08yy8jyqy/JSHUpcv3JWOc1d/exr6k51KDLLKUGIjNM/NMxrB1q54uyKVIdyUh9cMQ+Al/epHULCpQQhMs7mQ230Dg7zkQgniOqyPJaW5/HiXiUICZcShMg4L+5pJjPdItn+MN6Vq+bz8r5j9A4MpzoUmcWUIETG+fmeZtYtLSM/OxrzL53M1efMp39oRN1dJVShJggzu87MdpvZPjO7J87xbDN7LDj+upktm3B8iZl1m9nnw4xTBKC2tYddDV1ctXp+qkOZ1PrlZeRnpfPcrqZUhyKzWGgJwszSgQeIrSOxBrjdzNZMOO1OoM3dzwLuB+6bcPxrwFNhxSgy3tM7GgC4dm1liiOZXHZGOh9eWcELu5rU3VVCE2YJYj2wz90PuPsA8CiwYcI5G4gtSgTwOHC1mRmAmX0cOAjsCDFGkTE/3dHI6spClpRHs3vrRFedM5/6jj62H+lIdSgyS4WZIBYBtePu1wX74p4TrDfRAZSbWQHwn4CvnOoJzOwuM9tsZpubm5uTFrjMPce6+9l0uHVGlB5GXbNmARlpxk+21ac6FJmlotpI/WXgfnc/5Uggd3/Q3de5+7qKiuh2S5Toe+bdRtzhmrULUh1KwkrysvjQynn8eFu9qpkkFGEmiCNA9bj7i4N9cc8xswygGGgBLgX+2swOAZ8D/rOZ3R1irDLH/fObR6ipyGdN1cxa5uTfnFvFkfZe3q5TNZMkX5gJYhOw0syWm1kWcBuwccI5G4E7gu1bgOc95sPuvszdlwFfB/6bu/9DiLHKHHakvZfXD7Zy8wWLCJrAZoxr1laSmW78y9tHUx2KzEKhJYigTeFu4GlgJ/ADd99hZvea2U3BaQ8Ra3PYB/wx8L6usCJh2/hW7Mt1wwUTm8iirzg3kytXzeeJt44wODyS6nBklgl1NJC7Pwk8OWHfF8dt9wG3TvIYXw4lOBHA3fnRm3VcvLR0xvRemujWddX89N1Gfra7mV9bM3PaUCT6otpILTIttr7Xzp7Gbj5x0eJUh3LGPrqqgnkFWfxgc+3kJ4ucBiUImdO+//p75Gelc9MFC1MdyhnLTE/j1y9azAu7mmjq6kt1ODKLKEHInNXRM8iPtx3l4xcuoiDicy9N5pOXVDM04jz6hkoRkjxKEDJn/d8ttfQPjfCbly5JdShTtqKigCvOruB7rx1mYEiN1ZIcShAyJw0Nj/Dtlw+xflkZaxcWpzqcpPj0B5fR1NXPU+9oZLUkhxKEzEk/2V7PkfZe7rqiJtWhJM1Hzq5g+bx8vvXSQY2slqRQgpA5x9355ksHWFGRPyOm9k5UWprxmStq2H6kQ6vNSVIoQcic8/yuJt450slnrlhBWtrMGjk9mV+/aDFVxTk88Py+VIcis4AShMwp7s7XntnDkrI8br5o5o2cnkxWRhp3XVHDG4daeWWfShEyNUoQMqc8vaOBHUc7+ezVK8lMn51v/9vXL2FRSS5/+eRORkbUFiFnbnZ+QkTiGBga4a+e2sWKinw+PoMHxk0mJzOdL1y3ih1HO3ni7YkTKIskTglC5ozvvHqIQy09/PmNa8iYpaWHUR87byHnLirmq/+6m77B4VSHIzPU7P6UiASauvr4u+f28pGzK7hy1ezpuXQyaWnGf77hHI529PHgiwdSHY7MUEoQMid86Ykd9A+N8MWPrUl1KNPm8hXl3HheFf/w/D72NXWlOhyZgZQgZNZ7ans9T73TwOd+dSUrKgpSHc60+vJNa8nLTucLj29jWA3WcpqUIGRWa+8Z4Isbd/CBRUXc9eHZM2o6UfMKsvnSx9aw9b12HnnlUKrDkRlGCUJmLXfnTx7fRnvPAPd94rxZ3zB9Mh+/YBFXrqrgr/51FzuOau1qSdzc/MTInPDQLw7yzLuN/On158yaCfnOhJnx1VvPpzQvkz/4p6109Q2mOiSZIZQgZFbacriVv3pqF9etreTTv7Is1eGk3LyCbP7+9ouobevlnh9u12R+khAlCJl1alt7+Mx3t7CoNJf7bjkPs9k139KZWr+8jD+5dhU/2VbP3z23N9XhyAwws5fREpmgvWeAO779BoPDzsOfuoTi3MxUhxQpn7mihn1N3Xz92b0sLs3jlotn7lrcEj4lCJk1egaG+N3vbKautZfv/c6lc65LayLMjP9287k0dPRxz//bRnlB1pwYOChnRlVMMiv0DAzx6W9vYsvhNr72yfNZv7ws1SFFVlZGGt/4dxexuqqQz3x3Cz/f05zqkCSiQk0QZnadme02s31mdk+c49lm9lhw/HUzWxbs/zUz22Jm24O/V4UZp8xs3f1D/Pb/3sSmQ63c/8kLuPG82TsRX7IU5WTyvTsvZeX8An73O5t5YXdTqkOSCAotQZhZOvAAcD2wBrjdzCbOc3An0ObuZwH3A/cF+48BH3P3c4E7gO+GFafMbE2dfXzyH19l06E27v/kBWy4YPat8RCWkrws/ul3giTxyGb+7+baVIckERNmCWI9sM/dD7j7APAosGHCORuAR4Ltx4Grzczc/U13Pxrs3wHkmll2iLHKDLS3sYubv/EKB48d51t3rFNyOAMleVk8etdlXL6inD95fBtff3aP1pCQMWEmiEXA+J8kdcG+uOe4+xDQAZRPOOcTwFZ375/4BGZ2l5ltNrPNzc2qR51LfrztKBseeJn+oREevesyNbROQWFOJg/dcQmfuGgxX392L7/7nc109GgwnUS8kdrM1hKrdvpMvOPu/qC7r3P3dRUVFdMbnKRE3+AwX/mXHdz9/TdZXVnIj//wQ5y3uCTVYc14WRlp/M2t5/GVm9by4t5mbvyHl3jniKblmOvCTBBHgOpx9xcH++KeY2YZQDHQEtxfDPwI+C133x9inDJDvHOkg4/9/S/49suH+NQHl/HoXZdTWZyT6rBmDTPjjg8u47HPXM7wsHPzN17mH57fy+DwSKpDkxQJM0FsAlaa2XIzywJuAzZOOGcjsUZogFuA593dzawE+Alwj7u/HGKMMgP0DQ7z9Wf38PEHXqajd5BHfns9X75pLVkZkS4Az1gXLSnlx5/9MNeureRvfrqHm7/xMrsaOlMdlqSAhTkni5ndAHwdSAcedve/NLN7gc3uvtHMcoj1ULoQaAVuc/cDZvbnwJ8C4+cDuMbdT9oXb926db558+aw/imSIi/sbuLLG3dwuKWHm85fyL0b1lKSl5XqsOaMp7bX81+eeIf2nkHu+OAyPnv1So1On2XMbIu7r4t7bLZM2qUEMbvsbujiq0/v4tmdTdTMy+crG9by4ZVqZ0qF1uMDfPXpXTy6qZayvCw+f+0qbr148ZydPn22UYKQGePQseN8/dk9PPH2UQqyMvj9K1dw54eWk52RnurQ5rx3jnTwlX/ZwaZDbSwrz+OzV69kwwWLSE/TZIgzmRKERN6exi6+9dIBfrj1CBnpxqc+uJzf+0iNqpMixt15dmcTX3tmDzvrO6mpyOcPrzqLG89bSKZKFDOSEoREkrvz6v4WvvnSAV7Y3UxOZhq3XbKE//DRFcwvUu+kKBsZcX76bgP3P7OX3Y1dVBbl8FsfXMpvrl+ipD7DKEFIpLQeH+CHW+t4bFMte5u6mVeQxR2XL+PfXraUsnx9ucwkIyPOz/c089AvDvKLfcfIyUzj4xcs4jcuqebC6hKtxTEDKEFIyg2PxEoLj256j5/uaGRgeIQLqku4fX01Gy5YRE6m2hhmul0NnXz7F4fY+PZRegeHOWt+Ab+xbjE3X7iYikLNlBNVShCSEiMjzubDbfx421Ge3N7Ase5+inMzufnCRdy2vprVlUWpDlFC0N0/xE+2HeUHm+vYcriN9DTjgyvKueHcKq5dW6lSYsQoQci06Rsc5rUDLbywq4mndzTS0NlHTmYaV62ez43nLeSq1fNVWphD9jV188Otdfxkez2HW3pITzMurynnug9UcuXq+SwqyU11iHOeEoSE6mh7Lz/b3czzu5p4ed8xegeHyclM48MrK7jxvCp+9ZwF5Gdr8cK5zN15t76TJ7fX8+T2Bg4eOw7A2QsKuHLVfD6yqoJ1S8s0Oj4FlCAkqRo7+3jtQAuvHWjh1f0tHGrpAWBxaS5XrZ7PVavnc1lNuUoKEpe7s7/5OD/b3cQLu5t442Arg8NOflY6lywv49Ll5VxaU8a5i4rVdXYaKEHIGRsecfY3d/PWe+28WdvO6wdaOBD8+ivMyeDS5eVcVlPGR86u4Kz5Beq1Iqetu3+IV/Yd4+d7mnn9YCv7mroByMtK5+KlpVy6vIyLlpTygcXFFOVomo9kU4KQhLg7jZ39vF3Xzlu17bz1Xjvbj3TQ3T8ExBLCJcvKuLymnMtXlHNOVZFG0UrSHevu542Drbx2oIXXD7Syu7ELADOomZfP+dUlXFBdwvmLS1hdVahR9lOkBCHv0zc4zN7GbnbWd7KzoZNd9V3sauikLVgoJiPNOKeqiAtGP4zVJdTMyydNCUGmWXvPAG/XdfB2bTvbgh8vx7oHAMhMN1ZUFLC6spBVlUWsripkdWUhlUU5Ks0mSAliDmvvGWB/83EOHjvOgeZuDjQfZ29TFwePHWd0ZcnczHRWVRZyTlUhqxYUcu7iEtYuLFIbgkSSu3O0oy9IGB3sbuhkV0MX9R19Y+cU52ayqrKQsxcUUDOvgOUV+dTMy2dRSa4mGZzgVAlCXUtmOHenubufurZe6tp6qW3t4dCx4xw4FksKrccHxs7NSDOWluexoqKAf3PeQs6pLOScqiKWlOWpZCAzhpmxqCSXRSW53HBu1dj+jp5BdjV0sruxi10NXeyq7+SJt47S1Tc0dk5murGkLI/l8wqoqchnWXk+1WW5LC7NY2FJjqqrJlCCiLiBoRGau/tp6OjjaHtvkAh6fulv/9Avr/hVUZhNzbx8rl1bSc28fGoq8qmpKKC6VL+eZPYqzsvk0ppyLq05say9u9N6fCBWgg5+NB0MStQv7m1mYMJnZ0FRNotL81hcmhvc8lhUkktVcQ7zi3IoysmYU1VXShApMjzitPUM0NjZF9xiSaCp65e3R+taxyvNy2RxaR5nLyjkqtXzqS7L+6U3s8YciMSYGeUF2ZQXZLNuWdkvHRseceo7ejnS1jtWAh/90bX1vTZ+vK2e4ZFfroLPzUxnQVE2C4pyqCzOYUHR6C2byqIc5hfmMK8wi7ys2fEZnB3/iggYGXHaewdpPd7Pse4BWo8P0NLdT8vxAVpG7x/vH9tu6xlgJE7zz7yCrLE33/nVJWNvxgVF2SwsiSWBAiUAkSlLT7OgtJDHpXGODw2P0NjVT11rDw2dfTR19tPQ2Rds97H1vTYaO/vfVwoByMlMozw/m3kFWbEElT/+74nteQXZlOZnRrZqS980E/QNDtPZO0h77yAdvYN09Izb7h2ko2dgbHv8OSf7wgcoycukLD+LefnZrKgoYP3yrLE3zIKibOYX5VBZlENFYbYGBolEREZ62lhbx8m4O+09gzR29cVK/Z2jPwr7aT0+wLHjsVqCd4920nK8n8Hh+F8SeVnplORmUpSbSUleJiW5WZTkZVIcbBeP7Q/25WVRkptJXlZ6qFVecz5BtHT3c/s3Xxv70u8bfP+vgVFmUJidQUle7D+sODeThSW5FOdmUp6fRVlwm1eQTVnwS6E0L0tf+iKzlJlRmp9FaX7WpJNPujtd/UO0dMcSyLHuWK1Ca/fA2A/O9p5BOnsHOXCsm/ae2P2B4ZN/J2WmG4U5mVy0pJRv3RG3I9KUzPkEkZeVwfJ5+bEsnZc59sU/eisJ9pXkZlGQk6GBYSJyRsyMopxMinIyWT4vP6Fr3J2+wZEggQyMJY2O0e3eQbr6BllQGM4CW3M+QeRmpfOP/z75mVdEZKrMjNysdHKz0qksnv5VFlX3ISIicYWaIMzsOjPbbWb7zOyeOMezzeyx4PjrZrZs3LE/DfbvNrNrw4xTRETeL7QEYWbpwAPA9cAa4HYzWzPhtDuBNnc/C7gfuC+4dg1wG7AWuA74RvB4IiIyTcIsQawH9rn7AXcfAB4FNkw4ZwPwSLD9OHC1xfpsbQAedfd+dz8I7AseT0REpkmYCWIRUDvufl2wL+457j4EdADlCV6Lmd1lZpvNbHNzc3MSQxcRkRndSO3uD7r7OndfV1FRkepwRERmlTATxBGgetz9xcG+uOeYWQZQDLQkeK2IiIQozASxCVhpZsvNLItYo/PGCedsBO4Itm8BnvfYAhUbgduCXk7LgZXAGyHGKiIiE4Q2UM7dh8zsbuBpIB142N13mNm9wGZ33wg8BHzXzPYBrcSSCMF5PwDeBYaAP3D34VM935YtW46Z2eEphDwPODaF68OiuE6P4jo9UY0LohvbbItr6ckOzJoV5abKzDafbFWlVFJcp0dxnZ6oxgXRjW0uxTWjG6lFRCQ8ShAiIhKXEsQJD6Y6gJNQXKdHcZ2eqMYF0Y1tzsSlNggREYlLJQgREYlLCUJEROKa0wnCzL5qZrvMbJuZ/cjMSsYdS9l042Z2q5ntMLMRM1s3bv8yM+s1s7eC2/+azrhOFVtwLBJTtJvZl83syLjX6YZUxRLEc8pp71PFzA6Z2fbgNdqcwjgeNrMmM3tn3L4yM3vGzPYGf0sjElfK31tmVm1mL5jZu8Fn8Y+C/cl/zdx9zt6Aa4CMYPs+4L5gew3wNpANLAf2A+nTGNc5wCrgZ8C6cfuXAe+k+DU7WWwpfc0mxPhl4POpfn8FsaQHr0UNkBW8RmtSHVcQ2yFgXgTiuAK4aPx7G/hr4J5g+57Rz2YE4kr5ewuoAi4KtguBPcHnL+mv2ZwuQbj7Tz02iyzAa8TmfIIUTzfu7jvdffd0Pd/pOEVsmqI9vkSmvZ/T3P1FYjMpjDd+KYBHgI9PZ0xw0rhSzt3r3X1rsN0F7CQ223XSX7M5nSAm+G3gqWA7oenGU2S5mb1pZj83sw+nOphxovaa3R1UHT6ciuqJcaL2uoznwE/NbIuZ3ZXqYCZY4O71wXYDsCCVwUwQlfcWwSqcFwKvE8JrFtpcTFFhZs8ClXEO/Zm7PxGc82fE5nz6pyjFFUc9sMTdW8zsYuCfzWytu3dGILZpdaoYgf8J/AWxL8C/AP6W2A8A+WUfcvcjZjYfeMbMdgW/miPF3d3MotIfPzLvLTMrAP4f8Dl374yttRaTrNds1icId//VUx03s08BNwJXe1B5xzRMNz5ZXCe5ph/oD7a3mNl+4GwgqQ2MZxIb0zxFe6Ixmtk3gR+HFUcCIjt1vbsfCf42mdmPiFWHRSVBNJpZlbvXm1kV0JTqgADcvXF0O5XvLTPLJJYc/sndfxjsTvprNqermMzsOuALwE3u3jPuUCSnGzezCgvW5jazGmJxHUhtVGMi85oFH45RNwPvnOzcaZDItPfTzszyzaxwdJtYh41Uvk4TjV8K4A4gKiXXlL+3LFZUeAjY6e5fG3co+a9ZKlvjU30j1pBaC7wV3P7XuGN/Rqz3yW7g+mmO62ZiddX9QCPwdLD/E8COINatwMdS8JrFjS3Vr9mEGL8LbAe2BR+aqhS/z24g1tNkP7FqupTFMi6mGmI9qt4O3lMpiwv4P8SqTweD99adxJYefg7YCzwLlEUkrpS/t4APEavi2jbuu+uGMF4zTbUhIiJxzekqJhEROTklCBERiUsJQkRE4lKCEBGRuJQgREQkLiUIERGJSwlCRETiUoIQCYGZ/YmZfTbYvt/Mng+2rzKzaZvzS2QqlCBEwvESMDrb7jqgIJg/58NEZ74jkVNSghAJxxbgYjMrIjYtyavEEsWHiSUPkcib9bO5iqSCuw+a2UHgU8ArxObNuRI4i9gCLyKRpxKESHheAj5PrErpJeD3gDddE6DJDKEEIRKel4itH/yqx9YR6EPVSzKDaDZXERGJSyUIERGJSwlCRETiUoIQEZG4lCBERCQuJQgREYlLCUJEROJSghARkbj+P0dp9qkzFRx7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Z,_ = integrate.quad(f_prior, -100, 100)\n",
    "yx = [f_prior(w)/Z for w in wx]\n",
    "plt.plot(wx, yx)\n",
    "plt.xlabel('w')\n",
    "plt.ylabel('Probability')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "w_true = -2.5529898158340787\n",
      "w_map =  -3.139874975387916\n"
     ]
    }
   ],
   "source": [
    "model = LogisticRegressionMAP(fit_intercept=False)\n",
    "model.fit(X, y)\n",
    "w_map = model.coef_[0]\n",
    "print('w_true =', w_true[0])\n",
    "print('w_map = ', w_map[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_posterior(w):\n",
    "    return np.exp(compute_log_posterior(X, y, np.array([w])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoI0lEQVR4nO3deXxcZ33v8c9P0mixNtuS91UmdhIHJ3HiLEDshkAWeEEWGiC0oaGlzS00N0AvpaHQkBsKJaUFGqBA2uSSmwJJWxJwqWlisnMJiWUnsWMn3jfJtizLsvZdv/vHzDgTZWyNrJk5ozPf9+ull86cZfTzeDRfned5znPM3RERERmpIOgCREQkNykgREQkKQWEiIgkpYAQEZGkFBAiIpJUUdAFpEttba0vXLgw6DJkItq6Nfr99NODrUMkAOvXrz/i7tOSbQtNQCxcuJD6+vqgy5CJ6NJLo9+feirIKkQCYWZ7T7RNTUwiIpKUAkJERJJSQIiISFKh6YMQEQEYGBigoaGB3t7eoEvJKaWlpcydO5dIJJLyMQoIEQmVhoYGKisrWbhwIWYWdDk5wd1paWmhoaGBurq6lI9TE5OIhEpvby81NTUKhwRmRk1NzZjPqhQQIhI6Coc3O5XXRE1Mkvdau/vp6R8i0tHHtMqSoMsRyRk6g5C8du+vd7P1UAf7jnZz7Xf/H23dA0GXJJIzFBCSt/a1dHPXf7/GlEnFnDmrigNtPXznye1BlyWSMxQQkre+9/RODKirLae6LMK1587hx8/vo6d/KOjSRHKC+iAkL7X1DPCzFxu59tw5FP86+nfSB8+fyyMvNvL4a0287+zZAVcoafHpT8NLL6X3Oc89F771rZPu8vWvf52SkhJuvfVWPvOZz/Dyyy/zxBNP8MQTT3Dvvffyox/96E3HVFRU8IlPfII1a9Ywa9YsvvrVr/K5z32Offv28a1vfYurr76aPXv28NGPfpSuri4AvvOd7/D2t7+dp556ittvv53Kykp27NjBO9/5Tv7pn/6JgoLxnQPoDELy0qOvHKJnYIiPXDT/+LqLFtUwrbKE/37lUICVSRisXLmSZ599FoD6+no6OzsZGBjg2WefZdWqVUmP6erq4rLLLmPz5s1UVlbyxS9+kbVr1/LII49w++23AzB9+nTWrl3Lhg0beOihh7j11luPH//CCy/w7W9/my1btrBz504efvjhcf87dAYheekXmw4yd0oZ58ytPr6usMB4x1tq+PWOFtxdQyXDYJS/9DPl/PPPZ/369bS3t1NSUsJ5551HfX09zz77LHfffXfSY4qLi7nqqqsAWLZsGSUlJUQiEZYtW8aePXuA6FXit9xyCy+99BKFhYVs27bt+PEXXnghixYtAuAjH/kIv/71r7n++uvH9e/QGYTknY7eAX6z4wjvXTbrTSFw8aIajnT2sbO5K6DqJAwikQh1dXX88Ic/5O1vfzsrV67kySefZMeOHZx55pknPCb+fiwoKKCkpOT48uDgIADf/OY3mTFjBi+//DL19fX09/cfP37kezkdf+AoICTv/GZnC4PDzjtPn/6mbRctqgHghd1Hs12WhMzKlSv5+7//e1atWsXKlSv5/ve/z/Lly8f1wd3W1sasWbMoKCjggQceYGjo9QEVL7zwArt372Z4eJiHHnqISy65ZNz/BgWE5J1ntjVTXlzI+QumvGnbwppJVJUW8cqBtgAqkzBZuXIlBw8e5G1vexszZsygtLSUlStXjus5P/nJT3L//fdzzjnn8Nprr1FeXn582wUXXMAtt9zCmWeeSV1dHdddd914/wnqg5D889zOFi5eVENx0Zv/PjIzzppdzeYD7QFUJmHyrne9i4GB1y+8TOwvSKazs/P48h133JF02+LFi9m4cePx9Xfdddfx5aqqKn7xi1+Mp+Q30RmE5JUjnX3sOtLFBXVTT7jPW+dU8erBdgaGhrNYmUju0RmE5JX1e1sBuGDhm5uX4pbOrqJ/cJjdR7pYMqMyW6VJnrjooovo6+t7w7oHHniAZcuWnfJzXnrppVwav7d6GikgJK/U7zlKcVEBb51TfcJ9Fk+PhsLOw50KiAkql4cpP//884H8XHcf8zEZbWIys6vMbKuZ7TCz25Js/3Mz22JmG83scTNbkLDtJjPbHvu6KZN1Sv5Yt6eVc+dOpqSo8IT71NVGO/52NneecB/JXaWlpbS0tJzSB2JYxW8YVFpaOqbjMnYGYWaFwHeBy4EGYJ2ZrXb3LQm7vQiscPduM/sE8HfAh81sKvAlYAXgwPrYsa2ZqlfCr6d/iFca27h51aKT7ldeUsTs6lJdCzFBzZ07l4aGBpqbm4MuJafEbzk6FplsYroQ2OHuuwDM7EHgGuB4QLj7kwn7/xa4MbZ8JbDW3Y/Gjl0LXAX8JIP1Ssi93HCMwWFnxUn6H+LeMr1CZxATVPwiNRm/TDYxzQH2JzxuiK07kY8DvxzLsWZ2s5nVm1m9/lqQ0WxsOAbAOXMnj7rvotpydjV3qZlC8lpODHM1sxuJNid9fSzHufs97r7C3VdMmzYtM8VJaGxqbGd2dSk1FaPfNW5+TTmdfYMc0w2EJI9lMiAagXkJj+fG1r2Bmb0b+AJwtbv3jeVYkbF4pbHtpKOXEs2bUgbAvqPdmSxJJKdlMiDWAYvNrM7MioEbgNWJO5jZcuAHRMPhcMKmR4ErzGyKmU0BroitEzkl7b0D7D7SxdlzUwuI+TWTANjfqoCQ/JWxTmp3HzSzW4h+sBcC97n7ZjO7E6h399VEm5QqgH+PjVne5+5Xu/tRM/sy0ZABuDPeYS1yKjY3RqfOSP0MIhoQOoOQfJbRC+XcfQ2wZsS62xOW332SY+8D7stcdZJPXmmMTr63LMWAKC8poqa8mP0KCMljOdFJLZJpmxrbUu6gjps7pYzGY70ZrEoktykgJC9sPtDG0tmpnT3EzagqpalNASH5SwEhodc3OMSelm5On1kxpuNmVpdyqF0BIflLASGht/tIF0PDPuaJ92ZUldLWM0DvwNDoO4uEkAJCQm97U+xmK9PHFhAzq6ITmx1SM5PkKQWEhN72pg4KDBZNKx995wQzq6MBcVABIXlKASGht62pk4U15ZRGTjzFdzIzYmcQTeqHkDylgJDQ23a4g8UzxtZBDa+fQaijWvKVAkJCrW9wiL0t3WPufwCoKCmioqRIfRCStxQQEmq7mqMjmE7lDAJgRlWJmpgkbykgJNS2H46OYDrVe0vrWgjJZwoICbXtTR0UFtiYRzDF6WpqyWcKCAm1bU0dLKiZREnR2EYwxc2sKuVwRx/Dw7qznOQfBYSE2vamTpacQgd13MzqUgaHnSNdfaPvLBIyCggJrd6BIfa0dJ1yBzW8fi3E4XYFhOQfBYSE1q7mLoYdFp9iBzVAbUUxAEc6FRCSfxQQElrbD3cAsGQcZxA15dH7R7R09qelJpGJRAEhobW9qZPCAqOu9tRGMAHUxM4gjnYpICT/KCAktLY1dbBwHCOYIHo1dXFhgTqpJS8pICS0th/uPOUL5OLMjJqKYo6qiUnykAJCQql3YIi9LV3j6qCOq6kopkVNTJKHFBASSjubO6MjmKafegd13NTyElo0iknykAJCQil+F7nxNjEB1JbrDELykwJCQmlbUwdF4xzBFDe1vFjDXCUvKSAklLYf7mRhbTnFReN/i9dUlNAzMER3/2AaKhOZOBQQEkrbmzrGdYFcovi1EDqLkHyjgJDQ6R0YYu/RU7uLXDI15bGAUD+E5BkFhITOjsOduKengxqiTUyARjJJ3lFASOjE52AazyyuiXQGIflKASGhs62pk6ICY2HN+EcwgfogJH8pICR0tjd1UJemEUwAk4qLKIsUqolJ8o4CQkInHXMwjTS1vFgzukreUUBIqPT0D7HvaHfa+h/iaiuKOaKAkDyjgJBQ2dmc3hFMcZMnFXOsWwEh+UUBIaGyrSk2gikNk/QlmjwpwrHugbQ+p0iuU0BIqGxr6iRSaCxMwxxMiaZMKqZVZxCSZzIaEGZ2lZltNbMdZnZbku2rzGyDmQ2a2fUjtg2Z2Uuxr9WZrFPCIz6CKVKY3rd2dVmEjt5BBoeG0/q8IrmsKFNPbGaFwHeBy4EGYJ2ZrXb3LQm77QM+Bnw2yVP0uPu5mapPwmnb4Q7Onjs57c87ZVIEgLaegeNXVouEXSbPIC4Edrj7LnfvBx4Erkncwd33uPtGQH+Wybh19w/S0NrDkjTNwZRoSuxq6mM96oeQ/JHJgJgD7E943BBbl6pSM6s3s9+a2bXJdjCzm2P71Dc3N4+jVAmDnYe7YiOY0ttBDdEmJkAjmSSv5HIn9QJ3XwH8HvAtM3vLyB3c/R53X+HuK6ZNm5b9CiWnHB/BlOYhrhDtpAY0kknySiYDohGYl/B4bmxdSty9MfZ9F/AUsDydxUn4bDvcQaTQWFAzKe3PPTnWB9GqgJA8ksmAWAcsNrM6MysGbgBSGo1kZlPMrCS2XAu8A9hy8qMk321v6mRRbUXaRzBB9EI5UBOT5JeMBYS7DwK3AI8CrwL/5u6bzexOM7sawMwuMLMG4IPAD8xsc+zwM4F6M3sZeBL42ojRTyJvsq2pI+1TbMRVlhRRYGpikvySsWGuAO6+BlgzYt3tCcvriDY9jTzuN8CyTNYm4RIfwfShFfNG3/kUFBRYdLqNHp1BSP7I5U5qkZTtONwJZGYEU9zksoj6ICSvKCAkFLY1RQMiEyOY4iZPitCmgJA8ooCQUNje1EFxYQELpqZ/BFPcZM3HJHlGASGhsK2pg0XTyinKwAimOM3oKvlGASGhsK2pM6PNSwCTy3RPCMkvCgiZ8Lr6Bmk81sOSNN8DYqQpkyJ09Q/RP6ipwyQ/KCBkwouPYMr4GcTxCft0FiH5QQEhE158DqZMDnGF6DBX0MVykj9SCggze7+ZKUwkJ20/3ElxUQHzMziCCV6fsK+1S2cQkh9S/dD/MLDdzP7OzM7IZEEiY7WtqYNFtZkdwQSvT/ndpntCSJ5I6TfK3W8kOpvqTuCHZvZc7F4MmW30FUnB9qZOlmS4/wEUEJJ/Uv6Ty93bgf8geme4WcB1wAYz+58Zqk1kVB29A9ERTBnufwConqSAkPySah/ENWb2CNH7MkSAC939PcA5wP/KXHkiJxefYuP0mVUZ/1mVJUWYQbsCQvJEqrO5fgD4prs/k7jS3bvN7OPpL0skNfERTGfMzHwTU0GBUVUa0X2pJW+k2sR0aGQ4mNldAO7+eNqrEknR1kMdTCouZM7ksqz8vOqyiJqYJG+kGhCXJ1n3nnQWInIqth7qYMmMSgoKLCs/TwEh+eSkAWFmnzCzTcAZZrYx4Ws3sDE7JYok5+5sberg9CyMYIqbPEkBIfljtD6IHwO/BP4WuC1hfYe7H81YVSIpONLZz9Gufk7PQv9DXFVZhMbWnqz9PJEgjRYQ7u57zOzPRm4ws6kKCQnS1kPRDupsBoSamCSfpHIG8T5gPeBAYkOvA4syVJfIqLY2BRcQ7o5Zdvo9RIJy0oBw9/fFvtdlpxyR1G091E5NeTG1FSVZ+5mTyyIMDjvd/UOUl6Q6SlxkYjrpO9zMzjvZdnffkN5yRFK3takzq2cP8MbpNhQQEnajvcP/4STbHLgsjbWIpGx42Nne1MGHVszL6s+tTpjye3aWrr0QCcpoTUzvzFYhImPR0NpDd/9QVq6gTqQJ+ySfjNbEdJm7P2FmH0i23d0fzkxZIif32qF2AJZkOyA0YZ/kkdGamH4HeAJ4f5JtDiggJBCv30UumDMITdgn+WC0JqYvxb7/YXbKEUnNa4c6mDuljIosdxQf74PQfaklD6Q63XeNmd1tZhvMbL2Z/aOZ1WS6OJET2dbUkfX+B4CKkiIKC0xNTJIXUp2s70GgGfhd4PrY8kOZKkrkZPoHh9nV3JX15iUAM9PV1JI3Uj0/n+XuX054/Ddm9uFMFCQymu2HOxgcdpbOzvxNgpKJBsRgID9bJJtSPYN4zMxuMLOC2NeHgEczWZjIiWw5EB3BtHRWMAFRVRbhWLf6ICT8Rhvm2sHrczB9GvjX2KYCoBP4bCaLE0lmy8F2JhUXsqCmPJCfX10WoU0BIXlgtFFM2W/kFRnFlgPtnDGzksIs3SRopOqyCPtaugL52SLZlPIYQTObAiwGSuPrRt6GVCTT3J0tB9u55tzZgdUwWZ3UkidSCggz+2PgU8Bc4CXgYuA5NBeTZFlDaw8dvYOcGVD/A7w+5ffwsGftVqciQUi1k/pTwAXA3tj8TMuBY5kqSuREthwMtoMaogEx7NDZr5FMEm6pBkSvu/cCmFmJu78GnJ65skSS23KgnQKDM2YGGxAAbd1qZpJwSzUgGsxsMvAzYK2Z/RzYO9pBZnaVmW01sx1mdluS7atiV2cPmtn1I7bdZGbbY183pVinhNyWg+3U1ZZTVlwYWA2asE/yRUp9EO5+XWzxDjN7EqgG/vtkx5hZIfBd4HKgAVhnZqvdfUvCbvuAjzFiuKyZTQW+BKwgOsx2fezY1lTqlfDacqCd8xZMCbQGTfkt+SLVMwjM7DwzuxU4G2hw99EGgl8I7HD3XbF9HwSuSdzB3fe4+0ZgeMSxVwJr3f1oLBTWAlelWquEU1v3AI3HegLtfwAFhOSPVCfrux24H6gBaoH/Y2ZfHOWwOcD+hMcNsXWpSOlYM7vZzOrNrL65uTnFp5aJavPBNgDOnBXs5TkKCMkXqV4H8fvAOQkd1V8jOtz1bzJUV0rc/R7gHoAVK1Z4kLVI5r3SGA2IZXOqA61jsvogJE+k2sR0gIQL5IASoHGUYxqBxBsGz03hmHQcKyG1saGNOZPLqKkoCbSOskghkULjmEYxSciNNhfTt4l2ErcBm81sbezx5cALozz3OmCxmdUR/XC/Afi9FOt6FPhq7OptgCuAz6d4rITUpsY2zp4b7NkDaMpvyR+jNTHVx76vBx5JWP/UaE/s7oNmdgvRD/tC4D5332xmdwL17r7azC6IPe8U4P1m9r/d/Sx3P2pmXyYaMgB3uvvR1P9ZEjZt3QPsbenmwxfMG33nLKgqi+i2oxJ6o03Wd3982cyKgSWxh1vdfdTfDndfA6wZse72hOV1RJuPkh17H3DfaD9D8sOmWP/D2XMmB1tIjOZjknyQ6lxMlxIdxbSH6NTf88zsJk3WJ9mysfEYEHwHdVx1WYQjnZryW8It1VFM/wBc4e5bAcxsCfAT4PxMFSaSaFNDGwtqJh2/ijlo1WURdjR3Bl2GSEalOoopEg8HAHffBuTGb6rkhY0NbZw9d3LQZRwXvWmQmpgk3FINiPVm9i9mdmns6595vQNbJKNaOvtoPNbD2TnSvARQPamYjr5Bhod1+Y2EV6oB8afAFuDW2NcW4BOZKkokUbyDelkODHGNqy6L4A4dvZryW8Jr1D6I2KR7L7v7GcA3Ml+SyBttbGjDDM6aHewcTIni020c6+nPmX4RkXQb9QzC3YeArWY2Pwv1iLzJS/uP8ZZpFVSW5s4HseZjknyQ6iimKUSvpH4BOH63dne/OiNVicS4Oxv2tXLl0plBl/IGCgjJB6kGxF9ntAqRE9h1pItj3QOct2By0KW8gSbsk3ww2lxMpUQ7qE8DNgH3urt65SRrNuyN3iPqvPnB3iRopON9EBrqKiE2Wh/E/UTv6rYJeA/RC+ZEsmbDvmNUlRbxlmkVQZfyBmpiknwwWhPTUndfBmBm9zL6DK4iabVhbyvnzp9CQYEFXcoblEYKKS4q0IR9EmqjnUEcf/eraUmyrb13gG2HOzg/x5qX4jRhn4TdaGcQ55hZe2zZgLLYYwPc3XNnYLqEzsv7j+FOznVQx1WXRdQHIaE22nTfhdkqRGSkDXuPYQbnzpscdClJ6aZBEnapTrUhknXr97WyZHplTl0gl0gBIWGngJCcNDTsvLi3lfMW5Gb/A0D1JAWEhJsCQnLSlgPtdPQNcvGiqUGXckI6g5CwU0BITvrtrhYALl5UE3AlJ1ZdFqGzb5DBoeGgSxHJCAWE5KTnd7dQV1vOjKrSoEs5ofjFcu2a8ltCSgEhOWdo2Hl+99Gcbl4Czcck4aeAkJzz6sF2OnoHuagud5uXQNNtSPgpICTnxPsfLsrxM4jXJ+zrD7gSkcxQQEjO+e2uoyyomcSs6rKgSzkpnUFI2CkgJKcMDTsv7G7h4hxvXgKoindSKyAkpBQQklO2HGinvXcw55uXQGcQEn4KCMkpz2xvBuCSxbUBVzK6kqJCyiKFmrBPQksBITnl6W3NLJ1VxfTK3L3+IZGuppYwU0BIzujoHWDD3lZWLZkWdCkpU0BImCkgJGc8t7OFwWHndyZSQGjCPgkxBYTkjKe3NVNeXMj5OTyD60g6g5AwU0BITnB3ntnezNveUktx0cR5WyogJMwmzm+ihNqelm72H+3hd5bk/uilRAoICTMFhOSEp7YeBphQHdQAk8sidPcP0T+oKb8lfBQQkhMe29zEadMrWFBTHnQpY1KtGV0lxBQQErjWrn5e2HOUK8+aEXQpY6arqSXMMhoQZnaVmW01sx1mdluS7SVm9lBs+/NmtjC2fqGZ9ZjZS7Gv72eyTgnW468dZmjYufKsmUGXMmZVCggJsaJMPbGZFQLfBS4HGoB1Zrba3bck7PZxoNXdTzOzG4C7gA/Htu1093MzVZ/kjsc2H2JWdSnL5lQHXcqYTdaEfRJimTyDuBDY4e673L0feBC4ZsQ+1wD3x5b/A3iXmVkGa5Ic09M/xDPbm7li6Qwm4n/98XtC9OieEBI+mQyIOcD+hMcNsXVJ93H3QaANiM/zXGdmL5rZ02a2MtkPMLObzazezOqbm5vTW71kxdPbmukdGOaKCdi8BAl9EJqwT0IoVzupDwLz3X058OfAj82sauRO7n6Pu69w9xXTpk2s4ZES9djmQ1SXRbiwLven907m9T6IwYArEUm/TAZEIzAv4fHc2Lqk+5hZEVANtLh7n7u3ALj7emAnsCSDtUoAegeGeGxLE1csnUGkMFf/Vjm5SGEB5cWFamKSUMrkb+U6YLGZ1ZlZMXADsHrEPquBm2LL1wNPuLub2bRYJzdmtghYDOzKYK0SgF+92kRn3yDXLR/Z8jixTCkvprVLASHhk7FRTO4+aGa3AI8ChcB97r7ZzO4E6t19NXAv8ICZ7QCOEg0RgFXAnWY2AAwDf+ruRzNVqwTjZy82MrOqlIsW5f7tRU+mtqKEFgWEhFDGAgLA3dcAa0asuz1huRf4YJLjfgr8NJO1SbCOdvXz1NZm/uiSOgoLJt7opUS1FSU0tHYHXYZI2k3Mhl+Z8P5r00EGh51rz53YzUsA0yqLOdKpMwgJHwWEBOJnLzayZEYFZ86qDLqUcaspL+FoVx9Dwx50KSJppYCQrNtxuJP1e1u5dvmcCXlx3Ei1FcUMOxzr1lmEhIsCQrLuR8/vJVJofPD8eaPvPAHUVpYAqJlJQkcBIVnV0z/ET9c3cNVbZzEt9sE60dVWxAOiL+BKRNJLASFZ9Z8vH6C9d5AbL5ofdClpU1tRDCggJHwUEJJV//r8XpbMqJiwU2skEz+DaO5QQEi4KCAka17ef4yNDW38/kULQtE5HVddFiFSaLpYTkJHASFZ84NndlJZWsQHzpv41z4kMjNqyks4ojMICRkFhGTFzuZOfvnKIf7gbQuoLI0EXU7a1VYWqw9CQkcBIVnxg6d3UlxYwB++oy7oUjKiprxEw1wldBQQknEH23p45MVGbrhg3vEO3bCprSihRWcQEjIKCMm4Hzy9C3f4k1WLgi4lY2pj8zG5a7oNCQ8FhGTUniNd/Otv9/LBFfOYO2VS0OVkzLSKEvqHhmnv1Z3lJDwUEJJRX390K8VFBXzm8sVBl5JR8avCD7f3BlyJSPooICRjXtzXyn9tOsifrFzE9MrSoMvJqFnVZQAcbFNASHgoICQjhoedr/zXq9RWlIS67yFuVnU0AA+29QRciUj6KCAkIx6q30/93lb+4solVJRk9MaFOWFGVSlmcOCYziAkPBQQknaH23v56ppXuXjRVD60IhxTeo+muKiA2ooSnUFIqCggJO3u+M/N9A0O87cfODtUcy6NZnZ1qfogJFQUEJJWP3+pkTWbDvGpdy2mrrY86HKyalZ1GQeO6QxCwkMBIWmzq7mTv3p4EysWTOF/5EHH9Ehzp5TR0NrDsO5NLSGhgJC06B0Y4pYfv0ikqIC7P7KcosL8e2stqJlE3+AwhzWrq4RE/v0WS9q5O3/1yCa2HGznGx86h9mTy4IuKRALaqJNantbugKuRCQ9FBAybt/81XYe3tDIZ969hMvOmBF0OYFZeDwgugOuRCQ9FBAyLg+t28fdj2/nQyvmcuu7Tgu6nEDNnlxKUYGx96jOICQcFBByyh5at4/bHt7EqiXT+Mp1y/JqSGsyRYUFzJs6id1HFBASDgoIOSUPPLeHv/zpJlYtnsY9Hz2fSB52SiezeHoFrx3qCLoMkbTQb7WMydCw89U1r/LXP9/Mu8+czj1/cD6lkcKgy8oZZ8ysZM+RLnoHhoIuRWTcFBCSsrbuAf74/nXc88wu/uBtC/jejedTUqRwSHTGrCqGHbY3dQZdisi4hX8WNUmLZ7c38xf/vpEjnX185bq38vsXLQi6pJx0+sxKAF492M6yudUBVyMyPjqDkJM62tXP5x/exEfvfYGK0iIe+eQ7FA4nUVdTTnVZhPV7W4MuRWTcdAYhSfUODPGj5/fxj7/aRlf/EH98SR2fvfJ09TeMoqDAWLFgCuv2HA26FJFxU0DIG7T1DPDv9fv5wTO7aO7o45LTavnS+5eyeEZl0KVNGBfUTeXx1w5zuKM39HfSk3BTQAjuzoZ9rfzkhf38YuMBegeGeduiGu6+YTkXL5qa99c3jNXKxbV87Zfw+KuH+ciF84MuR+SUKSDyVN/gEOv3tPLYliYe23yIA229lBcXct3yufzehfPVwToOS2dVUVdbzi82HlBAyISmgMgD7s6h9l62HGhnY0Mbz+9u4cV9x+gbHKakqIBVS6bxmcuX8J5ls/Li9qCZZmZcfc5s7n5iO1sPdRwf2SQy0WT008DMrgL+ESgE/sXdvzZiewnwf4HzgRbgw+6+J7bt88DHgSHgVnd/NJO1TnRDw05rdz8Hj/Wy72j38a89R7p47VA7rd0DABQYLJ1dxY0XL+DiRTW847QaJhUrFNLtY29fyL2/3s1X1rzKDz92AQUFaqaTiSdjnwxmVgh8F7gcaADWmdlqd9+SsNvHgVZ3P83MbgDuAj5sZkuBG4CzgNnAr8xsibtP6MtT3Z2hYWdgyOkfGmYg/jU44vGQMzA0TM/AEJ29g3T1DdLZN0hX3xCdfQN09g3R2TfI0a4+Wjr7OdLZx9Gufkbep6amvJh5Uydx5VkzWTq7iqWzqjhjVpXOErJgSnkxf/meM/jrn73Cpx96ic9cvoSFNZPUnyMTSiY/KS4Edrj7LgAzexC4BkgMiGuAO2LL/wF8x6K/QdcAD7p7H7DbzHbEnu+5dBfZ2tXPB773G9wdB4bdcSf25Qw7ONF10Q/g2Dp//bvz5v2HHUhYjm8br7JIIeUlRVSUFDI1FgDL509hWkUxtZUlTK8sZUHNJOZNnaQgCNiNF82nrbufb6zdxuqXDzCpuJDqsghFhUakoEBnFQEI6yt+xqwqvv2R5Wl/3kx+gswB9ic8bgAuOtE+7j5oZm1ATWz9b0ccO2fkDzCzm4GbAebPP7XOwKJC461zqjGizS9mhhkYFnscWy4AEtYVmGEk378guuJN6wyIFBYQKSogUlhAcaFRVBhdjhQaxfHloujj0kghFSVFVJQUUV5SRHlxYV7eqW2iMjNuuWwx1y6fw1Nbm9lxuJOuvkEGh6NniK47k2ZV9E+5cJo3JTM36ZrQf2K6+z3APQArVqw4pf/9ytJIRpJXJG7ulEnceLGuPpeJJ5N/jjYC8xIez42tS7qPmRUB1UQ7q1M5VkREMiiTAbEOWGxmdWZWTLTTefWIfVYDN8WWrweecHePrb/BzErMrA5YDLyQwVpFRGSEjDUxxfoUbgEeJTrM9T5332xmdwL17r4auBd4INYJfZRoiBDb79+IdmgPAn820UcwiYhMNBntg3D3NcCaEetuT1juBT54gmO/Anwlk/WJiMiJaUiMiIgkpYAQEZGkFBAiIpKUAkJERJIyD8nlnGbWDOwdx1PUAkfSVE46qa6xUV1jk6t1Qe7WFra6Frj7tGQbQhMQ42Vm9e6+Iug6RlJdY6O6xiZX64LcrS2f6lITk4iIJKWAEBGRpBQQr7sn6AJOQHWNjeoam1ytC3K3trypS30QIiKSlM4gREQkKQWEiIgkldcBYWZfN7PXzGyjmT1iZpMTtn3ezHaY2VYzuzLLdX3QzDab2bCZrUhYv9DMeszspdjX97NZ18lqi20L7DUbUccdZtaY8Dq9N6haYvVcFXtNdpjZbUHWksjM9pjZpthrVB9gHfeZ2WEzeyVh3VQzW2tm22Pfp+RIXYG/t8xsnpk9aWZbYr+Ln4qtT/9r5u55+wVcARTFlu8C7ootLwVeBkqAOmAnUJjFus4ETgeeAlYkrF8IvBLwa3ai2gJ9zUbUeAfw2aDfX7FaCmOvxSKgOPYaLQ26rlhte4DaHKhjFXBe4nsb+DvgttjybfHfzRyoK/D3FjALOC+2XAlsi/3+pf01y+szCHd/zN0HYw9/S/TOdQDXAA+6e5+77wZ2ABdmsa5X3X1rtn7eWJyktkBfsxx2IbDD3Xe5ez/wINHXSmLc/Rmi94NJdA1wf2z5fuDabNYEJ6wrcO5+0N03xJY7gFeBOWTgNcvrgBjhj4BfxpbnAPsTtjXE1uWCOjN70cyeNrOVQReTINdes1tiTYf3BdE8kSDXXpdEDjxmZuvN7OagixlhhrsfjC0fAmYEWcwIufLewswWAsuB58nAa5bRGwblAjP7FTAzyaYvuPvPY/t8geid636US3UlcRCY7+4tZnY+8DMzO8vd23Ogtqw6WY3A94AvE/0A/DLwD0T/AJA3usTdG81sOrDWzF6L/dWcU9zdzSxXxuPnzHvLzCqAnwKfdvd2Mzu+LV2vWegDwt3ffbLtZvYx4H3AuzzWeAc0AvMSdpsbW5e1uk5wTB/QF1teb2Y7gSVAWjsYT6U2svCaJUq1RjP7Z+AXmaojBVl9XcbC3Rtj3w+b2SNEm8NyJSCazGyWux80s1nA4aALAnD3pvhykO8tM4sQDYcfufvDsdVpf83yuonJzK4CPgdc7e7dCZtWAzeYWYmZ1QGLgReCqDGRmU0zs8LY8iKide0KtqrjcuY1i/1yxF0HvHKifbNgHbDYzOrMrJjofddXB1gPAGZWbmaV8WWiAzaCfJ1GWg3cFFu+CciVM9fA31sWPVW4F3jV3b+RsCn9r1mQvfFBfxHtSN0PvBT7+n7Cti8QHX2yFXhPluu6jmhbdR/QBDwaW/+7wOZYrRuA9wfwmiWtLejXbESNDwCbgI2xX5pZAb/P3kt0pMlOos10gdWSUNMioiOqXo69pwKrC/gJ0ebTgdh76+NADfA4sB34FTA1R+oK/L0FXEK0iWtjwmfXezPxmmmqDRERSSqvm5hEROTEFBAiIpKUAkJERJJSQIiISFIKCBERSUoBISIiSSkgREQkKQWESAaY2V+Y2a2x5W+a2ROx5cvMLGtzfomMhwJCJDOeBeKz7a4AKmLz56wkd+Y7EjkpBYRIZqwHzjezKqLTkjxHNChWEg0PkZwX+tlcRYLg7gNmthv4GPAbovPmvBM4jegNXkRyns4gRDLnWeCzRJuUngX+FHjRNQGaTBAKCJHMeZbo/YOf8+h9BHpR85JMIJrNVUREktIZhIiIJKWAEBGRpBQQIiKSlAJCRESSUkCIiEhSCggREUlKASEiIkn9f1SpK52ApvWVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Z,_ = integrate.quad(f_posterior, -100, 100)\n",
    "yx = [f_posterior(w)/Z for w in wx]\n",
    "plt.plot(wx, yx)\n",
    "plt.xlabel('w')\n",
    "plt.axvline(x=w_map[0], color='red', label='w_map')\n",
    "plt.ylabel('Probability')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}