Unbiased inverse mean.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 442,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from scipy.linalg import sqrtm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 443,
   "metadata": {},
   "outputs": [],
   "source": [
    "def u():\n",
    "    a = 1\n",
    "    b = 3\n",
    "    return np.random.uniform(a, b, 1)[0]\n",
    "\n",
    "eps = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 444,
   "metadata": {},
   "outputs": [],
   "source": [
    "def u():\n",
    "    return np.random.normal(2, 5, 1)[0]\n",
    "\n",
    "eps = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = np.array([1, 2])\n",
    "cov = np.array([[1, 1], [1, 1]])\n",
    "\n",
    "def u():\n",
    "    return np.random.multivariate_normal(mean, cov, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [],
   "source": [
    "base = np.array([[10, 1], [3, 4]])\n",
    "\n",
    "def u():\n",
    "    return base + 2 * np.random.normal(np.array([0, 0]), np.ones((2, 2)), size=(1, 2, 2))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10.45909839,  5.04959902],\n",
       "       [ 7.49204888,  5.56011157]])"
      ]
     },
     "execution_count": 255,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.10810811, -0.02702703],\n",
       "       [-0.08108108,  0.27027027]])"
      ]
     },
     "execution_count": 256,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.inv(base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 445,
   "metadata": {},
   "outputs": [],
   "source": [
    "def m1(k):\n",
    "    return sum(u() for _ in range(k)) / k\n",
    "\n",
    "def m2(k):\n",
    "    return sum(u()**2 for _ in range(k)) / k\n",
    "\n",
    "def w(k):\n",
    "    m1k = m1(k)\n",
    "    return np.min([1 / k * m1k, m1k / m2(k), eps])\n",
    "\n",
    "\n",
    "def p(k):\n",
    "    wk = w(k)\n",
    "    return 1 - np.sqrt(1 / k * sum((1 - wk * u())**2 for _ in range(k)))\n",
    "\n",
    "def q(n, k):\n",
    "    pk = p(k)\n",
    "    return pk * (1 - pk) ** n\n",
    "\n",
    "def beta(k, N):\n",
    "    wk = w(k)\n",
    "    qN = q(N, k)\n",
    "    return wk / qN * np.prod([1 - wk * u() for _ in range(N)])\n",
    "\n",
    "def naive(N):\n",
    "    return 1 / np.mean([u() for _ in range(N)])\n",
    "\n",
    "def naive2(N):\n",
    "    return np.mean([1 / u() for _ in range(N)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [],
   "source": [
    "def m1(k):\n",
    "    samples = [u() for _ in range()]\n",
    "    return np.sum(samples, axis=0) / k\n",
    "\n",
    "def m2(k):\n",
    "    samples = [u() for _ in range(k)]\n",
    "    summed = np.zeros((2, 2))\n",
    "    for sample in samples:\n",
    "        summed += sample.dot(sample)\n",
    "    return summed / k\n",
    "\n",
    "def w(k):\n",
    "    return eps\n",
    "    m1k = m1(k)\n",
    "    return np.min([1 / k * m1k, m1k / m2(k), eps])\n",
    "\n",
    "\n",
    "def p(k):\n",
    "    wk = w(k)\n",
    "    samples = [np.eye(2) - wk * u() for _ in range(k)]\n",
    "    summed = np.zeros((2, 2))\n",
    "    for sample in samples:\n",
    "        summed += sample.dot(sample)\n",
    "    \n",
    "    return np.eye(2) - sqrtm(1 / k * summed)\n",
    "\n",
    "def q(n, k):\n",
    "    pk = p(k)\n",
    "    qn = pk\n",
    "    for _ in range(n):\n",
    "        qn = qn.dot(np.eye(2) - pk)\n",
    "    return qn\n",
    "\n",
    "def beta(k, N):\n",
    "    wk = w(k)\n",
    "    qN = q(N, k)\n",
    "    samples = [np.eye(2) - wk * u() for _ in range(N)]\n",
    "    prod = samples[0]\n",
    "    for sample in samples[1:]:\n",
    "        prod = prod.dot(sample)\n",
    "    return wk * np.linalg.inv(qN).dot(prod)\n",
    "\n",
    "def naive(N):\n",
    "    return np.linalg.inv(np.mean([u() for _ in range(N)], axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 447,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fb0f1346bd0>"
      ]
     },
     "execution_count": 447,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa0klEQVR4nO3df3RU1b338fc3yWA0KD+CSATaYK8VA8RiIyqg9UpJKb9ytUil0McWf7S21rR1+UjFYmzparsWqwgu216WVPSptmIsigZrFWhZlraYpBD5IVeKgIkBQoBouIkmYT9/zCFNYkJIZpJJZn9ea2XNOfvsOWfvYfJhZ58zZ8w5h4iIxL+EWDdARES6hwJfRMQTCnwREU8o8EVEPKHAFxHxRFKsG3A6gwYNcunp6bFuhohIr1FUVHTEOXd+a9t6dOCnp6dTWFgY62aIiPQaZra/rW2a0hER8YQCX0TEEwp8ERFPKPBFRDyhwBcR8UT8BX7Jalg6GvL6hx9LVse6RSIiPUKPviyzw0pWw4t3Q11NeL3q3fA6QObs2LVLRKQHiK8R/vof/TvsT6mrCZeLiHguvgK/qrRj5SIiHomvwO83rGPlIiIeia/An7QIQmc3LwudHS4XEfFcfAV+5myYsRz6DQcs/DhjuU7YiogQb1fpQDjcFfAiIh8TXyN8ERFpkwJfRMQTCnwREU8o8EVEPKHAFxHxhAJfRMQTCnwREU8o8EVEPKHAFxHxhAJfRMQTCnwREU8o8EVEPKHAjwZ9j66I9ALxd7fM7qbv0RWRXkIj/Ejpe3TjVsHeArLzs8l8IpPs/GwK9hbEukkiEfF2hF+wt4Blxcs4eOIgQ1KGkHt5LtMumtbxHfWG79EtWR3+D6iqNPx1j5MW6a+PdhTsLSBvcx61DbUAlJ8oJ29zHkDn3iciPYCXI/xTv8zlJ8pxuMZf5s6M4P43JbVD5d3u1JRT1buA+/eUk84znNay4mWNYX9KbUMty4qXxahF4oOu/qvSy8CP5i/zsgH9qDFrVlZjxrIB/SJqY9RoyqlTDp44yNQPqnnlQBnb3jnAKwfKmPpBNQdPHIx10yRORXMg2paoBL6ZTTGz3Wa2x8wWtLL9a2ZWYWZbg5/bonHczmrrl7Yzv8y/C9WTlzqA9xITOQm8l5hIXuoAfheqj7CVUdIbppx6oDl1SeRVHuPChgYSgAsbGsirPMacOm9nQaWLLStexvXHjzQbZFx//EhU/6qM+N1rZonAo8BkoBR4w8zWOud2tqj6jHPurkiPFw1DUoYw9uDb5B6rYkhDAwcTE1k2oB//HHJxp/a1Dse6c/s2K09LGRKt5kam37BgOqeVcmlT7rEqznauWdnZzpF7rCpGLZJ4N/bgHvIqjzW+704NMvLYE7VjmGvxpu7wDsyuBvKcc18I1n8A4Jz7aZM6XwOyOhr4Az95qZt8/28ial9rqo+9wznvl5HQpOsnDf73vKH0HTCiQ/s6WnuUfVX7OMnJxrIEEkjvl87A5IHRanLnnajAHXkbc/9un7MEbNDFkHJ+DBvWw+17ve1t6RO7rx3ijboDmwmdPPnx8oQEQp8Yf8b7Wf3N8UXOuazWtkVjSmco0HQIWRqUtfQlMysxs3wzG97WzszsDjMrNLPCurq6KDTv4/qeONIs7AESXLi8owYmDyS9Xzp9EvoA0CehT88Je+BoYiL7Qkl8FJxn+MiMfaEkjiYmxrhlPVzSWR0rF4lQa2F/uvLOiMYIfxYwxTl3W7D+VeDKpqN5M0sFqp1zH5rZN4AvO+eub2/fWVlZrrCwMKL2tSqvP9Bavw3yjkf/eDGUnZ9N+Ynyj5WnpaTxp1l/ikGLeomWH6gDCJ0NM5brklbpGktHtzH9Ohy+t/2Md2NmXTrCLwOajtiHBWWNnHOVzrkPg9XHgM9G4bid19b8dRzOa0fzBLVXMmeHw73fcMDCjwp76UqTFoUHFU2Fzg6XR0k0Ljl4A7jYzEYQDvqbga80rWBmac65U8PMmcCuKBy38yYtan30FsUXtqcYkjKk1RH+kJ5yUrkny5ytgJfuc+q91oUfkow48J1z9WZ2F/AKkAj8xjm3w8x+BBQ659YCd5vZTKAeOAp8LdLjRqQbXtieIvfy3GafGAVITkwm9/LcGLZKRFrVxYOMiOfwu1KXzeF7Jmq3kRCRHu90c/j6FIkHpl00TQEvIn7eWkFExEcKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERT0Ql8M1sipntNrM9Zragle1nmdkzwfZ/mFl6NI4rIiJnLuLAN7NE4FHgi0AGMMfMMlpUuxU45pz7D2Ap8PNIjysiIh0TjRH+OGCPc26vc+4j4PdATos6OcATwXI+MMnMLArHFhGRMxSNwB8KvNtkvTQoa7WOc64eqAJSW9uZmd1hZoVmVlhRURGF5omICPTAk7bOuRXOuSznXNb5558f6+aIiMSNaAR+GTC8yfqwoKzVOmaWBPQDKqNwbBEROUPRCPw3gIvNbISZ9QFuBta2qLMWuCVYngVscM65KBxbRETOUFKkO3DO1ZvZXcArQCLwG+fcDjP7EVDonFsLrAT+n5ntAY4S/k9BRES6UcSBD+CcWwesa1G2qMlyLXBTNI4lIiKd0+NO2oqISNdQ4IuIeEKBLyLiCQW+iIgnFPgiIp5Q4IuIeEKBLyLiCQW+iIgnFPgiIp5Q4IuIeEKBLyLiCQW+iIgnFPgiIp5Q4IuIeEKBLyLiCQW+iIgnFPgiIp5Q4IuIeEKBLyLiCQW+iIgnFPgiIp5Q4IuIeEKBLyLiCQW+iIgnFPgiIp5Q4IuIeEKBLyLiCQW+iIgnFPgiIp5Q4IuIeEKBLyLiiYgC38wGmtmrZvZ28DigjXoNZrY1+FkbyTFFRKRzIh3hLwDWO+cuBtYH662pcc59JviZGeExRUSkEyIN/BzgiWD5CeC/ItyfiIh0kUgD/wLnXHmwfBC4oI16yWZWaGZ/N7P/Ot0OzeyOoG5hRUVFhM0TEZFTktqrYGavAUNa2bSw6YpzzpmZa2M3n3TOlZnZRcAGM3vTOfev1io651YAKwCysrLa2p+IiHRQu4HvnPt8W9vM7JCZpTnnys0sDTjcxj7Kgse9ZvZnYCzQauCLiEjXiHRKZy1wS7B8C/BCywpmNsDMzgqWBwETgJ0RHldERDoo0sD/GTDZzN4GPh+sY2ZZZvZYUOdSoNDMtgEbgZ855xT4IiLdrN0pndNxzlUCk1opLwRuC5Y3A2MiOY6IiEROn7QVEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERT0QU+GZ2k5ntMLOTZpZ1mnpTzGy3me0xswWRHFNERDon0hH+duBGYFNbFcwsEXgU+CKQAcwxs4wIjysiIh2UFMmTnXO7AMzsdNXGAXucc3uDur8HcoCdkRxbREQ6pjvm8IcC7zZZLw3KWmVmd5hZoZkVVlRUdHnjRER80e4I38xeA4a0smmhc+6FaDfIObcCWAGQlZXlor1/ERFftRv4zrnPR3iMMmB4k/VhQZmIiHSj7pjSeQO42MxGmFkf4GZgbTccV0REmoj0sswbzKwUuBooMLNXgvILzWwdgHOuHrgLeAXYBax2zu2IrNkiItJRkV6lswZY00r5e8DUJuvrgHWRHEtERCKjT9qKiHhCgS8i4gkFvoiIJyKawxcR6ay6ujpKS0upra2NdVN6peTkZIYNG0YoFDrj5yjwRSQmSktLOffcc0lPT2/v9izSgnOOyspKSktLGTFixBk/T1M6IhITtbW1pKamKuw7wcxITU3t8F9HCnwRiRmFfed15rVT4IuIeEKBLyK9QsHeArLzs8l8IpPs/GwK9hZ0extWrVrFXXfd1eq28ePHd9lx9+3bx+jRoyPej07aikiPV7C3gLzNedQ2hOesy0+Uk7c5D4BpF02LYcv+bfPmzbFuQrsU+CIScz/f8nPeOvpWm9tLKkr46ORHzcpqG2pZ9NdF5P9PfqvPGTlwJPeNu++0x923bx/Tp09n+/btACxZsoTq6mr+/Oc/c+WVV7Jx40aOHz/OypUrueaaawB49913ue666ygrK2PevHk8+OCDAPTt25fq6mqqq6vJycnh2LFj1NXVsXjxYnJycjhx4gSzZ8+mtLSUhoYGfvjDH/LlL3+ZoqIivv/971NdXc2gQYNYtWoVaWlpFBUVMX/+fACys7PP7IVsh6Z0RKTHaxn27ZVHQ319PVu2bOHhhx/moYceaizfsmULzz33HCUlJTz77LMUFhY2e15ycjJr1qyhuLiYjRs3cs899+Cc449//CMXXngh27ZtY/v27UyZMoW6ujq+853vkJ+f3xjwCxcuBODrX/86jzzyCNu2bYtanzTCF5GYa28knp2fTfmJ8o+Vp6Wk8fiUx7ukTTfeeCMAn/3sZ9m3b19j+eTJk0lNTW2s8/rrr5OVldW43TnH/fffz6ZNm0hISKCsrIxDhw4xZswY7rnnHu677z6mT5/ONddcw/bt29m+fTuTJ08GoKGhgbS0NI4fP87x48e59tprAfjqV7/Kyy+/HHGfFPgi0uPlXp7bbA4fIDkxmdzLcyPab1JSEidPnmxcb3pd+1lnnQVAYmIi9fX1jeUtL4dsuf7UU09RUVFBUVERoVCI9PR0amtr+fSnP01xcTHr1q3jgQceYNKkSdxwww2MGjWKv/3tb832cfz48Yj61RZN6YhIjzftomnkjc8jLSUNw0hLSSNvfF7EJ2wvuOACDh8+TGVlJR9++CEvvfRSu8959dVXOXr0KDU1NTz//PNMmDCh2faqqioGDx5MKBRi48aN7N+/H4D33nuPc845h3nz5nHvvfdSXFzMJZdcQkVFRWPg19XVsWPHDvr370///v15/fXXgfB/ItGgEb6I9ArTLpoW9StyQqEQixYtYty4cQwdOpSRI0e2+5xx48bxpS99idLSUubNm9dsOgdg7ty5zJgxgzFjxpCVldW4zzfffJN7772XhIQEQqEQv/rVr+jTpw/5+fncfffdVFVVUV9fz3e/+11GjRrF448/zvz58zGzqJ20Ned67veEZ2VluZYnREQkPuzatYtLL7001s3o1Vp7Dc2syDmX1Vp9TemIiHhCgS8i4gkFvoiIJxT4IiKeUOCLiHhCgS8i4gkFvoj0DiWrYeloyOsffixZHfEu27rt8G233cbOnTsj3n9b0tPTOXLkSJftvy364JWI9Hwlq+HFu6GuJrxe9W54HSBzdtQP99hjj0V9nz2BAl9EYu/lBXDwzba3l74BDR82L6urgRfugqInWn/OkDHwxZ+1e+j6+nrmzp1LcXExo0aN4sknn2Tq1KksWbKErKws7rzzTt544w1qamqYNWtW450zFyxYwNq1a0lKSiI7O5slS5ZQUVHBN7/5TQ4cOADAww8/zIQJE6isrGTOnDmUlZVx9dVXE6sPvCrwRaTnaxn27ZV3wO7du1m5ciUTJkxg/vz5/PKXv2y2/Sc/+QkDBw6koaGBSZMmUVJSwtChQ1mzZg1vvfUWZtZ4s7Pc3Fy+973vMXHiRA4cOMAXvvAFdu3axUMPPcTEiRNZtGgRBQUFrFy5MuJ2d4YCX0Rir72R+NLR4WmclvoNh69H9lWHw4cPb7wB2rx581i+fHmz7atXr2bFihXU19dTXl7Ozp07ycjIIDk5mVtvvZXp06czffp0AF577bVmc//vv/8+1dXVbNq0iT/84Q8ATJs2jQEDBkTU5s7SSVsR6fkmLYLQ2c3LQmeHyyN0utsdv/POOyxZsoT169dTUlLCtGnTqK2tJSkpiS1btjBr1ixeeuklpkyZAsDJkyf5+9//ztatW9m6dStlZWX07ds34jZGiwJfRHq+zNkwY3l4RI+FH2csj8oJ2wMHDjTenvjpp59m4sSJjdvef/99UlJS6NevH4cOHWr8EpLq6mqqqqqYOnUqS5cubfxWquzsbB555JHG52/duhWAa6+9lqeffhqAl19+mWPHjkXc7s7QlI6I9A6Zs7vkipxLLrmERx99lPnz55ORkcGdd97Jiy++CMBll13G2LFjGTlyZLOpnw8++ICcnBxqa2txzvGLX/wCgOXLl/Ptb3+bzMxM6uvrufbaa/n1r3/Ngw8+yJw5cxg1ahTjx4/nE5/4RNT7cSYiuj2ymd0E5AGXAuOcc63ey9jM9gEfAA1AfVu37mxJt0cWiV+6PXLkOnp75EhH+NuBG4H/PoO6/+mc6/5PGoiICBBh4DvndsHHT3qIiEjP010nbR3wJzMrMrM7TlfRzO4ws0IzK6yoqOim5omIxL92R/hm9howpJVNC51zL5zhcSY658rMbDDwqpm95Zzb1FpF59wKYAWE5/DPcP8iItKOdgPfOff5SA/inCsLHg+b2RpgHNBq4IuISNfo8ikdM0sxs3NPLQPZhE/2iohIN4oo8M3sBjMrBa4GCszslaD8QjNbF1S7AHjdzLYBW4AC59wfIzmuiEhP9vzzz3fp7ZU7K9KrdNYAa1opfw+YGizvBS6L5DgiIr3J888/z/Tp08nIyIh1U5rRJ21FJOYeenEHO997P6r7zLjwPB6cMeq0dX7729+yfPlyPvroI6688krmz5/P7bffzpYtW2hoaGDcuHE888wzpKenk5OTw7Fjx6irq2Px4sXk5OQA8OSTT7JkyRLMjMzMTO68807Wrl3LX/7yFxYvXsxzzz3Hpz71qaj2rbMU+CLipV27dvHMM8/w17/+lVAoxLe+9S12797NzJkzeeCBB6ipqWHevHmMHj2a+vp61qxZw3nnnceRI0e46qqrmDlzJjt37mTx4sVs3ryZQYMGcfToUQYOHMjMmTOZPn06s2bNinU3m1Hgi0jMtTcS7wrr16+nqKiIK664AoCamhoGDx7MokWLuOKKK0hOTm68VbJzjvvvv59NmzaRkJBAWVkZhw4dYsOGDdx0000MGjQIgIEDB3Z7Pzoi7gK/YG8By4qXcfDEQYakDCH38lymXTQt1s0SkR7GOcctt9zCT3/602bl5eXlVFdXU1dXR21tLSkpKTz11FNUVFRQVFREKBQiPT2d2traGLW88+Lq9sgFewvI25xH+YlyHI7yE+Xkbc6jYG9kX5AgIvFn0qRJ5Ofnc/jwYQCOHj3K/v37+cY3vsGPf/xj5s6dy3333QdAVVUVgwcPJhQKsXHjRvbv3w/A9ddfz7PPPktlZWXjPgDOPfdcPvjggxj06vTiaoS/rHgZtQ3N/9etbahlWfEyjfJFpJmMjAwWL15MdnY2J0+eJBQKkZOTQygU4itf+QoNDQ2MHz+eDRs2MHfuXGbMmMGYMWPIyspi5MiRAIwaNYqFCxfyuc99jsTERMaOHcuqVau4+eabuf3221m+fDn5+fk95qRtRLdH7modvT1y5hOZOD7eH8MouaUkmk0TkQjp9siR6+jtkeNqSmdISmu3/Gm7XETEJ3EV+LmX55KcmNysLDkxmdzLc2PUIhGRniOu5vBPzdPrKh2R3sE5p+/T6KTOTMfHVeBDOPQV8CI9X3JyMpWVlaSmpir0O8g5R2VlJcnJye1XbiLuAl9Eeodhw4ZRWlqKvuioc5KTkxk2bFiHnqPAF5GYCIVCjBgxItbN8EpcnbQVEZG2KfBFRDyhwBcR8USP/qStmVUA+1sUDwKOxKA5seRjn8HPfqvP/uiqfn/SOXd+axt6dOC3xswK2/rYcLzysc/gZ7/VZ3/Eot+a0hER8YQCX0TEE70x8FfEugEx4GOfwc9+q8/+6PZ+97o5fBER6ZzeOMIXEZFOUOCLiHii1wS+mU0xs91mtsfMFsS6PdFkZr8xs8Nmtr1J2UAze9XM3g4eBwTlZmbLg9ehxMwuj13LO8/MhpvZRjPbaWY7zCw3KI/bfptZspltMbNtQZ8fCspHmNk/gr49Y2Z9gvKzgvU9wfb0mHYgAmaWaGb/NLOXgnUf+rzPzN40s61mVhiUxfT93SsC38wSgUeBLwIZwBwzy4htq6JqFTClRdkCYL1z7mJgfbAO4dfg4uDnDuBX3dTGaKsH7nHOZQBXAd8O/k3jud8fAtc75y4DPgNMMbOrgJ8DS51z/wEcA24N6t8KHAvKlwb1eqtcYFeTdR/6DPCfzrnPNLnePrbvb+dcj/8BrgZeabL+A+AHsW5XlPuYDmxvsr4bSAuW04DdwfJ/A3Naq9ebf4AXgMm+9Bs4BygGriT8acukoLzxvQ68AlwdLCcF9SzWbe9EX4cRDrfrgZcAi/c+B+3fBwxqURbT93evGOEDQ4F3m6yXBmXx7ALnXHmwfBC4IFiOu9ci+LN9LPAP4rzfwdTGVuAw8CrwL+C4c64+qNK0X419DrZXAand2uDoeBj4v8DJYD2V+O8zgAP+ZGZFZnZHUBbT97fuh98LOOecmcXl9bNm1hd4Dviuc+79pt98FI/9ds41AJ8xs/7AGmBkbFvUtcxsOnDYOVdkZtfFuDndbaJzrszMBgOvmtlbTTfG4v3dW0b4ZcDwJuvDgrJ4dsjM0gCCx8NBedy8FmYWIhz2Tznn/hAUx32/AZxzx4GNhKcz+pvZqcFX03419jnY3g+o7N6WRmwCMNPM9gG/Jzyts4z47jMAzrmy4PEw4f/cxxHj93dvCfw3gIuDM/t9gJuBtTFuU1dbC9wSLN9CeI77VPn/Cc7qXwVUNfkTsdew8FB+JbDLOfeLJpvitt9mdn4wssfMziZ8zmIX4eCfFVRr2edTr8UsYIMLJnh7C+fcD5xzw5xz6YR/bzc45+YSx30GMLMUMzv31DKQDWwn1u/vWJ/Y6MAJkKnA/xCe81wY6/ZEuW+/A8qBOsJzd7cSnrdcD7wNvAYMDOoa4SuW/gW8CWTFuv2d7PNEwnOcJcDW4GdqPPcbyAT+GfR5O7AoKL8I2ALsAZ4FzgrKk4P1PcH2i2Ldhwj7fx3wkg99Dvq3LfjZcSqzYv3+1q0VREQ80VumdEREJEIKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ88f8BJhJFASA+hLwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "exact = np.linalg.inv(base)\n",
    "exact = 1/2\n",
    "\n",
    "for n in 2 ** np.arange(4, 10):\n",
    "    beta_sample = np.mean([beta(n, n) for _ in range(1000)])\n",
    "    naive_sample = np.mean([naive(n) for _ in range(1000)])\n",
    "    \n",
    "    if n == 2 ** 4:\n",
    "        plt.plot(n, beta_sample, marker='o', color='tab:green', label='unbiased')\n",
    "        plt.plot(n, naive_sample, marker='o', color='tab:orange', label='biased')\n",
    "    else:\n",
    "        plt.plot(n, beta_sample, marker='o', color='tab:green')\n",
    "        plt.plot(n, naive_sample, marker='o', color='tab:orange')\n",
    "    \n",
    "plt.axhline(exact, label='exact')\n",
    "plt.legend(loc='best')\n",
    "    \n",
    "#     plt.errorbar(n, np.mean(beta_errs), yerr=np.std(beta_errs), marker='o', color='tab:green')\n",
    "#     plt.errorbar(n, np.mean(naive_errs), yerr=np.std(naive_errs), marker='o', color='tab:blue')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 276,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.10779232, -0.02863408],\n",
       "       [-0.08016343,  0.26781071]])"
      ]
     },
     "execution_count": 276,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta(1000, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.10826602, -0.02592476],\n",
       "       [-0.08116969,  0.26595243]])"
      ]
     },
     "execution_count": 275,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "naive(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 7,  6],\n",
       "       [18, 19]])"
      ]
     },
     "execution_count": 179,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.dot(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 4,  1],\n",
       "       [ 9, 16]])"
      ]
     },
     "execution_count": 177,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.prod([a, a], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 456,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 3\n",
    "var = 2\n",
    "\n",
    "def g(x):\n",
    "    return x ** 2\n",
    "\n",
    "def estimate(n):\n",
    "    samples = np.random.normal(mu, var, n)\n",
    "    return np.mean(samples)\n",
    "\n",
    "def bestimate(n):\n",
    "    samples = np.random.normal(mu, var, n)\n",
    "    return np.mean(samples) + 1 / n \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 460,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7fb0d946a790>"
      ]
     },
     "execution_count": 460,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAD8CAYAAAC7IukgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWKElEQVR4nO3df6xc5Z3f8fcnxiQXiq43BFbBdgMNyKrlOHZ1RdOSbLRQgtnUEJHihY0aiSDTdJdkV6qooFoRiCKFBalpSQkV0SK2Uha4sOCNyQ8HBUcoFQXMAgaWuDI0WXzTrbub2A2Jg3/w7R8zhuvrufYdz7l35nreL2k0c5555jzPmXs0Hz3nOfecVBWSJDXpHf3ugCTp+GO4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhp3Qr87MCje85731JlnntnvbkjSvPLMM8/8XVWdNrXccGk788wz2bJlS7+7IUnzSpKfdCr3sJgkqXGGiySpcYaLJKlxhoskqXGGSy+2jsNXVsBNi1rPW8f73SNJGgieLXasto7Dxs/Dvj2t5d2vtZYBVq7rX78kaQA4cjlW3//i28Fy0L49rXJJGnKGy7HavaO7ckkaIobLsRpd0l25JA0Rw+VYXXAjLBw5tGzhSKtckoac4XKsVq6DtbfD6FIgree1tzuZL0kc52eLJTkZ+BqwF/hBVX2j0QZWrjNMJKmDnkYuSd6V5Kkkzyd5KcnNHeosTbI5yV+36/xhD+3dnWRnkhc7vLcmybYk25Nc3y6+DHiwqtYDlxxru5Kk7vR6WOwN4Pyq+iCwCliT5ENT6uwH/l1VLQc+BPxBkuWTKyQ5PckpU8rO7tDePcCaqYVJFgB3ABcDy4Er220sAV5rVzvQ3aZJko5VT+FSLa+3Fxe2HzWlzv+uqr9qv/4F8DKweMqqPgpsSPJOgCTrga92aO9x4GcdunIusL2qXq2qvcB9wKXADloBA9Nsa5K1Se7avXv30TZXkjRDPU/oJ1mQ5DlgJ/BoVT15hLpnAquBQ+pU1QPAJuD+JJ8CPgNc3kU3FvP2CAVaobIYeAj4ZJI7gY2dPlhVG6vqmtHR0S6akyQdSc8T+lV1AFiVZBHwcJIVVdVpTuQfAH8B/FFV/b8O67k1yX3AncD7J42IeunbL4Grel2PJKk7jZ2KXFW7gM10nhNZSCtYvlFVD3X6fJKPACuAh4EvdNn8BLB00vKSdpkkqQ96PVvstPaIhSQjwIXAj6bUCfCnwMtV9R+nWc9q4C5a8yRXAacm+VIXXXkaOCfJWUlOBK4Avtnl5kiSGtLryOW9wOYkW2n9wD9aVY8AJPl2kjOA84B/DZyf5Ln243emrOckYF1VvVJVbwKfBg67L3OSe4EngGVJdiS5GqCq9gPX0pq3eRkYr6qXetw2SdIxSlUdvdYQGBsbqy1btvS7G5I0ryR5pqrGppZ7+RdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjDBdJUuMMF0lS4wwXSVLjTuh3B2ZTkpOBrwF7gR9U1Tf63CVJGgpHHbkkeVeSp5I8n+SlJDdPU+/uJDuTvDil/MdJXkjyXJItvXT2CG2sSbItyfYk10966zLgwapaD1zSS9uSpJmbyWGxN4Dzq+qDwCpgTZIPdah3D7BmmnX8dlWtqqqxTm8mOT3JKVPKzp5JG0kWAHcAFwPLgSuTLG+/vQR4rf36wDR9kyQ17KjhUi2vtxcXth/Vod7jwM+OsR8fBTYkeSdAkvXAV2fYxrnA9qp6tar2AvcBl7bf20ErYMD5JUmaMzP6wU2yIMlzwE7g0ap6sos2CvhekmeSXNOxQtUDwCbg/iSfAj4DXD7D9S/m7dEJtAJlcfv1Q8Ank9wJbOz04SRrk9y1e/fuGTYnSTqaGYVLVR2oqlW0RgHnJlnRRRsfrqp/Quuw1R8k+a1p2rgV+DVwJ3DJpNHSMauqX1bVVVX1b6ebzK+qjVV1zejoaK/NSZLaujpUVFW7gM1MP7fS6TMT7eedwMO0DmMdJslHgBXtOl/oolsTwNJJy0vaZZKkPpnJ2WKnJVnUfj0CXAj8aCYrT3LywYn69mnBHwNe7FBvNXAXrbmSq4BTk3xphtvwNHBOkrOSnAhcAXxzhp+VJM2CmYxc3gtsTrKV1g/5o1X1CECSbyc5o/36XuAJYFmSHUmuBn4T+GGS54GngG9V1Xc7tHESsK6qXqmqN4FPAz+ZWqlTG1W1H7iW1pzNy8B4Vb3UzZcgSWpWqg478WsojY2N1ZYtPf0bjiQNnSTPdPo3E0/PlSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znCRJDXOcJEkNc5wkSQ1znAZJlvH4Ssr4KZFreet4/3ukaTj1An97oDmyNZx2Ph52Lentbz7tdYywMp1/euXpOOSI5dh8f0vvh0sB+3b0yqXpIYZLsNi947uyiWpB4bLsBhd0l25JPXAcBkWF9wIC0cOLVs40iqXpIYZLsNi5TpYezuMLgXSel57u5P5kmaFZ4sNk5XrDBNJc8KRiySpcYaLJKlxx/1hsSQnA18D9gI/qKpv9LlLknTc62nkkuRdSZ5K8nySl5LcPE29u5PsTPJiL+0daV1J1iTZlmR7kusnvXUZ8GBVrQcu6bV9SdLR9XpY7A3g/Kr6ILAKWJPkQx3q3QOsmW4lSU5PcsqUsrOnqX7YupIsAO4ALgaWA1cmWd5+ewnwWvv1gSNsiySpIT2FS7W83l5c2H5Uh3qPAz87wqo+CmxI8k6AJOuBr07TZqd1nQtsr6pXq2ovcB9wafu9HbQCBjpsb5K1Se7avXv3EbonSepGzxP6SRYkeQ7YCTxaVU92u46qegDYBNyf5FPAZ4DLu1jFYt4enUArUBa3Xz8EfDLJncDGDm1vrKprRkdHu+22JGkaPU/oV9UBYFWSRcDDSVZUVddzK1V1a5L7gDuB908aEfXav18CVzWxLknSzDR2KnJV7QI2c4S5lSNJ8hFgBfAw8IUuPz4BLJ20vKRdJknqg17PFjutPWIhyQhwIfCjY1jPauAuWvMkVwGnJvlSF6t4GjgnyVlJTgSuAL7ZbT8kSc3odeTyXmBzkq20fuAfrapHAJJ8O8kZ7df3Ak8Ay5LsSHL1lPWcBKyrqleq6k3g08BPOjXYaV1VtR+4lta8zcvAeFW91OO2SZKOUaoOO7lrKI2NjdWWLVv63Q1JmleSPFNVY1PLvfyLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhoskqXGGiySpcYaLJKlxhosG09Zx+MoKuGlR63nreL97JKkLJ/S7A9Jhto7Dxs/Dvj2t5d2vtZYBVq7rX78kzZgjFw2e73/x7WA5aN+eVrmkecFw0eDZvaO7ckkD57gOlyQnJ/mzJF9P8ql+90czNLqku3JJA+eo4ZLkXUmeSvJ8kpeS3DxNvTVJtiXZnuT6SeU/TvJCkueSbOmls0nuTrIzyYszaRu4DHiwqtYDl/TStubQBTfCwpFDyxaOtMolzQszGbm8AZxfVR8EVgFrknxocoUkC4A7gIuB5cCVSZZPqvLbVbWqqsY6NZDk9CSnTCk7u0PVe4A1XbS9BHit/frAUbZTg2LlOlh7O4wuBdJ6Xnu7k/nSPHLUs8WqqoDX24sL24+aUu1cYHtVvQqQ5D7gUuCvZ9iPjwKfTfI7VfVGkvW0Rh0XT+nL40nO7KLtHbQC5jmmCdIka4G1Z5/dKcvUNyvXGSbSPDajOZckC5I8B+wEHq2qJ6dUWczbIwRo/agvbr8u4HtJnklyTaf1V9UDwCbg/vbcyGeAy2e4DUdq+yHgk0nuBDZO0/bGqrpmdHR0hs1Jko5mRv/nUlUHgFVJFgEPJ1lRVS8e5WMHfbiqJpKcDjya5EdV9XiHNm5tjzruBN5fVa8ftqYuVdUvgat6XY8kqTtdnS1WVbuAzUyZ9wAmgKWTlpe0y6iqg887gYdpHcY6TJKPACvadb7QRbembVuS1B8zOVvstPaIhSQjwIXAj6ZUexo4J8lZSU4ErgC+2T4V+JT2Z08GPgYcNuJJshq4i9ZcyVXAqUm+NMNt6Nj2DD/bkw3PTnDeLY9x1vXf4rxbHmPDs2aaJMHMRi7vBTYn2Urrh/zRqnoEIMm3k5xRVfuBa2nNm7wMjFfVS8BvAj9M8jzwFPCtqvpuhzZOAtZV1StV9SbwaeAnUysluRd4AliWZEeSq4/Q9qza8OwENzz0AhO79lDAxK493PDQCwaMJAFpnQymsbGx2rJl5v+Gc94tjzGxa89h5YsXjfDfrz+/ya5J0sBK8kynfzM5rv9Dfzb9tEOwHKlckoaJ4XKMzlg00lW5JA0Tw+UYXXfRMkYWLjikbGThAq67aFmfeiRJg8P7uRyjT6xu/Z/mbZu28dNdezhj0QjXXbTsrXINma3jrVsC7N7RusDmBTd6hQENNcOlB59YvdgwkTc3kzrwsJjUK29uJh3GcJF65c3NpMMYLlKvvLmZdBjDReqVNzeTDmO4SL3y5mbSYTxbTANpw7MT8+s0b29uJh3CcNHAOXhR0D37WnemPnhRUGCwA0bSWzwspoFz26ZtbwXLQXv2HeC2Tdv61CNJ3TJcNHC8KKg0/xkuGjheFFSa/wwXDRwvCirNf07oa+B4UVBp/jNcNJC8KKg0v3lYTJLUOMNFktQ4w0WS1DjDRZLUOMNFktQ4w0WS1DjDRZLUOMNFGjZbx+ErK+CmRa3nreP97pGOQ/4TpTRMto7Dxs/DvvZFQHe/1loG70ejRjlykYbJ97/4drActG9Pq1xqkOEiDZPdO7orHwQexpuXDBdpmIwu6a683w4extv9GlBvH8YzYAae4SINkwtuhIVT7ouzcKRVPog8jDdvGS7SMFm5DtbeDqNLgbSe194+uJP58/EwngDPFpOGz8p1gxsmU40uaR8S61Cu3m0db40Cd+9ofacX3NjYvuHIRdLgmm+H8eaTWZ7PMlwkDa6V63j6Azfzt5zGmxX+ltN4+gM3z5+R1yCb5fksD4tJGlgbnp3ghqffx559//mtspGnF/DlpRODeafSWTzM1LTavYN0Ud6t43rkkuTkJH+W5OtJPtXv/kiDYMOzE5x3y2Ocdf23OO+Wx9jw7ES/uzSt2zZtY8++A4eU7dl3gNs2betTj45g6zj7//Jzhxxm2v+XnxvY06b/D+/pqrxbsx4uSd6V5Kkkzyd5KcnNPazr7iQ7k7zY4b01SbYl2Z7k+nbxZcCDVbUeuORY25WOFxueneCGh15gYtceCpjYtYcbHnphYAPmp7v2dFXeT7/6zo2ccODXh5SdcODX/Oo7gzk/9OW9l/OrOvGQsl/ViXx57+WNrD9V1ciKpm0gCXByVb2eZCHwQ+APq+p/TKpzOrCnqn4xqezsqto+ZV2/BbwO/LeqWjGpfAHwP4ELgR3A08CVwKXAd6rquSR/XlW/N10/3/2+f1wX/oe7G9hiaXA9+ze72HvgzcPKT1zwDlb/w0Vz36GjmFf9/fEPp3/vzA/PXT9m6Nm/2cUpB37O0uzknezjDRbyWp3OLxb8Rlff7fhn//kzVTU2tXzWRy7V8np7cWH7MTXRPgpsSPJOgCTrga92WNfjwM86NHMusL2qXq2qvcB9tIJlB3DwnMWO25pkbZK79u/f392GSfNQpx/qI5X329J3j/COKRMA70irfNC8wcKuyvtt6btH+HlGea7O4claznN1Dj/PaGPf7ZxM6LdHFs8AZwN3VNWTk9+vqgeSnAXcn+QB4DO0RiEztRiYfDL8DuCfArcD/yXJx4GNnT5YVRuBjWNjY+vv/zf/rIsmpfnnvFseY6LDIaXFi0YY1P1/w7MT3LZpGz/dtYczFo1w3UXLBnIy/6YvfZd/v+9rnJS9b5X9qk7k1oW/z03H8Xc7/tnO5XMSLlV1AFiVZBHwcJIVVfXilDq3JrkPuBN4/6TRTi/t/hK4qtf1SMeL6y5axg0PvXDIJPnIwgVcd9GyPvbqyD6xevFAhslUqz5+DTc+vJ8/qvs4I3/PT+tU/hNX8OGPX9Pvrk1rNr/bOT0Vuap2JdkMrAEOCZckHwFWAA8DXwCu7WLVE8DSSctL2mWSJjn4QzIfRgLzTes7/H1+d9MFfrfMzYT+acC+drCMAN8D/qSqHplUZzXw58C/BP4X8A3glar64w7rOxN4ZMqE/gm0JvQvoBUqTwO/V1UvzbSfY2NjtWXLlmPYQkkaXkn6M6EPvBfYnGQrrR/9RycHS9tJwLqqeqWq3gQ+Dfxk6oqS3As8ASxLsiPJ1QBVtZ/WSGcT8DIw3k2wSJKaNesjl/nCkYskda+fIxdJ0pDx2mJSA+bL6bLSXDFchog/gLPj4CVVDp7ee/CSKoDfr4aWh8WGxHy7ptR8Mq8urijNEcNlSPgDOHvm08UVpbliuAwJfwBnzxmLOl+LabpyaRgYLkPCH8DZc91FyxhZuOCQskG/pIo02wyXIeEP4Oz5xOrFfPmyD7B40QihdRHIL1/2ASfzNdQ8W2xIeE2p2TVfLq4ozRXDZYj4AyhprnhYTJLUOMNFktQ4w0WS1DjDRZLUOMNFktQ47+fSluT/0uEGZTM0Cuyeg8/NpH6vdY703nuAvzvKugfFsf5N+tXOXOxDM617tHrDsP/A3OxD82n/ma7++6rqtMNqVpWPHh/AXXPxuZnU77XOUd7b0u/verb/Jv1qZy72oZnWPVq9Ydh/mvzbzkUbg/QbdPDhYbFmbJyjz82kfq91jnVbBs1cbUdT7czFPjTTukerNwz7D8zNtsyn/aer+h4W04wl2VIdbmcqzYT7z3Bx5KJu3NXvDmhec/8ZIo5cJEmNc+QiSWqc4SJJapzhIklqnOGiY5LkHyX50yQP9rsvmp+SfCLJ15Pcn+Rj/e6PmmW46C1J7k6yM8mLU8rXJNmWZHuS6wGq6tWquro/PdWg6nIf2lBV64HPAr/bj/5q9hgumuweYM3kgiQLgDuAi4HlwJVJls991zRP3EP3+9Aft9/XccRw0Vuq6nHgZ1OKzwW2t0cqe4H7gEvnvHOaF7rZh9LyJ8B3quqv5rqvml2Gi45mMfDapOUdwOIkpyb5r8DqJDf0p2uaJzruQ8DngH8B/Kskn+1HxzR7Tuh3BzQ/VdXf0zpWLh2TqroduL3f/dDscOSio5kAlk5aXtIuk2bKfWgIGS46mqeBc5KcleRE4Argm33uk+YX96EhZLjoLUnuBZ4AliXZkeTqqtoPXAtsAl4GxqvqpX72U4PLfUgHeeFKSVLjHLlIkhpnuEiSGme4SJIaZ7hIkhpnuEiSGme4SJIaZ7hIkhpnuEiSGme4SJIa9/8BFsKqyJXkXnYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for n in 2 ** np.arange(2, 10):\n",
    "    meang = np.mean([estimate(n) for _ in range(1000)])\n",
    "    meanb = np.mean([bestimate(n) for _ in range(1000)])\n",
    "\n",
    "    plt.loglog(n, meang, 'o', color='tab:blue')\n",
    "    plt.loglog(n, meanb, 'o', color='tab:orange')\n",
    "\n",
    "    \n",
    "plt.axhline(mu)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 466,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQH0lEQVR4nO3dUWhc153H8d8vrlzUh5Wg8S5EstcpCqIhERhE9iEvhTRrh40aE4qa7L7VxLiQBvbBEEMJIS/O4oeAaCDrbYLpwyarhmDsjYu2pKShkAcrNTgORsQYFkt58KaLpy+COOa/D6OJR7KlndHcuffMPd8PGGWOZu78rWv95uTcc89xRAgAUH/3VF0AAKAcBD4AZILAB4BMEPgAkAkCHwAyQeADQCa+VXUBW7n33ntj7969VZcBAAPlk08++TIidm1sTzrw9+7dq8XFxarLAICBYvu/79ae5JCO7RnbJxuNRtWlAEBtJBn4EXE2Ig6PjIxUXQoA1EaSgQ8AKB6BDwCZqF/gX5yXXntIenm0+fXifNUVAUASkp6l07WL89LZF6Sbq83HjWvNx5I0NVtdXQCQgHr18D945XbYt9xcbbYDQObqFfiN5e7aASAj9Qr8kfHu2gEgI/UK/MdekoaG17cNDTfbASBz9Qr8qVlpZk4a2S3Jza8zc1ywBQAlOkvH9oykmYmJie5fPDVLwAPAXSTZw2dpBQAoXpKBDwAoHoEPAJkg8AEgEwQ+AGSCwAeATBD4AJAJAh8AMkHgA0AmCHwAyASBDwCZIPABIBMEPgBkgsAHgEwQ+ACQidIC3/b3bL9p+92y3hMAcFtHgW/7LdvXbV/a0H7A9pLtK7Zf3OoYEXE1Ig71UiwAYPs63fHqlKRfSvp1q8H2DkmvS3pc0rKk87bPSNoh6fiG1/80Iq73XC0AYNs6CvyI+Mj23g3Nj0i6EhFXJcn2O5Keiojjkp7cbkG2D0s6LEl79uzZ7mEAABv0MoY/Jula2+Pltba7sv1d229I2mf72GbPi4iTETEdEdO7du3qoTwAQLvSNjGPiD9LOlLW+wEA1uulh78iaXfb4/G1tp7ZnrF9stFoFHE4AIB6C/zzkh6wfb/tnZKekXSmiKIi4mxEHB4ZGSnicAAAdT4t821JH0uatL1s+1BEfC3peUkLki5Lmo+Iz/pXKgCgF53O0nl2k/Zzks4VWpGaQzqSZiYmJoo+NABkK8mlFRjSAYDiJRn4tXBxXnrtIenl0ebXi/NVVwQgc6VNy8zKxXnp7AvSzdXm48a15mNJmpqtri4AWUuyhz/w0zI/eOV22LfcXG22A0BFkgz8gR/Dbyx31w4AJUgy8AfeyHh37QBQAgK/Hx57SRoaXt82NNxsB4CKJBn4Az+GPzUrzcxJI7slufl1Zo4LtgAq5YiouoZNTU9Px+LiYtVlAMBAsf1JRExvbE+yh48CcT8AgDXMw68z7gcA0IYefp1xPwCANkkG/sBftE0F9wMAaJNk4A/8jVep4H4AAG2SDHwUhPsBALQh8OuM+wEAtGGWTt1NzRLwACTRwweAbCQZ+MzSAYDiJRn4zNIBgOLVbgz/9IUVnVhY0hc3VnXf6LCO7p/UwX1jVZcFAJWrVeCfvrCiY+99qtWbtyRJKzdWdey9TyWJ0AeQvSSHdLbrxMLSN2Hfsnrzlk4sLFVUEQCko1aB/8WN1a7aASAntQr8+0aHu2pHSViiGUhCrQL/6P5JDQ/tWNc2PLRDR/dPVlQRvlmiuXFNUtxeopnQB0qXZOBvdx7+wX1jOv70wxobHZYljY0O6/jTD3PBtkos0Qwkgy0O0V8vj0q6278xSy/fKLcWIBObbXFYq2mZuFPl9yWMjK8N59ylHUCpkhzSQTFa9yWs3FhV6PZ9CacvrJRXBEs0A8kg8GssifsSWKIZSAZDOjWWzH0JLNEMJIEefo1xXwKAdgR+jXFfAoB2DOnUWGs2DquHApAI/No7uG+MgAcgiSEdAMhGkoHPFocAULwkA58tDgGgeEkGPgCgeAQ+8sG6/Mgcs3SQh9a6/K2lmlvr8kvcBYxs0MNHHliXHyDwkYnGcnftQA0R+MjDZuvvsy4/MkLg98npCyt69NXf6/4X39ejr/6+3DXocSfW5Qe4aNsPrY1HWmvRtzYekcQyB1VpXZj94JXmMM7IeDPsuWCLjBD4fbDVxiMEfoVYlx+ZY0inD5LZeAQA2hD4fcDGIwBSROD3ARuPAEgRY/h9wMYjAFJUWuDbPijpHyT9laQ3I+K/ynrvKrDxCIDUdDSkY/st29dtX9rQfsD2ku0rtl/c6hgRcToinpN0RNJPtl8yAGA7Ou3hn5L0S0m/bjXY3iHpdUmPS1qWdN72GUk7JB3f8PqfRsT1tf/+xdrrAAAl6ijwI+Ij23s3ND8i6UpEXJUk2+9Ieioijkt6cuMxbFvSq5J+GxF/2uy9bB+WdFiS9uzZ00l5AIAO9DJLZ0zStbbHy2ttm/m5pB9K+rHtI5s9KSJORsR0REzv2rWrh/IAAO1Ku2gbEXOS5sp6PwDAer0E/oqk3W2Px9faemZ7RtLMxMREEYdDxU5fWGGKKpCAXoZ0zkt6wPb9tndKekbSmSKKYhPz+mgtJLdyY1Wh2wvJsXpoAtjyMTudTst8W9LHkiZtL9s+FBFfS3pe0oKky5LmI+Kz/pWKQbTVQnKoUGvLx8Y1SXF7y0dCv9Y6naXz7Cbt5ySdK7QiMaRTJywkl6ittnxkRdHaSnItHYZ06oOF5BLFlo9ZSjLwUR8sJJcotnzMEoGPvjq4b0zHn35YY6PDsqSx0WEdf/phZulULaUtH7l4XBpHRNU13KFtDP+5zz//vOpygHq6OF/9lo+ti8ft1xOGhqWZuWpqqfrnURDbn0TE9B3tKQZ+y/T0dCwuLlZdBoB+ee2htZlCG4zslv750p3t/ZLSB08BNgt8hnQAVCeVi8dbzVqqEQIfQHVSuXicygdPnyUZ+LZnbJ9sNBpVlwKgn1K5eJzKB0+fJRn4zMMHMjE12xwnH9ktyc2vVYybp/LB02fsaQugWlOz1V8Ybb1/TWbpbIbABwApjQ+ePktySAcAULwkA5+LtgBQvCQDn4u2AFA8xvCBkrEDGKpC4AMlau0A1toUprUDmKRsQ58PwPIkOaQD1BU7gK3HFpjlIvCBErED2Hp8AJYrycBnlg7qih3A1uMDsFxJBj6zdFBX7AC2Hh+A5Uoy8IG6Ygew9fgALBezdICSHdw3lm3Ab9T6OTBLpxwEPoBK8QFYHoZ0ACATBD4AZILAB4BMJDmGb3tG0szExETVpQC1xZIG+Umyh888fKC/WNIgT0n28IF+oEd721ZLGuT6M8kBgY8ssErleixpkKckh3SAorFI13osaZAnAh9ZoEe7Hksa5InARxaS6tFenJdee0h6ebT59eJ86SWwpk+eGMNHFo7un1w3hi9V1KO9OC+dfUG6ufZ/Fo1rzceSNDVbaiksaZAfAh9ZSGaRrg9euR32LTdXm+0lBz7S0++ZZAQ+spFEj7ax3F07slHGTDLG8IEyjYx3145snFhY0uO3/qA/7nxBV7/9j/rjzhf0+K0/FDqTLMnAZ4tD1NZjL0lDGy4UDw0325G16b/8Tq8O/Urj93ypeyyN3/OlXh36lab/8rvC3iPJwGdpBdTW1Kw0MyeN7Jbk5teZOcbvoWM7f6Pv+Kt1bd/xVzq28zeFvQdj+EDZpmYJeNzhb/RlV+3bkWQPHwBy402u42zWvh0EPgCkoITrOwQ+AKSghOs7jOEDQCr6fH2HHj4AZILAB4BMEPgAkAkCHwAyQeADQCYIfADIBIEPAJkg8AEgE6UFvu3v237D9ru2f1bW+wIAmjoKfNtv2b5u+9KG9gO2l2xfsf3iVseIiMsRcUTSrKRHt18yAGA7Ou3hn5J0oL3B9g5Jr0t6QtKDkp61/aDth23/54Y/f732mh9Jel/SucL+BgCAjnS0lk5EfGR774bmRyRdiYirkmT7HUlPRcRxSU9ucpwzks7Yfl/Sv2+7agBA13pZPG1M0rW2x8uS/m6zJ9v+gaSnJX1bW/TwbR+WdFiS9uzZ00N5AIB2pa2WGREfSvqwg+edlHRSkqanp6O/VQFAPnqZpbMiaXfb4/G1NgBAgnoJ/POSHrB9v+2dkp6RdKaIomzP2D7ZaDSKOBwAQJ1Py3xb0seSJm0v2z4UEV9Lel7SgqTLkuYj4rMiioqIsxFxeGRkpIjDAQDU+SydZzdpPyemWALAQEhyaQWGdACgeEkGPkM6AFC8JAMfAFC8JAOfIR0AKF6Sgc+QDgAUL8nABwAUj8AHgEwQ+ACQiSQDn4u2AFC8JAOfi7YAULwkAx8AUDwCHwAyQeADQCaSDHwu2gJA8ZIMfC7aAkDxkgx8AEDxCHwAyASBDwCZIPABIBNJBj6zdACgeEkGPrN0AKB4SQY+AKB4BD4AZILAB4BMEPgAkAkCHwAykWTgMy0TAIqXZOAzLRMAipdk4AMAikfgA0AmCHwAyASBDwCZIPABIBMEPgBkgsAHgEwQ+ACQCQIfADKRZOCztAIAFC/JwGdpBQAoXpKBDwAoHoEPAJn4VtUFAEAKTl9Y0YmFJX1xY1X3jQ7r6P5JHdw3VnVZhSLwAWTv9IUVHXvvU63evCVJWrmxqmPvfSpJtQp9hnQAZO/EwtI3Yd+yevOWTiwsVVRRfxD4ALL3xY3VrtoHFYEPIHv3jQ531T6oCHwA2Tu6f1LDQzvWtQ0P7dDR/ZMVVdQfXLQFkL3WhVlm6QBABg7uG6tdwG/EkA4AZILAB4BMEPgAkAkCHwAyQeADQCYcEVXXsCnb/yPphqTt7IQy0sXrOnnuVs/Zzvfu1n6vpC//nzrK0M3Prl/HKvP8bfX93M9fL8dL+Xew7ufvbyNi1x2tEZH0H0kn+/26Tp671XO28727tUtarPrn3cvPfFDPX7fnKafzNyjnsNvv5XT+2v8MwpDO2RJe18lzt3rOdr633b9XGYqsbRDO31bfz/389XK8lH8Hczp/30h6SCc3thcjYrrqOrA9nL/BlsP5G4Qefk5OVl0AesL5G2y1P3/08AEgE/TwASATBD4AZILAB4BMEPiJsv0922/afrfqWrA9tg/a/jfb/2H776uuB92x/X3bb9h+1/bPqq6nCAR+iWy/Zfu67Usb2g/YXrJ9xfaLkhQRVyPiUDWVYjNdnsPTEfGcpCOSflJFvVivy/N3OSKOSJqV9GgV9RaNwC/XKUkH2hts75D0uqQnJD0o6VnbD5ZfGjp0St2fw1+sfR/VO6Uuzp/tH0l6X9K5csvsDwK/RBHxkaT/3dD8iKQraz36ryS9I+mp0otDR7o5h276F0m/jYg/lV0r7tTt72BEnImIJyT9U7mV9geBX70xSdfaHi9LGrP9XdtvSNpn+1g1paFDdz2Hkn4u6YeSfmz7SBWFoSOb/Q7+wPac7X9VTXr47GmbqIj4s5pjvxhQETEnaa7qOrA9EfGhpA8rLqNQ9PCrtyJpd9vj8bU2DA7O4WDL5vwR+NU7L+kB2/fb3inpGUlnKq4J3eEcDrZszh+BXyLbb0v6WNKk7WXbhyLia0nPS1qQdFnSfER8VmWd2BzncLDlfv5YPA0AMkEPHwAyQeADQCYIfADIBIEPAJkg8AEgEwQ+AGSCwAeATBD4AJAJAh8AMvF/pwdZf51jxJ8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for n in 2 ** np.arange(2, 12):\n",
    "    meang = np.mean([estimate(n) for _ in range(1000)])\n",
    "    meanb = np.mean([bestimate(n) for _ in range(1000)])\n",
    "\n",
    "    plt.loglog(n, np.abs(mu - meang), 'o', color='tab:blue')\n",
    "    plt.loglog(n, np.abs(mu - meanb), 'o', color='tab:orange')\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 328,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7fb0f0b71990>"
      ]
     },
     "execution_count": 328,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAD8CAYAAAC2PJlnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYDUlEQVR4nO3df3DVdX7v8eeLABKoC12U4RJwgeqkMqw2bIbZ1ovOtB1+iCKDq9VraweZbN3RNt65Qwd22pE71+2q3OmU3VF6afHqzGX9FTEjG7vhD686d3QcA9glmGYHs7USYKDV0GtMFgjv/nFOMDmekHO+JPmek7weMwzJ5/s5n8/7yNfzyuf7+fJFEYGZmVkSk9IuwMzMypdDxMzMEnOImJlZYg4RMzNLzCFiZmaJOUTMzCyxyWkXMNauuuqqWLhwYdplmJmVjQMHDvxrRFyd79iEC5GFCxfS0tKSdhlmZmVD0sdDHfPlLDMzS8whYmZmiTlEzMwsMYeImZklNuE21pNoPNTJ9uZ2jnf1MG9WJZtXVbO+pirtsszMUucQGUbjoU627j1Mz7k+ADq7eti69zCAg8TMJryyvpwlabGk3ZIaRmuO7c3tFwOkX8+5PrY3t4/WlGZmZaPkQkTSM5JOSWrNaV8tqV3SUUlbACKiIyI2jWY9x7t6imo3M5tISi5EgGeB1QMbJFUATwFrgCXAvZKWjEUx82ZVFtVuZjaRlFyIRMTbwKc5zcuBo9mVx1ngBeCOQseU9F1JLZJaTp8+XVQ9m1dVUzmlYlBb5ZQKNq+qLmocM7PxqORCZAhVwCcDvj8GVEmaLelvgRpJW4d6cUTsiojaiKi9+uq8j38Z0vqaKn644ZtUzapEQNWsSn644ZveVDczo8zvzoqIfwMeHO151tdUOTTMzPIol5VIJ7BgwPfzs21mZpaicgmR94HrJC2SNBW4B3gt5ZrMzCa8kgsRSc8D7wLVko5J2hQR54GHgWagDXgpIo6kWaeZmZXgnkhE3DtE++vA62NcjpmZXULJrUTMzKx8OETMzCwxh4iZmSXmEDEzs8QcImZmlphDxMzMEnOImJlZYg4RMzNLzCFiZmaJOUTMzCwxh4iZmSXmEDEzs8QcImZmlphDxMzMEnOImJlZYg4RMzNLzCFiZmaJOUTMzCwxh4iZmSXmEDEzs8QcImZmlphDxMzMEnOImJlZYg4RMzNLzCFiZmaJOUTMzCwxh4iZmSU2Oe0CLpekGcDTwFngzYjYk3JJZmYTRkErEUn1klolHZH0yBB9/mv2eKuk5yVNS1KQpGcknZLUmufYakntko5K2pJt3gA0REQdsC7JnGZmlsywISJpKVAHLAduBG6TdG1Onyrgz4DaiFgKVAD35PSZI+nKnLZB42Q9C6zOU0cF8BSwBlgC3CtpCTAf+CTbrW+492NmZiOnkJXI9cB7EfFFRJwH3iLz03+uyUClpMnAdOB4zvFbgEZJVwBIqgN+nDtIRLwNfJpn/OXA0YjoiIizwAvAHcAxMkFS6PsxM7MRUsiHbiuwQtJsSdOBW4EFAztERCfwP4F/AU4AZyJif06fl4Fm4EVJ9wEPAHcVUWsVX644IBMeVcBe4E5JO4F9Q71Y0u2Sdp05c6aIKc3M7FKGDZGIaAOeAPYDPwM+IOeykaRfJ7MqWATMA2ZI+sM8Yz0J9AI7gXUR8fll1k9EdEfExoj43qU21SNiX0R8d+bMmZc7pZmZZRV0+ScidkfEtyLiZuAz4Bc5XX4f+GVEnI6Ic2RWB7+TO46kFcBS4FXg0SJr7WTwCmh+ts3MzFJS6N1Zc7K/X0NmP+QnOV3+Bfi2pOmSBPwe0JYzRg2wi8yKZSMwW9JjRdT6PnCdpEWSppLZuH+tiNebmdkIK3Qj+hVJH5LZc3goIroAJL0uaV5EvAc0AAeBw9lxd+WMMR24OyI+iogLwP3Ax7kTSXoeeBeolnRM0iaA7Kb+w2T2VdqAlyLiSFHv1szMRpQiIu0axlRtbW20tLSkXYaZWdmQdCAiavMd8y2xZmaWmEPEzMwSc4iYmVliDhEzM0vMIWJmZok5RMzMLDGHiJmZJeYQMTOzxBwiZmaWmEPEzMwSc4iYmVliDhEzM0vMIWJmZok5RMzMLDGHiJmZJeYQMTOzxBwiZmaWmEPEzMwSc4iYmVliDhEzM0vMIWJmZok5RMzMLDGHiJmZJeYQMTOzxBwiZmaWmEPEzMwSc4iYmVliDhEzM0vMIWJmZolNTruAkSBpBvA0cBZ4MyL2pFySmdmEUNBKRFK9pFZJRyQ9MkSfWZIaJP2TpDZJv520KEnPSDolqTWnfbWkdklHJW0ZcGgD0BARdcC6pPOamVlxhg0RSUuBOmA5cCNwm6Rr83TdAfwsIn4z268tZ5w5kq7Macs3DsCzwOqcvhXAU8AaYAlwr6Ql2cPzgU+yX/cN957MzGxkFLISuR54LyK+iIjzwFtkfvK/SNJM4GZgN0BEnI2IrpxxbgEaJV2RfU0d8ON8E0bE28CnOc3LgaMR0RERZ4EXgDuyx46RCZIh35Ok2yXtOnPmzDBv18zMClVIiLQCKyTNljQduBVYkNNnEXAa+N+SDkn6++w+xUUR8TLQDLwo6T7gAeCuImqt4svVBmSCoyr79V7gTkk7gX35XhwR+yLiuzNnzixiSjMzu5RhQyQi2oAngP3Az4AP+Oolo8nAMmBnRNQA3cCWnD5ExJNAL7ATWBcRn19O8QPG7Y6IjRHxPW+qm5mNnYI21iNid0R8KyJuBj4DfpHT5RhwLCLey37fQCZUBpG0AlgKvAo8WmStnQxeAc3PtpmZWUoKvTtrTvb3a8jsh/xk4PGIOAl8Iqk62/R7wIc5Y9QAu8jsY2wEZkt6rIha3weuk7RI0lTgHuC1Il5vZmYjrNC/bPiKpA/J7Dc81L9pLul1SfOyff4U2CPp58BvAX+VM8Z04O6I+CgiLgD3Ax/nm0zS88C7QLWkY5I2ZTf1Hyazr9IGvBQRRwqs38zMRoEiIu0axlRtbW20tLSkXYaZWdmQdCAiavMd82NPzMwsMYeImZkl5hAxM7PEHCJmZpaYQ8TMzBJziJiZWWIOETMzS8whYmZmiTlEzMwsMYeImZkl5hApQFNHEysbVnLDczewsmElTR1NaZdkZlYSJqddQKlr6mhi2zvb6O3rBeBE9wm2vbMNgLWL16ZYmZlZ+rwSGcaOgzsuBki/3r5edhzckVJFZmalwyEyjJPdJ4tqNzObSBwiw5g7Y25R7WZmE4lDZBj1y+qZVjFtUNu0imnUL6tPqaJLazzUyU2Pv8GiLU3c9PgbNB7yvyBsZqPHG+vD6N8833FwBye7TzJ3xlzql9WX5KZ646FOtu49TM+5PgA6u3rYuvcwAOtrqtIszczGKYdIAdYuXluSoZFre3P7xQDp13Ouj+3N7Q4RMxsVvpw1jhzv6imq3czscjlExpF5syqLajczu1wOkXFk86pqKqdUDGqrnFLB5lXVKVVkZuOd90TGkf59j+3N7Rzv6mHerEo2r6r2foiZjRqHyDizvqbKoWFmY8aXs8zMLDGHiJmZJeYQMTOzxBwiZmaWmEPEzMwSc4iYmVli4+IWX0kzgKeBs8CbEbEn5ZLMzCaEglYikuoltUo6IumRS/SrkHRI0k8vpyhJz0g6Jak1p321pHZJRyVtGXBoA9AQEXXAusuZ28zMCjdsiEhaCtQBy4EbgdskXTtE93qgbYhx5ki6MqdtqHGeBVbn9K0AngLWAEuAeyUtyR6eD3yS/XrwY2zNzGzUFLISuR54LyK+iIjzwFtkfvIfRNJ8YC3w90OMcwvQKOmKbP864Mf5OkbE28CnOc3LgaMR0RERZ4EXgDuyx46RCZJC35OZmY2AQj5wW4EVkmZLmg7cCizI0+9vgD8HLuQbJCJeBpqBFyXdBzwA3FVErVV8udqATHD0P99jL3CnpJ3AvnwvlnS7pF1nzpwpYkozM7uUYUMkItqAJ4D9wM+AD8i5ZCTpNuBURBwYZqwngV5gJ7AuIj5PVvZXxu2OiI0R8b2hNtUjYl9EfHfmzJkjMaWZmVHgpZ+I2B0R34qIm4HPgF/kdLkJWCfpn8lcZvpdSf8ndxxJK4ClwKvAo0XW2sngFdD8bJuZmaWk0Luz5mR/v4bMfshPBh6PiK0RMT8iFgL3AG9ExB/mjFED7CKzj7ERmC3psSJqfR+4TtIiSVOz87xWxOvNzGyEFboJ/YqkD8nsNzwUEV0Akl6XNK/AMaYDd0fERxFxAbgf+DhfR0nPA+8C1ZKOSdqU3dR/mMy+ShvwUkQcKXBuMzMbBYqItGsYU7W1tdHS0pJ2GWZmZUPSgYiozXfMt8OamVliDhEzM0vMIWJmZok5RMzMLDGHiJmZJeYQMTOzxBwiZmaWmEPEzMwSc4iYmVliDhEzM0vMIWJmZok5RMzMLDGHiJmZJeYQMTOzxBwiZmaWmEPEzMwSc4iYmVliDhEzM0vMIWJmZok5RMzMLDGHiJmZJeYQGWeaOppY2bCSG567gZUNK2nqaEq7JDMbxyanXYCNnKaOJra9s43evl4ATnSfYNs72wBYu3htipWZ2Xjllcg4suPgjosB0q+3r5cdB3ekVJGZjXcOkXHkZPfJotrNzC6XQ2QcmTtjblHtZmaXyyEyjtQvq2daxbRBbdMqplG/rD6lisxsvPPG+jjSv3m+4+AOTnafZO6MudQvq/emupmNGofIOLN28VqHhpmNGV/OMjOzxMbFSkTSDOBp4CzwZkTsSbkkM7MJoeCViKR6Sa2Sjkh6JM/xBZL+r6QPs30S7+ZKekbSKUmteY6tltQu6aikLdnmDUBDRNQB65LOa2ZmxSkoRCQtBeqA5cCNwG2Srs3pdh74bxGxBPg28JCkJTnjzJF0ZU5b7jgAzwKr89RRATwFrAGWAPdm55gPfJLt1lfIe7LS0Hiok5sef4NFW5q46fE3aDzUmXZJZlaEQlci1wPvRcQXEXEeeIvMT/8XRcSJiDiY/fr/A21AVc44twCNkq4AkFQH/Dh3soh4G/g0Tx3LgaMR0RERZ4EXgDuAY2SCpJj3ZClrPNTJ1r2H6ezqIYDOrh627j3sIDErI4V+4LYCKyTNljQduBVYMFRnSQuBGuC9ge0R8TLQDLwo6T7gAeCuIuqt4ssVB2TCowrYC9wpaSewb4iabpe068yZM0VMZ6Npe3M7PecGLxx7zvWxvbk9pYrMrFgFbaxHRJukJ4D9QDfwAUNcNpL0a8ArwCMR8e95xnpS0gvATuA3IuLzhLUPHLMb2DhMn33Avtra2rrLnc9GxvGunqLazaz0FHzpJyJ2R8S3IuJm4DPgF7l9JE0hEyB7ImJvvnEkrQCWAq8CjxZZbyeDV0Dzs21WhubNqiyq3cxKTzF3Z83J/n4Nmf2Qn+QcF7AbaIuIvx5ijBpgF5l9jI3AbEmPFVHv+8B1khZJmgrcA7xWxOuthGxeVU3llIpBbZVTKti8qjqlisysWMVsQr8i6UMyew4PRUQXgKTXJc0DbgL+CPhdSR9kf92aM8Z04O6I+CgiLgD3Ax/nTiTpeeBdoFrSMUmbALKb+g+T2VdpA16KiCNFvAcrIetrqvjhhm9SNasSAVWzKvnhhm+yvib3fgwzK1WKiLRrGFO1tbXR0tKSdhlWhhoPdbK9uZ3jXT3Mm1XJ5lXVDjybECQdiIjafMfGxd9YNxtt/bcj999N1n87MuAgsQnNf6fCrAC+HdksP4eIWQF8O7JZfg4RswL4dmSz/BwiZgXw7chm+Xlj3awA/ZvnvjvLbDCHiKWqqaOpbP453/U1VQ4NsxwOEUtNU0cT297ZRm9fLwAnuk+w7Z1tACUbJGY2mPdELDU7Du64GCD9evt62XFwR0oVmVmxHCKWmpPdJ4tqN7PS4xCx1MydMbeodjMrPQ4RS039snqmVUwb1DatYhr1y+pTqsjMiuWNdUtN/+Z5udydZWZf5RCxVK1dvNahYVbGfDnLzMwSc4iYmVliDhEzM0vMIWJmZok5RMzMLDGHiFmBmjqaWNmwkhueu4GVDStp6mhKuySz1PkWX7MC+GGRZvl5JWJWAD8s0iw/r0TMClCOD4tsPNTpf0TLRp1DxKwAc2fM5UT3ibztpajxUCdb9x6m51wfAJ1dPWzdexigJIPEgVe+fDnLrADl9rDI7c3tFwOkX8+5PrY3t6dU0dD6A6+zq4fgy8BrPNSZdmlWAK9EzApQbg+LPN7Vw+SvHeKKq5vRlC7i3Cx+dXoVx7tq0i7tKy4VeF6NlD6HiFmByulhkVfNPULPzL1o0jkANLWLaf9pL5XTpwKl9R6Od/UU1W6lxZezzMahK+Y0XwyQfpp0jivmNKdU0dDmzaosqt1Ki0PEbBz693Oni2pP0+ZV1VROqRjUVjmlgs2rqlOqyIrhy1lm41A53U22vqaKf/zsDV755d9xoeIzJvX9Oncuqivp/RDfTfYlr0TMxqFyupusqaOJnx7/ETH5MySIyZ/x0+M/KtnHyjQe6uT7+5+ja/ajzPjNLXTNfpTv73+uZO8mazzUyU2Pv8GiLU3c9PgbI17nuFiJSJoBPA2cBd6MiD0pl2SWqnK6m+xSTwMoxXp/8NYeJs1pGHTTQsxp4AdvTWZ9zZ+nXN1g/YGn2f/AjLlddJ2bxff3rwH+eMRWTiO6EpFUL6lV0hFJj1zGOM9IOiWpNc+x1ZLaJR2VtCXbvAFoiIg6YF3Sec3Gk7WL17L/O/v5+R//nP3f2V+SH8hQfk8D+GLGvrw3LXwxY19KFQ2tP/AmTe1CgklTu5g0p4EfvDVyP2ePWIhIWgrUAcuBG4HbJF2b02eOpCtz2gb1yXoWWJ1njgrgKWANsAS4V9ISYD7wSbZbX+7rzKx0DbVPU4r7NwCTpnQV1Z6msQi8kbycdT3wXkR8ASDpLTIrhCcH9LkFeFDSrRHxK0l12T5rBg4UEW9LWphnjuXA0YjoyM7xAnAHcIxMkHzAMMHYcbqbP/hf7xb/7sxsVERvPb1n/pkLXLjYNolJxMyFJfn/6q9OPURfnPtKe4WmlFy9vSfvBuU5EIxYrSN5OasVWCFptqTpwK3AgoEdIuJloBl4UdJ9wAPAXUXMUcWXKw7IhEcVsBe4U9JOIG/ESrpd0q7z588XMZ2ZjbavT/s6C2cuZOqkqQBMnTSVhTMX8vVpX0+5svy+8bUFKOeTWYhvfG3BEK9IT8WkKUW1JzFiK5GIaJP0BLAf6CazKvjKpaWIeDK7gtgJ/EZEfD4Cc3cDG4fpsw/YV1tbW/fin/z25U5pZhNYU0dTWdy00NTxKX/5/x7lXPzqYtsUXcH/+M//nbWLC/8cfOnBoY+N6N1ZEbEb2A0g6a/IrBQGkbQCWAq8CjwKPFzEFJ0MXt3Mz7aZmY2ZcnkEzljcpTeiISJpTkScknQNmb2Ob+ccrwF2AbcBvwT2SHosIv6iwCneB66TtIhMeNwD/JcRewNmZuPMaAfeSP9lw1ckfUhmX+KhiOjKOT4duDsiPoqIC8D9wMe5g0h6HngXqJZ0TNImgIg4T2bl0gy0AS9FxJERfg9mZlYgRUTaNYyp2traaGlpSbsMM7OyIelARNTmO+bHnpiZWWIOETMzS8whYmZmiU24PRFJp8mzmV+gmcCZUX5doX2H65f0+FXAvxYwf6lI+meSxhxjcf4U2v9y+1zqWDmdQ2Nx/ozkPGl9Bn0jIq7O2zsi/KvAX8Cu0X5doX2H65f0ONCS9n/nsfgzSWOOsTh/Cu1/uX2GOVY259BYnD8jOU8pfQb1//LlrOIkfWpZMa8rtO9w/S73eLkYi/cxUnOMxflTaP/L7ePzJ515SukzCJiAl7Ps0iS1xBC38pkVwufQxOKViOXalXYBVvZ8Dk0gXomYmVliXomYmVliDhEzM0vMIWJmZok5ROySJC2WtFtSQ9q1WPmRtF7S30l6UdLKtOuxkecQmYAkPSPplKTWnPbVktolHZW0BSAiOiJiUzqVWikq8vxpjIg64EHgD9Ko10aXQ2RiehZYPbBBUgXwFLAGWALcK2nJ2JdmZeBZij9//iJ73MYZh8gEFBFvA5/mNC8HjmZXHmeBF4A7xrw4K3nFnD/KeAL4h4g4ONa12uhziFi/KuCTAd8fA6okzZb0t0CNpK3plGZlIO/5A/wp8PvAdyQ9mEZhNrpG9N9Yt/EnIv6NzPVss6JFxI+AH6Vdh40er0SsXyewYMD387NtZoXw+TNBOUSs3/vAdZIWSZoK3AO8lnJNVj58/kxQDpEJSNLzwLtAtaRjkjZFxHngYaAZaANeiogjadZppcnnjw3kBzCamVliXomYmVliDhEzM0vMIWJmZok5RMzMLDGHiJmZJeYQMTOzxBwiZmaWmEPEzMwSc4iYmVli/wFSbO8a/JqSFwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for n in 2 ** np.arange(2, 10):\n",
    "    meang = np.mean([g(estimate(n)) for _ in range(10000)])\n",
    "    plt.loglog(n, meang, 'o', color='tab:blue')\n",
    "    plt.loglog(n, mu ** 2 + var / n, 'o', color='tab:green')\n",
    "#     plt.loglog(n, mu ** 2 + 2 * mu * var / n, 'o', color='tab:orange')\n",
    "\n",
    "    \n",
    "plt.axhline(g(mu))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 324,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAM1ElEQVR4nO3dT4ij9R3H8c+3cYRQZAbsXiZqq6wGpD0MBC9e22YLnTpYoXpsFxcP9hgwtHctObXUtmxbkULRioTgUksuUrx4cLY5jCKBRShu5qC2ZA7lAcfh28NMpjPbHc2fZ/L8ku/7BXuYX57k+e7uw5vs82SfmLsLALD8vlL0AACA+SD4ABAEwQeAIAg+AARB8AEgCIIPAEHcUfQAt2Nmm5I277rrrqcfeuihoscBgIVx/fr1T939wu0es5Q/h1+r1Xx7e7voMQBgYZjZdXev3e4xTukAQBAEHwCCIPgAEATBB4AgkvyUziw6vYFa3b52h5nW18pq1Kva2qgUPRYAFG6pgt/pDdRs7yjbP5AkDYaZmu0dSSL6AMJbqlM6rW7/OPYj2f6BWt1+QRMBQDqSDL6ZbZrZ1b29vYmetzvMJloHgEiSDL67X3P3K6urqxM9b32tPNE6AESSZPCn1ahXVV4pnVorr5TUqFcLmggA0rFUF21HF2b5lA4A/L+lCr50GH0CDwD/b6lO6QAAzkbwASAIgg8AQRB8AAiC4ANAEAQfAIJIMvjT3loBAHC2JIM/7a0VAABnSzL4AID8EXwACILgA0AQBB8AgiD4ABAEwQeAIAg+AARB8AEgCIIPAEEQfAAIguADQBBL9522i6TTG/CF6wDmhuAXpNMbqNneUbZ/IEkaDDM12zuSRPQBnIskT+lEuD1yq9s/jv1Itn+gVrdf0EQAll2SwY9we+TdYTbROgDMKsngR7C+Vp5oHQBmRfAL0qhXVV4pnVorr5TUqFcLmgjAsuOibUFGF2b5lA6AeSH4BdraqBB4AHPDKR0ACILgA0AQBB8AgiD4ABAEwQeAIAg+AARB8AEgCIIPAEEQfAAIguADQBAEHwCCIPgAEATBB4AgCD4ABEHwASAIgg8AQfAFKBhbpzfgG7qABTa34JvZVyX9RtJnkv7u7n+e174xu05voGZ7R9n+gSRpMMzUbO9IEtEHFsRMp3TM7CUz+9jM3rtl/ZKZ9c3shpk9d7T8uKTX3f1pST+YZb+Yv1a3fxz7kWz/QK1uv6CJAExq1nP4L0u6dHLBzEqSXpT0PUkPS3rKzB6WdI+kj442O10OJG93mE20DiA9MwXf3d+W9O9blh+RdMPdP3T3zyS9KukxSTd1GP0v3K+ZXTGzbTPb/uSTT2YZDzlaXytPtA4gPefxKZ2K/vdOXjoMfUVSW9IPzey3kq6d9WR3v+ruNXevXbhw4RzGwzQa9arKK6VTa+WVkhr1akETAZjU3C7auvt/JP14XvtDvkYXZvmUDrC4ziP4A0n3nvj5nqM1LLitjQqBBxbYeZzSeVfSg2Z2v5ndKelJSW9M8gJmtmlmV/f29s5hPACIadaPZb4i6R1JVTO7aWaX3f1zSc9K6kr6QNJr7v7+JK/r7tfc/crq6uos4wEATpjplI67P3XG+puS3pzltQEA+eJeOgAQBMEHgCCSDD4XbQEgf0kGn4u2AJC/JIMPAMgfwQeAIAg+AASRZPC5aAsA+Usy+Fy0BYD8JRl8AED+CD4ABEHwASCIJIPPRVsAyF+SweeiLWbV6Q306Atv6f7n/qpHX3hLnR7fwQPM7SsOgXnp9AZqtneU7R9IkgbDTM32jiTxjV0ILcl3+MAsWt3+cexHsv0Dtbr9giYC0kDwsXR2h9lE60AUBB9LZ32tPNE6EAXBx9Jp1Ksqr5ROrZVXSmrUqwVNBKQhyYu2ZrYpafPixYtFj4IFNLow2+r2tTvMtL5WVqNe5YItwjN3L3qGM9VqNd/e3i56DABYGGZ23d1rt3uMUzoAEATBB4AgCD4ABEHwASAIgg8AQRB8AAiC4ANAEEkGn/vhA0D+kgw+98MHgPwlGXwAQP4IPgAEQfABIAiCDwBBEHwACILgA0AQBB8AgiD4ABAEwQeAIJIMPrdWAID8JRl8bq0AAPlLMvgAgPzdUfQAAKROb6BWt6/dYab1tbIa9aq2NipFj4UlQ/CBgnV6AzXbO8r2DyRJg2GmZntHkog+csUpHaBgrW7/OPYj2f6BWt1+QRNhWRF8oGC7w2yidWBaBB8o2PpaeaJ1YFoEHyhYo15VeaV0aq28UlKjXi1oIiwrLtoCBRtdmOVTOjhvBB9IwNZGhcDj3HFKBwCCIPgAEATBB4AgCD4ABJFk8Lk9MgDkL8ngc3tkAMhfksEHAOSP4ANAEAQfAIIg+AAQBMEHgCAIPgAEQfABIAiCDwBBEHwACILgA0AQBB8AgiD4ABAEwQeAIAg+AARB8AEgiDuKHgDAYun0Bmp1+9odZlpfK6tRr2pro1L0WBgDwQcwtk5voGZ7R9n+gSRpMMzUbO9IEtFfAJzSATC2Vrd/HPuRbP9ArW6/oIkwCYIPYGy7w2yidaSF4AMY2/paeaJ1pIXgAxhbo15VeaV0aq28UlKjXi1oIkxibhdtzewBST+TtOruT8xrvwDyM7owy6d0FtNYwTezlyR9X9LH7v7NE+uXJP1SUknSH9z9hbNew90/lHTZzF6fbWQARdraqBD4BTXuO/yXJf1a0p9GC2ZWkvSipO9IuinpXTN7Q4fxf/6W5//E3T+eeVoAwNTGCr67v21m37hl+RFJN47eucvMXpX0mLs/r8N/DUzFzK5IuiJJ991337QvAwC4xSwXbSuSPjrx882jtdsys7vN7HeSNsysedZ27n7V3WvuXrtw4cIM4wEATprbRVt3/5ekZ+a1PwDAabO8wx9IuvfEz/ccrQEAEjRL8N+V9KCZ3W9md0p6UtIbeQxlZptmdnVvby+PlwMAaMzgm9krkt6RVDWzm2Z22d0/l/SspK6kDyS95u7v5zGUu19z9yurq6t5vBwAQON/SuepM9bflPRmrhMBAM4Ft1YAgCAIPgAEkWTwuWgLAPlLMvhctAWA/CUZfABA/gg+AATBl5gDWGqd3oD79x9JMvhmtilp8+LFi0WPAmCBdXoDNds7x1+8PhhmarZ3JClk9JM8pcNFWwB5aHX7x7EfyfYP1Or2C5qoWEkGHwDysDvMJlpfdgQfwNJaXytPtL7sCD6ApdWoV1VeKZ1aK6+U1KhXC5qoWFy0BbC0Rhdm+ZTOIXP3omc4U61W8+3t7aLHAICFYWbX3b12u8c4pQMAQRB8AAiC4ANAEAQfAIIg+AAQRJLB5wtQACB/SQafe+kAQP6SDD4AIH8EHwCCIPgAEATBB4AgCD4ABEHwASAIgg8AQSQZfP7jFQDkL8ng8x+vACB/SQYfAJA/gg8AQRB8AAiC4ANAEAQfAIIg+AAQBMEHgCAIPgAEQfABIIgkg8+tFQAgf0kGn1srAED+kgw+ACB/BB8AgiD4ABAEwQeAIAg+AARB8AEgCIIPAEEQfAAI4o6iBwAAHOr0Bmp1+9odZlpfK6tRr2pro5Lb6xN8AEhApzdQs72jbP9AkjQYZmq2dyQpt+hzSgcAEtDq9o9jP5LtH6jV7ee2D4IPAAnYHWYTrU+D4ANAAtbXyhOtTyPJ4HN7ZADRNOpVlVdKp9bKKyU16tXc9pFk8Lk9MoBotjYqev7xb6myVpZJqqyV9fzj3+JTOgCwjLY2KrkG/lZJvsMHAOSP4ANAEAQfAIIg+AAQBMEHgCDM3Yue4Uxm9omkf0759FVJ03yQf5Lnjbvtl2037eNfk/TpGPtPxbR/J0XsYx7Hz7jbz7rNFz22SMfQPI6fPPdTVIO+7u4Xbru1uy/lL0lXz/t54277ZdtN+7ik7aL/nOfxd1LEPuZx/Iy7/azbfMljC3MMzeP4yXM/KTVo9GuZT+lcm8Pzxt32y7ab9fFFMY/fR177mMfxM+72s27D8VPMflJqkKTET+ngi5nZtrvXip4Di4tjKJZlfocfwdWiB8DC4xgKhHf4ABAE7/ABIAiCDwBBEHwACILgLxEze8DM/mhmrxc9CxaTmW2Z2e/N7C9m9t2i50G+CH7izOwlM/vYzN67Zf2SmfXN7IaZPSdJ7v6hu18uZlKkasJjqOPuT0t6RtKPipgX54fgp+9lSZdOLphZSdKLkr4n6WFJT5nZw/MfDQviZU1+DP386HEsEYKfOHd/W9K/b1l+RNKNo3f0n0l6VdJjcx8OC2GSY8gO/ULS39z9H/OeFeeL4C+miqSPTvx8U1LFzO42s99J2jCzZjGjYUHc9hiS9FNJ35b0hJk9U8RgOD98p+0Scfd/6fDcKzAVd/+VpF8VPQfOB+/wF9NA0r0nfr7naA0YF8dQQAR/Mb0r6UEzu9/M7pT0pKQ3Cp4Ji4VjKCCCnzgze0XSO5KqZnbTzC67++eSnpXUlfSBpNfc/f0i50S6OIYwws3TACAI3uEDQBAEHwCCIPgAEATBB4AgCD4ABEHwASAIgg8AQRB8AAiC4ANAEP8Fq5gAGzKGzi4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for n in 2 ** np.arange(2, 10):\n",
    "    meang = np.mean([g(estimate(n)) for _ in range(10000)])\n",
    "    plt.loglog(n, np.abs(meang - mu ** 2 + 2 * mu * var / n), 'o', color='tab:blue')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 329,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7fb138d7e090>"
      ]
     },
     "execution_count": 329,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAD8CAYAAAC7IukgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZGUlEQVR4nO3df4zUdZ7n8eeLBqYbYrp3GBBo3AHE9C2BkXb7iDtzrjlmhXYQNd7oadwzpwSP5Nxj9w8ucreJmJhzVi6ZY+eMhhld9+buUGCAKLjTkkCO7OkpDaiAyAWZYaWF7Vmd5ga2+dW874/6whRFFV3V/e2uKur1SCpQn8+nvt93yTf18vP9fKu+igjMzMzSNKLcBZiZ2fXH4WJmZqlzuJiZWeocLmZmljqHi5mZpc7hYmZmqRtZ7gIqwTe+8Y2YOnVqucswM6squ3fv/oeIGJ+vz+ECTJ06lc7OznKXYWZWVSQdLdTn02JmZpY6h4uZmaXuug0XSdMlvSJpQ7lrMTOrNf2Gi6R6SR9I+kjSAUnPFhi3TNL+ZMyfDqYoSa9K6pa0P6e9XdIhSYclPX2tbUTEkYhYPJg6zMxsYIpZ0D8LzIuIU5JGAX8r6W8i4v9cGiBpFrAEmAucA34uaUtEHM4aMwHojYjfZLXNyB6T5TXgvwL/LWtsHfAicBdwDNgl6c2I+ETSbOD5nG08ERHdRby/Adm8t4tVHYf4oqeXyU0NLF/Qwv2tzUO1OzOzqtLvzCUyTiVPRyWP3J9S/j3g/Yj4x4i4APwv4IGcMXcCmyV9DUDSEuBHBfa5E/gqp3kucDiZkZwDXgfuS8bvi4h7ch79BoukRZLWnDx5sr+hV9i8t4sVG/fR1dNLAF09vazYuI/Ne7tK2o6Z2fWqqDUXSXWSPgS6gW0R8X7OkP3AHZLGSRoDfA+4KXtARKwHOoA3JD0KPAE8WEKtzcDnWc+PJW2Fah4n6WWgVdKKfGMi4q2IeLKxsbGEMmBVxyF6z/dd0dZ7vo9VHYdK2o6Z2fWqqO+5REQfMEdSE7BJ0qyI2J/Vf1DSXwDvAKeBD4G+PNt5QdLrwEvAzVkzotRFxJfA0qHY9hc9vSW1m5nVmpKuFouIHmAH0J6n75WI+P2I+EPg18D/zR0j6Q5gFrAJeKbEWru4cjY0JWkbdpObGkpqNzOrNcVcLTY+mbEgqYHMgvqnecZNSP78XTLrLf8zp78VWENmneRxYJyk50qodRdwi6RpkkYDDwNvlvD61Cxf0ELDqLor2hpG1bF8QUs5yjEzqzjFzFwmATskfUzmA35bRGwBkPS2pMnJuJ9J+gR4C/i3ySwn2xjgoYj4LCIuAo8BeX86QNJa4D2gRdIxSYuTCwWeIrNucxBYFxEHSnmzabm/tZnnH5hNc1MDApqbGnj+gdm+WszMLKGI3Au/ak9bW1v4t8XMzEojaXdEtOXru26/oW9mZuXjcDEzs9Q5XMzMLHUOFzMzS53DxczMUudwMTOz1DlczMwsdQ4XMzNLncPFzMxS53AxM7PUOVzMzCx1DhczM0udw8XMzFLncDEzs9Q5XMzMLHUOFzMzS53DxczMUudwMTOz1DlczMwsdQ4XMzNLncPFzMxS53AxM7PUOVzMzCx11224SJou6RVJG8pdi5lZrek3XCTVS/pA0keSDkh6tsC4P0v690taK6l+IAVJelVSt6T9efraJR2SdFjS09faTkQciYjFA6nBzMwGp5iZy1lgXkTcCswB2iXdnj1AUjPw74C2iJgF1AEP54yZIOmGnLYZefb3GtCe2yipDngRuBuYCTwiaaak2ZK25DwmFPG+zMxsiPQbLpFxKnk6KnlEnqEjgQZJI4ExwBc5/XcCmyV9DUDSEuBHefa3E/gqz/bnAoeTGck54HXgvojYFxH35Dy6+3tfSQ2LJK05efJkMcPNzKxIRa25SKqT9CHQDWyLiPez+yOiC/jPwN8Bx4GTEfFOzpj1QAfwhqRHgSeAB0uotRn4POv5saStUM3jJL0MtEpakW9MRLwVEU82NjaWUIaZmfWnqHCJiL6ImANMAeZKmpXdL+l3gPuAacBkYKykP86znReAM8BLwL1ZM6LURcSXEbE0Im6OiOeHaj9mZna1kq4Wi4geYAdXr4n8EfCLiPhVRJwHNgLfzn29pDuAWcAm4JkSa+0Cbsp6PiVpMzOzClPM1WLjJTUlf28A7gI+zRn2d8DtksZIEvBd4GDOdlqBNWRmOI8D4yQ9V0Ktu4BbJE2TNJrMBQNvlvB6MzMbJsXMXCYBOyR9TOYDfltEbAGQ9LakyckazAZgD7Av2e6anO2MAR6KiM8i4iLwGHA0d2eS1gLvAS2SjklaDBARF4CnyKzbHATWRcSBkt+xmZkNOUXku/CrtrS1tUVnZ2e5yzAzqyqSdkdEW76+6/Yb+mZmVj4OFzMzS53DxczMUudwMTOz1DlczMwsdQ4XMzNLncPFzMxS53AxM7PUOVzMzCx1DhczM0udw8XMzFLncDEzs9Q5XMzMLHUOFzMzS53DxczMUudwMTOz1DlczMwsdQ4XMzNLncPFzMxS53AxM7PUOVzMzCx1DhczM0udw8XMzFLncDEzs9Rd1+EiabqkVyRtKHctZma1pKhwkVQv6QNJH0k6IOnZnP4WSR9mPf6fpD8daFGSXpXULWl/Tnu7pEOSDkt6ur/tRMSRiFg80DrMzGxgip25nAXmRcStwBygXdLtlzoj4lBEzImIOcDvA/8IbMrdiKQJkm7IaZuRZ3+vAe054+qAF4G7gZnAI5JmJn2zJW3JeUwo8r2ZmVnKRhYzKCICOJU8HZU8osDw7wKfRcTRPH13AkslfS8izkpaAjxAJjCy97dT0tSc184FDkfEEQBJrwP3AZ9ExD7gnmLeSzZJi4BFM2bkyzczMxuootdcJNVJ+hDoBrZFxPsFhj4MrM3XERHrgQ7gDUmPAk8ADxZZQjPwedbzY0nbtWoeJ+lloFXSijz1vBURTzY2NhZZgpmZFaOomQtARPQBcyQ1AZskzYqI3DWR0cC9wFUf5FnbeSGZdbwE3BwRpwqNHayI+BJYOlTbNzOz/Eq+WiwieoAd5KyJJO4G9kTE3xd6vaQ7gFlk1mSeKWHXXcBNWc+nJG1mZlZhir1abHwyY0FSA3AX8GmeoY9Q4JRY8tpWYA2ZtZLHgXGSniuy1l3ALZKmJTOkh4E3i3ytmZkNo2JnLpOAHZI+JvMhvy0itgBIelvSZEljyYTOxmtsZwzwUER8FhEXgceAqxb+Ja0F3gNaJB2TtDgiLgBPkVmzOQisi4gDRdZvZmbDSJkLwWpbW1tbdHZ2lrsMM7OqIml3RLTl67uuv6FvZmbl4XAxM7PUOVzMzCx1DhczM0udw8XMzFLncDEzs9Q5XMzMLHUOFzMzS53DxczMUudwMTOz1DlczMwsdQ4XMzNLncPFzMxS53AxM7PUOVzMzCx1DhczM0udw8XMzFLncDEzs9Q5XMzMLHUOFzMzS53DxczMUudwMTOz1DlczMwsdQ4XMzNLncPFzMxSd92Gi6Tpkl6RtKHctZiZ1Zp+w0VSvaQPJH0k6YCkZwuMa5K0QdKnkg5K+oOBFiXpVUndkvbntLdLOiTpsKSnr7WNiDgSEYsHWoOZmQ1cMTOXs8C8iLgVmAO0S7o9z7jVwM8j4p8AtwIHszslTZB0Q07bjAL7fA1ozxlbB7wI3A3MBB6RNDPpmy1pS85jQhHvzczMhsDI/gZERACnkqejkkdkj5HUCPwh8K+T15wDzuVs6k5gqaTvRcRZSUuAB8iERe4+d0qamtM8FzgcEUeSfb4O3Ad8EhH7gHv6ey+5JC0CFs2YUSjjzMxsIIpac5FUJ+lDoBvYFhHv5wyZBvwK+CtJeyX9RNLY7AERsR7oAN6Q9CjwBPBgCbU2A59nPT+WtBWqeZykl4FWSSvyjYmItyLiycbGxhLKMDOz/hQVLhHRFxFzgCnAXEmzcoaMBG4DXoqIVuA0cNWaSES8AJwBXgLujYhTuWPSEhFfRsTSiLg5Ip4fqv2YmdnVSrpaLCJ6gB3krIeQmUUcy5rRbCATNleQdAcwC9gEPFNirV3ATVnPpyRtZmZWYYq5Wmy8pKbk7w3AXcCn2WMi4gTwuaSWpOm7wCc522kF1pBZJ3kcGCfpuRJq3QXcImmapNHAw8CbJbzezMyGSTEzl0nADkkfk/mA3xYRWwAkvS1pcjLuT4D/kYybA/ynnO2MAR6KiM8i4iLwGHA03w4lrQXeA1okHZO0OCIuAE+RWbc5CKyLiAMlvFczMxsmylwMVtva2tqis7Oz3GWYmVUVSbsjoi1f33X7DX0zMysfh4uZmaWu3y9R2vVh894uVnUc4oueXiY3NbB8QQv3txb8mpCZ2aA4XGrA5r1drNi4j97zfQB09fSyYuM+AAeMmQ0JnxarAas6Dl0Olkt6z/exquNQmSoys+udw6UGfNHTW1K7mdlgOVxqwOSmhpLazcwGy+FSA5YvaKFhVN0VbQ2j6li+oKXAK8zMBscL+jXg0qK9rxYzs+HicKkR97c2O0zMbNj4tJiZmaXO4WJmZqlzuJiZWeocLmZmljqHi5mZpc7hYmZmqXO4mJlZ6hwuZmaWOoeLmZmlzuFiZmapc7gM0NYjW5m/YT7f+utvMX/DfLYe2VrukszMKoZ/W2wAth7Zysp3V3Km7wwAx08fZ+W7KwFYOH1hGSszM6sMnrkMwOo9qy8HyyVn+s6wes/qMlVkZlZZHC4DcOL0iZLazcxqjcNlACaOnVhSu5lZrXG4DMCy25ZRX1d/RVt9XT3LbltWporMzCrLdbugL2k68B+Bxoj4fprbvrRov3rPak6cPsHEsRNZdtsyL+abmSX6DRdJ9cBO4GvJ+A0R8Uyecb8EfgP0ARciom2gRUl6FbgH6I6IWVnt7cBqoA74SUT8oNA2IuIIsFjShoHWcS0Lpy90mJiZFVDMabGzwLyIuBWYA7RLur3A2H8eEXPyBYukCZJuyGmbUWA7rwHtOWPrgBeBu4GZwCOSZiZ9syVtyXlMKOK9mZnZEOh35hIRAZxKno5KHjGAfd0JLJX0vYg4K2kJ8ACZsMjd505JU3Oa5wKHkxkJkl4H7gM+iYh9ZGY6JZG0CFg0Y0ahjDMzs4EoakFfUp2kD4FuYFtEvJ9nWADvSNot6cmrOiPWAx3AG5IeBZ4AHiyh1mbg86znx5K2QjWPk/Qy0CppRb4xEfFWRDzZ2NhYQhlmZtafohb0I6IPmCOpCdgkaVZE7M8Z9s8iois5HbVN0qcRsTNnOy8kM46XgJsj4hRDJCK+BJYO1fbNzKywki5FjogeYAc56yFJX1fyZzewicxprCtIugOYlfRfdVFAP7qAm7KeT0nazMyswvQbLpLGJzMWJDUAdwGf5owZe2mxXtJYYD6wP2dMK7CGzDrJ48A4Sc+VUOsu4BZJ0ySNBh4G3izh9WZmNkyKmblMAnZI+pjMB/y2iNgCIOltSZOBG4G/lfQR8AGwNSJ+nrOdMcBDEfFZRFwEHgOO5tuhpLXAe0CLpGOSFkfEBeApMus2B4F1EXGg1DdsZmZDT5mLwWpbW1tbdHZ2lrsMM7OqIml3oe80+udfzMwsdQ4XMzNLncPFzMxS53AxM7PUOVzMzCx1DhczM0udw8XMzFLncDEzs9Q5XMzMLHUOFzMzS53DxczMUudwMTOz1DlczMwsdQ4XMzNLncPFzMxS53AxM7PUjSx3AWa5Nu/tYlXHIb7o6WVyUwPLF7Rwf2tzucsysxI4XKyibN7bxYqN++g93wdAV08vKzbuA3DAmFURnxarEVuPbGX+hvl866+/xfwN89l6ZGu5S8prVcehy8FySe/5PlZ1HCpTRWY2EJ651ICtR7ay8t2VnOk7A8Dx08dZ+e5KABZOX1jGyq72RU9vSe1mVpk8c6kBq/esvhwsl5zpO8PqPavLVFFhk5saSmo3s8rkcKkBJ06fKKm9nJYvaKFhVN0VbQ2j6li+oKVMFZnZQDhcasDEsRNLai+n+1ubef6B2TQ3NSCguamB5x+Y7cV8syrjNZcasOy2ZVesuQDU19Wz7LZlZayqsPtbmx0mZlXO4VIDLi3ar96zmhOnTzBx7ESW3bas4hbzzez64XCpEQunL3SYmNmw8ZqLmZml7roNF0nTJb0iaUO5azEzqzX9houkekkfSPpI0gFJz15jbJ2kvZK2DKYoSa9K6pa0P6e9XdIhSYclPX2tbUTEkYhYPJg6zMxsYIqZuZwF5kXErcAcoF3S7QXGLgMO5uuQNEHSDTltMwps5zWgPWdsHfAicDcwE3hE0sykb7akLTmPCUW8NzMzGwL9hktknEqejkoekTtO0hRgIfCTApu6E9gs6WvJ+CXAjwrscyfwVU7zXOBwMiM5B7wO3JeM3xcR9+Q8uvt7b5IWSVpz8uTJ/oaamVkJilpzSU53fQh0A9si4v08w/4L8O+Bi/m2ERHrgQ7gDUmPAk8AD5ZQazPwedbzY0lboZrHSXoZaJW0okBNb0XEk42NjSWUYWZm/SkqXCKiLyLmAFOAuZJmZfdLugfojojd/WznBeAM8BJwb9aMKHUR8WVELI2ImyPi+aHaj5mZXa2kq8UiogfYQc56CPAd4F5JvyRzumqepP+e+3pJdwCzgE3AMyXW2gXclPV8StJmZmYVppirxcZLakr+3gDcBXyaPSYiVkTElIiYCjwMbI+IP87ZTiuwhsw6yePAOEnPlVDrLuAWSdMkjU7282YJrzcbEpv3dvGdH2xn2tNb+c4PtrN5r/+fx6yYmcskYIekj8l8wG+LiC0Akt6WNLnIfY0BHoqIzyLiIvAYcDTfQElrgfeAFknHJC2OiAvAU2TWbQ4C6yLiQJH7tipSLTc2g9/eObOrp5fgt3fOdMBYrVPEVRd+1Zy2trbo7OwsdxnG1Tc2g8yPbK789sqK/Pma7/xgO115bmTW3NTA/356XhkqMhs+knZHRFu+vuv2G/pWnarpxmbgO2eaFeJwsYpSTTc2A98506wQh4tVlGq6sRn4zplmhThcrKIsu20Z9XX1V7RV+o3NfOdMs6v5fi5WUarxxma+c6bZ1RwuVnF8YzOz6ufTYmZmljqHi9kgVdOXPs2Gi0+LmQ1C7pc+j58+zsp3VwL41J7VNM9czAah2r706d9Bs+HimYvZIFTTlz4v/Q5a7/k+4Le/gwb4ajdLncPFbBAmjp3I8dPH87ZXmlUdhzjf0MnY3+1Ao3qI802c/dUCVnWMrthw2by3i1Udh/iip5fJTQ0sX9BSsbXalXxazGwQqulLn90X36V+0kZGjO5BghGje6iftJHui++Wu7S8/IvT1c0zF7NBqKYvfTbc+A4x4vwVbRpxnoYb3wH+vDxFXUO1zbQ8y7qSw8VskKrlS58xsqek9nK7NNNSEoi6NNM6DlBZtzPwetbVHC5mNWJSgfWhSRW4PgTVNdOqtlkWwLPbf8rPfvFjLtb9mhF9v8O/mLaEZ+b9q9S27zUXsxpRTetDUF0zrWpbz3p2+09Zf/SHxMhfI0GM/DXrj/6QZ7f/NLV9OFzMasTC6QtZ+e2VTBo7CSEmjZ1UsXf4hMIzqkqcaTXc+M7l03eX/HaWVXl+9osf5633Z7/4cWr78G2Oga9/8/firv/warnLMLMsX535il+e/CUXuXi5bQQjmNo4la/Xf72MlV2t8+8L3ya97ca8dwEuq10nOpGubo+Afzqx+HrXLf12wdsce83FzCrSpQA59ptjnLt4jtEjRjPlhikVFywAo0eM5tzFc3nbK5EYCVwo0J4OhwswffxY3vg3f1DuMsysSm098tUVvzEHmfWszGnHyvtseXb7YdYf/eEVp8bi4ige/Oaf8cy84utdt7Rwn8PFzGyQqun7TkDmqrDtXHG12PdTvlrMay5AW1tbdHYWPmdqZmZXk1RwzcVXi5mZWeocLmZmljqHi5mZpc7hYmZmqXO4mJlZ6ny1GCDpV8DRAb68ETg5DK8rdnx/467Vf62+bwD/UMT+K8FA/03KtZ/hOIaG4/i5Vn81HT8wPMdQNR0/hcZ/MyLG5x0dEX4M4gGsGY7XFTu+v3HX6u+nr7Pc/62H+t+kXPsZjmNoOI6fa/VX0/GT5r/tcOyj0j6DLj18Wmzw3hqm1xU7vr9x1+of6HupNMP1PtLaz3AcQ8Nx/JSyn0o3HO+jmo6fksf7tJgVRVJnFPiylFl/fPzUHs9crFhryl2AVTUfPzXGMxczM0udZy5mZpY6h4uZmaXO4WJmZqlzuNiASJou6RVJG8pdi1UfSfdL+rGkNyTNL3c9lj6Hi10m6VVJ3ZL257S3Szok6bCkpwEi4khELC5PpVaJSjx+NkfEEmAp8C/LUa8NLYeLZXsNaM9ukFQHvAjcDcwEHpE0c/hLsyrwGqUfP3+e9Nt1xuFil0XETuCrnOa5wOFkpnIOeB24b9iLs4pXyvGjjL8A/iYi9gx3rTb0HC7Wn2bg86znx4BmSeMkvQy0SlpRntKsCuQ9foA/Af4I+L6kpeUozIbWyHIXYNUpIr4kc77crGQR8ZfAX5a7Dhs6nrlYf7qAm7KeT0nazIrh46dGOVysP7uAWyRNkzQaeBh4s8w1WfXw8VOjHC52maS1wHtAi6RjkhZHxAXgKaADOAisi4gD5azTKpOPH8vmH640M7PUeeZiZmapc7iYmVnqHC5mZpY6h4uZmaXO4WJmZqlzuJiZWeocLmZmljqHi5mZpc7hYmZmqfv/p03hf4mdrwkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for n in 2 ** np.arange(2, 10):\n",
    "    meang = np.mean([1 / estimate(n) for _ in range(10000)])\n",
    "    plt.loglog(n, meang, 'o', color='tab:blue')\n",
    "    plt.loglog(n, 1 / mu + var / (n * mu ** 3), 'o', color='tab:green')\n",
    "\n",
    "    \n",
    "plt.axhline(1 / mu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def u():\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 435,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N():\n",
    "    return np.random.randint(1, 5, 1)[0]\n",
    "\n",
    "def p(n):\n",
    "    return 1 / 5\n",
    "\n",
    "# def g(x, n):\n",
    "#     if n == 0:\n",
    "#         return 1 / x\n",
    "    \n",
    "#     return n * (-1) ** n / (x ** (n + 1))\n",
    "\n",
    "def g(x, n):\n",
    "    if n == 0:\n",
    "        return x ** 2\n",
    "    elif n == 1:\n",
    "        return 2 * x\n",
    "    elif n == 2:\n",
    "        return 2\n",
    "    return 0\n",
    "\n",
    "def V(x, n):\n",
    "    samples = np.random.normal(mu, var, n)\n",
    "    pn = p(n)\n",
    "    nfac = math.factorial(n)\n",
    "#     print(n, nfac)\n",
    "    return g(x, n) / nfac * np.prod(samples - x) / pn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 439,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-91.32752217673789"
      ]
     },
     "execution_count": 439,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.mean(np.random.normal(mu, var, n)) ** 2\n",
    "Ns = [N() for _  in range(100)]\n",
    "np.mean([V(x, int(Ns_i)) for Ns_i in Ns])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}