Trying the RMS alsorithm from the Uncertainty in Neural Networks: Approximately Bayesian Ensembling paper by Pearce et al (2020) on trivial example of estimating mean for normal distribution with known variance

Pearce_etal_2020_for_normal_mean.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trying the RMS alsorithm from the [*Uncertainty in Neural Networks: Approximately Bayesian Ensembling*](https://arxiv.org/abs/1810.05546v5) paper by Pearce et al (2020) on trivial example of estimating mean for normal distribution with known variance. In this case we know the exact solution for the problem, since it can bre calculated using [conjugacy](https://en.wikipedia.org/wiki/Conjugate_prior#When_likelihood_function_is_a_continuous_distribution)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as sp\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "\n",
    "n = 15\n",
    "\n",
    "# true parameters\n",
    "μ = 5\n",
    "σ = 2.7\n",
    "\n",
    "x = sp.norm(μ, σ).rvs(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior_μ = 0\n",
    "prior_σ = 10\n",
    "prior_dist = sp.norm(prior_μ, prior_σ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5.003623405947105, 0.6954491092861294)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def posterior(x, prior_μ, prior_σ):\n",
    "    n = len(x)\n",
    "    σ2 = σ ** 2\n",
    "    prior_σ2 = prior_σ ** 2\n",
    "    \n",
    "    # see: https://en.wikipedia.org/wiki/Conjugate_prior#When_likelihood_function_is_a_continuous_distribution\n",
    "    post_σ2 = 1 / (1/prior_σ2 + n/σ2) \n",
    "    post_μ = post_σ2 * (prior_μ/prior_σ2 + np.sum(x)/σ2)\n",
    "    \n",
    "    return float(post_μ), np.sqrt(post_σ2)\n",
    "\n",
    "post_μ, post_σ = posterior(x, prior_μ, prior_σ)\n",
    "post_dist = sp.norm(post_μ, post_σ)\n",
    "\n",
    "post_μ, post_σ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4.9764284 , 0.69544911],\n",
       "       [4.95463788, 0.69544911],\n",
       "       [5.01882196, 0.69544911],\n",
       "       ...,\n",
       "       [5.00404823, 0.69544911],\n",
       "       [4.99212787, 0.69544911],\n",
       "       [5.0073286 , 0.69544911]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = 5000\n",
    "results = []\n",
    "\n",
    "for _ in range(R):\n",
    "    acc_μ = prior_dist.rvs(1)\n",
    "    results.append(posterior(x, acc_μ, prior_σ))\n",
    "    \n",
    "results = np.vstack(results)\n",
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5.003886029404278, 0.048175803294772626)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(results[:, 0]), np.std(results[:, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEKCAYAAAAyx7/DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3xU9Z3/8dcnN0II4RLCHRtkWUQIRBoQQVAQK7Zqrdqq9bJeWtTHesO2yrY+fuvW1h9adavVlgeF+rPKCltbdNcVq1SpRa0YFCuKrNZCuaiEKOFOSPj8/jgnwySZZCYwQzjwfj4e85gzM99z5ntyMu/5zvd8zznm7oiISHRktXcFRESkbRTcIiIRo+AWEYkYBbeISMQouEVEIkbBLSISMQpuEZGIUXCLiESMgluOaGb2rpmdmsblDTGzt8xsm5ndmK7lirSF6chJSRczWwP0AuqBHcCzwA3uvv0glvctd1+crjoeLDObC2x19+ntXRc5eqnFLel2trsXAqOA0cDt7VEJM8vJ0PxfAN491PURiafgloxw9w3AImA4gJkNNbMlZrYl7L44p6Gsmd1mZhvC7ofVZnaamT0GHAP8t5ltN7Nbw7J9zey3ZlZlZn+L764wszXhsv4C7DCznPC5KXFlWqtHs/nj18nMXgQmAQ+FdfrHg1leWOYHZvaLuMfdzGyvmeUf8B9fjngKbskIMxsAfBl4y8xygf8Gngd6AjcA88L+4iHA9cBod+8MnAGscffLgL8TtuDd/R4zywqX8zbQDzgNuNnMzoh764uBrwBd3b2uSZ1arEcq87v7ZOBPwPXhr4q/HczyQmXAirjH5cBqd9+doKwIoOCW9HvKzLYAS4E/AncBY4FCYKa717r7i8AzBKFWD3QAjjezXHdf4+5/bWHZo4ESd/9huJyPgF8CF8WVedDd17n7rgTzt1aPVObPxPISBffbKby3HMXU7ybpdm7TnYlm1hdY5+774p5eC/Rz9w/N7GbgDmCYmf0euMXdNyZY9heAvuEXQ4NsglZwg3Wt1K3FeqQ4f1qXZ2Z5wCDgnbinR9I4yEWaUYtbDoWNwICwq6PBMcAGAHf/D3c/mSCYHbg7LNN0yNM64G/u3jXu1tndvxxXprVhUq3WI4X5072844EN7r4TwMwMOBW1uCUJBbccCq8TDA+81cxyw3HVZwPzw37uyWbWAdgN7CLoPgH4FDg2bjnLgK3hDr+OZpZtZsPNbPTB1iPd65Xi/GVATzMbZGYdgTsJvrzWHGB95Cih4JaMc/da4BzgTGAz8HPgcnd/n6B/e2b4/CcEO/m+H876f4HbwxEb33X3eoJgLCfYMbgZmAN0SUM90r1eqSgDfk8w+uZDgi+qj4AfHEh95OihA3BE2omZLQLmuPtv27suEi1qcYu0nzJgVXtXQqJHLW6RdmBm3Qi6Rjq5+972ro9Ei4JbRCRi1FUiIhIxGTkAp0ePHl5aWpqJRYuIHJGWL1++2d1LUimbkeAuLS2lsrIyE4sWETkimdnaVMuqq0REJGIU3CIiEaPgFhGJGJ0dUCTD9u7dy/r169m9W6fYFsjPz6d///7k5uYe8DIU3CIZtn79ejp37kxpaSnBCQDlaOXuVFdXs379egYOHHjAy1FXiUiG7d69m+LiYoW2YGYUFxcf9K8vBbfIIaDQlgbp+F9IKbjNbHp4IdSVZvaELmQqItJ+kga3mfUDbgQq3H04waWiLmp9LhE5XFRXV1NeXk55eTm9e/emX79+sce1tbXtXb1m7r///gPqSvjBD37ASy+9lIEaHX5S3TmZA3Q0s71AAcElm0Qio3TG/zR6vGbmV9qpJodecXExK1YEl7G84447KCws5Lvf/W6jMu6Ou5OV1f69p/fffz9XXXUV+fmp/7Cvq6vjxz/+cZvep66ujpycaI7PSLqV3H0DcC/wd+BjoMbdn29azsymmVmlmVVWVVWlv6YiklYffvghw4cP59prr2XUqFGsW7eOrl27xl6fP38+3/rWtwD49NNPOe+886ioqGDMmDH8+c9/bra8OXPm8LWvfY0zzjiDIUOG8KMf/Sj22j333MPw4cMZPnw4P/vZzwDYtm0bZ555JiNHjmT48OE8+eST/Pu//zubNm1iwoQJTJkyBYBFixZx0kknMWrUKC688EJ27NgBQP/+/bnzzjsZP348Cxcu5NJLL+Wpp54C4IUXXqC8vJyysjK+/e1vx35ZNJ0nqpJ+3YTnDf4qMBDYAvzGzC5198fjy7n7bGA2QEVFhc4VK5LIohnwyTvJy7VF7zI4c+YBzfree+/xyCOPMGvWLOrq6losd+ONN3LrrbcyduxY1qxZw1lnncXKlSublVu2bBkrV64kLy+P0aNHc9ZZZ1FbW8u8efNYtmwZ9fX1jBkzhlNOOYVVq1ZRWlrKokWLAKipqaFLly7cd999/OlPf6Jr165s2rSJmTNn8oc//IGCggJ+/OMf88ADD/D97wdXt+vUqROvvPIKAE8//TQAO3fu5KqrrmLJkiUMGjSISy65hNmzZ3P99dc3myeqUvmdMIXgytpVAGb2O2Ac8Hirc4nIYW/QoEGMHp38WsuLFy9m9erVsceff/45u3btomPHjo3KnXHGGXTr1g2Ac889l6VLl7Jnzx7OP/98CgoKGj0/adIkZsyYwYwZMzj77LMZP358s/d99dVXee+99xg3bhwAtbW1nHzyybHXL7zwwmbzrFq1isGDBzNo0CAALr/8cubOnRsL7kTzRE0qwf13YKyZFRBcgfs0QKf+EzkQB9gyzpROnTrFprOysoi/sEr8DkJ3Z9myZeTl5bW6vKZD3cyMli7WMnToUCorK3n22Wf53ve+x1lnnRVrSce/79SpU3nssceS1j9+ntYkmidqUunjfh14EngTeCecZ3aG6yUih1hWVhbdunXjgw8+YN++fY36gKdMmcLDDz8ce9yws7Op559/ni1btrBz506efvppxo8fz8SJE1m4cCG7du1i+/btPP3000yYMIENGzZQWFjIZZddxi233MKbb74JQOfOndm2bRsA48aN449//CMfffQRADt27OCDDz5odT2OP/54Pvjgg9g8jz/+OKeccsqB/2EOQyntUnX3fwX+NcN1EZF2dvfddzN16lSOOeYYjj/+ePbs2QPAww8/zHXXXccjjzxCXV0dkyZNahTkDU4++WS++c1v8te//pXLLruM8vJyAC6++OJYl8x1111HWVkZzz77LDNmzCArK4u8vDxmzZoFwLRp05gyZQoDBgxg8eLFzJ07lwsvvDC2g/Guu+5i8ODBLa5DQUEBc+fO5bzzzqO+vp4TTzyRb3/722n9O7W3jFxzsqKiwnUhBTmctOdwwFWrVjF06NBD9n7tZc6cOaxcuZKf/vSn7V2Vw16i/wkzW+7uFanM3/6DNkVEpE2iOfpcRA47DWO+JfPU4hYRiRgFt4hIxCi4RUQiRsEtIhIx2jkpcog1HZp4sFIZ2pidnU1ZWRl1dXUMHDiQxx57jK5du7JmzRoGDhzI7bffzp133gnA5s2b6dOnD9dccw0PPfQQq1ev5pprrmHLli3s2bOHCRMmMHv24XsM3qxZsygoKODyyy8/6GXdddddjY7mHDduHK+++upBL/dgqcUtchTo2LEjK1asYOXKlXTv3r3RwTPHHnsszzzzTOzxb37zG4YNGxZ7fOONNzJ9+nRWrFjBqlWruOGGGw5ZvVs78VVLrr322rSENgTBHe9wCG1QcIscdU466SQ2bNgQe9yxY8fYeUMAFixYwDe+8Y3Y6x9//DH9+/ePPS4rK2u2zO3bt3PaaacxatQoysrKYmfqW7NmDccddxz/9E//xIgRI7jgggvYuXMnAKWlpdx2222MGTOGMWPG8OGHHwJwxRVXcMsttzBp0iRuu+02PvvsM84991xGjBjB2LFj+ctf/gIEXyg//OEPAfj973/PxIkT2bdvH3fccQf33nsvAKeeeirTp09n4sSJDB06lDfeeIPzzjuPwYMHc/vtt8fqf+655/LFL36RYcOGxX5NzJgxg127dlFeXs4ll1wCQGFhIRCcD+V73/sew4cPp6ysjAULFgCwZMkSTj31VC644AKOO+44LrnkkqTnTjkQ6ioROYrU19fzhz/8gauvvrrR8xdddBHz58+nd+/eZGdn07dvXzZuDK6XMn36dCZPnsy4ceP40pe+xJVXXtnovN0A+fn5LFy4kKKiIjZv3szYsWM555xzAFi9ejVz585l/PjxXHXVVfz85z+PXcihqKiIZcuW8etf/5qbb7451vL/3//9XxYvXkx2djY33HADJ5xwAk899RQvvvgil19+OStWrGDmzJmMHj2aCRMmcOONN/Lss88mvBBEXl4eL7/8Mg888ABf/epXWb58Od27d2fQoEFMnz6d4uJifvWrX9G9e3d27drF6NGjOf/885k5cyYPPfRQwvOy/O53v2PFihW8/fbbbN68mdGjRzNx4kQA3nrrLd5991369u3L+PHjeeWVVxqd0TAd1OIWOQo0tByLi4v57LPPOP300xu9PnXqVF544QWeeOKJZqc9vfLKK1m1ahVf//rXWbJkCWPHjo2dw6SBu/P973+fESNGMGXKFDZs2MCnn34KwIABA2KnbL300ktZunRpbL6LL744dv/aa6/Fnv/6179OdnY2AEuXLuWyyy4DYPLkyVRXV1NTU0NBQQG//OUvOf3007n++utjp3FtquELpKysjGHDhtGnTx86dOjAsccey7p16wB48MEHGTlyJGPHjmXdunVJT2S1dOlSLr74YrKzs+nVqxennHIKb7zxBgBjxoyhf//+ZGVlUV5ezpo1a1pd1oFQcIscBRr6uNeuXUttbW2zE0Tl5eXxxS9+kfvuu4/zzz+/2fx9+/blqquu4umnnyYnJ6fZRRTmzZtHVVUVy5cvZ8WKFfTq1St2WthEp3pNNh1/6tVEXQ0NZd955x2Ki4tjvw4S6dChAxCc/bBhuuFxXV0dS5YsYfHixbz22mu8/fbbnHDCCUmvedla90f8e2RnZx9QP30yCm6Ro0iXLl148MEHuffee9m7d2+j177zne9w9913U1xc3Oj55557Llb2k08+obq6mn79+jUqU1NTQ8+ePcnNzeWll15i7dq1sdf+/ve/x1rTTzzxRKNug4a+4QULFnDSSSclrPPEiROZN28eEPQh9+jRg6KiItauXct9993HW2+9xaJFi3j99dcP5E9CTU0N3bp1o6CggPfff7/RZdlyc3Ob/Z0a6rRgwQLq6+upqqri5ZdfZsyYMQf0/gdCfdwih1h7X6j4hBNOYOTIkcyfP58JEybEnh82bFij0SQNnn/+eW666abYxXt/8pOf0Lt370ZlLrnkEs4++2wqKiooLy/nuOOOi702dOhQHn30Ua655hoGDx7MddddF3ttz549nHjiiezbt48nnngiYX3vuOMOrrzySkaMGEFBQQGPPvoo7s7VV1/NvffeS9++fZk7dy5XXHFFrLuiLaZOncqsWbMYMWIEQ4YMYezYsbHXpk2bxogRIxg1alTsywPga1/7Gq+99hojR47EzLjnnnvo3bs377//fpvf/0DotK5yVNBpXdtHa9enLC0tpbKykh49erRDzdqXTusqInKUSRrcZjbEzFbE3baa2c2HonIiEm2lpaUJW9sQtMaPxtZ2OiTt43b31UA5gJllAxuAha3OJCKNuHuz0RVydEpH93Rbu0pOA/7q7muTlhQRIDg4pbq6OiNH0Em0uDvV1dWxHb0Hqq2jSi4CEu76NbNpwDSAY4455qAqJXIk6d+/P+vXr6eqqqq9qyKHgfz8/EanEDgQKQe3meUB5wD/kuh1d58NzIZgVMlB1UrkCJKbm8vAgQPbuxpyBGlLV8mZwJvu/mmmKiMiIsm1JbgvpoVuEhEROXRSCm4zKwBOB36X2eqIiEgyKfVxu/tOoDhpQRERyTgdOSkiEjEKbhGRiFFwi4hEjIJbRCRiFNwiIhGj4BYRiRgFt4hIxOjSZXJUas8r4ogcLLW4RUQiRsEtIhIxCm4RkYhRcIuIRIyCW0QkYjSqRI5ITUeNiBxJ1OIWEYkYBbeISMQouEVEIibVS5d1NbMnzex9M1tlZidlumIiIpJYqjsnHwCec/cLzCwPKMhgnUREpBVJg9vMioCJwBUA7l4L1Ga2WiIi0pJUukqOBaqAR8zsLTObY2admhYys2lmVmlmlVVVVWmvqIiIBFIJ7hxgFPALdz8B2AHMaFrI3We7e4W7V5SUlKS5miIi0iCV4F4PrHf318PHTxIEuYiItIOkwe3unwDrzGxI+NRpwHsZrZWIiLQo1VElNwDzwhElHwFXZq5KIiLSmpSC291XABUZrouIiKRAR06KiESMgltEJGIU3CIiEaPgFhGJGAW3iEjEKLhFRCJGwS0iEjEKbhGRiFFwi4hEjIJbRCRiFNwiIhGj4BYRiRgFt4hIxCi4RUQiRsEtIhIxCm4RkYhJ6UIKZrYG2AbUA3XurosqiIi0k1QvXQYwyd03Z6wmIiKSEnWViIhETKrB7cDzZrbczKYlKmBm08ys0swqq6qq0ldDERFpJNXgHu/uo4AzgX82s4lNC7j7bHevcPeKkpKStFZSRET2Sym43X1jeL8JWAiMyWSlRESkZUmD28w6mVnnhmngS8DKTFdMREQSS2VUSS9goZk1lP8Pd38uo7USEZEWJQ1ud/8IGHkI6iIiIinQcEARkYhRcIuIRIyCW0QkYhTcIiIRo+AWEYkYBbeISMQouEVEIkbBLSISMQpuEZGIUXCLiESMgltEJGIU3CIiEaPgFhGJGAW3iEjEKLhFRCJGwS0iEjEKbhGRiEk5uM0s28zeMrNnMlkhERFpXVta3DcBqzJVERERSU1KwW1m/YGvAHMyWx0REUkm1Rb3T4FbgX0tFTCzaWZWaWaVVVVVaamciIg0lzS4zewsYJO7L2+tnLvPdvcKd68oKSlJWwVFRKSxVFrc44FzzGwNMB+YbGaPZ7RWIiLSoqTB7e7/4u793b0UuAh40d0vzXjNREQkIY3jFhGJmJy2FHb3JcCSjNRERERSoha3iEjEKLhFRCJGwS0iEjEKbhGRiFFwi4hEjIJbRCRiFNwiIhGj4BYRiRgFt4hIxCi4RUQiRsEtIhIxCm4RkYhRcIuIRIyCW0QkYhTcIiIRo+AWEYkYBbeISMQkvQKOmeUDLwMdwvJPuvu/ZrpiIodS6Yz/iU2vmfmVdqyJSHKpXLpsDzDZ3bebWS6w1MwWufufM1w3ERFJIGlwu7sD28OHueHNM1kpERFpWUp93GaWbWYrgE3AC+7+eoIy08ys0swqq6qq0l1PEREJpRTc7l7v7uVAf2CMmQ1PUGa2u1e4e0VJSUm66ykiIqE2jSpx9y3AEmBqRmojIiJJJQ1uMysxs67hdEdgCvB+pismIiKJpTKqpA/wqJllEwT9f7r7M5mtlkjbxA/nEznSpTKq5C/ACYegLiIikgIdOSkiEjEKbhGRiFFwi4hEjIJbRCRiFNwiIhGj4BYRiRgFt4hIxCi4RUQiRsEtIhIxCm4RkYhRcIuIRIyCW0QkYhTcIiIRo+AWEYkYBbeISMQouEVEIiaVS5cNMLOXzGyVmb1rZjcdioqJiEhiqVy6rA74jru/aWadgeVm9oK7v5fhuomISAJJW9zu/rG7vxlObwNWAf0yXTEREUmsTX3cZlZKcP3J1xO8Ns3MKs2ssqqqKj21ExGRZlIObjMrBH4L3OzuW5u+7u6z3b3C3StKSkrSWUcREYmTUnCbWS5BaM9z999ltkoiItKaVEaVGDAXWOXu92e+SiIi0ppUWtzjgcuAyWa2Irx9OcP1EhGRFiQdDujuSwE7BHUREZEU6MhJEZGIUXCLiESMgltEJGIU3CIiEaPgFhGJGAW3iEjEKLhFRCJGwS0iEjEKbhGRiFFwi4hEjIJbRCRiFNwiIhGj4BYRiRgFt4hIxCi4RUQiRsEtIhIxCm4RkYhJ5ZqTvzKzTWa28lBUSEREWpdKi/v/AVMzXA8REUlR0uB295eBzw5BXUREJAVp6+M2s2lmVmlmlVVVVelarIiINJG24Hb32e5e4e4VJSUl6VqsiIg0oVElIiIRk9PeFRA53JTO+J9Gj9fM/Eo71UQksVSGAz4BvAYMMbP1ZnZ15qslIiItSdridveLD0VFRNqqactY5GihPm4RkYhRcIuIRIyCW0QkYhTcIiIRo+GAcnhxh7rdULsDarfDnu3h9LbwPrzV7+Xa7LfJpY48qyOXOjqwl1zqyKWeXKsji31Y3KINb3QPsI8s6sim1nPYS04wTTC914PHvPoR5BZAXiHkFUBeJ8jtFNzH33ILwAyRTFNwS2bs2we7t8D2TbDr87jbZ00eh7edn8PumiCsvT6lt5iRG76VG7XkUEtueJ9DnWdTH/eD0sMIdxoHaxb7yLV6cqgPQ7+OPOrIoZ4c2xcUen5BautsWdChCDp2g45dIb9rcN+x2/7p/PBxpx7QqWdwn99FgS9touCW1LkHIbv90yCQd1Q1md4UPN5RFdz21SVejmUH4VXQPbjv3Ad6Hh8EWF5h2IIthA6FcS3azo1bt9l5HPdvL7KXHOrJzsjqZrGPXOpYfcdk2LurcYu/dgfsbfJ4z7bgy2rXlv33NevCL6ctLX8hZedBp5K4MC+BwpLwuZ7BdOe+0Lm3Ql4ABbc02FcfhO7WjbB1Q9z9x/unt30M9bXN583KhcKewa1zH+gzAgp77W9RFhSHrdDw1qFzWsJnNx0Oehmt2UcWe8gLwjK/y8EtzD34NbFrS/gLYzNsr9r/Jddw274JNq2CHZsS/61zC4K/cVEY5LHpPuF0HyjsDTl5B1dfOawpuI8GdbVB6DYK5SbT2z9t3iLM7hCEQlE/GHDi/rAo7Bm2BMNbfle1ApMxC76wOnSGrgOSl3eHPVthx2bY9kmw/bZ9HHyRNkyvW9byl2mnknDb9YcuDbd+0GVAMF3YC7Iy80tFMk/BHXXuwYd763qoSXDbuiH44MftkAOCroiifsGHe9Ck8EPed/9zRf2C1rECuX2Y7W/pFw9quZw77PwsLtg3xn1Jb4TPPoK/vRzs3I2XlRN0v8QCPQz3+KBXt8xhS8F9uKvd0TyIa9YHfac1G4LHdbsbz5PTcf8H8h9OC1pZ8YFc1Bfyi9pnfSS9zKBTcXDrPbzlcrtrwv+bDcH/Tuz/aH3Qcn93YfN9EnmFcYEe11pvCPqifpCT2e4qSUzB3Z7q68LWUZMwjoX0+qA/NJ5lBX2YXfpDn5Fw3Jf3f6AaPlwF3dVSksYaWu+9hiV+fV990L++dUOT/8Mw5D9+O+iDbyq+S6aobxDqTRsJCve0U3BnSu3O5j9dE937vsbz5XfZH8QDxoQtnAH7Wzmd+0B2bvuskxy5srKDHZtFfaB/ReIye3c3bqlv3Rg0LrZuhM/XwNpXgtE0TSnc007B3Vb1dbCzOtiZ19Ba3voxbNu4f8fR1o2J/4HzOgcfjM59YODEBP2K/YKdVyKHo9z8oL+9tT73PduDz0As2MPuvJoNYbgvDbptmmoI9859oXOv4FdlYc9wZ3jv4LlOPTVaJqTghuBn4s7qxuOQG01/Ggzd2v5pUK7pjj4s2Etf1Ae6DYQvjAuHZvULg7pvcK9QPig6jWsEdCiEDoOhx+CWy+zZ3jjUt26M23+zDjZUBjvcm33OgI7dwzDvtX+EU0OwF/aCgnD4aUH3I3rUzJEX3Pv2wZ6aYE/7rs+DoN35WXC/67P9jxte27E5GFPbtMsCICc//MfoBd1KYcDocHxySRjU4fjZwl6QfeT9KUUyokMhlPxjcGtJ/d6g0bTtk7AR9Qls+zS4374peL76w6AxlWg4JBYe5FW8/1iChvuCHuF09/1B37FrsDM2IvuGDr+0aThQYXcN7N4ajGWNTde08PzW4MCGndVBILd0hJplh0frdQ/uux8b9Oc1HCzSMC65IZzTdKCIRJsuZdYOsnP3D1FtTaOjecNfxDuqg8ZYrGFWDdV/DUbP7KxuPR/yu+w/NUH8dMfwcdPpDkX7x+fndjxkeZFScJvZVOABIBuY4+4zM1Kbn/xD+IdN0PqNl5Ub/uGKgvsORVAyJPw2LY4L53C64XGHIsjSCRFFjhhm+z/jPYcmL99wDp2d1Y2DfXdN49MV7K6JO2VB+FxLp3CI1SUbuh4DN61Iz7q1Imlwm1k28DBwOrAeeMPM/svd30t7bcq/GYZyXCDnF+3/ZmsI65x8tYRFpO2ysvYHPa30wzflDnt3Nj4PTcNJ0fZsDc5Ts2dbcGDTIZDKu4wBPnT3jwDMbD7wVSD9wX36D9O+SJF0U9fJUchs/wnOuvRr79qkFNz9gHVxj9cDJzYtZGbTgGnhw+1mthroAWw+2EoeprRu0ZT2dbO707m0A6ZtFk3x6/aFVGdKJbgT9Uk0G6fj7rOB2Y1mNKt09xZG80eb1i2ajtR1O1LXC7RuiaSyp249EH86s/7Axra+kYiIpEcqwf0GMNjMBppZHnAR8F+ZrZaIiLQkaVeJu9eZ2fXA7wmGA/7K3d9NcfmzkxeJLK1bNB2p63akrhdo3Zox9wSHlYqIyGFLR6OIiESMgltEJGIOOrjNLN/MlpnZ22b2rpn9W4Iyp5pZjZmtCG//52Df91Ays2wze8vMnknwmpnZg2b2oZn9xcxGtUcdD0SS9YrsNjOzNWb2TljvygSvR3mbJVu3KG+3rmb2pJm9b2arzOykJq9HcrulsF5t3mbpOD5zDzDZ3bebWS6w1MwWufufm5T7k7uflYb3aw83AauARNf7OpPg2NnBBAcm/YIEBygdplpbL4j2Npvk7i0dtBHlbQatrxtEd7s9ADzn7heEI9gKmrwe1e2WbL2gjdvsoFvcHtgePswNb0fMHk8z6w98BZjTQpGvAr8O/w5/BrqaWZ9DVsEDlMJ6Hckiuc2OZGZWBEwE5gK4e627N70aSeS2W4rr1WZp6eMOf3KvADYBL7j76wmKnRR2pywysxYufHdY+ilwK9DSKQsTnRKg/U9mkFyy9YLobjMHnjez5eGpGJqK6jaD5OsG0dxuxwJVwCNh990cM+vUpEwUt1sq6wVt3OoZoMEAAAHLSURBVGZpCW53r3f3coKjKseYWdPLTb8JfMHdRwI/A55Kx/tmmpmdBWxy9+WtFUvw3GH9iyPF9YrkNguNd/dRBD+t/9nMJjZ5PXLbLE6ydYvqdssBRgG/cPcTgB3AjCZlorjdUlmvNm+ztI4qCX8CLAGmNnl+a0N3irs/C+SaWY90vneGjAfOMbM1wHxgspk93qRMFE8JkHS9IrzNcPeN4f0mYCHBGS7jRXGbAcnXLcLbbT2wPu7X+pMEgde0TNS2W9L1OpBtlo5RJSVm1jWc7ghMAd5vUqa3WXACbTMbE75v9cG+d6a5+7+4e393LyU41P9Fd7+0SbH/Ai4P93iPBWrc/eNDXde2SGW9orrNzKyTmXVumAa+BKxsUixy2wxSW7eobjd3/wRYZ2ZDwqdOo/mpoyO33VJZrwPZZukYVdIHeNSCCy5kAf/p7s+Y2bVhxWcBFwDXmVkdsAu4yCN8yGaTdXsW+DLwIbATuLIdq3ZQjpBt1gtYGH4OcoD/cPfnjpBtlsq6RXW7AdwAzAtHXnwEXHmEbLdk69XmbaZD3kVEIkZHToqIRIyCW0QkYhTcIiIRo+AWEYkYBbeISMQouEVEIkbBLSISMf8fS4tWul9TExQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "grid = np.linspace(post_μ - 2*post_σ, post_μ + 2*post_σ, 1000)\n",
    "\n",
    "plt.hist(results[:, 0], density=True, label=\"RMS approximation\")\n",
    "plt.plot(grid, post_dist.pdf(grid), label=\"True posterior\")\n",
    "plt.title('Posterior for $\\mu$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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
}