poisson_estimator.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://lips.cs.princeton.edu/the-poisson-estimator/\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "# say we want to calculate e^E(X), where X ~ N(0, 1). That should be 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [],
   "source": [
    "# naive method\n",
    "def naive_method(N):\n",
    "    x = np.random.normal(0, 1, N)\n",
    "    x_hat = x.mean()\n",
    "    return np.exp(x_hat)\n",
    "\n",
    "def possion_estimator(delta):\n",
    "    J = np.random.poisson(delta)\n",
    "    x = np.random.normal(0, 1, J)\n",
    "    return np.prod(x / delta) * np.exp(delta)\n",
    "\n",
    "def k_trial(k, method, N=10):\n",
    "    return np.array([method(N) for _ in range(k)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_mean = []\n",
    "res_std = []\n",
    "for k in [10, 100, 1000, 10000, 100000, 1000000]:\n",
    "    res = k_trial(k, naive_method)\n",
    "    res_mean.append(res.mean())\n",
    "    res_std.append(res.std()/np.sqrt(k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7f1074d0a910>"
      ]
     },
     "execution_count": 186,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGiCAYAAADEJZ3cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqJElEQVR4nO3df1RU953/8deAMGhWJuIPfqyANDUqMRKCP/gRbE1TkEZWtz0re7oSdRNTW7eRspvTZU1MzI9S28Yk/iI1jRJOEjSpRc2pJuJ3q8RIYiDiNtpjtMEDy86UaiIj+hUV7/cPv85m5EcYFObj+Hyc8zkn9zOfe3nfe047Lz/3zufaLMuyBAAAYLAgfxcAAADwVQgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4PgeWqqoq5ebmKiYmRjabTVu2bOl2/O9+9zt9+9vf1vDhwxUeHq60tDS9++67HcZt3rxZiYmJstvtSkxMVEVFha+lAQCAAOVzYDlz5oySkpK0evXqHo2vqqrSt7/9bW3fvl21tbWaNm2acnNzdeDAAc+Y6upq5eXlKT8/XwcPHlR+fr5mz56tDz/80NfyAABAALJdy8sPbTabKioqNGvWLJ/2u+OOO5SXl6elS5dKkvLy8uR2u7Vjxw7PmOnTp2vIkCEqLy/vbXkAACBADOjvP3jp0iWdPn1aERERnr7q6mr95Cc/8RqXnZ2tF154ocvjtLW1qa2tzeu4n3/+uYYOHSqbzXbd6wYAANefZVk6ffq0YmJiFBTU9Y2ffg8szz33nM6cOaPZs2d7+lwulyIjI73GRUZGyuVydXmc4uJiLVu2rM/qBAAA/aexsVEjR47s8vN+DSzl5eV68skntXXrVo0YMcLrs6tnRSzL6nampKioSIWFhZ7tlpYWxcXFqbGxUeHh4del3rPnL2rys/9HkrR/ybc0KLTf8x0AAAHN7XYrNjZWgwcP7nZcv30Db9q0SQ8++KDeeust3XfffV6fRUVFdZhNaW5u7jDr8mV2u112u71Df3h4+HULLAPOX1SQfZDnuAQWAAD6xlc9ztEv67CUl5dr3rx5euONN3T//fd3+DwtLU2VlZVefTt37lR6enp/lAcAAAzn85RBa2urjh075tmur69XXV2dIiIiFBcXp6KiIjU1NamsrEzS5bDywAMP6MUXX1RqaqpnJmXgwIFyOBySpMWLF2vq1Klavny5Zs6cqa1bt2rXrl3au3fv9ThHAABwg/N5hqWmpkbJyclKTk6WJBUWFio5OdnzE2Wn06mGhgbP+F//+te6ePGiFi1apOjoaE9bvHixZ0x6ero2btyoDRs2aMKECSotLdWmTZs0ZcqUaz0/AAAQAK5pHRaTuN1uORwOtbS0XNeHbhOXXl6V9/BT2TzDAgDAddbT72/eJQQAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4/kcWKqqqpSbm6uYmBjZbDZt2bKl2/FOp1Pf//73NWbMGAUFBamgoKDDmNLSUtlstg7t3LlzvpYHAAACkM+B5cyZM0pKStLq1at7NL6trU3Dhw/XkiVLlJSU1OW48PBwOZ1OrxYWFuZreQAAIAAN8HWHnJwc5eTk9Hj8qFGj9OKLL0qS1q9f3+U4m82mqKgoX8sBAAA3AWOeYWltbVV8fLxGjhypGTNm6MCBA92Ob2trk9vt9moAACAwGRFYxo4dq9LSUm3btk3l5eUKCwtTRkaGjh492uU+xcXFcjgcnhYbG9uPFQMAgP5kRGBJTU3VnDlzlJSUpMzMTL355pu6/fbbtWrVqi73KSoqUktLi6c1Njb2Y8UAAKA/+fwMS38ICgrSpEmTup1hsdvtstvt/VgVAADwFyNmWK5mWZbq6uoUHR3t71IAAIABfJ5haW1t1bFjxzzb9fX1qqurU0REhOLi4lRUVKSmpiaVlZV5xtTV1Xn2/etf/6q6ujqFhoYqMTFRkrRs2TKlpqZq9OjRcrvdWrlyperq6rRmzZprPD0AABAIfA4sNTU1mjZtmme7sLBQkjR37lyVlpbK6XSqoaHBa5/k5GTPf9fW1uqNN95QfHy8jh8/Lkk6deqUHn74YblcLjkcDiUnJ6uqqkqTJ0/uzTkBAIAAY7Msy/J3EdeD2+2Ww+FQS0uLwsPDr8sxz56/qMSl70qSDj+VrUGhRj7yAwDADaun399GPsMCAADwZQQWAABgPAILAAAwHoGlG+2X/vfxnv31n3ttAwCA/kNg6cI7nzh134o9nu15Gz7SPcv/U+984vRjVQAA3JwILJ145xOnfvjax/qLu82r39VyTj987WNCCwAA/YzAcpX2S5aWvX1Ynd38udK37O3D3B4CAKAfEViusr/+czlbznX5uSXJ2XJO++s/77+iAAC4yRFYrtJ8uuuw0ptxAADg2hFYrjJicNh1HQcAAK4dgeUqkxMiFO0Ik62Lz22Soh1hmpwQ0Z9lAQBwUyOwXCU4yKYnci+/Rfrq0HJl+4ncRAUHdRVpAADA9UZg6cT08dEqmXO3RoTbvfqjHGEqmXO3po+P9lNlAADcnHj9cBemj49WxteH6c4nd0qSSudPUubo4cysAADgB8ywdOPL4WRyQgRhBQAAPyGwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDyfA0tVVZVyc3MVExMjm82mLVu2dDve6XTq+9//vsaMGaOgoCAVFBR0Om7z5s1KTEyU3W5XYmKiKioqfC0NAAAEKJ8Dy5kzZ5SUlKTVq1f3aHxbW5uGDx+uJUuWKCkpqdMx1dXVysvLU35+vg4ePKj8/HzNnj1bH374oa/lAQCAAGSzLMvq9c42myoqKjRr1qwejf/mN7+pu+66Sy+88IJXf15entxut3bs2OHpmz59uoYMGaLy8vIeHdvtdsvhcKilpUXh4eE9PYVunT1/UYlL35UkHX4qW4NCB1yX4wIAgMt6+v1txDMs1dXVysrK8urLzs7Wvn37utynra1NbrfbqwEAgMBkRGBxuVyKjIz06ouMjJTL5epyn+LiYjkcDk+LjY3t6zIBAICfGBFYpMu3l77MsqwOfV9WVFSklpYWT2tsbOzrEgEAgJ8Y8VBGVFRUh9mU5ubmDrMuX2a322W32/u6NAAAYAAjZljS0tJUWVnp1bdz506lp6f7qSIAAGASn2dYWltbdezYMc92fX296urqFBERobi4OBUVFampqUllZWWeMXV1dZ59//rXv6qurk6hoaFKTEyUJC1evFhTp07V8uXLNXPmTG3dulW7du3S3r17r/H0AABAIPA5sNTU1GjatGme7cLCQknS3LlzVVpaKqfTqYaGBq99kpOTPf9dW1urN954Q/Hx8Tp+/LgkKT09XRs3btRjjz2mxx9/XLfddps2bdqkKVOm9OacAABAgLmmdVhMwjosAADceG6odVgAAAC6Q2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIw3wN8FAADQG+2XLO2v/1zNp89pxOAwTU6IUHCQzd9lBRxTrjOBBQBww3nnE6eWvX1YzpZznr5oR5ieyE3U9PHRfqwssJh0nbklBAC4obzziVM/fO1jry9RSXK1nNMPX/tY73zi9FNlgcW06+xzYKmqqlJubq5iYmJks9m0ZcuWr9xnz549SklJUVhYmL72ta/ppZde8vq8tLRUNputQzt37lwXRwQAc7VfslT955PaWtek6j+fVPsly98lBYz2S5aWvX1YnV3RK33L3j7MNb9GJl5nn28JnTlzRklJSZo/f76+973vfeX4+vp6fec739GCBQv02muv6f3339ePfvQjDR8+3Gv/8PBwHTlyxGvfsLAwX8sDAL8yaQo9EO2v/7zDv/i/zJLkbDmn/fWfK+22of1XWIAx8Tr7HFhycnKUk5PT4/EvvfSS4uLi9MILL0iSxo0bp5qaGv3qV7/yCiw2m01RUVE9Pm5bW5va2to82263u8f7AkBfuDKFfvW/Oa9MoZfMuZvQco2aT/ds5r2n49A5E69znz/DUl1draysLK++7Oxs1dTU6MKFC56+1tZWxcfHa+TIkZoxY4YOHDjQ7XGLi4vlcDg8LTY2tk/qB4CeMHEKPRCNGNyzmfeejkPnTLzOfR5YXC6XIiMjvfoiIyN18eJFnThxQpI0duxYlZaWatu2bSovL1dYWJgyMjJ09OjRLo9bVFSklpYWT2tsbOzT80DfOXv+okb9++816t9/r7PnL/q7HKBXfJlCR+9NTohQtCNMXf2o1qbLt+AmJ0T0Z1kBx8Tr3C+/ErLZvE/Zsiyv/tTUVM2ZM0dJSUnKzMzUm2++qdtvv12rVq3q8ph2u13h4eFeDQD8xcQp9EAUHGTTE7mJktThy/TK9hO5iazHco1MvM59HliioqLkcrm8+pqbmzVgwAANHdr5gzpBQUGaNGlStzMsAGASE6fQA9X08dEqmXO3ohze1zLKEcZzQteRade5zxeOS0tL09tvv+3Vt3PnTk2cOFEhISGd7mNZlurq6nTnnXf2dXkAcF1cmUJ3tZzr9DkWmy7/Hz23Kq6P6eOj9e3EKCNWYA1kJl1nnwNLa2urjh075tmur69XXV2dIiIiFBcXp6KiIjU1NamsrEyStHDhQq1evVqFhYVasGCBqqur9corr6i8vNxzjGXLlik1NVWjR4+W2+3WypUrVVdXpzVr1lyHUwSAvndlCv2Hr30sm+QVWrhV0TeCg2z8dLkfmHKdfb4lVFNTo+TkZCUnJ0uSCgsLlZycrKVLl0qSnE6nGhoaPOMTEhK0fft27d69W3fddZeefvpprVy50usnzadOndLDDz+scePGKSsrS01NTaqqqtLkyZOv9fwAfAkPOPct06bQgUDi8wzLN7/5Tc9Ds50pLS3t0PeNb3xDH3/8cZf7PP/883r++ed9LQUAjGPSFDoQSHj5IQBcZ6ZMoQOBhJcfAgAA4xFYgJvIl1dZ3V//OauuArhhEFiAm8Q7nzh134o9nu15Gz7SPcv/s99fEQ8AvUFgAW4CV17K9xd3m1f/lZfyEVoAmI7AAgQ4XsoHIBAQWIAAx0v5AAQCAgsQ4HgpH4BAQGABAhwv5QMQCAgsQIC78lK+rtZZtUmK5qV8AAxHYAEC3JWX8knqEFp4KR+AGwWBBbgJXHkp34hwu1c/L+UDcKPgXULwu6tXX80cPZx/7feB6eOjlfH1YbrzyZ2SpNL5k7jWAG4YzLDAr1h9tX99OZzwBmEANxICC/yG1VcBAD1FYIFfsPoqAMAXPMMCv/Bl9dW024b2X2EBblDoAB3/+f3+LgMAfMYMC/yC1VcBAL4gsMAvWH0VAOALAgv8gtVXAQC+ILDAL1h9FQDgCwIL/IbVVwEAPcWvhOBXrL4KAOgJZljgd6y+CgD4KgQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABjP58BSVVWl3NxcxcTEyGazacuWLV+5z549e5SSkqKwsDB97Wtf00svvdRhzObNm5WYmCi73a7ExERVVFT4WhoAAAhQPr+t+cyZM0pKStL8+fP1ve997yvH19fX6zvf+Y4WLFig1157Te+//75+9KMfafjw4Z79q6urlZeXp6efflp///d/r4qKCs2ePVt79+7VlClTfC1QCg729bQ6d/6iBp4/97/HvcDLrfsE1xkAbl5nzvRomM2yLKu3f8Nms6miokKzZs3qcsxPf/pTbdu2TX/60588fQsXLtTBgwdVXV0tScrLy5Pb7daOHTs8Y6ZPn64hQ4aovLy80+O2tbWpra3Ns+12uxUbG6sWSeG9PSEAANCv3JIcklpaWhQe3vU3eJ8/w1JdXa2srCyvvuzsbNXU1OjChQvdjtm3b1+Xxy0uLpbD4fC02NjY6188AAAwQp/PvbtcLkVGRnr1RUZG6uLFizpx4oSio6O7HONyubo8blFRkQoLCz3bV2ZY9D//I3WT0Hxx9vxFpTy9S5JU+/h9GhTKrYq+wHUGgJuY2y3FxHzlsH75ZrDZbF7bV+5Cfbm/szFX932Z3W6X3W7v+MEtt1xu10PIRf3f0LD/PS5fpH2D6wwAN6/29h4N6/NvhqioqA4zJc3NzRowYICGDh3a7ZirZ10QmAaFDtDxn9/v7zIAAAbr82dY0tLSVFlZ6dW3c+dOTZw4USEhId2OSU9P7+vyAADADcDnGZbW1lYdO3bMs11fX6+6ujpFREQoLi5ORUVFampqUllZmaTLvwhavXq1CgsLtWDBAlVXV+uVV17x+vXP4sWLNXXqVC1fvlwzZ87U1q1btWvXLu3du/c6nCIAALjR+TzDUlNTo+TkZCUnJ0uSCgsLlZycrKVLl0qSnE6nGhoaPOMTEhK0fft27d69W3fddZeefvpprVy50msNl/T0dG3cuFEbNmzQhAkTVFpaqk2bNvm+BgsAAAhI17QOi0ncbrccDsdX/o7bF2fPX1Ti0nclSYefyubXKwAAXGc9/f7mXUIAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgs3doN3iIMAIAZmGEBAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADG61VgWbt2rRISEhQWFqaUlBS999573Y5fs2aNxo0bp4EDB2rMmDEqKyvz+ry0tFQ2m61DO3fuXG/KAwAAAWaArzts2rRJBQUFWrt2rTIyMvTrX/9aOTk5Onz4sOLi4jqMLykpUVFRkV5++WVNmjRJ+/fv14IFCzRkyBDl5uZ6xoWHh+vIkSNe+4aFhfXilAAAQKCxWZZl+bLDlClTdPfdd6ukpMTTN27cOM2aNUvFxcUdxqenpysjI0O//OUvPX0FBQWqqanR3r17JV2eYSkoKNCpU6d6eRqS2+2Ww+FQS0uLwsPDe30cAADQf3r6/e3TLaHz58+rtrZWWVlZXv1ZWVnat29fp/u0tbV1mCkZOHCg9u/frwsXLnj6WltbFR8fr5EjR2rGjBk6cOBAt7W0tbXJ7XZ7NQAAEJh8CiwnTpxQe3u7IiMjvfojIyPlcrk63Sc7O1u/+c1vVFtbK8uyVFNTo/Xr1+vChQs6ceKEJGns2LEqLS3Vtm3bVF5errCwMGVkZOjo0aNd1lJcXCyHw+FpsbGxvpwKAAC4gfTqoVubzea1bVlWh74rHn/8ceXk5Cg1NVUhISGaOXOm5s2bJ0kKDg6WJKWmpmrOnDlKSkpSZmam3nzzTd1+++1atWpVlzUUFRWppaXF0xobG3tzKgAA4AbgU2AZNmyYgoODO8ymNDc3d5h1uWLgwIFav369zp49q+PHj6uhoUGjRo3S4MGDNWzYsM6LCgrSpEmTup1hsdvtCg8P92oAACAw+RRYQkNDlZKSosrKSq/+yspKpaend7tvSEiIRo4cqeDgYG3cuFEzZsxQUFDnf96yLNXV1Sk6OtqX8gAAQIDy+WfNhYWFys/P18SJE5WWlqZ169apoaFBCxculHT5Vk1TU5NnrZVPP/1U+/fv15QpU/TFF19oxYoV+uSTT/Tqq696jrls2TKlpqZq9OjRcrvdWrlyperq6rRmzZrrdJoAAOBG5nNgycvL08mTJ/XUU0/J6XRq/Pjx2r59u+Lj4yVJTqdTDQ0NnvHt7e167rnndOTIEYWEhGjatGnat2+fRo0a5Rlz6tQpPfzww3K5XHI4HEpOTlZVVZUmT5587WcIAABueD6vw2Iq1mEBAODG0yfrsAAAAPgDgQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjNerwLJ27VolJCQoLCxMKSkpeu+997odv2bNGo0bN04DBw7UmDFjVFZW1mHM5s2blZiYKLvdrsTERFVUVPSmNAAAEIB8DiybNm1SQUGBlixZogMHDigzM1M5OTlqaGjodHxJSYmKior05JNP6tChQ1q2bJkWLVqkt99+2zOmurpaeXl5ys/P18GDB5Wfn6/Zs2frww8/7P2ZAQCAgGGzLMvyZYcpU6bo7rvvVklJiadv3LhxmjVrloqLizuMT09PV0ZGhn75y196+goKClRTU6O9e/dKkvLy8uR2u7Vjxw7PmOnTp2vIkCEqLy/vUV1ut1sOh0MtLS0KDw/35ZQAAICf9PT726cZlvPnz6u2tlZZWVle/VlZWdq3b1+n+7S1tSksLMyrb+DAgdq/f78uXLgg6fIMy9XHzM7O7vKYV47rdru9GgAACEw+BZYTJ06ovb1dkZGRXv2RkZFyuVyd7pOdna3f/OY3qq2tlWVZqqmp0fr163XhwgWdOHFCkuRyuXw6piQVFxfL4XB4WmxsrC+nAgAAbiC9eujWZrN5bVuW1aHviscff1w5OTlKTU1VSEiIZs6cqXnz5kmSgoODe3VMSSoqKlJLS4unNTY29uZUAADADcCnwDJs2DAFBwd3mPlobm7uMENyxcCBA7V+/XqdPXtWx48fV0NDg0aNGqXBgwdr2LBhkqSoqCifjilJdrtd4eHhXg0AAAQmnwJLaGioUlJSVFlZ6dVfWVmp9PT0bvcNCQnRyJEjFRwcrI0bN2rGjBkKCrr859PS0jocc+fOnV95TAAAcHMY4OsOhYWFys/P18SJE5WWlqZ169apoaFBCxculHT5Vk1TU5NnrZVPP/1U+/fv15QpU/TFF19oxYoV+uSTT/Tqq696jrl48WJNnTpVy5cv18yZM7V161bt2rXL8ysiAABwc/M5sOTl5enkyZN66qmn5HQ6NX78eG3fvl3x8fGSJKfT6bUmS3t7u5577jkdOXJEISEhmjZtmvbt26dRo0Z5xqSnp2vjxo167LHH9Pjjj+u2227Tpk2bNGXKlGs/QwAAcMPzeR0WU7EOCwAAN54+WYcFAADAHwgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGC8XgWWtWvXKiEhQWFhYUpJSdF7773X7fjXX39dSUlJGjRokKKjozV//nydPHnS83lpaalsNluHdu7cud6UBwAAAozPgWXTpk0qKCjQkiVLdODAAWVmZionJ0cNDQ2djt+7d68eeOABPfjggzp06JDeeustffTRR3rooYe8xoWHh8vpdHq1sLCw3p0VAAAIKD4HlhUrVujBBx/UQw89pHHjxumFF15QbGysSkpKOh3/wQcfaNSoUXrkkUeUkJCge+65Rz/4wQ9UU1PjNc5msykqKsqrAQAASD4GlvPnz6u2tlZZWVle/VlZWdq3b1+n+6Snp+u///u/tX37dlmWpb/85S/67W9/q/vvv99rXGtrq+Lj4zVy5EjNmDFDBw4c6LaWtrY2ud1urwYAAAKTT4HlxIkTam9vV2RkpFd/ZGSkXC5Xp/ukp6fr9ddfV15enkJDQxUVFaVbb71Vq1at8owZO3asSktLtW3bNpWXlyssLEwZGRk6evRol7UUFxfL4XB4WmxsrC+nAgAAbiC9eujWZrN5bVuW1aHvisOHD+uRRx7R0qVLVVtbq3feeUf19fVauHChZ0xqaqrmzJmjpKQkZWZm6s0339Ttt9/uFWquVlRUpJaWFk9rbGzszakAAIAbwABfBg8bNkzBwcEdZlOam5s7zLpcUVxcrIyMDD366KOSpAkTJuiWW25RZmamnnnmGUVHR3fYJygoSJMmTep2hsVut8tut/tSPgAAuEH5NMMSGhqqlJQUVVZWevVXVlYqPT29033Onj2roCDvPxMcHCzp8sxMZyzLUl1dXadhBgAA3Hx8mmGRpMLCQuXn52vixIlKS0vTunXr1NDQ4LnFU1RUpKamJpWVlUmScnNztWDBApWUlCg7O1tOp1MFBQWaPHmyYmJiJEnLli1TamqqRo8eLbfbrZUrV6qurk5r1qy5jqcKAABuVD4Hlry8PJ08eVJPPfWUnE6nxo8fr+3btys+Pl6S5HQ6vdZkmTdvnk6fPq3Vq1frX//1X3Xrrbfq3nvv1fLlyz1jTp06pYcfflgul0sOh0PJycmqqqrS5MmTr8MpAgCAG53N6uq+zA3G7XbL4XCopaVF4eHh/i4HAAD0QE+/v3mXEAAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADG61VgWbt2rRISEhQWFqaUlBS999573Y5//fXXlZSUpEGDBik6Olrz58/XyZMnvcZs3rxZiYmJstvtSkxMVEVFRW9KAwAAAcjnwLJp0yYVFBRoyZIlOnDggDIzM5WTk6OGhoZOx+/du1cPPPCAHnzwQR06dEhvvfWWPvroIz300EOeMdXV1crLy1N+fr4OHjyo/Px8zZ49Wx9++GHvzwwAAAQMm2VZli87TJkyRXfffbdKSko8fePGjdOsWbNUXFzcYfyvfvUrlZSU6M9//rOnb9WqVfrFL36hxsZGSVJeXp7cbrd27NjhGTN9+nQNGTJE5eXlndbR1tamtrY2z3ZLS4vi4uLU2Nio8PBwX04JAAD4idvtVmxsrE6dOiWHw9H1QMsHbW1tVnBwsPW73/3Oq/+RRx6xpk6d2uk+77//vhUaGmr9/ve/ty5dumS5XC5r6tSp1g9+8APPmNjYWGvFihVe+61YscKKi4vrspYnnnjCkkSj0Wg0Gi0AWmNjY7cZZIB8cOLECbW3tysyMtKrPzIyUi6Xq9N90tPT9frrrysvL0/nzp3TxYsX9Xd/93datWqVZ4zL5fLpmJJUVFSkwsJCz/alS5f0+eefa+jQobLZbL6cVreuJD9mbvoW17n/cK37B9e5f3Cd+0dfXmfLsnT69GnFxMR0O86nwHLF1YHAsqwuQ8Lhw4f1yCOPaOnSpcrOzpbT6dSjjz6qhQsX6pVXXunVMSXJbrfLbrd79d16660+nknPhYeH8z+GfsB17j9c6/7Bde4fXOf+0VfXudtbQf+fT4Fl2LBhCg4O7jDz0dzc3GGG5Iri4mJlZGTo0UcflSRNmDBBt9xyizIzM/XMM88oOjpaUVFRPh0TAADcXHz6lVBoaKhSUlJUWVnp1V9ZWan09PRO9zl79qyCgrz/THBwsKTLsyiSlJaW1uGYO3fu7PKYAADg5uLzLaHCwkLl5+dr4sSJSktL07p169TQ0KCFCxdKuvxsSVNTk8rKyiRJubm5WrBggUpKSjy3hAoKCjR58mTP/arFixdr6tSpWr58uWbOnKmtW7dq165d2rt373U81d6x2+164oknOtx+wvXFde4/XOv+wXXuH1zn/mHCdfb5Z83S5YXjfvGLX8jpdGr8+PF6/vnnNXXqVEnSvHnzdPz4ce3evdszftWqVXrppZdUX1+vW2+9Vffee6+WL1+uv/3bv/WM+e1vf6vHHntMn332mW677TY9++yz+u53v3vtZwgAAG54vQosAAAA/Yl3CQEAAOMRWAAAgPEILAAAwHgEFgAAYDwCSxeqqqqUm5urmJgY2Ww2bdmyxd8lBaTi4mJNmjRJgwcP1ogRIzRr1iwdOXLE32UFnJKSEk2YMMGzSmVaWprXy0bRN4qLi2Wz2VRQUODvUgLOk08+KZvN5tWioqL8XVZAampq0pw5czR06FANGjRId911l2pra/u9DgJLF86cOaOkpCStXr3a36UEtD179mjRokX64IMPVFlZqYsXLyorK0tnzpzxd2kBZeTIkfr5z3+umpoa1dTU6N5779XMmTN16NAhf5cWsD766COtW7dOEyZM8HcpAeuOO+6Q0+n0tD/+8Y/+LingfPHFF8rIyFBISIh27Nihw4cP67nnnuvTV+F0pVfvEroZ5OTkKCcnx99lBLx33nnHa3vDhg0aMWKEamtrPWv74Nrl5uZ6bT/77LMqKSnRBx98oDvuuMNPVQWu1tZW/dM//ZNefvllPfPMM/4uJ2ANGDCAWZU+tnz5csXGxmrDhg2evlGjRvmlFmZYYJSWlhZJUkREhJ8rCVzt7e3auHGjzpw5o7S0NH+XE5AWLVqk+++/X/fdd5+/SwloR48eVUxMjBISEvSP//iP+uyzz/xdUsDZtm2bJk6cqH/4h3/QiBEjlJycrJdfftkvtRBYYAzLslRYWKh77rlH48eP93c5AeePf/yj/uZv/kZ2u10LFy5URUWFEhMT/V1WwNm4caM+/vhjFRcX+7uUgDZlyhSVlZXp3Xff1csvvyyXy6X09HSdPHnS36UFlM8++0wlJSUaPXq03n33XS1cuFCPPPKI5/U7/YlbQjDGv/zLv+i//uu/jHiHVCAaM2aM6urqdOrUKW3evFlz587Vnj17CC3XUWNjoxYvXqydO3cqLCzM3+UEtC/fsr/zzjuVlpam2267Ta+++qoKCwv9WFlguXTpkiZOnKif/exnkqTk5GQdOnRIJSUleuCBB/q1FmZYYIQf//jH2rZtm/7whz9o5MiR/i4nIIWGhurrX/+6Jk6cqOLiYiUlJenFF1/0d1kBpba2Vs3NzUpJSdGAAQM0YMAA7dmzRytXrtSAAQPU3t7u7xID1i233KI777xTR48e9XcpASU6OrrDP2rGjRunhoaGfq+FGRb4lWVZ+vGPf6yKigrt3r1bCQkJ/i7ppmFZltra2vxdRkD51re+1eGXKvPnz9fYsWP105/+VMHBwX6qLPC1tbXpT3/6kzIzM/1dSkDJyMjosNTEp59+qvj4+H6vhcDShdbWVh07dsyzXV9fr7q6OkVERCguLs6PlQWWRYsW6Y033tDWrVs1ePBguVwuSZLD4dDAgQP9XF3g+I//+A/l5OQoNjZWp0+f1saNG7V79+4Ov9LCtRk8eHCH569uueUWDR06lOeyrrN/+7d/U25uruLi4tTc3KxnnnlGbrdbc+fO9XdpAeUnP/mJ0tPT9bOf/UyzZ8/W/v37tW7dOq1bt67/i7HQqT/84Q+WpA5t7ty5/i4toHR2jSVZGzZs8HdpAeWf//mfrfj4eCs0NNQaPny49a1vfcvauXOnv8u6KXzjG9+wFi9e7O8yAk5eXp4VHR1thYSEWDExMdZ3v/td69ChQ/4uKyC9/fbb1vjx4y273W6NHTvWWrdunV/qsFmWZfV/TAIAAOg5HroFAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPH+H+pz2SA8WU0IAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# lower variance but biased\n",
    "plt.errorbar(np.log10([10, 100, 1000, 10000, 100000, 1000000]), res_mean, yerr=res_std, fmt='o')\n",
    "plt.ylim(0.8, 1.2)\n",
    "plt.axhline(1, color='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_mean = []\n",
    "res_std = []\n",
    "for k in [10, 100, 1000, 10000, 100000, 1000000]:\n",
    "    res = k_trial(k, possion_estimator, N=1)\n",
    "    res_mean.append(res.mean())\n",
    "    res_std.append(res.std()/np.sqrt(k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7f107482ca00>"
      ]
     },
     "execution_count": 191,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGiCAYAAADEJZ3cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqFUlEQVR4nO3de3RV9Z338c/J7QQsOXLNZUhCpApEJBPDJZeGihdCKhmYliHzdAjgIEqHKWCmrj4pomK1KZ2KyC0WK0ZaDWhpAFdBCPNUIiWKiRymgktxDJMMc05TUHIIPAQI+/mDhzOe5iInJDk/Tt6vtfZa7t/57p3v3iuaj799OTbLsiwBAAAYLCTQDQAAAHwVAgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMJ7fgaWyslJ5eXmKi4uTzWbTtm3bOqz/7W9/q/vuu0+DBw9WVFSUMjIytHv37lZ1W7duVXJysux2u5KTk1VeXu5vawAAIEj5HVjOnj2rlJQUrV279prqKysrdd9992nnzp2qqanRpEmTlJeXp0OHDnlrqqqqlJ+fr4KCAh0+fFgFBQWaOXOm3nvvPX/bAwAAQch2PV9+aLPZVF5erunTp/u13e233678/Hw9/vjjkqT8/Hx5PB7t2rXLWzNlyhT1799fZWVlnW0PAAAEibCe/oGXL1/WmTNnNGDAAO9YVVWVHnnkEZ+6nJwcrVq1qt39NDc3q7m52We/n3/+uQYOHCibzdblfQMAgK5nWZbOnDmjuLg4hYS0f+GnxwPLs88+q7Nnz2rmzJneMbfbrejoaJ+66Ohoud3udvdTXFys5cuXd1ufAACg59TX12vo0KHtft6jgaWsrExPPvmktm/friFDhvh89pezIpZldThTUlRUpMLCQu96Y2OjEhISVF9fr6ioqC7p99yFSxr/zL9Jkg4uvUd9I3o83/UKnGcA6L08Ho/i4+PVr1+/Dut67C/Dli1bNG/ePL3xxhu69957fT6LiYlpNZvS0NDQatbly+x2u+x2e6vxqKioLgssYRcuKcTe17tf/pB2D84zAOCrbufokfewlJWVae7cuXrttdd0//33t/o8IyNDFRUVPmN79uxRZmZmT7QHAAAM5/f/yjY1NenTTz/1rtfW1srpdGrAgAFKSEhQUVGRTpw4oU2bNkm6ElZmz56t559/Xunp6d6ZlD59+sjhcEiSFi9erIkTJ2rFihWaNm2atm/frr1792r//v1dcYwAAOAG5/cMS3V1tVJTU5WamipJKiwsVGpqqvcRZZfLpbq6Om/9L37xC126dEkLFy5UbGysd1m8eLG3JjMzU5s3b9bLL7+sMWPGqLS0VFu2bNGECROu9/gAAEAQ8HuG5a677lJHr24pLS31WX/77bevab8zZszQjBkz/G0HAAD0AnyXEAAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACM53dgqaysVF5enuLi4mSz2bRt27YO610ul7773e9qxIgRCgkJ0ZIlS1rVlJaWymaztVrOnz/vb3sAACAI+R1Yzp49q5SUFK1du/aa6pubmzV48GAtXbpUKSkp7dZFRUXJ5XL5LJGRkf62BwAAglCYvxvk5uYqNzf3muuHDRum559/XpK0cePGdutsNptiYmL8bQcAAPQCxtzD0tTUpMTERA0dOlRTp07VoUOHOqxvbm6Wx+PxWQAAQHAyIrCMHDlSpaWl2rFjh8rKyhQZGamsrCwdO3as3W2Ki4vlcDi8S3x8fA92DAAAepIRgSU9PV2zZs1SSkqKsrOz9frrr+u2227TmjVr2t2mqKhIjY2N3qW+vr4HOwYAAD3J73tYekJISIjGjRvX4QyL3W6X3W7vwa4AAECgGDHD8pcsy5LT6VRsbGygWwEAAAbwe4alqalJn376qXe9trZWTqdTAwYMUEJCgoqKinTixAlt2rTJW+N0Or3b/vnPf5bT6VRERISSk5MlScuXL1d6erpuvfVWeTwerV69Wk6nU+vWrbvOwwMAAMHA78BSXV2tSZMmedcLCwslSXPmzFFpaalcLpfq6up8tklNTfX+c01NjV577TUlJibq+PHjkqTTp0/roYcektvtlsPhUGpqqiorKzV+/PjOHBMAAAgyfgeWu+66S5Zltft5aWlpq7GO6iXpueee03PPPedvKwAAoJcw8h4WAACALyOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDx/A4slZWVysvLU1xcnGw2m7Zt29Zhvcvl0ne/+12NGDFCISEhWrJkSZt1W7duVXJysux2u5KTk1VeXu5vawAAIEj5HVjOnj2rlJQUrV279prqm5ubNXjwYC1dulQpKSlt1lRVVSk/P18FBQU6fPiwCgoKNHPmTL333nv+tgcAAIJQmL8b5ObmKjc395rrhw0bpueff16StHHjxjZrVq1apfvuu09FRUWSpKKiIu3bt0+rVq1SWVmZvy0CAIAgY8Q9LFVVVZo8ebLPWE5Ojg4cONDuNs3NzfJ4PD4LAAAITkYEFrfbrejoaJ+x6Ohoud3udrcpLi6Ww+HwLvHx8d3dJgAACBAjAosk2Ww2n3XLslqNfVlRUZEaGxu9S319fXe3CAAAAsTve1i6Q0xMTKvZlIaGhlazLl9mt9tlt9u7uzUAAGAAI2ZYMjIyVFFR4TO2Z88eZWZmBqgjAABgEr9nWJqamvTpp59612tra+V0OjVgwAAlJCSoqKhIJ06c0KZNm7w1TqfTu+2f//xnOZ1ORUREKDk5WZK0ePFiTZw4UStWrNC0adO0fft27d27V/v377/OwwMAAMHA78BSXV2tSZMmedcLCwslSXPmzFFpaalcLpfq6up8tklNTfX+c01NjV577TUlJibq+PHjkqTMzExt3rxZjz32mJYtW6bhw4dry5YtmjBhQmeOCQAABBm/A8tdd90ly7La/by0tLTVWEf1V82YMUMzZszwtx0AANALGHHTLXq3vhFhOv7T+wPdBgDAYEbcdAsAANARAgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGM/vwFJZWam8vDzFxcXJZrNp27ZtX7nNvn37lJaWpsjISN1yyy164YUXfD4vLS2VzWZrtZw/f97f9gAAQBDyO7CcPXtWKSkpWrt27TXV19bW6lvf+pays7N16NAh/ehHP9KiRYu0detWn7qoqCi5XC6fJTIy0t/2AABAEArzd4Pc3Fzl5uZec/0LL7yghIQErVq1SpI0atQoVVdX6+c//7m+853veOtsNptiYmKueb/Nzc1qbm72rns8nmveFgAA3Fi6/R6WqqoqTZ482WcsJydH1dXVunjxonesqalJiYmJGjp0qKZOnapDhw51uN/i4mI5HA7vEh8f3y39AwCAwOv2wOJ2uxUdHe0zFh0drUuXLunkyZOSpJEjR6q0tFQ7duxQWVmZIiMjlZWVpWPHjrW736KiIjU2NnqX+vr6bj0OAAAQOH5fEuoMm83ms25Zls94enq60tPTvZ9nZWXpzjvv1Jo1a7R69eo292m322W327upYwAAYJJun2GJiYmR2+32GWtoaFBYWJgGDhzYdlMhIRo3blyHMywAAKD36PbAkpGRoYqKCp+xPXv2aOzYsQoPD29zG8uy5HQ6FRsb293tAQCAG4DfgaWpqUlOp1NOp1PSlceWnU6n6urqJF25t2T27Nne+gULFug///M/VVhYqI8++kgbN27USy+9pB/84AfemuXLl2v37t367LPP5HQ6NW/ePDmdTi1YsOA6Dw8AAAQDv+9hqa6u1qRJk7zrhYWFkqQ5c+aotLRULpfLG14kKSkpSTt37tQjjzyidevWKS4uTqtXr/Z5pPn06dN66KGH5Ha75XA4lJqaqsrKSo0fP/56jg0AAAQJm3X1DtgbnMfjkcPhUGNjo6Kiorpkn+cuXFLy47slSUefylHfiB65RxkAgF7jWv9+811CAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAYAudu7CJQ3737/TsP/9O527cCnQ7QBBgcACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBagF+HpFQA3KgILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjhQW6AQAAYK6Wy5YO1n6uhjPnNaRfpMYnDVBoiK3H+yCwAACANr31oUvL3zwqV+N571isI1JP5CVryujYHu3F70tClZWVysvLU1xcnGw2m7Zt2/aV2+zbt09paWmKjIzULbfcohdeeKFVzdatW5WcnCy73a7k5GSVl5f72xoAAOgib33o0vd+/YFPWJEkd+N5fe/XH+itD1092o/fMyxnz55VSkqKHnjgAX3nO9/5yvra2lp961vf0vz58/XrX/9af/jDH/RP//RPGjx4sHf7qqoq5efn68c//rH+9m//VuXl5Zo5c6b279+vCRMm+NugFBrq72G17cIl9blw/n/2e5EJKdzg+J3uGZxn3OBaLltavv2IrDY+syTZJC3fcUT3Jfa7/stDZ89eU5nNsqy2+rm2jW02lZeXa/r06e3W/PCHP9SOHTv00UcfeccWLFigw4cPq6qqSpKUn58vj8ejXbt2eWumTJmi/v37q6ysrM39Njc3q7m52bvu8XgUHx+vRklRnT0gAACgqvg79L++W/yVdWWvFSmj/o/X9bM8khySGhsbFRXV/l/wbn9KqKqqSpMnT/YZy8nJUXV1tS5evNhhzYEDB9rdb3FxsRwOh3eJj4/v+uYBAOiFGr7Wv0vrukK3z1O63W5FR0f7jEVHR+vSpUs6efKkYmNj261xu93t7reoqEiFhYXe9aszLPrv/5Y6SGj+OHfhktJ+vFeSVLPsXvWNYFoXNzZ+p3vGmfMXlfaT38uySS8W3KnsWwYG5KmKYMd57j5Djn8h/erwV9e98ao07DpDi8cjxcV9ZVmP/NfKZvP9Bbp6FerL423V/OXYl9ntdtnt9tYf3HTTlaUrhF/S/42I/J/98h933Oj4ne52b33o0hM7juhCeIQkaW7ZhwF7qiKYcZ671/hRfRXr+FjuxvNt3sdikxTjiNT4UX8lXW9IbGm5prJuvyQUExPTaqakoaFBYWFhGjhwYIc1fznrAgAmu/pUxZ88zT7jgXqqIlhxnrtfaIhNT+QlS7oSTr7s6voTeck9OqPV7YElIyNDFRUVPmN79uzR2LFjFR4e3mFNZmZmd7cHAF2i5bKl5W8ebfepCkla/uZRtVzu9HMOEOe5J00ZHauSWXcqxhHpMx7jiFTJrDt7fCbL7/ngpqYmffrpp9712tpaOZ1ODRgwQAkJCSoqKtKJEye0adMmSVeeCFq7dq0KCws1f/58VVVV6aWXXvJ5+mfx4sWaOHGiVqxYoWnTpmn79u3au3ev9u/f3wWHCADd72Dt563eV/FlliRX43kdrP1cGcMH9lxjQYbz3LOmjI7VfckxN+abbqurqzVp0iTv+tUbX+fMmaPS0lK5XC7V1dV5P09KStLOnTv1yCOPaN26dYqLi9Pq1at93uGSmZmpzZs367HHHtOyZcs0fPhwbdmyxf93sABAgDScaf+PaGfq0DbOc88LDbEZEf78Dix33XWXOnp1S2lpaauxb37zm/rggw863O+MGTM0Y8YMf9sBACMM6Rf51UV+1KFtnOfei29rBoAuMD5pgGIdka1uULzKpivfwTI+aUBPthV0OM+9F4EFALqAiU9VBCPOc+9FYAGALnL1qYohUb7viArUUxXBivPcO/HWKADoQlNGxyrr64N0x5N7JEmlD4xT9q2D+T/+LsZ57n2YYQGALvblP5qBegS0N+A89y4EFgAAYDwuCQEAbkh9I8J0/Kf3B7oN9BBmWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGCwt0AwAQbPpGhOn4T+8PdBtAUGGGBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjdSqwrF+/XklJSYqMjFRaWpreeeedDuvXrVunUaNGqU+fPhoxYoQ2bdrk83lpaalsNlur5fz5851pDwAABBm/v615y5YtWrJkidavX6+srCz94he/UG5uro4ePaqEhIRW9SUlJSoqKtKLL76ocePG6eDBg5o/f7769++vvLw8b11UVJQ+/vhjn20jIyM7cUgAACDY+B1YVq5cqXnz5unBBx+UJK1atUq7d+9WSUmJiouLW9X/6le/0sMPP6z8/HxJ0i233KJ3331XK1as8AksNptNMTExnT0OAAAQxPy6JHThwgXV1NRo8uTJPuOTJ0/WgQMH2tymubm51UxJnz59dPDgQV28eNE71tTUpMTERA0dOlRTp07VoUOHOuylublZHo/HZwEAAMHJr8By8uRJtbS0KDo62mc8Ojpabre7zW1ycnL0y1/+UjU1NbIsS9XV1dq4caMuXryokydPSpJGjhyp0tJS7dixQ2VlZYqMjFRWVpaOHTvWbi/FxcVyOBzeJT4+3p9DAQAAN5BO3XRrs9l81i3LajV21bJly5Sbm6v09HSFh4dr2rRpmjt3riQpNDRUkpSenq5Zs2YpJSVF2dnZev3113XbbbdpzZo17fZQVFSkxsZG71JfX9+ZQwEAADcAvwLLoEGDFBoa2mo2paGhodWsy1V9+vTRxo0bde7cOR0/flx1dXUaNmyY+vXrp0GDBrXdVEiIxo0b1+EMi91uV1RUlM8CAACCk1+BJSIiQmlpaaqoqPAZr6ioUGZmZofbhoeHa+jQoQoNDdXmzZs1depUhYS0/eMty5LT6VRsbKw/7QEAgCDl91NChYWFKigo0NixY5WRkaENGzaorq5OCxYskHTlUs2JEye871r55JNPdPDgQU2YMEFffPGFVq5cqQ8//FCvvPKKd5/Lly9Xenq6br31Vnk8Hq1evVpOp1Pr1q3rosMEAAA3Mr8DS35+vk6dOqWnnnpKLpdLo0eP1s6dO5WYmChJcrlcqqur89a3tLTo2Wef1ccff6zw8HBNmjRJBw4c0LBhw7w1p0+f1kMPPSS32y2Hw6HU1FRVVlZq/Pjx13+EAADghmezLMsKdBNdwePxyOFwqLGxscvuZzl34ZKSH98tSTr6VI76Rvid7wCj8DsNwDTX+veb7xICAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAL0Ii2X/+e1SwdrP/dZBwCTEViAXuKtD126d+U+7/rcl9/XN1b8H731oSuAXQHAtSGwAL3AWx+69L1ff6A/eZp9xt2N5/W9X39AaAFgPAILEORaLlta/uZRtXXx5+rY8jePcnkIgNEILECQO1j7uVyN59v93JLkajyvg7Wf91xTAOAnAgsQ5BrOtB9WOlMHAIFAYAGC3JB+kV1aBwCBQGABgtz4pAGKdUTK1s7nNkmxjkiNTxrQk20BgF8ILECQCw2x6Ym8ZElqFVqurj+Rl6zQkPYiDQAEHoEF6AWmjI5Vyaw7NSTK7jMe44hUyaw7NWV0bIA6A4BrExboBgD0jCmjY5X19UG648k9kqTSB8Yp+9bBzKwAuCEwwwL0Il8OJ+OTBhBWANwwCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8ToVWNavX6+kpCRFRkYqLS1N77zzTof169at06hRo9SnTx+NGDFCmzZtalWzdetWJScny263Kzk5WeXl5Z1pDQAABCG/A8uWLVu0ZMkSLV26VIcOHVJ2drZyc3NVV1fXZn1JSYmKior05JNP6siRI1q+fLkWLlyoN99801tTVVWl/Px8FRQU6PDhwyooKNDMmTP13nvvdf7IAABA0LBZlmX5s8GECRN05513qqSkxDs2atQoTZ8+XcXFxa3qMzMzlZWVpX/913/1ji1ZskTV1dXav3+/JCk/P18ej0e7du3y1kyZMkX9+/dXWVnZNfXl8XjkcDjU2NioqKgofw6pXecuXFLy47slSUefylHfiLAu2S8QKPxOAzDNtf799muG5cKFC6qpqdHkyZN9xidPnqwDBw60uU1zc7MiIyN9xvr06aODBw/q4sWLkq7MsPzlPnNyctrd59X9ejwenwUAAAQnvwLLyZMn1dLSoujoaJ/x6Ohoud3uNrfJycnRL3/5S9XU1MiyLFVXV2vjxo26ePGiTp48KUlyu91+7VOSiouL5XA4vEt8fLw/hwIAAG4gnbrp1maz+axbltVq7Kply5YpNzdX6enpCg8P17Rp0zR37lxJUmhoaKf2KUlFRUVqbGz0LvX19Z05FAAAcAPwK7AMGjRIoaGhrWY+GhoaWs2QXNWnTx9t3LhR586d0/Hjx1VXV6dhw4apX79+GjRokCQpJibGr31Kkt1uV1RUlM8CAACCk1+BJSIiQmlpaaqoqPAZr6ioUGZmZofbhoeHa+jQoQoNDdXmzZs1depUhYRc+fEZGRmt9rlnz56v3CcAAOgd/H5EoLCwUAUFBRo7dqwyMjK0YcMG1dXVacGCBZKuXKo5ceKE910rn3zyiQ4ePKgJEyboiy++0MqVK/Xhhx/qlVde8e5z8eLFmjhxolasWKFp06Zp+/bt2rt3r/cpIgAA0Lv5HVjy8/N16tQpPfXUU3K5XBo9erR27typxMRESZLL5fJ5J0tLS4ueffZZffzxxwoPD9ekSZN04MABDRs2zFuTmZmpzZs367HHHtOyZcs0fPhwbdmyRRMmTLj+IwQAADc8v9/DYirewwJ8NX6nAZimW97DAgAAEAgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwXqcCy/r165WUlKTIyEilpaXpnXfe6bD+1VdfVUpKivr27avY2Fg98MADOnXqlPfz0tJS2Wy2Vsv58+c70x4AAAgyfgeWLVu2aMmSJVq6dKkOHTqk7Oxs5ebmqq6urs36/fv3a/bs2Zo3b56OHDmiN954Q++//74efPBBn7qoqCi5XC6fJTIysnNHBQAAgorfgWXlypWaN2+eHnzwQY0aNUqrVq1SfHy8SkpK2qx/9913NWzYMC1atEhJSUn6xje+oYcffljV1dU+dTabTTExMT4LAACA5GdguXDhgmpqajR58mSf8cmTJ+vAgQNtbpOZman/+q//0s6dO2VZlv70pz/pN7/5je6//36fuqamJiUmJmro0KGaOnWqDh061GEvzc3N8ng8PgsAAAhOfgWWkydPqqWlRdHR0T7j0dHRcrvdbW6TmZmpV199Vfn5+YqIiFBMTIxuvvlmrVmzxlszcuRIlZaWaseOHSorK1NkZKSysrJ07NixdnspLi6Ww+HwLvHx8f4cCgAAuIF06qZbm83ms25ZVquxq44ePapFixbp8ccfV01Njd566y3V1tZqwYIF3pr09HTNmjVLKSkpys7O1uuvv67bbrvNJ9T8paKiIjU2NnqX+vr6zhwKAAC4AYT5Uzxo0CCFhoa2mk1paGhoNetyVXFxsbKysvToo49KksaMGaObbrpJ2dnZevrppxUbG9tqm5CQEI0bN67DGRa73S673e5P+wAA4Abl1wxLRESE0tLSVFFR4TNeUVGhzMzMNrc5d+6cQkJ8f0xoaKikKzMzbbEsS06ns80wAwAAeh+/ZlgkqbCwUAUFBRo7dqwyMjK0YcMG1dXVeS/xFBUV6cSJE9q0aZMkKS8vT/Pnz1dJSYlycnLkcrm0ZMkSjR8/XnFxcZKk5cuXKz09Xbfeeqs8Ho9Wr14tp9OpdevWdeGhAugbEabjP73/qwsBwDB+B5b8/HydOnVKTz31lFwul0aPHq2dO3cqMTFRkuRyuXzeyTJ37lydOXNGa9eu1b/8y7/o5ptv1t13360VK1Z4a06fPq2HHnpIbrdbDodDqampqqys1Pjx47vgEAEAwI3OZrV3XeYG4/F45HA41NjYqKioqC7Z57kLl5T8+G5J0tGnctQ3wu98BwAAOnCtf7/5LiEAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjNepwLJ+/XolJSUpMjJSaWlpeueddzqsf/XVV5WSkqK+ffsqNjZWDzzwgE6dOuVTs3XrViUnJ8tutys5OVnl5eWdaQ0AAAQhvwPLli1btGTJEi1dulSHDh1Sdna2cnNzVVdX12b9/v37NXv2bM2bN09HjhzRG2+8offff18PPvigt6aqqkr5+fkqKCjQ4cOHVVBQoJkzZ+q9997r/JEBAICgYbMsy/JngwkTJujOO+9USUmJd2zUqFGaPn26iouLW9X//Oc/V0lJif7jP/7DO7ZmzRr97Gc/U319vSQpPz9fHo9Hu3bt8tZMmTJF/fv3V1lZWZt9NDc3q7m52bve2NiohIQE1dfXKyoqyp9Date5C5c0/pl/kyQdXHqP+kaEdcl+AQDAFR6PR/Hx8Tp9+rQcDkf7hZYfmpubrdDQUOu3v/2tz/iiRYusiRMntrnNH/7wBysiIsL63e9+Z12+fNlyu93WxIkTrYcffthbEx8fb61cudJnu5UrV1oJCQnt9vLEE09YklhYWFhYWFiCYKmvr+8wg/g1ZXDy5Em1tLQoOjraZzw6Olput7vNbTIzM/Xqq68qPz9f58+f16VLl/Q3f/M3WrNmjbfG7Xb7tU9JKioqUmFhoXf98uXL+vzzzzVw4EDZbDZ/DqtDV5NfV87coDXOc8/hXPcMznPP4Dz3jO48z5Zl6cyZM4qLi+uwrlPXOP4yEFiW1W5IOHr0qBYtWqTHH39cOTk5crlcevTRR7VgwQK99NJLndqnJNntdtntdp+xm2++2c8juXZRUVH8y9ADOM89h3PdMzjPPYPz3DO66zx3eCno//MrsAwaNEihoaGtZj4aGhpazZBcVVxcrKysLD366KOSpDFjxuimm25Sdna2nn76acXGxiomJsavfQIAgN7Fr6eEIiIilJaWpoqKCp/xiooKZWZmtrnNuXPnFBLi+2NCQ0MlXZlFkaSMjIxW+9yzZ0+7+wQAAL2L35eECgsLVVBQoLFjxyojI0MbNmxQXV2dFixYIOnKvSUnTpzQpk2bJEl5eXmaP3++SkpKvJeElixZovHjx3uvVy1evFgTJ07UihUrNG3aNG3fvl179+7V/v37u/BQO8dut+uJJ55odfkJXYvz3HM41z2D89wzOM89w4Tz7PdjzdKVF8f97Gc/k8vl0ujRo/Xcc89p4sSJkqS5c+fq+PHjevvtt731a9as0QsvvKDa2lrdfPPNuvvuu7VixQr91V/9lbfmN7/5jR577DF99tlnGj58uJ555hl9+9vfvv4jBAAAN7xOBRYAAICexHcJAQAA4xFYAACA8QgsAADAeAQWAABgPAJLOyorK5WXl6e4uDjZbDZt27Yt0C0FpeLiYo0bN079+vXTkCFDNH36dH388ceBbivolJSUaMyYMd63VGZkZPh82Si6R3FxsWw2m5YsWRLoVoLOk08+KZvN5rPExMQEuq2gdOLECc2aNUsDBw5U37599dd//deqqanp8T4ILO04e/asUlJStHbt2kC3EtT27dunhQsX6t1331VFRYUuXbqkyZMn6+zZs4FuLagMHTpUP/3pT1VdXa3q6mrdfffdmjZtmo4cORLo1oLW+++/rw0bNmjMmDGBbiVo3X777XK5XN7lj3/8Y6BbCjpffPGFsrKyFB4erl27duno0aN69tlnu/WrcNrTqe8S6g1yc3OVm5sb6DaC3ltvveWz/vLLL2vIkCGqqanxvtsH1y8vL89n/ZlnnlFJSYneffdd3X777QHqKng1NTXpH/7hH/Tiiy/q6aefDnQ7QSssLIxZlW62YsUKxcfH6+WXX/aODRs2LCC9MMMCozQ2NkqSBgwYEOBOgldLS4s2b96ss2fPKiMjI9DtBKWFCxfq/vvv17333hvoVoLasWPHFBcXp6SkJP393/+9Pvvss0C3FHR27NihsWPH6u/+7u80ZMgQpaam6sUXXwxILwQWGMOyLBUWFuob3/iGRo8eHeh2gs4f//hHfe1rX5PdbteCBQtUXl6u5OTkQLcVdDZv3qwPPvhAxcXFgW4lqE2YMEGbNm3S7t279eKLL8rtdiszM1OnTp0KdGtB5bPPPlNJSYluvfVW7d69WwsWLNCiRYu8X7/Tk7gkBGP88z//s/793//diO+QCkYjRoyQ0+nU6dOntXXrVs2ZM0f79u0jtHSh+vp6LV68WHv27FFkZGSg2wlqX75kf8cddygjI0PDhw/XK6+8osLCwgB2FlwuX76ssWPH6ic/+YkkKTU1VUeOHFFJSYlmz57do70wwwIjfP/739eOHTv0+9//XkOHDg10O0EpIiJCX//61zV27FgVFxcrJSVFzz//fKDbCio1NTVqaGhQWlqawsLCFBYWpn379mn16tUKCwtTS0tLoFsMWjfddJPuuOMOHTt2LNCtBJXY2NhW/1MzatQo1dXV9XgvzLAgoCzL0ve//32Vl5fr7bffVlJSUqBb6jUsy1Jzc3Og2wgq99xzT6snVR544AGNHDlSP/zhDxUaGhqgzoJfc3OzPvroI2VnZwe6laCSlZXV6lUTn3zyiRITE3u8FwJLO5qamvTpp59612tra+V0OjVgwAAlJCQEsLPgsnDhQr322mvavn27+vXrJ7fbLUlyOBzq06dPgLsLHj/60Y+Um5ur+Ph4nTlzRps3b9bbb7/d6iktXJ9+/fq1uv/qpptu0sCBA7kvq4v94Ac/UF5enhISEtTQ0KCnn35aHo9Hc+bMCXRrQeWRRx5RZmamfvKTn2jmzJk6ePCgNmzYoA0bNvR8Mxba9Pvf/96S1GqZM2dOoFsLKm2dY0nWyy+/HOjWgso//uM/WomJiVZERIQ1ePBg65577rH27NkT6LZ6hW9+85vW4sWLA91G0MnPz7diY2Ot8PBwKy4uzvr2t79tHTlyJNBtBaU333zTGj16tGW3262RI0daGzZsCEgfNsuyrJ6PSQAAANeOm24BAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYLz/B50QNB4OfsO5AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# larger variance but unbiased\n",
    "plt.errorbar(np.log10([10, 100, 1000, 10000, 100000, 1000000]), res_mean, yerr=res_std, fmt='o')\n",
    "plt.ylim(0.8, 1.2)\n",
    "plt.axhline(1, color='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.16"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}