{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "baa0180c",
   "metadata": {},
   "source": [
    "### Exercise 28"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "365ef1b2",
   "metadata": {},
   "source": [
    "#### Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1852671d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import minimize\n",
    "from scipy.integrate import quad\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27f3e03a",
   "metadata": {},
   "source": [
    "#### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ed1103c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigma(x):\n",
    "    return 1/(1 + np.exp(-x))\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    return -np.log(1 + x**4) + np.log(sigma(2*x))\n",
    "\n",
    "\n",
    "def f1(x):\n",
    "    return -(4*x**3)/(1+x**4) + 2 - 2*sigma(2*x)\n",
    "\n",
    "\n",
    "def f2(x):\n",
    "    return (4*x**6 - 12*x**2)/(1 + x**4)**2 - 4*sigma(2*x)*(1 - sigma(2*x))\n",
    "\n",
    "\n",
    "def f22(x):\n",
    "    return (4*x**6 - 12*x**2)/(1 + x**4)**2 - 4*np.exp(2*x) / (np.exp(2*x) + 1)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5bf9e34",
   "metadata": {},
   "source": [
    "#### Find Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4f5971e5",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\denni\\AppData\\Local\\Temp\\ipykernel_12700\\2578787126.py:2: RuntimeWarning: overflow encountered in exp\n",
      "  return 1/(1 + np.exp(-x))\n"
     ]
    }
   ],
   "source": [
    "Z, error = quad(lambda x: sigma(2*x)/(1 + x**4), -np.inf, np.inf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "350f58c5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z=1.1107 with error=2.2388e-09\n"
     ]
    }
   ],
   "source": [
    "print(f'Z={np.round(Z, 4)} with error={np.round(error, 13)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "99cdefa2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def p(x):\n",
    "    return sigma(2*x)/(1 + x**4) / Z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "468d587c",
   "metadata": {},
   "source": [
    "### Maximise density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "4e9253a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successful optimisation: True\n"
     ]
    }
   ],
   "source": [
    "res = minimize(lambda x: -f(x), x0=0)\n",
    "print(f'Successful optimisation: {res.success}')\n",
    "xhat = res.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "ca059728",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = -f2(xhat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bff1b8e7",
   "metadata": {},
   "source": [
    "####  Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "b337f000",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6HUlEQVR4nO3dd3xUZdr/8c+VkBACoUcQiAgSkBB6IPQmIEWkKkVQLIsNXfdxXevPVVd83FXx0dVVQMGyILpBFCkLqAgogoSIkICEIlUwIZAAGyDt+v2RsgmpQGbOJHO9X695kTnnnjNflMw159z3uW9RVYwxxngvH6cDGGOMcZYVAmOM8XJWCIwxxstZITDGGC9nhcAYY7xcFacDXKz69evr1Vdf7XQMY4ypULZs2XJcVYOL2lfhCsHVV19NdHS00zGMMaZCEZEDxe2zS0PGGOPlrBAYY4yXs0JgjDFezgqBMcZ4OSsExhjj5awQGGOMl7NCYIwxXq7C3UdgTGV36tQp3n33XTIyMoiIiKB///5ORzKVnBUCYzzIhg0bmDRpEgcO/PfenxUrVjBkyBAHU5nKzi4NGeMhMjIymDJlCr6+vnz33XekpKTw97//nUGDBjkdzVRyUtFWKIuIiFCbYsJUVr/88guZmZm0aNGiwPYjR46QnJxMmzZtHEpmKjoR2aKqEUXts0tDxniApKQk6tatS7NmzQrtU1Wuv/56fHx82Lp1Kz4+diJvypf9izLGYarKgAEDuO2224rcLyI88cQTbN++nWXLlrk5nfEGVgiMcdjy5cvZtm0bAwcOLLbNzTffzFVXXcVLL73kxmTGW1ghMMZhL7/8MldddRUTJ04stk2VKlX4wx/+wPr169m0aZMb0xlvYIXAGAcdOHCAb775ht/97nf4+fmV2PbOO++kdu3abNiwwU3pjLewzmJjHPTRRx8BcMstt5TaNigoiCNHjhAYGOjqWMbLWCEwxkEPPPAAHTp0KHK0UFFyi0BmZia+vr6ujGa8iEsvDYnIEBHZJSJ7ROSxIva/KiJbcx7xIpLsyjzGeJrq1atf9F3D48ePZ8qUKS5KZLyRywqBiPgCbwJDgTBgooiE5W+jqn9Q1Q6q2gH4O/Cpq/IY42nmzZvHX//6Vy72ps569erx2WefcfbsWRclM97GlWcEXYE9qrpPVdOAhcDIEtpPBD5yYR5jPMobb7zB0qVLEZGLet3IkSM5e/Ys33zzjWuCGa/jykLQGDiU7/nhnG2FiEhToBnwdTH7p4lItIhEJyYmlntQY9ztyJEjxMTEcMMNN1z0a/v27Uu1atVYvny5C5IZb+Qpw0cnAFGqmlnUTlWdraoRqhoRHBzs5mjGlL/cD/FLKQQBAQFcd911LF++/KIvKxlTFFeOGjoChOR73iRnW1EmAPe7MIsxHmXp0qU0bdqUsLCw0hsX4e6772b37t1kZmZSpYoN/jOXx5X/gjYDoSLSjOwCMAGYdGEjEbkWqAN878IsxngUX19fRo8efdH9A7ku5UzCmOK4rBCoaoaITAdWAr7AXFWNE5HngGhVXZLTdAKwUO0c13iRTz+9/AFyCQkJxMfH06tXr3JIZLyZrUdgjJuV181g06ZN45NPPiEpKcluLjOlKmk9Ak/pLDbGa4wYMYJJkwpdJb1offv2JSUlhZ9++qkcUhlvZoXAGDdKT09n7dq11KtX77KP1bdvXwDWrl172ccy3s0KgTFutHnzZlJTU+nfv/9lH6tJkya0aNHCbiwzl80KgTFu9O233wLQu3fvcjlev379WLduHZmZRd6CY0yZWCEwxo02bNhAaGgo5XVj5OOPP84PP/xg6xiby2J3ohjjRjfeeCODBw8ut+M1b9683I5lvJcVAmPc6I477ij3Yy5YsIBz58655NjGO9j5pDFucuDAAY4cKW6WlUu3YMECXn755XI/rvEeVgiMcZMXX3yRsLCwcu/Y7d69Ozt37uTkyZPlelzjPawQGOMmGzZsIDIystzvAu7evTsAmzZtKtfjGu9hhcAYNzh79ixxcXF07dq13I/dpUsXfHx8+P57m7fRXBorBMa4wbZt28jMzKRz587lfuygoCDatGnjkv4H4x1s1JAxbrBlyxYAlxQCgOjoaPz9/V1ybFP5WSEwxg1GjhxJvXr1CAkJKb3xJbAiYC6HXRoyxg0aN27M+PHjy7QQzZ49e3jkkUf4wx/+wPvvv09WVlaprzl06BAjRozg66+LXPbbmBJZITDGxc6fP8+sWbM4ePBgqW1nzJhB69atef3113nnnXeYOnUqAwcO5LfffivxdbVr12bp0qV899135RXbeBErBMa42LZt27jnnnvYvHlzie3+9re/8dRTTzF27Fj279/PqVOneOedd9i4cSM333wzGRkZxb42KCiIli1bEhMTU97xjRdwaSEQkSEisktE9ojIY8W0uVlEdohInIgscGUeY5yQ++FcUkfxsmXLePTRR5kwYQLz58/nyiuvRES48847mT17NuvWreOZZ54p8X06depkhcBcEpcVAhHxBd4EhgJhwEQRCbugTSjwONBTVdsAD7kqjzFO2bJlC3Xr1qVp06ZF7k9JSeHuu+8mPDyc9957r9ANZ5MnT2bq1Km8+OKL7N27t9j36dy5MwcPHuT48ePlmt9Ufq48I+gK7FHVfaqaBiwERl7Q5nfAm6p6EkBVE1yYxxhHbNmyhU6dOhXbUfz4449z9OhR5s6dS9WqVYts88ILL+Dn58fzzz9f7PtERkbSu3dvkpKSyiW38R6uLASNgUP5nh/O2ZZfS6CliHwnIhtFZEhRBxKRaSISLSLRiYmJLoprTPlLT08nNja22MtCu3fvZvbs2dx333106dKl2ONceeWV3HPPPXz44Yfs2bOnyDa9e/dm3bp1tGrVqlyyG+/hdGdxFSAU6AdMBOaISO0LG6nqbFWNUNWI8lrQwxh38PPz47fffuPhhx8ucv+zzz5L1apVefLJJ0s91qOPPoqvry+vv/56ie3KMtzUmPxcWQiOAPnvnmmSsy2/w8ASVU1X1V+AeLILgzGVRu3atYtckWzHjh0sWLCABx54gIYNG5Z6nIYNGzJ69Gjmz5/P+fPni2zz6KOP0rZt28vObLyLKwvBZiBURJqJiD8wAVhyQZvPyD4bQETqk32paJ8LMxnjVm+99RavvPJKkftmzpxJQEAAf/zjH8t8vDvuuIMTJ06wZMmFv0rZateuzY4dO0hOTr6UuMZLuawQqGoGMB1YCewEPlHVOBF5TkRuzGm2EkgSkR3AGuARVbWeLlNpzJ07l2XLlhXanpCQwD//+U9uu+026tevX+bjXXfddTRp0oR58+YVub9Tp04AbN269ZLyGu/k0j4CVV2uqi1V9RpVnZGz7WlVXZLzs6rq/6hqmKq2VdWFrsxjjDtlZGSwfft2OnbsWGjf22+/zfnz53nooYcu6pi+vr7ceuutrFy5kqIGTuQWArufwFwMpzuLjam09u7dy/nz52nXrl2B7RkZGbz99tsMHTr0kkb4jB07lqysLJYuXVpoX3BwME2aNLFCYC6KFQJjXCQ2NhaA8PDwAttXrlzJ0aNH+d3vfndJx+3YsSMhISF89tlnRe6fPn06/fv3v6RjG+9k01Ab4yJJSUnUrFmTsLACN9Qzb948goODGT58+CUdV0QYNWoU77zzDqmpqQQGBhbY/+ijj15yZuOd7IzAGBeZNm0aycnJVKtWLW9bYmIiS5YsYcqUKZe1hsCoUaM4e/Ysq1atKnJ/QkICp0+fvuTjG+9ihcAYF7pwWon58+eTnp7O7bffflnH7d27N7Vq1Sqyn2DHjh00aNCAL7744rLew3gPKwTGuMDZs2fp3r17gev4qsrcuXPp0qVLoX6Di+Xn58eAAQNYvXo1qlpgX2hoKH5+fmzfvv2y3sN4DysExrjAzz//zMaNG0lLS8vbFhMTw/bt27njjjvK5T0GDRrEwYMHC8095Ofnx7XXXsu2bdvK5X1M5WeFwBgXKGrE0Lx58wgICGDChAnl8h4DBw4E4Msvvyy0r127dnZGYMrMCoExLhAbG4u/vz+hodlTZ6Wnp/PRRx8xatQoateuXS7v0aJFC5o2bcrq1asL7Wvbti2HDh3i5MmT5fJepnKz4aPGuMD27du59tpr8fPzA2Dt2rWcOHGC8ePHl9t7iAiDBg3iX//6FxkZGVSp8t9f5xEjRtCgQYO89zemJHZGYIwLNG3alMGDB+c9j4qKonr16lx//fXl+j4DBgwgJSWlUH9AWFgYU6dOpUaNGuX6fqZysjMCY1zgrbfeyvs5MzOTxYsXM3z48AL3FJSH3r17A7Bu3bq8eYZyxcbGcubMGbp161au72kqHzsjMKacXTic87vvviMhIYGxY8eW+3s1adKEZs2asX79+kL77rvvvmIXxDEmPysExpSzOXPmEBISQkJC9hLcUVFRBAQEMGzYMJe8X+/evVm/fn2hApQ7cujC7cZcyAqBMeUsNjaW5ORkgoODycrK4tNPP2XIkCEuu17fp08fEhMT2bVrV4Htbdu25fTp0xw4cMAl72sqDysExpSz2NhYwsPDERF++OEHjhw54pLLQrn69OkDUOjyUO7013Y/gSmNFQJjypGqsn379rwbyaKiovDz8+OGG25w2Xu2aNGCBg0aFCoEuRmsEJjS2KghY8pRQkICx48fJzw8HFVl0aJFDBo0qNxuIiuKiNCtWzc2bdpUYHtQUBBr1qy57HmNTOXn0jMCERkiIrtEZI+IPFbE/qkikigiW3Med7kyjzGulpmZyfTp0+nduzc//vgj+/fvd+lloVyRkZHEx8dz4sSJAtv79et3UWsiG+/kskIgIr7Am8BQIAyYKCJhRTT9WFU75DzecVUeY9yhUaNG/P3vf6dTp05ERUXh6+vLyJEjXf6+ufcK/PDDDwW279y5kxdffLHA5HfGXMiVZwRdgT2quk9V04CFgOt/I4xx0LFjx0hPT8+7LNS/f3/q1avn8veNiIhARApdHoqJieHxxx8nPj7e5RlMxeXKQtAYOJTv+eGcbRcaKyLbRCRKREKKOpCITBORaBGJTkxMdEVWY8rFqFGjGDZsGHFxccTHx7vlshBk9weEh4ezcePGAttz+wdyZ0M1pihOjxr6ArhaVdsBq4H3i2qkqrNVNUJVI4KDg90a0JiyysrKIi4ujrCwMKKiovLWFnaXyMhINm3aVOAGslatWuHr60tcXJzbcpiKx5WF4AiQ/xt+k5xteVQ1SVXP5zx9B+jswjzGuNSBAwc4c+YM4eHhLFq0iN69e9OwYUO3vX+3bt04efIku3fvztsWEBBAaGionRGYErmyEGwGQkWkmYj4AxOAJfkbiMiV+Z7eCOx0YR5jXCr3wzYoKIjY2Fi3XRbKFRkZCVDk5aGff/7ZrVlMxeKyQqCqGcB0YCXZH/CfqGqciDwnIjfmNHtQROJE5CfgQWCqq/IY42q5hWDHjh0AjBkzxq3v37p1a4KCggp1GL/99tu2bKUpkVS0CakiIiI0Ojra6RjGFBITE8OGDRuYN28e/v7+fP/9927PMHDgQE6ePMmWLVvc/t7Gs4nIFlWNKGqf053FxlQanTp1Yvjw4cTExLj9slCuyMhIfvrpJ1JTU/O2nT59mgceeICVK1c6ksl4PisExpSD9PR0Vq1axYcffgjgWCHo1q0bmZmZxMTE5G2rVq0ac+bMKXKRe2PACoEx5SI+Pp7rr7+eDz74gI4dO9KsWTNHchTVYVylShVat25tI4dMsawQGFMOcj9k9+7dy7hx4xzLccUVV3D11VcX6jAODw+3QmCKZYXAmHIQGxuLj0/2r5NTl4Vyde3alc2bNxfYFh4ezuHDh0lOTnYmlPFoVgiMKQexsbFUrVqVNm3a0KpVK0ezdOnShQMHDpB/Opbw8HCaNGnCr7/+6mAy46msEBhTDrZu3crZs2cdvSyUq0uXLgAFzgqGDRvGoUOHCAsragJg4+2sEBhTDiZPngw4f1kIsoexikiBQiAiDiYyns4KgTHlYOPGjYSGhnrEamBBQUG0bt26UD/BE088wdSpU50JZTyaFQJjLtPq1av56quvGDVqlMd88+7SpQubN28uMBNpYmIiy5cvdzCV8VRWCIy5TC+//DKq6hH9A7m6dOlCQkIChw79d0mQ8PBwEhMTSUhIcDCZ8URlKgQiMkJErGgYU4SYmBj8/f3zOmk9QVEdxrZIjSlOWT/cxwO7ReRvInKtKwMZU5GkpKRw/PhxWrRo4TGXhQDat2+Pn59fgTWM27RpA1ghMIWVqRCo6mSgI7AXeE9Evs9ZPjLIpemM8XALFy4EoH///g4nKahq1aq0a9euwBlBgwYNGDRoELVr13YumPFIZb7co6qngCiyF6G/EhgNxIjIAy7KZozHmz9/PgBDhw51OElhXbt2ZcuWLWRlZQHZQ0hXrVrFrbfe6nAy42nK2kcwUkQWA98AfkBXVR0KtAcedl08YzzXmTNn2Lx5M3feeScDBgxwOk4hXbp04dSpU8THxxfYnpWVRUVbh8S4VlnPCMYAr6pqW1V9SVUTAFQ1FbjTZemM8WArVqzg3LlzTJkyhWrVqjkdp5CiOowXLVpErVq1OHz4sFOxjAcqayE4pqrr8m8Qkb8CqOpXxb1IRIaIyC4R2SMij5XQbqyIqIgUuXqOMZ4oKiqKwMBA9u/f73SUIrVu3Zrq1asX6ic4c+YM27dvdzCZ8TRlLQSDithW4kVREfEF3sxpFwZMFJFCE53kdDj/Hth04T5jPNXZs2dZunQpmZmZhRaL9xS+vr506tSpQCGwkUOmKCUWAhG5V0S2A9eKyLZ8j1+A0lbD7grsUdV9qppGdifzyCLa/QX4K3DuEvIb44hVq1aRmprK+fPnPWJaieJ06dKFrVu3kp6eDkCdOnVo3LixFQJTQGlnBAuAEcDnOX/mPjrnDCktSWPgUL7nh3O25RGRTkCIqi4r6UA5Q1WjRSQ6/9S6xjglKiqKGjVqAHh8ITh37lyBD35bpMZcqLRCoKq6H7gfOJ3vgYjUvZw3zrlTeSZlGHWkqrNVNUJVI4KDgy/nbY25bGlpaXzxxRe0bt0a8PxCAAU7jCdNmuRR02EY51UpZf8C4AZgC6BA/lsnFWhewmuPACH5njfJ2ZYrCAgHvsm5I7MhsEREblTV6DKlN8YBX331FSkpKVxzzTUcP36cevXqOR2pWM2bN6du3bps3ryZadOmAdh9BKaQEguBqt6Q8+elrMS9GQgVkWZkF4AJwKR8x04B6uc+F5FvgD9aETCeLioqipo1a/Lee+/h7+/vdJwSiQgRERGFpqROTk4mKyuLunUv68TeVBJlvaGsp4hUz/l5sojMFJGrSnqNqmYA04GVwE7gE1WNE5HnROTGyw1ujBMyMjL4/PPPGTFiBFWrVvWo+YWK06VLF2JjY0lNTQUgNTWVunXr8sYbbziczHiKsg4ffQtIFZHcO4n3Ah+W9iJVXa6qLVX1GlWdkbPtaVVdUkTbfnY2YDzd2rVrSUpKomfPnvTp04fvv//e6Uil6tKlC5mZmWzduhWAwMBAmjdvTlxcnLPBjMcoayHI0Ox70kcCb6jqm2Rf4zfGq+TeRFa7dm3Wr1+Pn5+f05FKVdyU1DZyyOQqayE4LSKPA5OBZTkjfjz/N8CYcpSZmcnixYsZNmwYu3btwsfHp0IsBt+oUSMaN25cqBDEx8dz/vx5B5MZT3Ex6xGcB+5U1WNkjwB6yWWpjPFAGzZs4LfffmPcuHFs27aN0NBQAgMDnY5VJrlLV+YKDw8nIyOj0IR0xjuVdT2CY6o6U1XX5zw/qKofuDaaMZ4lKiqKqlWrMmzYMLZt20a7du2cjlRmXbp0IT4+nuTkZAB69erFrFmzaNCggbPBjEco66ihMSKyW0RSROSUiJwWkVOuDmeMp8jKyuLTTz/l+uuvJzAwkLCwMPr27et0rDLL7SeIjs4ej9GkSROmTZvGFVdc4WQs4yHKemnob8CNqlpLVWuqapCq1nRlMGM8yebNmzl8+DDjxo3D19eXJUuWcP/99zsdq8wiIrIn9s1/eWj37t2sX7/eqUjGg5R2Z3Gu31R1p0uTGOPBFi1ahJ+fHyNGjEBVK8T9A/nVqVOHFi1aFCgETz31FNHR0ezdu9fBZMYTlPWMIFpEPhaRiTmXicaIyBiXJjPGQ6gqUVFRXHfdddSuXZvf//73dOzY0elYFy0yMpKNGzfmrU4WHh7Ovn37+M9//uNwMuO0shaCmkAqMJj/zkB6g6tCGeNJtm7dyi+//JI3UdtPP/3kkSuSlaZHjx4cPXqUgwcPAv+dLG/Hjh1OxjIeoEyXhlT1dlcHMcZTLVq0CF9fX0aOHImqsn37dm6++WanY1207t27A9nDYJs2bZpXCLZt25bXmWy8U1lHDbUUka9EJDbneTsRecq10YxxXu5lob59+1K/fn2OHDnCyZMnadu2rdPRLlrbtm2pXr163rQY11xzDdWrV+enn35yOJlxWlkvDc0BHgfSAVR1G9mziRpTqe3YsYNdu3blXRbKXeu3It1DkKtKlSp07dqVDRs2AODj48OKFSt4/PHHHU5mnFbWQhCoqj9csC2jvMMY42kWLVqEiDB69Ggge/H3adOmVcgzAsjuJ9i6dWteB3Hv3r258sorHU5lnFbWQnBcRK4hezEaRGQccNRlqYzxEFFRUfTs2ZOGDRsC0KlTJ2bNmkXt2rWdDXaJunfvTmZmZt6NZYcPH2bmzJkcO3bM4WTGSWUtBPcDs8hexP4I8BBwj6tCGeMJdu/ezfbt2wss63jgwAEyMzMdTHV5unXrBpDXT/Drr7/y8MMPs3HjRidjGYeVWAhE5H9E5H+AUcByYAbwNvApMNbl6Yxx0KJFiwAYMyb7lpm0tDRatGjB008/7WSsy1KvXj1atWqV108QHh6Oj49P3loFxjuVdkYQlPOIAO4F6gC1yT4b6OTSZMY4LCoqiq5duxISkr309s6dO8nIyKiw/QO5evTowffff4+qEhgYSMuWLa0QeLkSC4GqPquqz5I97XQnVf2jqj4MdAZKXKoSQESGiMguEdkjIo8Vsf8eEdkuIltF5FsR8fzJ3Y1X2L9/P1u2bClwWejHH38EqJB3FefXo0cPjh8/zp49ewDo0KGDFQIvV9Y+ggZAWr7naTnbiiUivsCbwFAgDJhYxAf9AlVtq6odyJ7YbmYZ8xjjUp9++ikAY8f+9wpoTEwMNWrUIDQ01KlY5SL/jWUA7du35/Dhw5w6ZRMKe6uyFoIPgB9E5BkReQbYBLxXymu6AntUdZ+qpgELyV7qMo+q5v+XV52cUUnGOC0qKooOHTrQvHnzvG0//vgjHTp0wMenrL82nql169bUqlUrr8P4/vvv58yZM9SsaRMKe6uyTjExQ0RWAL1zNt2uqj+W8rLGwKF8zw8DkRc2EpH7gf8B/IEBZcljjCsdOXKE77//nueff77A9ieeeCJvwraKzMfHh+7du/Ptt98CEBRky497uzJ/tVHVGFV9LedRWhEoM1V9U1WvAR4Fipy2QkSmiUi0iEQnJiaW11sbU6TFixcDBS8LAQwdOpRhw4Y5Eanc9enTh7i4OHJ/n1544QVeeeUVh1MZp7jyHPcIEJLveZOcbcVZSPYw1UJUdbaqRqhqRHBwcPklNKYIUVFRhIWFce211+Zt27VrF2vXriUjo3LcUN+vXz8A1q1bl/fn/PnzHUxknOTKQrAZCBWRZiLiT/bcREvyNxCR/L1uw4HdLsxjTKkSEhJYv359gdFCAO+++y6DBw+uFJeGIHvFsurVq7NmzRogu8M4Li6OtLS0Ul5pKiOXFQJVzQCmAyuBncAnqhonIs+JyI05zaaLSJyIbCW7n+A2V+Uxpiw+++wzsrKyCl0WiomJoW3btvj5+TmUrHz5+fnRq1cvvvnmGyB7CGlaWho///yzs8GMI1w6/EFVl6tqS1W9RlVn5Gx7WlWX5Pz8e1Vto6odVLW/qsa5Mo8xpYmKiqJFixYFbhpTVX788Uc6dapc91D279+fuLg4EhIS6NChA/DfeyWMd6nY4+CMKUcnTpxgzZo1jB07tsCaxAcPHuTEiROVrhDk9hOsXbuWli1bcvXVV9uylV6qrIvXG1PpLVmyhIyMjEL9AzExMUDFv6P4Qp06daJGjRqsWbOGm266iV9++cXpSMYhdkZgTI6oqCiaNm1K586dC2wfMmQIGzdupH379g4lcw0/Pz969+6d109gvJcVAmOAU6dOsXr1asaMGVPgshBAtWrViIyMJCAgwKF0rtOvXz927tzJsWPH2LRpE2FhYTbvkBeyQmAMsHTpUtLS0gqNFgJ47rnnKu18/fn7CerXr8/OnTv54YcLFyM0lZ0VAmPIvizUqFGjvAnZch07dow///nPlbYQdOrUiaCgINasWUPz5s2pW7cumzdvdjqWcTMrBMbrnTlzhhUrVjB27NhCE8pt2bIFoNKNGMpVpUoVevfuzddff42IEBERkbeMpfEeVgiM11u+fDnnzp0rNFoIYOPGjfj6+hbqQK5MBg8ezO7du9m7dy8RERFs376ds2fPOh3LuJEVAuP1oqKiaNCgAT179iy0b+PGjbRr147q1as7kMw9hg8fDsCKFSsYOHAgkydP5vTp0w6nMu5khcB4tdTUVJYtW8aYMWPw9fUtsE9V+eWXX/IWfK+sWrRoQWhoKMuXL6d///689957XHHFFU7HMm5kN5QZr/bvf/+b1NTUIi8LiQi7d+/2isskw4YNY9asWaSmplKtWjVOnDhBvXr1nI5l3MTOCIxXi4qKon79+vTp06fI/SJCYGCgm1O537Bhwzh37hzffPMNkydPpkePHk5HMm5khcB4rXPnzvHFF18wevRoqlQpfHL8zDPP8NBDD7k/mAP69OlDYGAgy5cvp3Xr1sTHx5OSkuJ0LOMmVgiM11q1ahVnzpwp8rIQwKJFi4iPj3dzKmcEBAQwcOBAli1bRkREBPDfobOm8rNCYLxWVFQUderUoX///oX2nTp1iri4OCIjCy2zXWkNGzaM/fv3U6tWLQA2bdrkcCLjLlYIjFc6f/48S5YsYdSoUUUuNrN582ZUtdKPGMpv6NChAGzYsIFrr702b3F7U/nZqCHjlb766itSUlKKvSyUO6VE165d3RnLUVdddRXh4eEsX76cZ599lho1ajgdybiJnREYrxQVFUWtWrW47rrritwfFBTE8OHDqVOnjpuTOWv48OGsW7eOQYMGMWzYMKfjGDdxaSEQkSEisktE9ojIY0Xs/x8R2SEi20TkKxFp6so8xgCkp6fz2WefceONN1K1atUi2zz44IMsXbrUzcmcN3bsWDIyMvj888/ZtGmTLV3pJVxWCETEF3gTGAqEARNFJOyCZj8CEaraDogC/uaqPMbkWrNmDSdPniz2slBGRgaq6uZUniEiIoKrrrqKqKgoxo4dy9/+Zr+S3sCVZwRdgT2quk9V04CFwMj8DVR1jaqm5jzdCDRxYR5jgOzLQjVq1GDw4MFF7l+4cCHBwcFeuXSjiDB27FhWr15NZGQk69ev99qi6E1cWQgaA4fyPT+cs604dwIritohItNEJFpEohMTE8sxovE2GRkZLF68mBEjRhS74tjatWvJzMzkqquucnM6zzBu3DjS0tIICgriyJEjHDhwwOlIxsU8orNYRCYDEcBLRe1X1dmqGqGqEcHBwe4NZyqVdevWcfz48WIvC+W26dWrV6FJ6LxFt27dCAkJYd++fQCsX7/e4UTG1VxZCI4AIfmeN8nZVoCIDASeBG5U1fMuzGMMn3zyCYGBgQwZMqTI/ceOHSM+Pr7YuYe8gY+PDxMnTmTDhg3UrFnT7ifwAq4sBJuBUBFpJiL+wARgSf4GItIRmEV2EUhwYRZjSE9PJyoqipEjRxY7kVzut9++ffu6M5rHueWWW8jMzGT69Om8+uqrTscxLuayQqCqGcB0YCWwE/hEVeNE5DkRuTGn2UtADeBfIrJVRJYUczhjLttXX31FUlISEyZMKLZNq1ateOSRR+jYsaMbk3metm3b0qZNG7755huvmH3V20lFGxEQERGhtqaquRS33XYbn3/+Ob/99lux9w+Y/3rxxRd5/PHHuffee5k0aRK9evVyOpK5DCKyRVUjitrnEZ3FxrjauXPnWLx4MWPGjCm2CJw8eZLvvvuO9PR0N6fzTLfeeisiwuzZs1m0aJHTcYwLWSEwXmHFihWcPn2aiRMnltimV69ebNu2zY3JPFejRo0YPnw4vr6+fPnll07HMS5khcB4hdybxIqacjrXl19+SZ06dejQoYP7gnm4O++8k7S0NGJjYzl27JjTcYyLWCEwld6ZM2f44osvuOmmm4pciQyyF6r/8ssvGTBggNfeP1CU4cOH561d/NVXXzmcxriKFQJT6S1ZsoSzZ8+WOFpoz549HDp0iIEDB7oxmefz8/Pj/vvvB2Dnzp0OpzGuYoXAVHoffvghISEh9OzZs9g2q1evBrBCUIS7774bX19fUlNTS29sKiQrBKZS+/XXX1m1ahVTpkzBx6f4f+533XUX69evp0WLFm5MVzE0atSIcePGMXfuXFvQvpKyQmAqtfnz55OVlcWtt95aYjt/f38bJ1+Ce++9l5SUFG655RanoxgXsEJgKi1V5f3336dbt260atWq2Hbffvstf/rTnzhx4oQb01Usffv2pUaNGqxevdrus6iErBCYSismJoa4uDhuu+22Ett9/PHHvPHGGzaVQiluuOEG0tLSmD17ttNRTDmzQmAqrffee4+qVasyfvz4YtuoKkuXLmXAgAHFrk9gsj344IMAPP/882RkZDicxpQnKwSmUkpLS+Ojjz5i5MiRJS5A/9NPP7F//35GjRrlvnAVVNeuXalZsybHjh1jwYIFTscx5cgKgamUli1bRlJSUqmXhRYvXoyPjw833nhjie0M+Pr68sILL9CsWTOeffZZ6yuoRKwQmEppzpw5NGzYsNh1iXNlZGQwfPhwrrjiCjclq9juv/9+3nzzTfbt28esWbOcjmPKiRUCU+ns3buXf//730ybNq3YKSVyzZgxgyVLbBmMixEWFkaHDh145plnSE5OdjqOKQdWCEyl89Zbb+Hj48O0adNKbHfy5Ek3Japcnn76afbt20dSUhLPP/+803FMObBCYCqV1NRU5s6dy+jRo2ncuHGx7VSVXr16MXXqVPeFqyTGjRvHqVOnGDx4MK+99hqxsbFORzKXyaWFQESGiMguEdkjIo8Vsb+PiMSISIaIjHNlFuMdFi5cyMmTJ5k+fXqJ7bZt28aOHTuIjIx0U7LK4/rrr6devXpUr16dWrVqce+995KVleV0LHMZXFYIRMQXeBMYCoQBE0Uk7IJmB4GpgI1FM5dNVXnzzTdp06YNffr0KbHtggULqFKlCjfddJOb0lUe/v7+jB8/nhUrVvDss8/y7bffWsdxBefKM4KuwB5V3aeqacBCYGT+Bqq6X1W3AfZ1wly2TZs2ERMTw/3334+IFNsuMzOTBQsWMHjwYOrXr+/GhJXHlClTSEtLo3nz5gwaNIhHHnmEX375xelY5hK5shA0Bg7le344Z9tFE5FpIhItItGJiYnlEs5UPq+//jpBQUFMnjy5xHYrV67k8OHD3H777W5KVvlERkZy5MgRhg4dyrvvvouvry9TpkyxewsqqArRWayqs1U1QlUjgoODnY5jPNC+ffv4+OOPufvuuwkKCiqxbb9+/fjggw/sJrLLICI0bNgQgMaNGzN79my+++47HnusUFegqQBKHmR9eY4AIfmeN8nZZky5e/nll6lSpQp/+MMfSm0bGBjIlClT3JCqcsvKymLkyJGEhoYyc+ZMNmzYwMyZM+nevTvjxtnYj4rElWcEm4FQEWkmIv7ABMDu3DHl7ujRo8ydO5fbbruNRo0aldh2zpw5vPrqq6iqm9JVXj4+PtSoUYN58+aRmprKSy+9RPfu3bn99tvZtWuX0/HMRXBZIVDVDGA6sBLYCXyiqnEi8pyI3AggIl1E5DBwEzBLROJclcdUXv/7v/9LRkYGf/rTn0psd+7cOZ566ilWr15dYmeyKbv77ruP5ORkPvzwQ/z9/fnkk08ICAhgxIgRJCQkOB3PlJWqVqhH586d1ZhcBw4cUH9/f73rrrtKbfvuu+8qoF9++aUbknmHrKws7dKli4aGhmpGRoaqqn733XdarVo17dChgyYnJzuc0OQCorWYz9UK0VlsTHGeffZZAP7f//t/JbbLyspi5syZtGvXjgEDBrgjmlcQEf74xz+ye/fuvDmbevToweLFi4mLi2PEiBG26H0FYIXAVFhbt25l3rx53HfffVx11VUlts39YHrkkUfsslA5GzNmDP/3f/9Hv3798rZdf/31/POf/+Tbb79l9OjR/Oc//3EuoCmVaAXrNIuIiNDo6GinYxiHqSr9+/cnNjaW3bt3l7j4DMB3333HP/7xDz744AN8fX3dlNLMmzePu+66i86dO7N06VKb7ttBIrJFVSOK2mdnBKZCWrhwIWvXruUvf/lLqUUAoGfPnsyfP9+KgAutWrWK8ePHF5h36Pbbb+ezzz4jNjaWnj17snfvXgcTmuJYITAVTlJSEr///e/p2rVrqVNNnz17lqeeeooTJ064KZ33On78OJ988gkffvhhge0jRozgq6++4sSJE0RGRrJ8+XKHEpriWCEwFc5DDz3EyZMneeedd0r9hv/KK68wY8YMtm7d6p5wXmzChAl069aNRx99lJSUlAL7unfvzsaNG2nSpAnDhw/nT3/6k01H4UmKG07kqQ8bPurd/vnPfyqgzzzzTKlt9+3bp4GBgTpmzBg3JDOqqps3b1YR0YceeqjI/ampqXrPPfcooJGRkbp9+3Y3J/RelDB81PEP9ot9WCHwXvHx8RoUFKQ9e/bU9PT0EttmZmZqv379NCgoSA8cOOCmhEZV9d5771UR0a1btxbb5uOPP9Z69epplSpV9LHHHtP//Oc/bkzonawQmAovOTlZr732Wq1Xr57u37+/1PZvvfWWAjpnzhw3pDP5nTp1SufNm6dZWVkltktMTNSpU6cqoFdffbV+8MEHeTelmfJXUiGwPgLj8dLS0rj55pvZs2cPixYtomnTpqW+ZsSIETzxxBPceeedbkho8gsKCmLq1KmICEeOHMn+xlmE+vXrM2/ePL7++mvq1KnDrbfeSvv27Vm8eLGteOZmVgiMR8vIyGDSpEmsWrWKWbNm0bdv3xLbJyUlkZGRQePGjZkxY4bdPOag+Ph4wsLCmDlzZont+vfvT3R0NB9//DHp6emMGTOGa6+9lrfeestuRHOX4k4VPPVhl4a8x9mzZ3XUqFEK6Kuvvlpq+5SUFO3YsaPefPPNrg9nSpWVlaVjx45VEdGFCxeW6TXp6en60UcfaZcuXRTQ2rVr63333ac//PBDqZeaTMmwPgJT0fz222/au3dvBfS1114rtX1ycrL26tVLfX19dfny5W5IaMoiNTVV+/Tpo35+fvrZZ5+V+XVZWVn67bff6qRJkzQgIEABbd26tT755JNWFC6RFQJToaxbt04bN26sAQEB+tFHH5Xa/tdff9XOnTtrlSpV9OOPP3ZDQnMxkpOTNTIyUn19fXXt2rWX9PpZs2Zp37591cfHRwFt1KiR3nPPPbpkyRI9efJk+YeuhKwQmArh5MmT+uCDD6qI6DXXXKM//vhjqa/JysrSjh07amBgoH7xxReuD2kuyenTp/W5554rddhvaY4fP67vv/++jhkzRqtXr66Aioh26NBBH3zwQY2KitIDBw7YGUMRSioENumccVxqaipz5sxhxowZHD9+nPvuu4+//vWvVK9evcTXiAjVqlVjw4YN1KhRg3bt2rkxtblUv/32GxMmTOCll14iIqLIOdDK5Ny5c2zatIm1a9eybt06NmzYwNmzZwGoV68enTp1olOnTrRv355WrVrRsmVLatSoUV5/jQqnpEnnrBAYx8TGxvLuu+/y4YcfkpSURP/+/XnllVfo2LFjsa9JS0tj/vz5PP3009x11138+c9/dmNiUx5++OEHRo8ezdGjR5k8eTJ//vOfueaaay77uGlpacTExBR4xMbGFpjKonHjxrRq1YpWrVoRGhpK06ZNCQkJ4aqrriI4OBgfn8o7kNIKgfEIGRkZbNq0iZUrV7Js2TJiYmLw8/Nj5MiRPPDAA/Tp06fY1+7YsYMFCxbw/vvvc/jwYTp37syrr75K79693fg3MOUlJSWFF154gddeey3vPpGPPvqo3If7nj9/nvj4eHbt2sWuXbsK/JycnFygrb+/PyEhIYSEhNCgQQOuuOKKIh/16tWjZs2aFW4mW8cKgYgMAV4DfIF3VPXFC/ZXBT4AOgNJwHhV3V/SMa0QeD5VJSkpiT179hAfH8+PP/6Y9w3tzJkz+Pj40LVrV8aPH88tt9xCcHBwoWMkJiayadMmhg8fjohw8803s2jRIgYOHMhDDz3EkCFD7B6BSuDo0aO8/vrrJCYm8s477wDZq801bNiQHj160LZtW6pUqVLu76uqnDhxgkOHDnHw4MG8Pw8ePMjhw4dJSEggISGhULHIr2bNmtSqVYvatWvnPXKf16xZk+rVq1O9enUCAwNL/DP3Zz8/P5f+m3akEIiILxAPDAIOA5uBiaq6I1+b+4B2qnqPiEwARqvq+JKOa4Wg/KkqmZmZpKen5z0yMjLyfj5//jxnzpzhzJkznD59utDPSUlJHDt2jGPHjvHbb7/x66+/cvr06bzjV6tWjXbt2tGuXTsiIyOJiIjAz8+PpKQk2rdvT82aNdm4cSNz587N+9Z27NgxAPbu3Uvz5s3Zs2cPQUFBNGjQwKn/TMYNTp06RVhYGEeOHAEgICCA0NBQpk+fzrRp0zh37hz/+te/CA4OJjg4mNq1a1OtWjXq1KlDtWrVyj3P+fPnOX78eF5hSEhI4Pjx46SkpJCcnJz354U/nzp1iszMzIt6LxGhatWqBR4BAQEFnj/yyCOMGjXqkv4uJRWC8i+1/9UV2KOq+3JCLARGAjvytRkJPJPzcxTwhoiIuqA6zZ07l4cffrjQ+ql+fn40adIEVeXo0aOcO3euwH5/f38aNmwIZH97SUtLK7C/atWqed9ojx49mvc/P/evEBAQQN26dQE4duxYoVvnAwICqFWrFqpKQkICF/7VAwICCAoKIisri6SkpEJ/r4CAAKpVq0ZWVlaR315y/yFlZWUV+HDOzejj43PZt/PXrFmTmjVrcvToUXx8fBARAgIC8PHx4bXXXmPq1KksWbKEsWPHMmfOnAKvXb9+Pb169WLfvn0sXryYli1bMnToUFq3bk2PHj1o0qQJAC1atLisjKZiqFmzZt638w0bNhAdHU18fHzeh/yvv/7KrbfeWuh1r732Gg8++CCxsbF07doVf39/fH198fHxwcfHh7///e/cfPPNREdHM27cuLztvr6+iAivv/46gwcPZv369fzud78rdPx3332X66+/nn//+9/MmDGj0P6FCxfSoUMHPv30U5588kkaNGiAqpKVlYWq8uabbxIcHMynn37KBx98QFZWVt6+rKwsbrvtNqpWrcrGjRuJjo4mLS2N8+fPk5KSgqrSrVs3srKyXHY5ypWFoDFwKN/zw0BkcW1UNUNEUoB6wPH8jURkGjANKHVt2uLUr1+fJk2aFPiwzB110qFDByD7Q//UqVMFXhcUFET79u0REbZu3VrglncRoVatWrRt2xYRITo6ukAhERHq1q1LeHg4ABs3bszruMo9Bbziiito3bo1AN9++23eh3Lu/kaNGtGyZUtUlXXr1hXYBxASEkKLFi3IzMxk/fr1hfY3a9aMZs2acf78eb7//vsCfzcfHx9at25NaGgo586dY/369QV+OXx8fOjRowdt2rTh7NmzLFu2DH9/f/z8/PL+vOOOO+jatSvx8fH84x//KPB3h+yFzKtUqUK7du14+eWXqVatGtWqVSMwMLDAf5tJkyYxadKkUv4vGm8gIjRt2pSmTZsyceLEAvtCQkKIj48nMTGRxMREUlJSOHv2LD179gSgTp06TJ8+nfPnz+d92GZmZhISEgJkF5p+/frlbc/94la7dm0g+/c99/Mgv6CgoLx2Re3PHeFWp06dIkevtW7dmpCQEI4ePcquXbsK7f/jH/+YVygWLlxYaP/cuXNdOuLJlZeGxgFDVPWunOdTgEhVnZ6vTWxOm8M5z/fmtDle1DHBLg0ZY8ylcGrN4iNASL7nTXK2FdlGRKoAtcjuNDbGGOMmriwEm4FQEWkmIv7ABGDJBW2WALfl/DwO+NoV/QPGGGOK57I+gpxr/tOBlWQPH52rqnEi8hzZtzovAd4FPhSRPcAJsouFMcYYN3JlZzGquhxYfsG2p/P9fA64yZUZjDHGlKzy3k9tjDGmTKwQGGOMl7NCYIwxXs4KgTHGeLkKN/uoiCQCBy7x5fW54K5lD2G5Lo7luniems1yXZzLydVUVQvP8EgFLASXQ0Sii7uzzkmW6+JYrovnqdks18VxVS67NGSMMV7OCoExxng5bysEs50OUAzLdXEs18Xz1GyW6+K4JJdX9REYY4wpzNvOCIwxxlzACoExxng5ry0EIvKwiKiI1Hc6C4CI/EVEtonIVhFZJSKNnM4EICIvicjPOdkWi0htpzMBiMhNIhInIlki4vgwPxEZIiK7RGSPiDzmdB4AEZkrIgk5C0B5DBEJEZE1IrIj5//h753OBCAiASLyg4j8lJPrWacz5SciviLyo4gsLe9je2UhEJEQYDBw0Oks+bykqu1UtQOwFHi6lPbushoIV9V2QDzwuMN5csUCY4B1TgcREV/gTWAoEAZMFJEwZ1MB8B4wxOkQRcgAHlbVMKAbcL+H/Pc6DwxQ1fZAB2CIiHRzNlIBvwd2uuLAXlkIgFeBPwEe01OuqvkXS66Oh2RT1VWqmpHzdCPZK805TlV3qmrhxV+d0RXYo6r7VDUNWAiMdDgTqrqO7HU+PIqqHlXVmJyfT5P94dbY2VSg2c7kPPXLeXjE76GINAGGA++44vheVwhEZCRwRFV/cjrLhURkhogcAm7Bc84I8rsDWOF0CA/UGDiU7/lhPOCDrSIQkauBjsAmh6MAeZdftgIJwGpV9YhcwP+R/eU1yxUHd+nCNE4RkS+BhkXsehJ4guzLQm5XUi5V/VxVnwSeFJHHgenAnz0hV06bJ8k+pZ/vjkxlzWUqLhGpASwCHrrgjNgxqpoJdMjpC1ssIuGq6mgfi4jcACSo6hYR6eeK96iUhUBVBxa1XUTaAs2An0QEsi9zxIhIV1U95lSuIswne2U3txSC0nKJyFTgBuA6d64pfRH/vZx2BAjJ97xJzjZTDBHxI7sIzFfVT53OcyFVTRaRNWT3sTjd2d4TuFFEhgEBQE0R+aeqTi6vN/CqS0Oqul1Vr1DVq1X1arJP4Tu5owiURkRC8z0dCfzsVJb8RGQI2aekN6pqqtN5PNRmIFREmomIP9lrby9xOJPHkuxvYe8CO1V1ptN5colIcO6oOBGpBgzCA34PVfVxVW2S85k1Afi6PIsAeFkh8HAvikisiGwj+9KVRwypA94AgoDVOUNb33Y6EICIjBaRw0B3YJmIrHQqS05n+nRgJdkdn5+oapxTeXKJyEfA90ArETksInc6nSlHT2AKMCDn39TWnG+7TrsSWJPzO7iZ7D6Cch+q6YlsigljjPFydkZgjDFezgqBMcZ4OSsExhjj5awQGGOMl7NCYIwxXs4KgTHGeDkrBMYY4+WsEBhzmUSkS856DQEiUj1nLvtwp3MZU1Z2Q5kx5UBEnid7HphqwGFV/V+HIxlTZlYIjCkHOXMMbQbOAT1yZrE0pkKwS0PGlI96QA2y52UKcDiLMRfFzgiMKQcisoTslcmaAVeq6nSHIxlTZpVyPQJj3ElEbgXSVXVBzvrFG0RkgKp+7XQ2Y8rCzgiMMcbLWR+BMcZ4OSsExhjj5awQGGOMl7NCYIwxXs4KgTHGeDkrBMYY4+WsEBhjjJf7/1aTC5UePigIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mplot = 200\n",
    "space = np.linspace(-4, 4, mplot)\n",
    "\n",
    "plt.plot(space, p(space), color='black', label='p')\n",
    "plt.plot(space, norm(loc=xhat, scale=A**(-1/2)).pdf(space), color='black', linestyle='dashed', label='Laplace approximation')\n",
    "#plt.plot(space, norm(loc=xhat, scale=1/2.0).pdf(space), color='black', linestyle='dashed', label='Laplace approximation')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('density')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}