{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cfb1c1d1",
   "metadata": {},
   "source": [
    "A toy program to generate a sequence of random\n",
    "numbers distributed according to ANY probability\n",
    "distribution using the Metropolis-Hastings algorithm.\n",
    "\n",
    "The probability distribution of interest is myPDF(x)\n",
    "It does not need to be normalized.\n",
    "It is hardcoded.  Easy to change.\n",
    "\n",
    "The user needs to supply a few things.\n",
    "\n",
    "First:\n",
    "A pdf that is used to move in the\n",
    "sequence of random numbers from x to y.\n",
    "This is the so called \"proposal distribution\" \n",
    "P(y|x).\n",
    "Note: this is NOT the probability of the \n",
    "nth random number being y if the nth-1 is x.\n",
    "But it is related to it.....\n",
    "In principle P(y|x) can be ANY function as long as it\n",
    "is non-zero in the domain of validity of myPDF. \n",
    "The proposal function is of the Metropolis type, ie,\n",
    "P(y|x) = P(x|y)\n",
    "It is hardcoded as proposal(x) and returns y\n",
    "Note that the returned value is still just a \"proposal\"\n",
    "so it can be accepted or rejected\n",
    "\n",
    "Second:\n",
    "An initial value of x.  \n",
    "It can be anything but it must have myPDF(x) > 0\n",
    "It probably makes sense for this to be somewhere in the bulk\n",
    "of the distribution of interest.\n",
    "\n",
    "Third: \n",
    "The number of members of the chain to generate  (n)\n",
    "\n",
    "Fourth:\n",
    "The number of initial members of the chain to \n",
    "ignore (the so called burn-in, nBurn)\n",
    "\n",
    "To visualize the chain, we fill a histogram of members of the chain.\n",
    "The parameters of the histogram are hardcoded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1d345a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f882bdd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PDF\n",
    "# Harcode a gaussian \n",
    "# No need to normalize!!!!\n",
    "def myPDF(x, mean, sigma):\n",
    "    return np.exp( -(x-mean)**2 / (2*sigma**2) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "380ee45f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Proposal\n",
    "# Uniform random step between -1 and 1 (simple, a bit stupid..)\n",
    "# It can be any function, does not need to be flat/\n",
    "def proposal(x):\n",
    "    step = 1.\n",
    "    return x - step + 2*step*np.random.rand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ba7261d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize random number\n",
    "np.random.seed(1629871)\n",
    "\n",
    "# Parameters of the gaussian\n",
    "mean  = 3\n",
    "sigma = 1.5\n",
    "\n",
    "# Parameters of the chain\n",
    "xstart = 0.\n",
    "n      = 10000\n",
    "nBurn  = 1000\n",
    "\n",
    "# initialize chain as list\n",
    "xlist = [xstart]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d89d2fe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# do n+nBurn-1 steps  (we alread have 1 step...)\n",
    "for i in range(n+nBurn-1):\n",
    "    xp    = proposal(xlist[-1])\n",
    "    fnow  = myPDF(xlist[-1], mean, sigma)\n",
    "    fnext = myPDF(xp,        mean, sigma)\n",
    "    if np.random.rand() < fnext/fnow:\n",
    "        xlist.append(xp)\n",
    "    else:\n",
    "        xlist.append(xlist[-1])\n",
    "        \n",
    "# put list into array removing the burn in part\n",
    "xr = np.array( xlist[nBurn : ] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e04f1e8b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAD1CAYAAACiJBXjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApGUlEQVR4nO3deZyNdf/H8ddnBpMlScauGQmlhTSV7nJXtMitUFqVZdSQSqko1V13q1YVlUzZQpZocaOVlp/c1kJJZWkme2TIzpjP74/rwmCYMzPnzPcsn+fjcR5zzvdcc877lPnMNd/ru4iqYowxJrLFuQ5gjDGm6KyYG2NMFLBibowxUcCKuTHGRAEr5sYYEwWsmBtjTBQo4eqNK1WqpMnJya7e3hhjItK8efM2qGrioe3OinlycjJz58519fbGGBORRCQzr3brZjHGmChgxdwYY6KAFXNjjIkCVsyNMSYKWDE3xpgoYMXcGGOigBVzY4yJAlbMjTEmClgxN+YokpNB5MDNJi2bcOVsBqgxkSAzE3JvxiXiLosxR2Nn5sYYEwXyLeYicoyIzBaRBSKySESe8NuHicjvIjLfvzXy20VE+ovIUhFZKCKNQ/wZjDEm5gXSzbILaKaqW0WkJDBdRD7xn+ulquMPOf5KoK5/Ow8Y6H81xhgTIvmematnq/+wpH/To3xLa+Bd//tmAhVEpFrRoxpjjDmSgPrMRSReROYDfwJfqOos/6ln/K6UV0QkwW+rAazI9e0r/TZjjDEhEtBoFlXdCzQSkQrAhyJyOtAHWAuUAtKBB4EnA33j9evXk5KSku9xaWlppKWlBfqyxhgTkdLT00lPTw/k0Ep5NYrq0XpM8vgGkceA7ar6Uq62i4EHVLWViAwCvlbV0f5zvwIXq+qa3K+TkpKitjmFCXcihw9NLOCPjDFBJSLzVPWwM+FARrMk+mfkiEhp4DLgl3394CIiQBvgJ/9bJgId/FEtTYDNhxZyY4wxwRVIN0s1YLiIxOMV/3GqOklEpolIIiDAfKCbf/wUoCWwFNgOdA56amOMMQfJt5ir6kLgrDzamx3heAXuLHo0Y4wxgbIZoMYYEwWsmBtjTBSwYm6MMVHAirkxxkQBK+bGFFLutc5tnXPjmq1nbkwh5V7r3NY5N67ZmbkxxkQBK+bGGBMFrJgbY0wUsGJujDFRwIq5McZEASvmxuSmamvcmohkxdzEtu3bYcgQuOYaOPFEiI+HuDioXBkuvZSe9IO1a12nNCZfNs7cxKT6STu5/o8X6ckrVCSLlfEnUvPGf3qzf+LiYNUqmD2bftwPSX2gUyd4+mkg0XFyY/JmxdzEnhkzmPxHB05mGbRpA/ffT62mF6AjD5/5U09+47fbXoP0dHj/fW5mANC+2CMbkx/rZjFR7bAp9wMGwEUXISh8+SV8+CFceCHeHiuHW0I9eOMNWLAAGjRgFLdA9+6we3cxfgpj8mfF3ES1fVPuNUfpndUHevTg4+yWtK01D5o3D/yFGjSAr7/mrfK9YOBAJiVcQ71aO0IX3JgCsmJuop8q3HMP3f9+Drp1o/XeD1n4R4WDDklKCmDRrBIl6Lb5BRg0iFYyhV9P/hfs3Bnq9MYExIq5iX5PPeV1r9x/P7z5pneB8xAZGQdGJWZm5vN6aWkwYgR89RW0bw9794YktjEFYRdATVTrxFB4/HHo2BFefDF4yxu2bw8bNsC998I99wCvB+d1jSmkfM/MReQYEZktIgtEZJGIPOG31xaRWSKyVETGikgpvz3Bf7zUfz45xJ/BmLzNmsUgusKll3qjUYK9Tu0998ADD8Abb9CRYcF9bWMKKJBull1AM1VtCDQCWohIE+B54BVVPRnIArr4x3cBsvz2V/zjjClef/0F11/PaqrD2LFQqlRo3qdvX2jWjIHcAfPnh+Y9jAlAvsVcPVv9hyX9mwLNgPF++3CgjX+/tf8Y//nmIrZ0vylGqt4kn7Vracd4qFgxdO9VogSMHs1fnADXXQfbtoXuvYw5ioAugIpIvIjMB/4EvgCWAZtUNds/ZCVQw79fA1gB4D+/GTghiJmNOaLkZEiNGwqTJtFj94tsSEoJ/ZtWrswtjISlS+Ghh0L/fsbkIaALoKq6F2gkIhWAD4FTivrG69evJyUl/x+0tLQ00tLSivp2JkbszVzBkPI9ofHF9J96F/2LabzWN1zs9aG/9hq0bQvNmhXPG5uokZ6eTnp6eiCHVsqrUbSAK8SJyGPADuBBoKqqZovI+cB/VPUKEfnMv/8/ESkBrAUS9ZA3SklJ0blz5xbovY05KlU+jbuSFmWnw8KFcNJJhXoZkYP39gzkR0QEdNt2OOssb+z5okVQrlyh3t+YoxGReap62JlwIKNZEv0zckSkNHAZsBj4CmjnH9YR+Ni/P9F/jP/8tEMLuTEhMWYMLfgMnnuu0IW8SMqUgaFD4Y8/vLHtxhSjfM/MReRMvAua8XjFf5yqPikiJwFjgIrAD8AtqrpLRI4BRgBnARuBG1V1+aGva2fmJqi2boX69Zm3uipnZ8/2lrItpEKfme87rksXePfd/eu5GBNMRzozL3A3S7BYMTdB9dBD8PzznM8M/qfnF+mlilzM16+H+vWhUSOYOjX449tNTCt0N4sxYe+336BfP+jUiZkUrZAHRWIiPPusN91//Pj8jzcmCKyYm8j3wANQurTXVx4ubr8dTj8dHn4Y9uxxncbEACvmJrJNnw7//a/XzVKlSlBeMvcKiklJhXyR+Hh4/nlv7Hlgw82MKRLrMzeRSxWaNoXly72iWaZMwH3cwZbn+6pCs2Zs+HYRJ+UsZQvlAe8XREZGsUc0UcL6zE30mTwZvvsOHnvMGxYYbkTghReolLOevx99MfAldo0pBCvmJjLt3ev1R598sjcUMFydcw5juR5efdVb/MuYELFibiLTBx/Ajz96k3NKlnSd5qie5DFvAa5+/VxHMVHMirmJPDk58PTTcMop3kqFYe5nToN27bzdjjZudB3HRCkr5ibyTJrkrb3y8MNFmulZrP79b9iyxetuMSYErJibyKLqnZWfdBLcdJPrNEeUnHzI8MYzzoBrroHXXqMCWa7jmShkxdxEli++gDlzoE8fb2OIMJWZeWCD6P3DEB97DP7+mx70dxnNRCkr5iZyqHoXPGvWhA4dXKcpuIYNoVUr7uJ12L7ddRoTZayYm8jx3XfejM/evUO3p2ch5Z41etSZo716kcgGGD78CAcYUzhWzE3Yyt3vnJyMN7SvYsWwHFeekXGgW+WgrpVDNW3KLM6Fl1/2xsobEyRWzE3Yyt3vHJ+5DD76CO64IzxnewZKhBfpBcv8z2NMkFgxNxHhHl7zLnjeeafrKEX2IW2hTh148UU3C8mYqGTF3IS/rCxSGQI33wzVqrlOU2S1kuLpvuw+mDWLdtW+cx3HRAkr5ib8vf025dgGPXu6ThIUGRnw5rZOcMIJtF/3sus4JkpYMTfhbc8e6N+fL2nuDe2LFmXKQFoaVzPRllE0QWHF3IS399+HVavox31HPOSw2ZaR4o47vK9vvuk2h4kK+RZzEaklIl+JyM8iskhE7vHb/yMiq0Rkvn9rmet7+ojIUhH5VUSuCOUHMFFuwACoV49PaXHEQ/KcbRkJatXiI9rAO+/Ajh2u05gIF8iZeTZwv6o2AJoAd4pIA/+5V1S1kX+bAuA/dyNwGtACeFNEImQ1JBNOzuJ7mDkTundHo/SPyAHc7a2kOHq06ygmwuX7E6Kqa1T1e//+FmAxUOMo39IaGKOqu1T1d2ApcG4wwprYcgcDvY2aO3Z0HSVkvuEib+PnAQNsmKIpkgKtVCQiycBZwCzgAuAuEekAzMU7e8/CK/Qzc33bSvIo/uvXrycl5bBt7A6TlpZGWlpaQWKaaJCVRXtGQfv2UKGC6zQhk5QkpP10N+l0pV217xi/9kLXkYwj6enppAe2+XelvBoDLuYiUg6YANyrqn+LyEDgKUD9ry8DqYG+XmJiIrahszmi4cMpw479k4T2rX2yT7RsipyRAWxrDzUf5Lp1AwAr5rEq0BNXEdmQV3tAHZEiUhKvkI9S1Q8AVHWdqu5V1RzgbQ50pawCauX69pp+mzGBycmBN99kBudDo0bA4WufRNVovrJlITWVa5kAa9a4TmMiVCCjWQQYDCxW1X652nNPxWsL/OTfnwjcKCIJIlIbqAvMDl5kE/WmTYMlS3iT7kc8JPcqhRE1HPFIunalBHth6FDXSUyECuTM/ALgVqDZIcMQXxCRH0VkIXAJ0BNAVRcB44CfgU+BO1XVloczgXvjDahUifG0O+Ihuc/Uo6G7hXr1mMYl8Pbb3l8mxhRQvn3mqjodkDyemnKU73kGeKYIuUysWrECJk6EXr3Y9fwxrtMUq0F0pVnGjd5uSlfY9AxTMNE5eNdErnfe8U63u3Z1naTYfUQbqFQJBg1yHcVEICvmJnzs9fuML7sMatd2nabY7SYBOnf2/jJZvdp1HBNhrJib8DFtmtfNkhrwCNfoc/vtB36pGVMAVsxN+Bg82NsWrk0b10mcSEoCqVeXqTQj49G3OSnJxg2YwFkxN+Fh40b48ENvxmdCgus0TuwbodN8bFeSyaT+H5+7jmQiiBVzEx7eew92747tLpZ92rSBxETSCGhqtzGAFXMTLoYMgbPO2j/jM6aVKgWdO3MV/7ULoSZgVsyNez/84N26dIncjSaC7bbbvBmh777rOomJEFbMjXtDhnj95DfdFLkbTQRb3bp8S1Pvv40tjWsCYMXcuLVzJ4waBW3beiNZzH5DSIUlS+C771xHMRHAirlx6+OPISvLLnzmYTztoFw57+zcmHxYMTdOfdtpMJmcSNzlza2f/BDbKAc33ADjxsGWLa7jmDBnxdy4k5nJhTu/JOnxzuRonPWT56VLF9i2zSvoxhyFFXPjzvDh3tdOnZzGCGtNmsApp1hXi8mXFXPjRk4ODB3KVJpDcrLrNGEpKQkkTuj1SyrMmAG//OI6kgljVsyNG199BRkZ3ogNk6d90/tfXHMr2cTz/KlDEbHffSZvVsyNG0OGQIUK3hre5uiqVqXE1f/iwSrD0d17omv/UxM0VsxN8cvKggkT4Oab2Ulp12kiQ2oqrFsHn37qOokJU1bMTfEbMwZ27fJGapjAtGwJVarYhVBzRPkWcxGpJSJficjPIrJIRO7x2yuKyBcissT/erzfLiLSX0SWishCEWkc6g9hwl/uNVcW3jsYGjb0FtYygSlZEjp0gEmTqMw612lMGArkzDwbuF9VGwBNgDtFpAHwEDBVVesCU/3HAFcCdf1bGjAw6KlNxNm/5sr8BZy5e57XbSB57RNujqhzZ8jO5hZGuk5iwlC+xVxV16jq9/79LcBioAbQGvAHCjMc9l/Jag28q56ZQAURqRbs4CZCDR3KLkp5m1CYgjn1VDj/fLow2BbfMocpUJ+5iCQDZwGzgCqqusZ/ai1Qxb9fA1iR69tW+m0m1u3aBSNGeCNYTjjBdZrIlJpKAxbDrFmuk5gwUyLQA0WkHDABuFdV/5ZcfyKrqopIgU4V1q9fT0pKSr7HpaWlkZaWVpCXNuFq4kTYuJEhpHKD6yyR6vrr2Xb7PZQdOtSbHWqiRnp6OunpAe0uVSmvRtEA/lwTkZLAJOAzVe3nt/0KXKyqa/xulK9Vtb6IDPLvjz70uNyvmZKSonPnzg0kuIkCIqAtroSffiJ+ZQZ7Nf5Au/UYFMgw6USnYz+AtWuhTBnXcUwxE5F5qnrYmXAgo1kEGAws3lfIfROBjv79jsDHudo7+KNamgCbDy3kJvbUZAV89hl06kQO8a7jRLQpVVJhyxZuLTvBZoSa/QLpM78AuBVoJiLz/VtL4DngMhFZAlzqPwaYAiwHlgJvA92DH9tEmo4M907BO3d2HSXijVvTFE4+mREXD0EVmxFqgAD6zFV1OnCkMWTN8zhegTuLmMtEk5wcOjMULrkETjrJdZrIJ+L9UnzkEVi2DKjjOpEJAzYD1ITet99Sh+W2m1AwdegAcXEwbJjrJCZMWDE3oTdkCJs4Dq69FvCXdvVng9rOQoVUsyZccQUMG0Yce12nMWHAirkJrc2bYfx4RnMTlPYW1dq3tKvtLFREqamwciWX8qXrJCYMWDE3oTVmDOzYYeuWh8JVV0HFiqRii28ZK+Ym1IYMgdNPZy75TxAzBZSQALfcQhs+gr/+cp3GOGbF3ITOTz/B7Nn+Ure2qFZIpKaSwG547z3XSYxjVsxNSCQnQ78zhrCbkiT2bG8XOkOlYUPm0djWOTdWzE1orM7czX2VRlDq2qtZr4l2oTOEhpAK8+fDDz+4jmIcsmJuQqIVk2DDBhtbXgxGc5PXfz50qOsoxiEr5iYkujAYatSAyy93HSXqZVER2raFkSOpn7Rz/xh+W7MltlgxN8G3ciUt+BQ6dYISAa+ybAopKQkuG5MKWVlcun3i/jH8tmZLbLFiboJv2DDiybEulmKSkQFfZDeDE0/kjbPtQmissmJugisnB4YMYSrNbFGt4hQf7/0l9PnnsGJFvoeb6GPF3ATX11/D778zmC6uk8SeTp28/pXhw/M91EQfK+YmuN55BypU4EPauk4Se2rXhmbNvFEtOTmu05hiZsXcBE9WFnzwAbRvz05Ku04Tm1JTYfly+PZb10lMMbNiboJn1CjYtQtuu811kth1zTVw3HE2IzQGWTE3waHqdbE0bgyNGrlOE7tKl4abboLx4ynPZtdpTDGyYm4KJDn5wMYSB01M+f57WLCA7t93sU0nXOvcGXbs4AbGuk5iipEVc1MgmZkHNpY4aGLK4MHs4BjezLrZNp1w7Zxz4LTTvH1XTczIt5iLyBAR+VNEfsrV9h8RWSUi8/1by1zP9RGRpSLyq4hcEargJoxs3w7vvcd42kGFCq7TGBFITeV8ZsLPP7tOY4pJIGfmw4AWebS/oqqN/NsUABFpANwInOZ/z5siEh+ssCZMTZgAmzfb2PJwcsst7KEEDB7sOokpJvkWc1X9FtgY4Ou1Bsao6i5V/R1YCpxbhHwmEgweDHXq8A0XuU5i9qlcmY9p7U0g2rnTdRpTDIqyCtJdItIBmAvcr6pZQA1gZq5jVvpth1m/fj0pKflvJZaWlkZaWloRYppQqsev8M038Mwz8IjtJhRO3qIb7f6a4I39v/lm13FMPtLT00lPTw/k0Ep5NYqq5vudIpIMTFLV0/3HVYANgAJPAdVUNVVEXgdmqupI/7jBwCeqOv7Q10xJSdG5c+cGEtyEERHvwuc+L8v93F+iP6xYgVSrSgD/nEwxiZMccurUg+rVbRJRFBGReap62JlwoUazqOo6Vd2rqjnA2xzoSlkF1Mp1aE2/zUSjHTvoxDBvokrVqq7TmEMocdC1K/zf/x12ITT3EFNb9zw6FKqYi0i1XA/bAvtGukwEbhSRBBGpDdQFZhctoglb48dzAhuhWzfXScyRdOoEpUrBoEEHFXDIY3ipiWj59pmLyGjgYqCSiKwEHgcuFpFGeN0sGUBXAFVdJCLjgJ+BbOBOVd0bkuTGvbfe4lfqUf/ii10nMUeSmAjXXgvvvsufm/qiWsZ1IhMi+RZzVb0pj+YjjndS1WeAZ4oSykSAhQthxgzeoh+viF34DGvdusHo0VzPOKCT6zQmRGwGqCmcQYMgIYHhdHSdxOSnaVM49VS68ZbrJCaErJibgtu6FUaMgBtu8DYTNmEpKcnvI48T7lnclSbMgvnzXccyIWLF3BTc6NGwZYtd+AxzGRkHLnK+trEDHHOM9xeViUpWzE0BKQwcCGeeCU2auA5jAnX88XDDDTBypPeL2EQdK+amQFKYCz/84J2V24XPyNKtm9dFNmqU6yQmBKyYmwK5i9ehXDlo3951FFNQ550HZ58NAwZgU3WjjxVzE7h167iRMd5ElPLlXacxBSUCPXp4s0GnTXOdxgSZFXMTuEGDSGA33HWX6ySmsG64ASpXhv79XScxQWbF3ARm924YOJApXAn16+9v3j/8zbaKiwwJCd56Lf/9Lyxf7jqNCSIr5iYw48fD2rX0p8dBzbmHv9lWcRGiWzeIj4c33nCdxASRFXMTmNdeg3r1+JzLXScxRVW9Olx3nbepyNathz1tKypGJivmJn+zZsHs2dCjh7esqol8PXrA5s3w7ruHPZV7025bUTFy2E+myV///t7olQ4dXCcxwXLeeXDOOTBgAEKO6zQmCKyYm6NbvRrGjYPUVDj2WNdpTLDsG6b4yy+0r/zl/m4Vu5AduayYm6N7/XXYu9eGI0aj666DKlUYcVa//d0qdiE7clkxN0e2ZYu3Dss110CdOq7TmGBLSPDOzj/7zFuf3kQ0K+bmyN55BzZtgl69XCcxoXLHHVC2LLz4ouskpoismJu87dkDr7wC//wnyTecZ/2p0er44+H222HMGPjjD9dpTBFYMTd5GzsWVqyA3r0PGqpm/alRqGdP73/uq6+6TmKKIN9iLiJDRORPEfkpV1tFEflCRJb4X4/320VE+ovIUhFZKCKNQxnehIgqvPACNGgAV17pOo0JtRNPhBtvhPR0yMpyncYUUiBn5sOAFoe0PQRMVdW6wFT/McCVQF3/lgYMDE5MU6w+/xx+/NHrK4+zP95iQq9esG0bvGX7hEaqfH9SVfVbYOMhza2B4f794UCbXO3vqmcmUEFEqgUpqykm37V5gVVUp1Tnm62fPFY0bAhXXOEt27Bzp+s0phAKe9pVRVXX+PfXAlX8+zWAFbmOW+m3mUgxcyYX7JxGjRfuZbeWsn7yWNK7N6xbB0OGuE5iCqFEUV9AVVVECrxtyfr160lJScn3uLS0NNLS0gqVzRTCk0+ynkok3nGH6ySmuF1yCfzjH9C3L3Tp4o1DN8UmPT2d9PT0QA6tlFdjYYv5OhGppqpr/G6UP/32VUCtXMfV9NsOk5iYyNy5cwv59iYk5syBTz7hZfryXLlyrtOY4iYCjz/udbcMHw52ElWsAj1xFZENebUXtptlItDRv98R+DhXewd/VEsTYHOu7hgT7p58EipW5A3udJ3EuHLZZd4iXM8+621IYiJGIEMTRwP/A+qLyEoR6QI8B1wmIkuAS/3HAFOA5cBS4G2ge0hSm+D7/nuYNAl69mQrtqBWzBKBxx7z1r4dMcJ1GlMAoo526U5JSVHrZgkjbdvC119DRgZS4TjbvD2WqcK558Jff1Hy91/ZoyUBb6OK3OubJyXZxXEXRGSeqh52wdEGERtYsAA++gjuvReOO851GuPavrPz33/nFkbub849E9g2rgg/VswNPPIIVKjgraBnDECrVnD22TzOE7Brl+s0JgBWzGNYcjJcKNNh8mSe4yFv0SVjwDs779uXZDJtVmiEsD7zGCai6AVNYflyyqxZynYt47djfeYGVJkadynNKy2E5cuR8sce9O/C/p24YX3m5jAtmQLffQePPUblpDK2zK05mAivVX0ONmzg8fL97N9FmLMz81iVk8PC+EacWWc7LF4MJUu6TmTCVbt23m5Ey5ZB5cr7m+3M3A07MzcHGzWKM/kRnn7aCrk5uqefhu3b4ZlnXCcxR2HFPBZt2wZ9+jCHFLj+etdpTLg75RRITfX2g/3tN9dpzBFYMY9Fzz8Pq1ZxL6/aeuUmME89BcccA/ff7zqJOQL7SY41mZne5r033sgMLnCdxkSKqlXh3//2lnz49FPXaUwerJjHmgcf9K5cPf+86yQm0vToASef7O0ZumeP6zTmEFbMY8n06d5Gzb16efs+GlMQCQnQrx/88gu88YbrNOYQNjQxVuzZA40bw+bN3lDEsmVtaJkpOFVo0QJmzaLa5sWsUdsVsrjZ0MRY9/LL8NNP8PrrULas6zQmUonAgAGwcyevcq/rNCYXK+axYPlyeOIJuOYauPpq12lMpKtXDx59lBsYB5Mnu05jfFbMo50qdO8OJUty3qz++6fs27R9UyS9e/NbyQZktupOOdlKcrLrQMaKeRRKTj5QsHtUHuNNxX7mGWavqnHQetS2sYAptFKlqPdVOkn8wdaej9na5mHAinkU2r+JwOo1PL7hLm9Px+62g58JsgsugK5d4bXXOJ8ZrtPEPCvm0UoVbr+d0uzwdlqPj3edyESjF16AE0/kvRIdKCdbEcG6XByxYh6thgyByZN5iOegfn3XaUy0Kl8ehg8nee9ytnZ9wLaTc6hIxVxEMkTkRxGZLyJz/baKIvKFiCzxv9r2NcUsiQxvlt4ll/A6d7mOY6LdP/8JDzwAgwbBlCmu08SsYJyZX6KqjXINYn8ImKqqdYGp/mNTXHbvZiw3eFc/hw7lxKQ4G71iQu+pp+CMM6BLFxrX/POgUVPW7VI8QtHN0hoY7t8fDrQJwXuYI+nTh/OY7XWzJCWRkWGjV0wxSEiAUaNg0ybm1b8Zzd67/9+ddbsUj6IWcwU+F5F5IpLmt1VR1TX+/bVAlSK+h8nHvqGIV8tE6NePYcfeBdde6zqWiTVnnOGt2TJ1Kjz5pOs0MadIa7OISA1VXSUilYEvgLuBiapaIdcxWap6WL95UlKSJiYm5vseaWlppKWl5XtcLBMBXbYcUlKgdm2YMcM7UzLGhc6dvRFUn34Kl19uawAFKD09nfT09HyPmzdvXqaqJh/aHrSFtkTkP8BW4HbgYlVdIyLVgK9V9bDhFLbQVvAcK1vYctr5sHo1zJkDdeq4jmRi2fbt3tyG1ath9mzk5DpWzIMo6AttiUhZETl2333gcuAnYCLQ0T+sI/BxYd/DBCAnh5Hc4i1L+v77VsiNe2XKwIcfevdbteI4NjmNEyuK0mdeBZguIguA2cBkVf0UeA64TESWAJf6j02oPPoorZkIr74KzZu7TmOM5+STYcIEWLrUG12Vne06UdQrdDFX1eWq2tC/naaqz/jtf6lqc1Wtq6qXqurG4MWNLbnXWDn0lpyMt8Fu376kczvceafjtMYc4uKL4a23uILP4e67reM8xGwGaBjbv8ZKHreUzPFeAW/Viu686VV4Y8JNly4MLP8gvPUWT8U9ZmPOQ6iE6wCmEKZNYxTt4R//gLFj2VvW/jea8HXHpr5w+wb+PfhpNmZWBHq6jhSV7Mw80nz9NVx1FUuoCxMnehebjAlnIt5U/3bteIX7IIDhd6bgrJhHkqlToWVLqF2b5kyFihVdJzImMPHxMHIkk2npLZv7+uuuE0UdK+aR4tNPoVUrb5TAtGn8aRNrTaRJSKDniR/wEa3h7rt59vgXXSeKKlbMI8HQoV4hP/VUmDYNKld2nciYQvktM4E2u9+H66/n4U29oVcvyMlxHSsqWDEPa+qtcZGaCs2aef3llSoB3gqIthqiiUglS8KoUbzJHfDSS3D99bBjx0FDcW3US8FZMQ9XW7cympvg8cehY0dvF/Ty5fc/bashmohWogQvnPgG9/EyORM+YGaZS6iavdJWWiwCK+bhaMkSaNKE63gf+vb1ullKlnSdypigysgU+ul9xH0wgSZlf2Lm7sbw5ZeuY0UsK+bhRJX7Kr3L5nopbFi0lg6VP4OHHrIJQSa6tW3rLRCXmAiXXw7/+Q8l2OM6VcSxYl5M8p2av349XHst/f7qyHFNG1IpYx6j1l3qOLUxxeTUU2H2bGjfHp54gv9xPixa5DpVRLFiXkyOODU/R7kwcyScfjpMnkwvXoCvvrKrmib2lC0LI0bAuHEkkQmNG3vdjLt2uU4WEayYu/Tjj3DRRYzkVu/0fM4cXqKXN8HCmFh13XWcxiK46ip4+GFvByPbKDpfMVPMA9nBo9hkZHjDDRs1gp9/5jbehv/9D848s1AvF1afLQTs80W2wny+9VSG8ePhk08gLg7+9S+48kqYNy8ECYsmbP7/qaqT29lnn63Fqbjf71DeosHLVbt3Vy1ZUjUhQbVnT9W//vKey31cAbn+bKFmny+yFebzHfRzsGuX6ksvqR5/vPfE1VerzpsXvIBFVNz//4C5mkdNjZkzc2dUYfp0xnOtNxU/Pd07K1+6FPr1g4oVbQKQMUdTqhTcf7/3F+1TT8G338LZZ8NFF9EtcQIlJNsmGxFD3SzFbu1anq74MgvjGkLTpjSP+wp694bff4e33oKaNfcfahOAjAlA+fLw6KPeD8lLL8Eff/DWhnZk1zoJffTf6OJfYnqykRXzYFq1ylvq88oroWZNHs16gDPPLQ0DB1Lh7xXelflcRdwYUwjHHeedqS9dSms+8oY1PvssnHoqcznbK/S//hpzOxvZrgZFsW0bzJgB33wDn30Gc+cCsIyTGEcvvqzekamzTnEc0pjIs6/r8ejiSUpqDZ+1hjVrYOxYcnqO8hbv6tXL29y8ZUu47DJvI5cTTiiO6M5YMQ9UdjYsXgw//ODdZs3yZq1lZ5NNPHM4h/+r8Cy9v2tNnVNPpY8IfVxnNiZCFbi7sVo1uPdernv1XjQzk5ZM4V/LJtNswDuUGTDAO6ZBA7jwQkhJgYYNvbkdUbS5S8iKuYi0AF4D4oF3VPW5UL1X0OTkeDMxly/31kdZssS7UPnbb/Dzz7BzJwA7pDTf61l8ywP8Uvkihi+9gPOPPZbzHcc3JtZ5vwSSgDu8286dMGcOff45nb7J02Hs2AM7HcXFQd263jj2OnUO3E46CWrVirj5HiEp5iISD7wBXAasBOaIyERV/TkU75en7GzYsQM2b4asLBpt2eJts7ZpE2RlwaZNjHhpHeW2rqE6q6nOaqqylpJkH3iNuDhITubb1SczZ2d3vqcxP3AWO2vVZ3lmPBcU24cxxhTKMcdA06Y8R1P6TsY7Yfv9d1i4EBYsgAULWP7RAmpmf0yp3OvBxMd7+wZUreqd9Vet6t0qV4YKFbx+e/9r9V27YONGOPZYpwviherM/FxgqaouBxCRMUBr4EAxX7XKG92xd++BW3b2wY/zatv3eNcu77fujh0H3/a1ZWcfFOgdgNatD2prSUVOOL06VK8O1RvwzLDqPDKgGtSuzSVpdfludTJ7lpciKclGmRgTyQ70wccBdfxb2/3PZSzbCytX0ix5GdPSl3k/8OvWeX3xa9fC/Pne4717D3vtiXCgPz4+HkqX9n6JlC594LbvcUIClCjh3eLjD9w/9HHu+3FxB25H+WUhGoIrviLSDmihqrf5j28FzlPVu3Ids4XARtOsBzYEIValIL1OOIrmzwb2+SKdfb7AXycxgONyVPXYQxudXQDNK4wxxpjCCdU481VArVyPa/ptxhhjQiBUxXwOUFdEaotIKeBG/K4lY4wxwReSYq6q2cBdwGfAYmCcqjpdaV5EXhSRX0RkoYh8KCIVXOYJFhFpISK/ishSEXnIdZ5gEpFaIvKViPwsIotE5B7XmYJNROJF5AcRmeQ6S7CJSAURGe//3C0WkagavSsiPf1/lz+JyGgROcZpnlBcAA1HInI5ME1Vs0XkeQBVfdBxrCLxh4D+Rq4hoMBNxToENIREpBpQTVW/F5FjgXlAm2j5fAAich+QApRX1Vau8wSTiAwH/k9V3/H/Qi+jqpscxwoKEakBTAcaqOoOERkHTFHVYa4yxczaLKr6uf8XA8BMvH78SLd/CKiq7gb2DQGNCqq6RlW/9+9vwfsrr4bbVMEjIjWBf+GPnI0mInIc8E9gMICq7o6WQp5LCaC0iJQAygCrXYaJmWJ+iFTgE9chgqAGsCLX45VEUbHLTUSSgbOAWY6jBNOrQG8gx3GOUKiNN6x4qN+N9I6IlHUdKlhUdRXwEvAHsAbYrKqfu8wUVcVcRL70+68OvbXOdcwjQDYwyl1SUxAiUg6YANyrqn+7zhMMItIK+FNVw2/rnOAoATQGBqrqWcA2IGqu6YjI8Xh/BdcGqgNlReQWl5miaqEtVT3qdvYi0gloBTTX6LhYEPVDQEWkJF4hH6WqH7jOE0QXAFeLSEvgGKC8iIxUVacFIYhWAitVdd9fUuOJomIOXAr8rqrrAUTkA+AfwEhXgaLqzPxo/IW/egNXq+p213mCJKqHgIqI4PW5LlbVfq7zBJOq9lHVmqqajPf/bVoUFXJUdS2wQkTq+03Nyb2cR+T7A2giImX8f6fN8a7pOBNVZ+b5eB1IAL7w/tszU1W7uY1UNP7InH1DQOOBIa6HgAbZBcCtwI8iMt9ve1hVbav2yHA3MMo/0VgOdHacJ2hUdZaIjAe+x+u2/QFwurNzzAxNNMaYaBYz3SzGGBPNrJgbY0wUsGJujDFRwIq5McZEASvmxhgTBayYG2NMFLBibowxUcCKuTHGRIH/B08qvZ+yJxT4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# now plot it...go out to 4 sigma\n",
    "low   = mean - 4*sigma\n",
    "high  = mean + 4*sigma\n",
    "nbins = 100\n",
    "bins  = np.linspace(low, high, nbins+1)\n",
    "ax    = plt.subplot(111)\n",
    "contents, bins, _ = ax.hist(xr, bins, histtype='step', color='blue')\n",
    "ax.set_xlim(bins[0], bins[-1])\n",
    "ax.tick_params(\"both\", direction='in', length=10, right=True)\n",
    "# Also show a gaussian curve of correct mean,sigma.\n",
    "# Properly normalized given number of entries and binsize.\n",
    "binsize = bins[1]-bins[0]\n",
    "xx      = np.arange(bins[0], bins[-1], 0.01)\n",
    "yy      = binsize * len(xr) * stats.norm.pdf(xx, mean, sigma)\n",
    "_ = ax.plot(xx, yy, color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "711a9a51",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}