delta_bulk_size.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-09-18 00:54:53.843618: W external/xla/xla/service/gpu/nvptx_compiler.cc:836] The NVIDIA driver's CUDA version is 12.2 which is older than the PTX compiler version (12.6.68). Because the driver is older than the PTX compiler version, XLA is disabling parallel compilation, which may slow down compilation. You should update your NVIDIA driver or use the NVIDIA-provided CUDA forward compatibility packages.\n"
     ]
    }
   ],
   "source": [
    "from delta_bulk import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'symmetric': True, 'canonicalize': False, 'aligned': True, 'delta_bulks': Array([[-1.1960341e-03,  4.7091188e-04,  9.5689128e-04, ...,\n",
      "        -6.1697447e-06, -1.1675453e-05,  5.8049827e-06],\n",
      "       [ 1.3619699e-03, -8.5998839e-04, -3.0566793e-04, ...,\n",
      "         2.8423116e-05,  2.6712756e-05, -6.9486164e-06],\n",
      "       [-8.3080743e-04, -4.0815197e-04, -7.6364668e-05, ...,\n",
      "         4.4501525e-05,  4.3573760e-05, -1.9222136e-05],\n",
      "       ...,\n",
      "       [ 5.5929285e-04,  8.0412795e-04,  8.8316377e-04, ...,\n",
      "        -1.1356422e-05, -9.1611510e-06,  1.1044276e-05],\n",
      "       [ 1.1526232e-03,  1.8599018e-04,  1.7479589e-04, ...,\n",
      "         5.2014344e-05, -7.7330697e-06,  3.8690858e-05],\n",
      "       [ 5.4587610e-04,  6.2748813e-04,  9.7122852e-04, ...,\n",
      "         1.1506224e-05,  3.2652792e-05, -9.9428726e-06]], dtype=float32)}\n"
     ]
    }
   ],
   "source": [
    "with open(\"delta_bulk_1726584467.5956047.pkl\", \"rb\") as f:\n",
    "    data = pickle.load(f)\n",
    "\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "db = data['delta_bulks']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "n = 16\n",
      "[]\n",
      "0.0015814677\n",
      "\n",
      "n = 16\n",
      "['aligned']\n",
      "0.0036289096\n",
      "\n",
      "n = 16\n",
      "['canonicalize', 'aligned']\n",
      "0.0040753856\n",
      "\n",
      "n = 400\n",
      "[]\n",
      "0.0013575461\n",
      "\n",
      "n = 400\n",
      "['aligned']\n",
      "0.0022507259\n",
      "\n",
      "n = 400\n",
      "['canonicalize', 'aligned']\n",
      "0.002232391\n",
      "\n",
      "n = 400\n",
      "['symmetric', 'canonicalize', 'aligned']\n",
      "0.0025880416\n",
      "\n",
      "n = 400\n",
      "['symmetric', 'aligned']\n",
      "0.0025422224\n"
     ]
    }
   ],
   "source": [
    "# open every pkl file in the directory\n",
    "import glob\n",
    "files = glob.glob(\"*.pkl\")\n",
    "for file in files:\n",
    "    print()\n",
    "    #print(file)\n",
    "    with open(file, \"rb\") as f:\n",
    "        data = pickle.load(f)\n",
    "    db = data['delta_bulks']\n",
    "    n = db.shape[0]\n",
    "    # if n > 100:\n",
    "    #     db = db[:100]\n",
    "    #     n = 100\n",
    "    print(f'n = {n}')\n",
    "    print([k for k in data.keys() if k != 'delta_bulks' and data[k]])\n",
    "\n",
    "    mean = np.mean(db, axis=0)\n",
    "    print(np.sum(mean**2) / np.mean(np.sum(db**2, axis=1)) - 1/(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bootstrap estimate of the mean\n",
    "mean = np.mean(db, axis=0)\n",
    "\n",
    "# boots = np.random.choice(db, size=(1000, db.shape[1]), replace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(2.1635191e-05, dtype=float32)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(np.mean(db, axis=0)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(0.00253596, dtype=float32)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(mean**2) / np.mean(np.sum(db**2, axis=1)) - 1/399"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(0.9974514, dtype=float32)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# covariance matrix\n",
    "cov = np.cov(db, rowvar=False)\n",
    "np.sum(np.diag(cov)) / np.mean(np.sum(db**2, axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(0.00269872, dtype=float32)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean100 = np.mean(db[:100], axis=0)\n",
    "np.sum(mean100**2)\n",
    "np.sum(mean100**2) / np.mean(np.sum(db[:100]**2, axis=1)) - 1/99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 400, 9610)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# numpy.random.Generator.choice\n",
    "\n",
    "g = np.random.default_rng(12345)\n",
    "boots = g.choice(db, size=(1000, 400), replace=True)\n",
    "\n",
    "boots.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(0.00504222, dtype=float32)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(mean**2) / np.mean(np.sum(db**2, axis=1)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(1.4895207, dtype=float32)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(means) / np.sum(mean**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(2.1635191e-05, dtype=float32)"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(mean**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiEAAAGdCAYAAADE96MUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhLUlEQVR4nO3de1TUZeLH8Q+IDKgwIMqAhkp7TE3dMkwlzS6ycoxTutK2tVrqmnaKLmpbya6Xtc0wa83No3ZZQ63UdM8p85JWmLZrqIhpoomWuloGdpNx1UDh+f3Rz9FRLAYGnwHfr3PmtHznmfk+3+ds8O4735kJMsYYAQAAXGTBticAAAAuTUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArAixPYFzVVRU6NChQ4qIiFBQUJDt6QAAgCowxujo0aNq0aKFgoOrdo4j4CLk0KFDSkhIsD0NAABQDQcPHtRll11WpbEBFyERERGSfjqIyMhIy7MBAABV4Xa7lZCQ4Pk7XhUBFyGnX4KJjIwkQgAAqGN8uZSCC1MBAIAVRAgAALCCCAEAAFYQIQAAwAoiBAAAWEGEAAAAK4gQAABgBRECAACsIEIAAIAVRAgAALCCCAEAAFYQIQAAwAoiBAAAWEGEAAAAK0JsTwBA3dNm7IpqP3b/lDQ/zgRAXcaZEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFjhU4SUl5dr/PjxSkxMVHh4uH71q1/pb3/7m4wxnjHGGE2YMEHx8fEKDw9XSkqK9uzZ4/eJAwCAui3El8HPPPOMZs+erXnz5qljx47avHmzhg0bJqfTqYcffliSNHXqVL3wwguaN2+eEhMTNX78eKWmpmrnzp0KCwurlYMALlVtxq6o9mP3T0nz40wAwHc+RcjHH3+s/v37Ky3tp19ebdq00cKFC7Vp0yZJP50FmT59usaNG6f+/ftLkubPny+Xy6W3335bd955p5+nDwAA6iqfXo657rrrlJOTo927d0uStm3bpv/85z/q16+fJGnfvn0qKipSSkqK5zFOp1Pdu3dXbm6uH6cNAADqOp/OhIwdO1Zut1vt27dXgwYNVF5ersmTJ2vQoEGSpKKiIkmSy+XyepzL5fLcd67S0lKVlpZ6fna73T4dAAAAqJt8ipDFixfrjTfe0IIFC9SxY0dt3bpVo0aNUosWLTRkyJBqTSArK0uTJk2q1mOB+qAm13XUxf0CwGk+vRzz2GOPaezYsbrzzjvVuXNn3X333Ro9erSysrIkSXFxcZKk4uJir8cVFxd77jtXZmamSkpKPLeDBw9W5zgAAEAd41OEHD9+XMHB3g9p0KCBKioqJEmJiYmKi4tTTk6O5363262NGzcqOTm50ud0OByKjIz0ugEAgPrPp5djbr31Vk2ePFmtWrVSx44d9cknn2jatGn64x//KEkKCgrSqFGj9NRTT6lt27aet+i2aNFCAwYMqI35AwCAOsqnCJkxY4bGjx+vBx54QIcPH1aLFi103333acKECZ4xjz/+uI4dO6aRI0fqyJEj6tWrl1atWsVnhAAAAC9B5uyPOw0AbrdbTqdTJSUlvDSDS8KldoEoH5IG1E/V+fvNd8cAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVhAhAADACiIEAABYQYQAAAAriBAAAGCFT98dAwA1VZOPqecj34H6hTMhAADACiIEAABYQYQAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVhAhAADACiIEAABYQYQAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVhAhAADACiIEAABYQYQAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVhAhAADACiIEAABYQYQAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVhAhAADACiIEAABYQYQAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVhAhAADACiIEAABYQYQAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVhAhAADACiIEAABYQYQAAAArQmxPAKgP2oxdYXsKAFDncCYEAABYQYQAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVvgcIV999ZUGDx6smJgYhYeHq3Pnztq8ebPnfmOMJkyYoPj4eIWHhyslJUV79uzx66QBAEDd51OE/PDDD+rZs6caNmyod999Vzt37tTf//53RUdHe8ZMnTpVL7zwgl588UVt3LhRjRs3Vmpqqn788Ue/Tx4AANRdPn1i6jPPPKOEhARlZ2d7tiUmJnr+tzFG06dP17hx49S/f39J0vz58+VyufT222/rzjvv9NO0AQBAXefTmZB33nlHXbt21e9+9zvFxsaqS5cueuWVVzz379u3T0VFRUpJSfFsczqd6t69u3Jzc/03awAAUOf5FCF79+7V7Nmz1bZtW61evVr333+/Hn74Yc2bN0+SVFRUJElyuVxej3O5XJ77zlVaWiq32+11AwAA9Z9PL8dUVFSoa9euevrppyVJXbp0UUFBgV588UUNGTKkWhPIysrSpEmTqvVYAJeWmnxR4P4paX6cCQB/8OlMSHx8vK688kqvbR06dNCBAwckSXFxcZKk4uJirzHFxcWe+86VmZmpkpISz+3gwYO+TAkAANRRPkVIz549VVhY6LVt9+7dat26taSfLlKNi4tTTk6O5363262NGzcqOTm50ud0OByKjIz0ugEAgPrPp5djRo8ereuuu05PP/207rjjDm3atEkvv/yyXn75ZUlSUFCQRo0apaeeekpt27ZVYmKixo8frxYtWmjAgAG1MX8AAFBH+RQh1157rd566y1lZmbqySefVGJioqZPn65BgwZ5xjz++OM6duyYRo4cqSNHjqhXr15atWqVwsLC/D55AABQdwUZY4ztSZzN7XbL6XSqpKSEl2ZQZ9TkgklcHFyYCtSu6vz95rtjAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACt8+pwQAKir+N4ZIPBwJgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWhNieAAAEujZjV1T7sfunpPlxJkD9wpkQAABgBRECAACsIEIAAIAVRAgAALCCCAEAAFYQIQAAwAoiBAAAWEGEAAAAK4gQAABgBRECAACsIEIAAIAVfHcM8P9q8v0gAADfcSYEAABYQYQAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVvA5IQBQi2ry+TP7p6T5cSZA4OFMCAAAsIIIAQAAVhAhAADACiIEAABYQYQAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVhAhAADACiIEAABYQYQAAAAriBAAAGAFEQIAAKwgQgAAgBVECAAAsIIIAQAAVhAhAADAihpFyJQpUxQUFKRRo0Z5tv3444/KyMhQTEyMmjRpovT0dBUXF9d0ngAAoJ6pdoTk5eXppZde0q9//Wuv7aNHj9ayZcu0ZMkSrVu3TocOHdLAgQNrPFEAAFC/VCtC/ve//2nQoEF65ZVXFB0d7dleUlKiOXPmaNq0abr55puVlJSk7Oxsffzxx9qwYYPfJg0AAOq+akVIRkaG0tLSlJKS4rU9Pz9fJ0+e9Nrevn17tWrVSrm5uZU+V2lpqdxut9cNAADUfyG+PmDRokXasmWL8vLyzruvqKhIoaGhioqK8trucrlUVFRU6fNlZWVp0qRJvk4DAOq9NmNX1Ojx+6ek+WkmQO3w6UzIwYMH9cgjj+iNN95QWFiYXyaQmZmpkpISz+3gwYN+eV4AABDYfIqQ/Px8HT58WNdcc41CQkIUEhKidevW6YUXXlBISIhcLpfKysp05MgRr8cVFxcrLi6u0ud0OByKjIz0ugEAgPrPp5dj+vTpo+3bt3ttGzZsmNq3b68nnnhCCQkJatiwoXJycpSeni5JKiws1IEDB5ScnOy/WQMAgDrPpwiJiIhQp06dvLY1btxYMTExnu3Dhw/XmDFj1LRpU0VGRuqhhx5ScnKyevTo4b9ZAwCAOs/nC1N/yfPPP6/g4GClp6ertLRUqampmjVrlr93AwAA6rggY4yxPYmzud1uOZ1OlZSUcH0ILqqavhMBCDS8OwYXU3X+fvPdMQAAwAoiBAAAWEGEAAAAK4gQAABgBRECAACsIEIAAIAVRAgAALCCCAEAAFYQIQAAwAoiBAAAWEGEAAAAK4gQAABgBRECAACsIEIAAIAVIbYnAPhTm7ErbE8BAFBFnAkBAABWECEAAMAKIgQAAFjBNSEAgPPU5Pqq/VPS/DgT1GecCQEAAFYQIQAAwAoiBAAAWEGEAAAAK4gQAABgBRECAACsIEIAAIAVfE4IAg7f/wL4B/8uIdBxJgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYIVPEZKVlaVrr71WERERio2N1YABA1RYWOg15scff1RGRoZiYmLUpEkTpaenq7i42K+TBgAAdZ9PEbJu3TplZGRow4YNev/993Xy5En17dtXx44d84wZPXq0li1bpiVLlmjdunU6dOiQBg4c6PeJAwCAui3El8GrVq3y+nnu3LmKjY1Vfn6+evfurZKSEs2ZM0cLFizQzTffLEnKzs5Whw4dtGHDBvXo0cN/MwcAAHVaja4JKSkpkSQ1bdpUkpSfn6+TJ08qJSXFM6Z9+/Zq1aqVcnNzK32O0tJSud1urxsAAKj/qh0hFRUVGjVqlHr27KlOnTpJkoqKihQaGqqoqCivsS6XS0VFRZU+T1ZWlpxOp+eWkJBQ3SkBAIA6pNoRkpGRoYKCAi1atKhGE8jMzFRJSYnndvDgwRo9HwAAqBt8uibktAcffFDLly/XRx99pMsuu8yzPS4uTmVlZTpy5IjX2ZDi4mLFxcVV+lwOh0MOh6M60wAAAHWYT2dCjDF68MEH9dZbb2nNmjVKTEz0uj8pKUkNGzZUTk6OZ1thYaEOHDig5ORk/8wYAADUCz6dCcnIyNCCBQu0dOlSRUREeK7zcDqdCg8Pl9Pp1PDhwzVmzBg1bdpUkZGReuihh5ScnMw7YwAAgBefImT27NmSpBtvvNFre3Z2toYOHSpJev755xUcHKz09HSVlpYqNTVVs2bN8stkAQBA/eFThBhjfnFMWFiYZs6cqZkzZ1Z7UgAAoP7ju2MAAIAVRAgAALCCCAEAAFYQIQAAwAoiBAAAWEGEAAAAK4gQAABgBRECAACsIEIAAIAVRAgAALDCp49tB6qqzdgVtqcAAAhwnAkBAABWECEAAMAKIgQAAFjBNSG4IK7rAADUJs6EAAAAK4gQAABgBRECAACs4JoQAEC9UJPr2PZPSfPjTFBVnAkBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFUQIAACwgggBAABWECEAAMAKIgQAAFhBhAAAACuIEAAAYAURAgAArCBCAACAFSG2J4Da1WbsCttTAIAq43fWpYUzIQAAwAoiBAAAWEGEAAAAK4gQAABgBRECAACsIEIAAIAVRAgAALCCzwkBAPgVn/WBquJMCAAAsIIIAQAAVhAhAADACq4JAQBc8mpyHcv+KWl+nMmlhTMhAADACiIEAABYQYQAAAAruCYEAIAa4HqS6uNMCAAAsIIIAQAAVhAhAADACiIEAABYQYQAAAAriBAAAGAFEQIAAKzgc0LqgJq8Bx0AgEDFmRAAAGAFEQIAAKwgQgAAgBVcEwIAQB1UH76zhjMhAADACiIEAABYQYQAAAAruCbkIuGzPgAA57rU/zZwJgQAAFhBhAAAACuIEAAAYEWtXRMyc+ZMPfvssyoqKtJVV12lGTNmqFu3brW1uyqrD++rBgCgPqiVMyFvvvmmxowZo4kTJ2rLli266qqrlJqaqsOHD9fG7gAAQB1UKxEybdo0jRgxQsOGDdOVV16pF198UY0aNdKrr75aG7sDAAB1kN9fjikrK1N+fr4yMzM924KDg5WSkqLc3NzzxpeWlqq0tNTzc0lJiSTJ7Xb7e2qSpIrS49V+bE3mVJP9AgDgT7XxN/b0cxpjqvwYv0fIt99+q/LycrlcLq/tLpdLu3btOm98VlaWJk2adN72hIQEf0+txpzTbc8AAICaq82/Z0ePHpXT6azSWOsfVpaZmakxY8Z4fq6oqND333+vmJgYBQUFeY11u91KSEjQwYMHFRkZebGnGlBYizNYizNYi5+wDmewFmewFmfUxloYY3T06FG1aNGiyo/xe4Q0a9ZMDRo0UHFxsdf24uJixcXFnTfe4XDI4XB4bYuKivrZfURGRl7y/wc6jbU4g7U4g7X4CetwBmtxBmtxhr/XoqpnQE7z+4WpoaGhSkpKUk5OjmdbRUWFcnJylJyc7O/dAQCAOqpWXo4ZM2aMhgwZoq5du6pbt26aPn26jh07pmHDhtXG7gAAQB1UKxHy+9//Xt98840mTJigoqIiXX311Vq1atV5F6v6yuFwaOLEiee9fHMpYi3OYC3OYC1+wjqcwVqcwVqcEShrEWR8eS8NAACAn/DdMQAAwAoiBAAAWEGEAAAAK4gQAABgRa1GyMyZM9WmTRuFhYWpe/fu2rRp08+OX7Jkidq3b6+wsDB17txZK1eu9LrfGKMJEyYoPj5e4eHhSklJ0Z49e857nhUrVqh79+4KDw9XdHS0BgwY4HV/Xl6e+vTpo6ioKEVHRys1NVXbtm2r8fH+nIu9FmvXrlVQUFClt7y8PM+4Tz/9VNdff73CwsKUkJCgqVOn+vfAKxGIa7F27Vr1799f8fHxaty4sa6++mq98cYb/j/4swTiOpzt888/V0RExC9+eKA/BOpaGGP03HPP6YorrpDD4VDLli01efJk/x78OQJ1LVavXq0ePXooIiJCzZs3V3p6uvbv3+/XYz+Xjb8hu3fvVv/+/dWsWTNFRkaqV69e+vDDD73GHDhwQGlpaWrUqJFiY2P12GOP6dSpU/456EoE4jps27ZNd911lxISEhQeHq4OHTroH//4h+8HZ2rJokWLTGhoqHn11VfNjh07zIgRI0xUVJQpLi6udPz69etNgwYNzNSpU83OnTvNuHHjTMOGDc327ds9Y6ZMmWKcTqd5++23zbZt28xtt91mEhMTzYkTJzxj/vWvf5no6Ggze/ZsU1hYaHbs2GHefPNNz/1Hjx41TZs2NUOHDjW7du0yBQUFJj093bhcLlNWVlZv1qK0tNR8/fXXXrd7773XJCYmmoqKCmOMMSUlJcblcplBgwaZgoICs3DhQhMeHm5eeumlWlmHQF6LyZMnm3Hjxpn169ebzz//3EyfPt0EBwebZcuWXVLrcFpZWZnp2rWr6devn3E6nbWyBqcF8lo89NBDpl27dmbp0qVm7969ZvPmzea999675NZi7969xuFwmMzMTPP555+b/Px807t3b9OlS5d6tRbGGNO2bVtzyy23mG3btpndu3ebBx54wDRq1Mh8/fXXxhhjTp06ZTp16mRSUlLMJ598YlauXGmaNWtmMjMzL6l1mDNnjnn44YfN2rVrzRdffGFee+01Ex4ebmbMmOHT8dVahHTr1s1kZGR4fi4vLzctWrQwWVlZlY6/4447TFpamte27t27m/vuu88YY0xFRYWJi4szzz77rOf+I0eOGIfDYRYuXGiMMebkyZOmZcuW5p///OcF55WXl2ckmQMHDni2ffrpp0aS2bNnj+8HWgU21uJcZWVlpnnz5ubJJ5/0bJs1a5aJjo42paWlnm1PPPGEadeune8HWUWBuhaVueWWW8ywYcOqdFy+CvR1ePzxx83gwYNNdnZ2rUdIoK7Fzp07TUhIiNm1a1e1j81XgboWS5YsMSEhIaa8vNyz7Z133jFBQUG19h9vNtbim2++MZLMRx995BnjdruNJPP+++8bY4xZuXKlCQ4ONkVFRZ4xs2fPNpGRkV6/S/0lUNehMg888IC56aabfDq+Wnk5pqysTPn5+UpJSfFsCw4OVkpKinJzcyt9TG5urtd4SUpNTfWM37dvn4qKirzGOJ1Ode/e3TNmy5Yt+uqrrxQcHKwuXbooPj5e/fr1U0FBgecx7dq1U0xMjObMmaOysjKdOHFCc+bMUYcOHdSmTRt/LYGHrbU41zvvvKPvvvvO61Nrc3Nz1bt3b4WGhnrtp7CwUD/88IPvB/sLAnktKlNSUqKmTZtW6dh8EejrsGbNGi1ZskQzZ86s1vH5IpDXYtmyZbr88su1fPlyJSYmqk2bNrr33nv1/fffV/t4f04gr0VSUpKCg4OVnZ2t8vJylZSU6LXXXlNKSooaNmxY7WO+EFtrERMTo3bt2mn+/Pk6duyYTp06pZdeekmxsbFKSkry7Kdz585eH76Zmpoqt9utHTt2+GcB/l8gr0NlqvM7s1Yi5Ntvv1V5efl5n5DqcrlUVFRU6WOKiop+dvzpf/7cmL1790qS/vrXv2rcuHFavny5oqOjdeONN3p+cURERGjt2rV6/fXXFR4eriZNmmjVqlV69913FRLi/w+QtbUW55ozZ45SU1N12WWX/eJ+zt6HPwXyWpxr8eLFysvLq5WvGgjkdfjuu+80dOhQzZ0796J8wVcgr8XevXv13//+V0uWLNH8+fM1d+5c5efn6/bbb/ftIKsokNciMTFR7733nv785z/L4XAoKipKX375pRYvXuzbQVaRrbUICgrSBx98oE8++UQREREKCwvTtGnTtGrVKkVHR//sfs7eh78E8jqc6+OPP9abb76pkSNH+nSM9erdMRUVFZKkv/zlL0pPT1dSUpKys7MVFBSkJUuWSJJOnDih4cOHq2fPntqwYYPWr1+vTp06KS0tTSdOnLA5/Vrz5ZdfavXq1Ro+fLjtqVhXlbX48MMPNWzYML3yyivq2LHjRZzdxXOhdRgxYoT+8Ic/qHfv3pZmdvFdaC0qKipUWlqq+fPn6/rrr9eNN96oOXPm6MMPP1RhYaGl2dauC61FUVGRRowYoSFDhigvL0/r1q1TaGiobr/9dpl69KHbxhhlZGQoNjZW//73v7Vp0yYNGDBAt956q77++mvb07tofF2HgoIC9e/fXxMnTlTfvn192letREizZs3UoEEDFRcXe20vLi5WXFxcpY+Ji4v72fGn//lzY+Lj4yVJV155ped+h8Ohyy+/XAcOHJAkLViwQPv371d2drauvfZa9ejRQwsWLNC+ffu0dOnS6h7yBdlai7NlZ2crJiZGt912W5X2c/Y+/CmQ1+K0devW6dZbb9Xzzz+ve+65p2oH5qNAXoc1a9boueeeU0hIiEJCQjR8+HCVlJQoJCREr776qm8HWgWBvBbx8fEKCQnRFVdc4dnWoUMHSfL8PvGnQF6LmTNnyul0aurUqerSpYt69+6t119/XTk5Odq4caNvB1oFttZizZo1Wr58uRYtWqSePXvqmmuu0axZsxQeHq558+b97H7O3oe/BPI6nLZz50716dNHI0eO1Lhx43w+xlqJkNDQUCUlJSknJ8ezraKiQjk5OUpOTq70McnJyV7jJen999/3jE9MTFRcXJzXGLfbrY0bN3rGJCUlyeFweP1XysmTJ7V//361bt1aknT8+HEFBwcrKCjIM+b0z6fPpPiTrbU4zRij7Oxs3XPPPee9dpucnKyPPvpIJ0+e9NpPu3btLnjKrSYCeS2kn96qmJaWpmeeecbnU4q+COR1yM3N1datWz23J598UhEREdq6dat++9vf1ui4KxPIa9GzZ0+dOnVKX3zxhWfb7t27Jcnz+8SfAnktTv/ePFuDBg08c/Q3W2tx/PhxSTrvWIODgz3HmZycrO3bt+vw4cNe+4mMjPT6D2B/COR1kKQdO3bopptu0pAhQ6r/1nWfLmP1waJFi4zD4TBz5841O3fuNCNHjjRRUVGeK4rvvvtuM3bsWM/49evXm5CQEPPcc8+Zzz77zEycOLHStxVFRUWZpUuXmk8//dT079//vLcVPfLII6Zly5Zm9erVZteuXWb48OEmNjbWfP/998YYYz777DPjcDjM/fffb3bu3GkKCgrM4MGDjdPpNIcOHapXa2GMMR988IGRZD777LPz5nXkyBHjcrnM3XffbQoKCsyiRYtMo0aNav0tuoG4FmvWrDGNGjUymZmZXm9V/O677y6pdTjXxXh3TKCuRXl5ubnmmmtM7969zZYtW8zmzZtN9+7dzW9+85taWonAXYucnBwTFBRkJk2aZHbv3m3y8/NNamqqad26tTl+/Hi9WYtvvvnGxMTEmIEDB5qtW7eawsJC86c//ck0bNjQbN261Rhz5i26ffv2NVu3bjWrVq0yzZs3r9W36AbiOmzfvt00b97cDB482Ot35uHDh306vlqLEGOMmTFjhmnVqpUJDQ013bp1Mxs2bPDcd8MNN5ghQ4Z4jV+8eLG54oorTGhoqOnYsaNZsWKF1/0VFRVm/PjxxuVyGYfDYfr06WMKCwu9xpSVlZlHH33UxMbGmoiICJOSkmIKCgq8xrz33numZ8+exul0mujoaHPzzTeb3Nxc/x78OWyshTHG3HXXXea666674Ly2bdtmevXqZRwOh2nZsqWZMmVKzQ60CgJxLYYMGWIknXe74YYbany8FxKI63CuixEhxgTuWnz11Vdm4MCBpkmTJsblcpmhQ4fWWpieFqhrsXDhQtOlSxfTuHFj07x5c3PbbbdVKWRrwsZa5OXlmb59+5qmTZuaiIgI06NHD7Ny5UqvMfv37zf9+vUz4eHhplmzZubRRx81J0+e9O/BnyUQ12HixImV/s5s3bq1T8cWZEw9uqoIAADUGfXq3TEAAKDuIEIAAIAVRAgAALCCCAEAAFYQIQAAwAoiBAAAWEGEAAAAK4gQAABgBRECAACsIEIAAIAVRAgAALCCCAEAAFb8H4wkSOENhOAIAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "var = np.mean(np.sum(db**2, axis=1))  # np.mean(np.sum(boots**2, axis=2), axis=1)\n",
    "means = np.sum(np.mean(boots, axis=1)**2, axis=1)\n",
    "plt.hist(means / var, bins=30)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAdU0lEQVR4nO3df3DX9X3A8Vf4kaBCEn6UBNYAsaWltYqKBaPtuqu5gkdbXLlVO26HztN2pu2Uugo7gWl/8EOnnAxh5SrqXZXp7rTtmNxcnGy1AS2yzlVG6YYDaxNaXRLEEih57w/P7zWC1iTfkDfh8bj7nsnn+/l8vu/P+8I3Tz/5fr+fkpRSCgCADA3q7wEAALwVoQIAZEuoAADZEioAQLaECgCQLaECAGRLqAAA2RIqAEC2hvT3AN6ss7MzXnrppRgxYkSUlJT093AAgHcgpRQHDhyI8ePHx6BBxTsPkl2ovPTSS1FTU9PfwwAAemDfvn3x7ne/u2j7yy5URowYERGvH2h5eXk/jwYAeCfa29ujpqam8Hu8WLILlTf+3FNeXi5UAOAkU+yXbXgxLQCQLaECAGRLqAAA2RIqAEC2hAoAkC2hAgBkS6gAANkSKgBAtoQKAJAtoQIAZEuoAADZEioAQLaECgCQLaECAGRrSH8PADi1TFq4qcfbvrB8dhFHApwMnFEBALIlVACAbAkVACBbQgUAyJZQAQCyJVQAgGwJFQAgW0IFAMiWUAEAsiVUAIBsCRUAIFtCBQDIllABALIlVACAbAkVACBbQgUAyJZQAQCyJVQAgGwJFQAgW0IFAMiWUAEAsiVUAIBsCRUAIFtCBQDIllABALIlVACAbAkVACBbQgUAyJZQAQCyJVQAgGwN6e8BALxTkxZu6vG2LyyfXcSRACeKMyoAQLaECgCQLaECAGRLqAAA2RIqAEC2hAoAkC2hAgBkq1uhcvTo0Vi8eHHU1tbGaaedFu95z3via1/7WqSUCuuklGLJkiUxbty4OO2006K+vj52795d9IEDAANft0JlxYoVsXbt2vibv/mb2LlzZ6xYsSJWrlwZq1evLqyzcuXKuOuuu2LdunWxbdu2OOOMM2LmzJlx6NChog8eABjYuvXJtD/84Q9jzpw5MXv265/wOGnSpHjwwQfj6aefjojXz6asWrUqbr755pgzZ05ERNx///1RVVUVjz76aFxxxRVFHj4AMJB164zKRRddFI2NjfHTn/40IiJ+/OMfxw9+8IO49NJLIyJiz5490dzcHPX19YVtKioqYsaMGdHU1HTcfXZ0dER7e3uXGwBARDfPqCxcuDDa29tjypQpMXjw4Dh69Gh84xvfiHnz5kVERHNzc0REVFVVddmuqqqqcN+bLVu2LG655ZaejB0AGOC6dUbloYceiu985zvxwAMPxLPPPhv33Xdf3H777XHffff1eACLFi2Ktra2wm3fvn093hcAMLB064zKX/zFX8TChQsLrzU5++yz43//939j2bJlMX/+/Kiuro6IiJaWlhg3blxhu5aWljj33HOPu8+ysrIoKyvr4fABgIGsW2dUXnvttRg0qOsmgwcPjs7OzoiIqK2tjerq6mhsbCzc397eHtu2bYu6uroiDBcAOJV064zKpz71qfjGN74REyZMiLPOOit27NgRd9xxR/zpn/5pRESUlJTE9ddfH1//+tdj8uTJUVtbG4sXL47x48fHZZdd1hfjBwAGsG6FyurVq2Px4sVx3XXXxf79+2P8+PHx+c9/PpYsWVJY56tf/WocPHgwrr322mhtbY2PfOQjsXnz5hg2bFjRBw8ADGwl6bc/VjYD7e3tUVFREW1tbVFeXt7fwwGKbNLCTf3yuC8sn90vjwunir76/e1aPwBAtoQKAJAtoQIAZEuoAADZEioAQLaECgCQLaECAGRLqAAA2RIqAEC2hAoAkC2hAgBkS6gAANnq1tWTgYGjNxcHdIE/4ERxRgUAyJZQAQCyJVQAgGwJFQAgW0IFAMiWUAEAsiVUAIBsCRUAIFtCBQDIllABALIlVACAbLnWD8Dv4LpI0H+cUQEAsiVUAIBsCRUAIFtCBQDIllABALIlVACAbAkVACBbQgUAyJZQAQCyJVQAgGwJFQAgW0IFAMiWUAEAsiVUAIBsCRUAIFtCBQDIllABALIlVACAbAkVACBbQgUAyJZQAQCyJVQAgGwN6e8BAJwIkxZu6u8hAD3gjAoAkC2hAgBkS6gAANkSKgBAtoQKAJAtoQIAZEuoAADZEioAQLaECgCQLaECAGRLqAAA2RIqAEC2hAoAkC2hAgBkS6gAANkSKgBAtoQKAJAtoQIAZGtIdzf4+c9/HjfddFM89thj8dprr8V73/ve2LBhQ1xwwQUREZFSiqVLl8b69eujtbU1Lr744li7dm1Mnjy56IOHU92khZv6ewgAfapbZ1T+7//+Ly6++OIYOnRoPPbYY/H888/HX//1X8fIkSML66xcuTLuuuuuWLduXWzbti3OOOOMmDlzZhw6dKjogwcABrZunVFZsWJF1NTUxIYNGwrLamtrC1+nlGLVqlVx8803x5w5cyIi4v7774+qqqp49NFH44orrijSsAGAU0G3zqh873vfiwsuuCD+6I/+KMaOHRvnnXderF+/vnD/nj17orm5Oerr6wvLKioqYsaMGdHU1HTcfXZ0dER7e3uXGwBARDfPqPzP//xPrF27NhYsWBB/+Zd/Gc8880x8+ctfjtLS0pg/f340NzdHRERVVVWX7aqqqgr3vdmyZcvilltu6eHw4eTndSYAb61bZ1Q6Ozvj/PPPj29+85tx3nnnxbXXXhvXXHNNrFu3rscDWLRoUbS1tRVu+/bt6/G+AICBpVuhMm7cuPjgBz/YZdkHPvCB2Lt3b0REVFdXR0RES0tLl3VaWloK971ZWVlZlJeXd7kBAER0M1Quvvji2LVrV5dlP/3pT2PixIkR8foLa6urq6OxsbFwf3t7e2zbti3q6uqKMFwA4FTSrdeo3HDDDXHRRRfFN7/5zfjsZz8bTz/9dHzrW9+Kb33rWxERUVJSEtdff318/etfj8mTJ0dtbW0sXrw4xo8fH5dddllfjB8AGMC6FSof/vCH45FHHolFixbFrbfeGrW1tbFq1aqYN29eYZ2vfvWrcfDgwbj22mujtbU1PvKRj8TmzZtj2LBhRR88ADCwdfuTaT/5yU/GJz/5ybe8v6SkJG699da49dZbezUwAADX+gEAsiVUAIBsCRUAIFtCBQDIllABALIlVACAbAkVACBbQgUAyJZQAQCyJVQAgGwJFQAgW0IFAMiWUAEAsiVUAIBsCRUAIFtCBQDI1pD+HgBw8pm0cFN/D+Gk0Zu5emH57JP2saFYnFEBALIlVACAbAkVACBbQgUAyJZQAQCyJVQAgGwJFQAgW0IFAMiWUAEAsiVUAIBsCRUAIFuu9QOQKddUAmdUAICMCRUAIFtCBQDIllABALIlVACAbAkVACBbQgUAyJZQAQCyJVQAgGwJFQAgW0IFAMiWUAEAsiVUAIBsCRUAIFtCBQDIllABALIlVACAbAkVACBbQgUAyJZQAQCyJVQAgGwJFQAgW0P6ewAwEExauKm/hwAwIDmjAgBkS6gAANkSKgBAtoQKAJAtoQIAZEuoAADZEioAQLaECgCQLaECAGRLqAAA2RIqAEC2hAoAkC2hAgBkS6gAANnqVagsX748SkpK4vrrry8sO3ToUDQ0NMTo0aNj+PDhMXfu3GhpaentOAGAU1CPQ+WZZ56Jv/3bv41zzjmny/Ibbrghvv/978fDDz8cW7ZsiZdeeik+85nP9HqgAMCpp0eh8uqrr8a8efNi/fr1MXLkyMLytra2+Pa3vx133HFHfPzjH49p06bFhg0b4oc//GFs3bq1aIMGAE4NPQqVhoaGmD17dtTX13dZvn379jhy5EiX5VOmTIkJEyZEU1NT70YKAJxyhnR3g40bN8azzz4bzzzzzDH3NTc3R2lpaVRWVnZZXlVVFc3NzcfdX0dHR3R0dBS+b29v7+6QAIABqltnVPbt2xd//ud/Ht/5zndi2LBhRRnAsmXLoqKionCrqakpyn4BgJNft0Jl+/btsX///jj//PNjyJAhMWTIkNiyZUvcddddMWTIkKiqqorDhw9Ha2trl+1aWlqiurr6uPtctGhRtLW1FW779u3r8cEAAANLt/70c8kll8Rzzz3XZdlVV10VU6ZMiZtuuilqampi6NCh0djYGHPnzo2IiF27dsXevXujrq7uuPssKyuLsrKyHg4fABjIuhUqI0aMiA996ENdlp1xxhkxevTowvKrr746FixYEKNGjYry8vL40pe+FHV1dXHhhRcWb9QAwCmh2y+m/V3uvPPOGDRoUMydOzc6Ojpi5syZcffddxf7YQCAU0BJSin19yB+W3t7e1RUVERbW1uUl5f393DgHZm0cFN/DwGK6oXls/t7CJxk+ur3t2v9AADZEioAQLaECgCQLaECAGRLqAAA2Sr625MBOPn15p1s3jFEMTmjAgBkS6gAANkSKgBAtoQKAJAtoQIAZEuoAADZEioAQLaECgCQLaECAGRLqAAA2RIqAEC2hAoAkC2hAgBkS6gAANkSKgBAtoQKAJAtoQIAZEuoAADZEioAQLaECgCQLaECAGRLqAAA2RrS3wOAXExauKm/hwDAmzijAgBkS6gAANkSKgBAtoQKAJAtoQIAZEuoAADZEioAQLaECgCQLaECAGRLqAAA2RIqAEC2XOuHAcX1egAGFmdUAIBsCRUAIFtCBQDIllABALIlVACAbAkVACBbQgUAyJZQAQCyJVQAgGwJFQAgW0IFAMiWa/2QHdfrAeANzqgAANkSKgBAtoQKAJAtoQIAZEuoAADZEioAQLa8PRmAourNRwy8sHx2EUfCQOCMCgCQLaECAGRLqAAA2RIqAEC2hAoAkC2hAgBkS6gAANnqVqgsW7YsPvzhD8eIESNi7Nixcdlll8WuXbu6rHPo0KFoaGiI0aNHx/Dhw2Pu3LnR0tJS1EEDAKeGboXKli1boqGhIbZu3RqPP/54HDlyJD7xiU/EwYMHC+vccMMN8f3vfz8efvjh2LJlS7z00kvxmc98pugDBwAGvm59Mu3mzZu7fH/vvffG2LFjY/v27fH7v//70dbWFt/+9rfjgQceiI9//OMREbFhw4b4wAc+EFu3bo0LL7yweCMHAAa8Xr1Gpa2tLSIiRo0aFRER27dvjyNHjkR9fX1hnSlTpsSECROiqanpuPvo6OiI9vb2LjcAgIhehEpnZ2dcf/31cfHFF8eHPvShiIhobm6O0tLSqKys7LJuVVVVNDc3H3c/y5Yti4qKisKtpqamp0MCAAaYHodKQ0ND/Od//mds3LixVwNYtGhRtLW1FW779u3r1f4AgIGjR1dP/uIXvxj/8A//EP/6r/8a7373uwvLq6ur4/Dhw9Ha2trlrEpLS0tUV1cfd19lZWVRVlbWk2EAAANct86opJTii1/8YjzyyCPxxBNPRG1tbZf7p02bFkOHDo3GxsbCsl27dsXevXujrq6uOCMGAE4Z3Tqj0tDQEA888EB897vfjREjRhRed1JRURGnnXZaVFRUxNVXXx0LFiyIUaNGRXl5eXzpS1+Kuro67/gBALqtW6Gydu3aiIj4gz/4gy7LN2zYEFdeeWVERNx5550xaNCgmDt3bnR0dMTMmTPj7rvvLspgAYBTS7dCJaX0O9cZNmxYrFmzJtasWdPjQQEARLjWDwCQMaECAGRLqAAA2RIqAEC2hAoAkC2hAgBkS6gAANkSKgBAtnp0UUIA6AuTFm7q8bYvLJ9dxJGQC2dUAIBsCRUAIFtCBQDIllABALIlVACAbAkVACBbQgUAyJbPUaFP9OazEADgDc6oAADZEioAQLaECgCQLaECAGRLqAAA2RIqAEC2hAoAkC2fo8Jb8lkoAPQ3Z1QAgGwJFQAgW0IFAMiW16gMcF5nApwqevN898Ly2UUcCcXkjAoAkC2hAgBkS6gAANkSKgBAtoQKAJAtoQIAZEuoAADZEioAQLaECgCQLaECAGRLqAAA2RIqAEC2hAoAkC2hAgBka0h/DwAA+tukhZt6vO0Ly2cXcSS8mTMqAEC2hAoAkC2hAgBkS6gAANkSKgBAtoQKAJAtoQIAZMvnqABAL/gMlr7ljAoAkC2hAgBkS6gAANnyGpUTpDd/wwSAU5UzKgBAtoQKAJAtoQIAZMtrVACgn/gMlt/NGRUAIFtCBQDIllABALIlVACAbAkVACBbQgUAyFafvT15zZo1cdttt0Vzc3NMnTo1Vq9eHdOnT++rh3vHfJQ9AAPBqfLW5j45o/J3f/d3sWDBgli6dGk8++yzMXXq1Jg5c2bs37+/Lx4OABig+iRU7rjjjrjmmmviqquuig9+8IOxbt26OP300+Oee+7pi4cDAAaoov/p5/Dhw7F9+/ZYtGhRYdmgQYOivr4+mpqajlm/o6MjOjo6Ct+3tbVFRER7e3uxhxYREZ0dr/XJfgHgZNEXv2Pf2GdKqaj7LXqo/OpXv4qjR49GVVVVl+VVVVXxX//1X8esv2zZsrjllluOWV5TU1PsoQEAEVGxqu/2feDAgaioqCja/vr9Wj+LFi2KBQsWFL7v7OyMV155JUaPHh0lJSVFeYz29vaoqamJffv2RXl5eVH2eaowdz1j3nrO3PWMees5c9czb563lFIcOHAgxo8fX9THKXqojBkzJgYPHhwtLS1dlre0tER1dfUx65eVlUVZWVmXZZWVlcUeVkRElJeX+yHsIXPXM+at58xdz5i3njN3PfPb81bMMylvKPqLaUtLS2PatGnR2NhYWNbZ2RmNjY1RV1dX7IcDAAawPvnTz4IFC2L+/PlxwQUXxPTp02PVqlVx8ODBuOqqq/ri4QCAAapPQuXyyy+PX/7yl7FkyZJobm6Oc889NzZv3nzMC2xPlLKysli6dOkxf2LidzN3PWPees7c9Yx56zlz1zMnat5KUrHfRwQAUCSu9QMAZEuoAADZEioAQLaECgCQrQETKq+88krMmzcvysvLo7KyMq6++up49dVX33abQ4cORUNDQ4wePTqGDx8ec+fOPeaD6kpKSo65bdy4sS8PpU+tWbMmJk2aFMOGDYsZM2bE008//bbrP/zwwzFlypQYNmxYnH322fGP//iPXe5PKcWSJUti3Lhxcdppp0V9fX3s3r27Lw+h3xR77q688spjfrZmzZrVl4fQL7ozbz/5yU9i7ty5MWnSpCgpKYlVq1b1ep8ns2LP3V/91V8d8zM3ZcqUPjyC/tGdeVu/fn189KMfjZEjR8bIkSOjvr7+mPU9zx3fO5m7ojzPpQFi1qxZaerUqWnr1q3p3/7t39J73/ve9LnPfe5tt/nCF76QampqUmNjY/rRj36ULrzwwnTRRRd1WSci0oYNG9IvfvGLwu3Xv/51Xx5Kn9m4cWMqLS1N99xzT/rJT36SrrnmmlRZWZlaWlqOu/5TTz2VBg8enFauXJmef/75dPPNN6ehQ4em5557rrDO8uXLU0VFRXr00UfTj3/84/TpT3861dbWnrRz9Fb6Yu7mz5+fZs2a1eVn65VXXjlRh3RCdHfenn766XTjjTemBx98MFVXV6c777yz1/s8WfXF3C1dujSdddZZXX7mfvnLX/bxkZxY3Z23P/7jP05r1qxJO3bsSDt37kxXXnllqqioSC+++GJhHc9zPZ+7YjzPDYhQef7551NEpGeeeaaw7LHHHkslJSXp5z//+XG3aW1tTUOHDk0PP/xwYdnOnTtTRKSmpqbCsohIjzzySJ+N/USaPn16amhoKHx/9OjRNH78+LRs2bLjrv/Zz342zZ49u8uyGTNmpM9//vMppZQ6OztTdXV1uu222wr3t7a2prKysvTggw/2wRH0n2LPXUqv/wOeM2dOn4w3F92dt982ceLE4/6y7c0+TyZ9MXdLly5NU6dOLeIo89Pbn4/f/OY3acSIEem+++5LKXme683cpVSc57kB8aefpqamqKysjAsuuKCwrL6+PgYNGhTbtm077jbbt2+PI0eORH19fWHZlClTYsKECdHU1NRl3YaGhhgzZkxMnz497rnnnqJfwvpEOHz4cGzfvr3L8Q4aNCjq6+uPOd43NDU1dVk/ImLmzJmF9ffs2RPNzc1d1qmoqIgZM2a85T5PRn0xd2948sknY+zYsfH+978//uzP/ixefvnl4h9AP+nJvPXHPnPUl8e5e/fuGD9+fJx55pkxb9682Lt3b2+Hm41izNtrr70WR44ciVGjRkWE57nezN0bevs8NyBCpbm5OcaOHdtl2ZAhQ2LUqFHR3Nz8ltuUlpYecwHEqqqqLtvceuut8dBDD8Xjjz8ec+fOjeuuuy5Wr15d9GPoa7/61a/i6NGjx3w68JuP97c1Nze/7fpv/Lc7+zwZ9cXcRUTMmjUr7r///mhsbIwVK1bEli1b4tJLL42jR48W/yD6QU/mrT/2maO+Os4ZM2bEvffeG5s3b461a9fGnj174qMf/WgcOHCgt0POQjHm7aabborx48cXfmF7nuv53EUU53muTz5Cv1gWLlwYK1aseNt1du7c2adjWLx4ceHr8847Lw4ePBi33XZbfPnLX+7Tx2Xgu+KKKwpfn3322XHOOefEe97znnjyySfjkksu6ceRMVBdeumlha/POeecmDFjRkycODEeeuihuPrqq/txZHlYvnx5bNy4MZ588skYNmxYfw/npPJWc1eM57msz6h85StfiZ07d77t7cwzz4zq6urYv39/l21/85vfxCuvvBLV1dXH3Xd1dXUcPnw4WltbuyxvaWl5y20iXv8/khdffDE6Ojp6fXwn0pgxY2Lw4MHHvKvp7Y63urr6bdd/47/d2efJqC/m7njOPPPMGDNmTPzsZz/r/aAz0JN564995uhEHWdlZWW8733v8zMXEbfffnssX748/umf/inOOeecwnLPcz2fu+PpyfNc1qHyrne9K6ZMmfK2t9LS0qirq4vW1tbYvn17YdsnnngiOjs7Y8aMGcfd97Rp02Lo0KHR2NhYWLZr167Yu3dv1NXVveWY/v3f/z1Gjhx50l28qrS0NKZNm9bleDs7O6OxsfEtj7eurq7L+hERjz/+eGH92traqK6u7rJOe3t7bNu27W3n8GTTF3N3PC+++GK8/PLLMW7cuOIMvJ/1ZN76Y585OlHH+eqrr8Z///d/n/I/cytXroyvfe1rsXnz5i6vdYzwPNebuTueHj3P9eqluBmZNWtWOu+889K2bdvSD37wgzR58uQub09+8cUX0/vf//60bdu2wrIvfOELacKECemJJ55IP/rRj1JdXV2qq6sr3P+9730vrV+/Pj333HNp9+7d6e67706nn356WrJkyQk9tmLZuHFjKisrS/fee296/vnn07XXXpsqKytTc3NzSimlP/mTP0kLFy4srP/UU0+lIUOGpNtvvz3t3LkzLV269LhvT66srEzf/e5303/8x3+kOXPmDNi37RVz7g4cOJBuvPHG1NTUlPbs2ZP++Z//OZ1//vlp8uTJ6dChQ/1yjH2hu/PW0dGRduzYkXbs2JHGjRuXbrzxxrRjx460e/fud7zPgaIv5u4rX/lKevLJJ9OePXvSU089lerr69OYMWPS/v37T/jx9ZXuztvy5ctTaWlp+vu///sub6E9cOBAl3U8z3V/7or1PDdgQuXll19On/vc59Lw4cNTeXl5uuqqq7r8oO3ZsydFRPqXf/mXwrJf//rX6brrrksjR45Mp59+evrDP/zD9Itf/KJw/2OPPZbOPffcNHz48HTGGWekqVOnpnXr1qWjR4+eyEMrqtWrV6cJEyak0tLSNH369LR169bCfR/72MfS/Pnzu6z/0EMPpfe9732ptLQ0nXXWWWnTpk1d7u/s7EyLFy9OVVVVqaysLF1yySVp165dJ+JQTrhizt1rr72WPvGJT6R3vetdaejQoWnixInpmmuuGXC/bFPq3ry98e/0zbePfexj73ifA0mx5+7yyy9P48aNS6Wlpen3fu/30uWXX55+9rOfncAjOjG6M28TJ0487rwtXbq0sI7nudd1d+6K9TxXktJJ+F5bAOCUkPVrVACAU5tQAQCyJVQAgGwJFQAgW0IFAMiWUAEAsiVUAIBsCRUAIFtCBQDIllABALIlVACAbAkVACBb/w/XmFdQy1hWbAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "var = np.mean((db[:, 0]**2))  # np.mean(np.sum(boots**2, axis=2), axis=1)\n",
    "means = (np.mean(boots[:, :, 0], axis=1))\n",
    "plt.hist(means / np.sqrt(var), bins=30)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 400 * np.sum(mean**2) - np.trace(cov)\n",
    "T = T / 400"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "VarT = (2/400) * np.sum(cov**2)\n",
    "Z = T / np.sqrt(VarT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(0.6771072, dtype=float32)"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Array(0.00437421, dtype=float32),\n",
       " np.float64(1.0433353106782822e-07),\n",
       " Array(13.542145, dtype=float32))"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T, VarT, Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "jax_py311",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}