zenn/20230924_imbalanced_ab_test/不均等に割り振ったABテストのシミュレーション.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyMRbafZRdFGd6owh89IZ9iR",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/tatamiya/blog_artifacts/blob/main/zenn/20230924_imbalanced_ab_test/%E4%B8%8D%E5%9D%87%E7%AD%89%E3%81%AB%E5%89%B2%E3%82%8A%E6%8C%AF%E3%81%A3%E3%81%9FAB%E3%83%86%E3%82%B9%E3%83%88%E3%81%AE%E3%82%B7%E3%83%9F%E3%83%A5%E3%83%AC%E3%83%BC%E3%82%B7%E3%83%A7%E3%83%B3.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "2q345y7UCb4V"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "from scipy.stats import bernoulli, norm, skew\n",
        "\n",
        "import matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "p = 0.001\n",
        "sample_size_a = 1000\n",
        "sample_size_b = 100\n",
        "num_iteration = 10000\n",
        "\n",
        "rs = np.random.RandomState(11)\n",
        "\n",
        "x_a = rs.binomial(sample_size_a, p, size=num_iteration)\n",
        "x_b = rs.binomial(sample_size_b, p, size=num_iteration)"
      ],
      "metadata": {
        "id": "0bt16P8PCjA2"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 分布の skewness と経験的に必要なサンプルサイズ\n",
        "skewness = bernoulli.stats(p, moments=\"s\")\n",
        "skewness, skewness ** 2 * 355"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "4NxWw9xHehGq",
        "outputId": "f43d4dea-307e-4202-de3b-e6425e7bd4c9"
      },
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(31.575322658699797, 353935.35535535525)"
            ]
          },
          "metadata": {},
          "execution_count": 3
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# 2標本の母比率の差の z 検定\n",
        "p_hat_a = x_a / sample_size_a\n",
        "p_hat_b = x_b / sample_size_b\n",
        "\n",
        "sum_of_sample_size = sample_size_a + sample_size_b\n",
        "\n",
        "p_hat_pooled = (x_a + x_b) / sum_of_sample_size\n",
        "\n",
        "z_values = (p_hat_a - p_hat_b) / np.sqrt(p_hat_pooled * (1 - p_hat_pooled) * (1/sample_size_a + 1/sample_size_b))\n",
        "\n",
        "z_alpha = norm.ppf(0.975)\n",
        "(z_values< -z_alpha).sum() / num_iteration, (z_values>z_alpha).sum() / num_iteration, (np.abs(z_values)>z_alpha).sum() / num_iteration"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "TIK4p4pADZlA",
        "outputId": "7e1801a5-76c4-42d7-a4ed-53cfae0fa5d5"
      },
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "<ipython-input-4-72c49d04bb60>:9: RuntimeWarning: invalid value encountered in divide\n",
            "  z_values = (p_hat_a - p_hat_b) / np.sqrt(p_hat_pooled * (1 - p_hat_pooled) * (1/sample_size_a + 1/sample_size_b))\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(0.0713, 0.0, 0.0713)"
            ]
          },
          "metadata": {},
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# 参考： 標本 CV 率の分布\n",
        "fig, ax = plt.subplots()\n",
        "\n",
        "ax.hist(p_hat_a, density=True, alpha=0.7, label=\"pattern A\")\n",
        "ax.hist(p_hat_b, density=True, alpha=0.7, label=\"pattern B\")\n",
        "\n",
        "ax.legend()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 447
        },
        "id": "ZLoHXpmNJj4M",
        "outputId": "82756e00-c5c7-4289-ec0b-f6186e6d9d58"
      },
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7f50f6f4a230>"
            ]
          },
          "metadata": {},
          "execution_count": 5
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtTUlEQVR4nO3de3QUVYLH8V/n0cEQumOApJMxAVxRiIIiSGh3VkUzZDA+RnEYkeE1HFA3uGLOIOQMooIjiA9EBmVXQRgHBuXo6IqKg0FwF0Ig0TA8JCILhh3SCYpJCJAHSe0fntTaJmqadOibzPdzTh27q27duvemqf5ZXQ+HZVmWAAAADBIW6gYAAAB8FwEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGCciFA34Gw0Njbq6NGj6tatmxwOR6ibAwAAWsGyLJ04cUJJSUkKC/vhYyQdMqAcPXpUycnJoW4GAAA4C0eOHNEFF1zwg2U6ZEDp1q2bpG866HK5QtwaAADQGlVVVUpOTra/x39IhwwoTT/ruFwuAgoAAB1Ma07P4CRZAABgHAIKAAAwDgEFAAAYp0OegwIA6Jwsy9KZM2fU0NAQ6qbgLISHhysiIiIotwAhoAAAjFBXV6fS0lKdOnUq1E1BG0RHRysxMVFOp7NN9RBQAAAh19jYqEOHDik8PFxJSUlyOp3ciLODsSxLdXV1OnbsmA4dOqS+ffv+6M3YfggBBQAQcnV1dWpsbFRycrKio6ND3RycpfPOO0+RkZH64osvVFdXpy5dupx1XZwkCwAwRlv+jxtmCNbfkE8CAAAwDgEFAAAYh3NQAADGmrxy5znd3vKJV53T7UnS5s2bNXz4cH399deKjY2151933XW64oor9Oyzz57zNn3X3XffrZdeeklr167VL3/5y3OyTY6gAADQidXV1bVp/VOnTmnt2rV68MEHtWLFiiC16scRUAAAOEvXXXedpk2bpmnTpsntdqtHjx566KGHZFmWXeaVV17RkCFD1K1bN3k8Ht11110qLy+XJB0+fFjDhw+XJJ1//vlyOByaOHGiJk6cqC1btmjx4sVyOBxyOBw6fPiwJGnPnj0aOXKkYmJilJCQoHHjxunLL79s1qbp06erR48eysjI0ObNm+VwOJSbm6shQ4YoOjpaV199tYqLi3+0j+vWrVNqaqpmzZqljz76SEeOHAniCH4/AgoAAG2watUqRUREaMeOHVq8eLGeeeYZvfTSS/by+vp6zZs3T7t27dKbb76pw4cPa+LEiZKk5ORkvf7665Kk4uJilZaWavHixVq8eLG8Xq+mTJmi0tJSlZaWKjk5WRUVFbr++us1aNAgFRQUaMOGDSorK9Po0aObtcnpdGrr1q1atmyZPf93v/udnn76aRUUFCgiIkK/+c1vfrR/y5cv169//Wu53W6NHDlSK1eubPugtQLnoLSgtb95huK3SgCAWZKTk7Vo0SI5HA5dcskl2r17txYtWqQpU6ZIkl8IuPDCC/Xcc8/pqquuUnV1tWJiYhQXFydJio+P9zsHxel0Kjo6Wh6Px573hz/8QYMGDdLjjz9uz1uxYoWSk5P12Wef6eKLL5Yk9e3bVwsXLrTLlJaWSpJ+//vf69prr5UkzZo1S5mZmaqpqfne+5UcOHBA27dv1xtvvCFJ+vWvf63s7GzNnj273W+kR0Bpg0BO3iLMAEDnNGzYML8va6/Xq6effloNDQ0KDw9XYWGhHnnkEe3atUtff/21GhsbJUklJSVKTU0NaFu7du3Shx9+qJiYmGbLDh48aAeUwYMHt7j+wIED7deJiYmSpPLycqWkpLRYfsWKFcrIyFCPHj0kSTfeeKMmT56sTZs26YYbbgio7YEioAAA0E5OnjypjIwMZWRkaPXq1erZs6dKSkqUkZFxVievVldX6+abb9YTTzzRbFlT4JCkrl27trh+ZGSk/bopVDUFpu9qaGjQqlWr5PP5FBER4Td/xYoVBBQAAEyWn5/v93779u3q27evwsPDtX//fn311VdasGCBkpOTJUkFBQV+5ZseqvfdJzg7nc5m86688kq9/vrr6t27t19oaA/vvvuuTpw4oU8++UTh4eH2/D179mjSpEmqqKjw+0kq2DhJFgCANigpKVF2draKi4v15z//WUuWLNH9998vSUpJSZHT6dSSJUv0P//zP/rP//xPzZs3z2/9Xr16yeFwaP369Tp27Jiqq6slSb1791Z+fr4OHz6sL7/8Uo2NjcrKytLx48c1ZswY7dy5UwcPHtT777+vSZMmNQszbbV8+XJlZmbq8ssv12WXXWZPo0ePVmxsrFavXh3U7X0XAQUAgDYYP368Tp8+raFDhyorK0v333+/pk6dKknq2bOnVq5caV+qu2DBAj311FN+6//kJz/Ro48+qlmzZikhIUHTpk2TJP32t79VeHi4UlNT7Z+GkpKStHXrVjU0NGjEiBEaMGCApk+frtjY2KA+x6isrEzvvPOORo0a1WxZWFiYbrvtNi1fvjxo22uJw/r2xdodRFVVldxutyorK+VyuYJef3vcuZCTZAHg+9XU1OjQoUPq06dPm56Ae66ZdLdXU/zQ3zKQ72+OoAAAAOMQUAAAgHG4igcAgLO0efPmUDeh0+IICgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4wQUUB555BE5HA6/qV+/fvbympoaZWVlqXv37oqJidGoUaNUVlbmV0dJSYkyMzMVHR2t+Ph4zZgxQ2fOnAlObwAAQKcQ8GXGl156qT744IP/r+BbDyt64IEH9M4772jdunVyu92aNm2abr/9dm3dulXSNw9CyszMlMfj0bZt21RaWqrx48crMjJSjz/+eBC6AwAAOoOAA0pERIQ8Hk+z+ZWVlVq+fLnWrFmj66+/XpL08ssvq3///tq+fbuGDRumv/71r9q3b58++OADJSQk6IorrtC8efM0c+ZMPfLII/YTHQEAkCSt+dW53d5dr57b7embe6kMHz5cX3/9td/TgU24jb7D4bBfh4eHKykpSXfccYfmz5+vqKiodt12wOegHDhwQElJSbrwwgs1duxYlZSUSJIKCwtVX1+v9PR0u2y/fv2UkpKivLw8SVJeXp4GDBighIQEu0xGRoaqqqq0d+/etvYFAAB8R11dXZvWf/nll1VaWqpDhw7p+eef1yuvvKLHHnssSK37fgEFlLS0NK1cuVIbNmzQCy+8oEOHDulf/uVfdOLECfl8PjmdTr/0J0kJCQny+XySJJ/P5xdOmpY3Lfs+tbW1qqqq8psAAAi16667TtOmTdO0adPkdrvVo0cPPfTQQ/r2c3hfeeUVDRkyRN26dZPH49Fdd92l8vJySdLhw4c1fPhwSdL5558vh8OhiRMnauLEidqyZYsWL15sn/N5+PBhSdKePXs0cuRIxcTEKCEhQePGjdOXX37ZrE3Tp09Xjx49lJGRoc2bN8vhcCg3N1dDhgxRdHS0rr76ahUXF/9oH2NjY+XxeJScnKybbrpJt956qz7++OMgjmLLAgooI0eO1C9/+UsNHDhQGRkZevfdd1VRUaHXXnutvdonSZo/f77cbrc9JScnt+v2AABorVWrVikiIkI7duzQ4sWL9cwzz+ill16yl9fX12vevHnatWuX3nzzTR0+fFgTJ06UJCUnJ+v111+XJBUXF6u0tFSLFy/W4sWL5fV6NWXKFJWWlqq0tFTJycmqqKjQ9ddfr0GDBqmgoEAbNmxQWVmZRo8e3axNTqdTW7du1bJly+z5v/vd7/T000+roKBAERER+s1vfhNQXz/77DNt2rRJaWlpZzlardemZ/HExsbq4osv1ueff66f/exnqqurU0VFhd9RlLKyMvucFY/Hox07dvjV0XSVT0vntTTJyclRdna2/b6qqoqQAgAwQnJyshYtWiSHw6FLLrlEu3fv1qJFizRlyhRJ8gsBF154oZ577jldddVVqq6uVkxMjOLi4iRJ8fHxft+fTqdT0dHRft+Pf/jDHzRo0CC/C0tWrFih5ORkffbZZ7r44oslSX379tXChQvtMqWlpZKk3//+97r22mslSbNmzVJmZqZqamrUpUuX7+3fmDFjFB4erjNnzqi2tlY33XSTcnJyzna4Wq1N90Gprq7WwYMHlZiYqMGDBysyMlK5ubn28uLiYpWUlMjr9UqSvF6vdu/ebR/akqSNGzfK5XIpNTX1e7cTFRUll8vlNwEAYIJhw4b5nUzq9Xp14MABNTQ0SPrmHM2bb75ZKSkp6tatmx0Qms7hDMSuXbv04YcfKiYmxp6abvdx8OBBu9zgwYNbXH/gwIH268TEREny+05uyaJFi1RUVKRdu3Zp/fr1+uyzzzRu3LiA2x6ogI6g/Pa3v9XNN9+sXr166ejRo3r44YcVHh6uMWPGyO12a/LkycrOzlZcXJxcLpfuu+8+eb1eDRs2TJI0YsQIpaamaty4cVq4cKF8Pp9mz56trKysdj8bGACAc+3kyZPKyMhQRkaGVq9erZ49e6qkpEQZGRlndfJqdXW1br75Zj3xxBPNljUFDknq2rVri+tHRkbar5tCVWNj4w9u0+Px6KKLLpIkXXLJJTpx4oTGjBmjxx57zJ7fHgIKKP/7v/+rMWPG6KuvvlLPnj3105/+VNu3b1fPnj0lfZOywsLCNGrUKNXW1iojI0PPP/+8vX54eLjWr1+ve++9V16vV127dtWECRM0d+7c4PYKAIBzJD8/3+/99u3b1bdvX4WHh2v//v366quvtGDBAvvUhIKCAr/yTbfYaDri8u3535135ZVX6vXXX1fv3r397kN2LoWHh0uSTp8+3a7bCah3a9eu/cHlXbp00dKlS7V06dLvLdOrVy+9++67gWwWAABjlZSUKDs7W3fffbc+/vhjLVmyRE8//bQkKSUlRU6nU0uWLNE999yjPXv2aN68eX7r9+rVSw6HQ+vXr9eNN96o8847TzExMerdu7fy8/N1+PBh+1yVrKwsvfjiixozZowefPBBxcXF6fPPP9fatWv10ksv2eEhmCoqKuTz+dTY2KgDBw5o7ty5uvjii9W/f/+gb+vbeBYPAABtMH78eJ0+fVpDhw5VVlaW7r//fk2dOlWS1LNnT61cuVLr1q1TamqqFixYoKeeespv/Z/85Cd69NFHNWvWLCUkJGjatGmSvjmtIjw8XKmpqfZPQ0lJSdq6dasaGho0YsQIDRgwQNOnT1dsbKzCwtrnK33SpElKTEzUBRdcoDFjxujSSy/Ve++91+5HcBzWty/W7iCqqqrkdrtVWVnZLifMTl65M+h1Lp94VdDrBIDOoqamRocOHVKfPn1+8IoS05hwt1fT/NDfMpDvb46gAAAA4xBQAACAcUJzCjAAAJ3A5s2bQ92ETosjKAAAwDgEFAAAYBwCCgDAGB3wwlJ8R7D+hgQUAEDINd2C/dSpUyFuCdqq6W/47dvqnw1OkgUAhFx4eLhiY2PtB9dFR0f7PYAP5rMsS6dOnVJ5ebliY2PbfFdbAgoAwAgej0fSjz9dF2aLjY21/5ZtQUABABjB4XAoMTFR8fHxqq+vD3VzcBYiIyOD9jwgAgoAwCjh4eHt8tA7dCycJAsAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjtCmgLFiwQA6HQ9OnT7fn1dTUKCsrS927d1dMTIxGjRqlsrIyv/VKSkqUmZmp6OhoxcfHa8aMGTpz5kxbmgIAADqRsw4oO3fu1L//+79r4MCBfvMfeOABvf3221q3bp22bNmio0eP6vbbb7eXNzQ0KDMzU3V1ddq2bZtWrVqllStXas6cOWffCwAA0KmcVUCprq7W2LFj9eKLL+r888+351dWVmr58uV65plndP3112vw4MF6+eWXtW3bNm3fvl2S9Ne//lX79u3Tn/70J11xxRUaOXKk5s2bp6VLl6quri44vQIAAB3aWQWUrKwsZWZmKj093W9+YWGh6uvr/eb369dPKSkpysvLkyTl5eVpwIABSkhIsMtkZGSoqqpKe/fubXF7tbW1qqqq8psAAEDnFRHoCmvXrtXHH3+snTt3Nlvm8/nkdDoVGxvrNz8hIUE+n88u8+1w0rS8aVlL5s+fr0cffTTQpgIAgA4qoCMoR44c0f3336/Vq1erS5cu7dWmZnJyclRZWWlPR44cOWfbBgAA515AAaWwsFDl5eW68sorFRERoYiICG3ZskXPPfecIiIilJCQoLq6OlVUVPitV1ZWJo/HI0nyeDzNruppet9U5ruioqLkcrn8JgAA0HkFFFBuuOEG7d69W0VFRfY0ZMgQjR071n4dGRmp3Nxce53i4mKVlJTI6/VKkrxer3bv3q3y8nK7zMaNG+VyuZSamhqkbgEAgI4soHNQunXrpssuu8xvXteuXdW9e3d7/uTJk5Wdna24uDi5XC7dd9998nq9GjZsmCRpxIgRSk1N1bhx47Rw4UL5fD7Nnj1bWVlZioqKClK3AABARxbwSbI/ZtGiRQoLC9OoUaNUW1urjIwMPf/88/by8PBwrV+/Xvfee6+8Xq+6du2qCRMmaO7cucFuCgAA6KAclmVZoW5EoKqqquR2u1VZWdku56NMXtn8CqW2Wj7xqqDXCQBARxLI9zfP4gEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGCeggPLCCy9o4MCBcrlccrlc8nq9eu+99+zlNTU1ysrKUvfu3RUTE6NRo0aprKzMr46SkhJlZmYqOjpa8fHxmjFjhs6cOROc3gAAgE4hoIBywQUXaMGCBSosLFRBQYGuv/563Xrrrdq7d68k6YEHHtDbb7+tdevWacuWLTp69Khuv/12e/2GhgZlZmaqrq5O27Zt06pVq7Ry5UrNmTMnuL0CAAAdmsOyLKstFcTFxenJJ5/UHXfcoZ49e2rNmjW64447JEn79+9X//79lZeXp2HDhum9997TTTfdpKNHjyohIUGStGzZMs2cOVPHjh2T0+ls1TarqqrkdrtVWVkpl8vVlua3aPLKnUGvc/nEq4JeJwAAHUkg399nfQ5KQ0OD1q5dq5MnT8rr9aqwsFD19fVKT0+3y/Tr108pKSnKy8uTJOXl5WnAgAF2OJGkjIwMVVVV2UdhAAAAIgJdYffu3fJ6vaqpqVFMTIz+8pe/KDU1VUVFRXI6nYqNjfUrn5CQIJ/PJ0ny+Xx+4aRpedOy71NbW6va2lr7fVVVVaDNDsh9ZbODX+ma2ODX+W13vdq+9QMAcA4FfATlkksuUVFRkfLz83XvvfdqwoQJ2rdvX3u0zTZ//ny53W57Sk5ObtftAQCA0Ao4oDidTl100UUaPHiw5s+fr8svv1yLFy+Wx+NRXV2dKioq/MqXlZXJ4/FIkjweT7OrepreN5VpSU5OjiorK+3pyJEjgTYbAAB0IG2+D0pjY6Nqa2s1ePBgRUZGKjc3115WXFyskpISeb1eSZLX69Xu3btVXl5ul9m4caNcLpdSU1O/dxtRUVH2pc1NEwAA6LwCOgclJydHI0eOVEpKik6cOKE1a9Zo8+bNev/99+V2uzV58mRlZ2crLi5OLpdL9913n7xer4YNGyZJGjFihFJTUzVu3DgtXLhQPp9Ps2fPVlZWlqKiotqlgwAAoOMJKKCUl5dr/PjxKi0tldvt1sCBA/X+++/rZz/7mSRp0aJFCgsL06hRo1RbW6uMjAw9//zz9vrh4eFav3697r33Xnm9XnXt2lUTJkzQ3Llzg9srAADQobX5Piih0N73QSl6IiPodV6RHBv0Ov1wFQ8AwHDn5D4oAAAA7YWAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYJyIUDfgH0XRkYpWlbsiObZd2wEAQEfAERQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgBBZT58+frqquuUrdu3RQfH69f/OIXKi4u9itTU1OjrKwsde/eXTExMRo1apTKysr8ypSUlCgzM1PR0dGKj4/XjBkzdObMmbb3BgAAdAoBBZQtW7YoKytL27dv18aNG1VfX68RI0bo5MmTdpkHHnhAb7/9ttatW6ctW7bo6NGjuv322+3lDQ0NyszMVF1dnbZt26ZVq1Zp5cqVmjNnTvB6BQAAOjSHZVnW2a587NgxxcfHa8uWLbrmmmtUWVmpnj17as2aNbrjjjskSfv371f//v2Vl5enYcOG6b333tNNN92ko0ePKiEhQZK0bNkyzZw5U8eOHZPT6fzR7VZVVcntdquyslIul+tsm/+9ip7ICHqdrXVFcuzZrXjXq0FtBwAAwRbI93ebzkGprKyUJMXFxUmSCgsLVV9fr/T0dLtMv379lJKSory8PElSXl6eBgwYYIcTScrIyFBVVZX27t3b4nZqa2tVVVXlNwEAgM7rrANKY2Ojpk+frn/+53/WZZddJkny+XxyOp2KjY31K5uQkCCfz2eX+XY4aVretKwl8+fPl9vttqfk5OSzbTYAAOgAzjqgZGVlac+ePVq7dm0w29OinJwcVVZW2tORI0fafZsAACB0Is5mpWnTpmn9+vX66KOPdMEFF9jzPR6P6urqVFFR4XcUpaysTB6Pxy6zY8cOv/qarvJpKvNdUVFRioqKOpumAgCADiigIyiWZWnatGn6y1/+ok2bNqlPnz5+ywcPHqzIyEjl5uba84qLi1VSUiKv1ytJ8nq92r17t8rLy+0yGzdulMvlUmpqalv6AgAAOomAjqBkZWVpzZo1euutt9StWzf7nBG3263zzjtPbrdbkydPVnZ2tuLi4uRyuXTffffJ6/Vq2LBhkqQRI0YoNTVV48aN08KFC+Xz+TR79mxlZWVxlAQAAEgKMKC88MILkqTrrrvOb/7LL7+siRMnSpIWLVqksLAwjRo1SrW1tcrIyNDzzz9vlw0PD9f69et17733yuv1qmvXrpowYYLmzp3btp4AAIBOo033QQkV7oPSAu6DAgAw3Dm7DwoAAEB7IKAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIwTcED56KOPdPPNNyspKUkOh0Nvvvmm33LLsjRnzhwlJibqvPPOU3p6ug4cOOBX5vjx4xo7dqxcLpdiY2M1efJkVVdXt6kjAACg8wg4oJw8eVKXX365li5d2uLyhQsX6rnnntOyZcuUn5+vrl27KiMjQzU1NXaZsWPHau/evdq4caPWr1+vjz76SFOnTj37XgAAgE4lItAVRo4cqZEjR7a4zLIsPfvss5o9e7ZuvfVWSdIf//hHJSQk6M0339Sdd96pTz/9VBs2bNDOnTs1ZMgQSdKSJUt044036qmnnlJSUlIbugMAADqDoJ6DcujQIfl8PqWnp9vz3G630tLSlJeXJ0nKy8tTbGysHU4kKT09XWFhYcrPz2+x3traWlVVVflNAACg8wpqQPH5fJKkhIQEv/kJCQn2Mp/Pp/j4eL/lERERiouLs8t81/z58+V2u+0pOTk5mM0GAACG6RBX8eTk5KiystKejhw5EuomAQCAdhTUgOLxeCRJZWVlfvPLysrsZR6PR+Xl5X7Lz5w5o+PHj9tlvisqKkoul8tvAgAAnVdQA0qfPn3k8XiUm5trz6uqqlJ+fr68Xq8kyev1qqKiQoWFhXaZTZs2qbGxUWlpacFsDgAA6KACvoqnurpan3/+uf3+0KFDKioqUlxcnFJSUjR9+nQ99thj6tu3r/r06aOHHnpISUlJ+sUvfiFJ6t+/v37+859rypQpWrZsmerr6zVt2jTdeeedXMEDAAAknUVAKSgo0PDhw+332dnZkqQJEyZo5cqVevDBB3Xy5ElNnTpVFRUV+ulPf6oNGzaoS5cu9jqrV6/WtGnTdMMNNygsLEyjRo3Sc889F4TuAACAzsBhWZYV6kYEqqqqSm63W5WVle1yPkrRExlBr7O1rkiObXXZoiMV9uslCY99b7nlE69qQ4sAAAiOQL6/O8RVPAAA4B8LAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOBGhbgDQoaz5VahbELi7Xg11CwAgYBxBAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh/ug/AOYvHJnq8sun3hVO7YEAIDW4QgKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4IQ0oS5cuVe/evdWlSxelpaVpx44doWwOAAAwRMhu1Pbqq68qOztby5YtU1pamp599lllZGSouLhY8fHxoWoWWombvwEA2lPIjqA888wzmjJliiZNmqTU1FQtW7ZM0dHRWrFiRaiaBAAADBGSIyh1dXUqLCxUTk6OPS8sLEzp6enKy8trVr62tla1tbX2+8rKSklSVVVVu7SvuuZMu9TbGlWn6ltd9tvtrDtdHZztt3JMA9lee/2dQiKAv48xOtP4m+y1iaFuQeBGrwx1C/APpun7wLKsHy9shcDf//53S5K1bds2v/kzZsywhg4d2qz8ww8/bEliYmJiYmJi6gTTkSNHfjQrdIiHBebk5Cg7O9t+39jYqOPHj6t79+5yOBxB3VZVVZWSk5N15MgRuVyuoNbd2TBWrcdYtR5j1XqMVesxVoFpr/GyLEsnTpxQUlLSj5YNSUDp0aOHwsPDVVZW5je/rKxMHo+nWfmoqChFRUX5zYuNjW3PJsrlcvEhbiXGqvUYq9ZjrFqPsWo9xiow7TFebre7VeVCcpKs0+nU4MGDlZuba89rbGxUbm6uvF5vKJoEAAAMErKfeLKzszVhwgQNGTJEQ4cO1bPPPquTJ09q0qRJoWoSAAAwRMgCyq9+9SsdO3ZMc+bMkc/n0xVXXKENGzYoISEhVE2S9M3PSQ8//HCzn5TQHGPVeoxV6zFWrcdYtR5jFRgTxsthWa251gcAAODc4Vk8AADAOAQUAABgHAIKAAAwDgEFAAAYp9MFlKVLl6p3797q0qWL0tLStGPHjh8sv27dOvXr109dunTRgAED9O677/ottyxLc+bMUWJios477zylp6frwIEDfmWOHz+usWPHyuVyKTY2VpMnT1Z1dXCejdOeQjFWvXv3lsPh8JsWLFgQ9L4FW7DH6o033tCIESPsuyEXFRU1q6OmpkZZWVnq3r27YmJiNGrUqGY3NzRRKMbquuuua/a5uueee4LZrXYRzLGqr6/XzJkzNWDAAHXt2lVJSUkaP368jh496ldHR91fSaEZL/ZZ33jkkUfUr18/de3aVeeff77S09OVn5/vVybon61gPFvHFGvXrrWcTqe1YsUKa+/evdaUKVOs2NhYq6ysrMXyW7dutcLDw62FCxda+/bts2bPnm1FRkZau3fvtsssWLDAcrvd1ptvvmnt2rXLuuWWW6w+ffpYp0+ftsv8/Oc/ty6//HJr+/bt1n/9139ZF110kTVmzJh2729bhGqsevXqZc2dO9cqLS21p+rq6nbvb1u0x1j98Y9/tB599FHrxRdftCRZn3zySbN67rnnHis5OdnKzc21CgoKrGHDhllXX311e3UzKEI1Vtdee601ZcoUv89VZWVle3UzKII9VhUVFVZ6err16quvWvv377fy8vKsoUOHWoMHD/arpyPurywrdOPFPusbq1evtjZu3GgdPHjQ2rNnjzV58mTL5XJZ5eXldplgf7Y6VUAZOnSolZWVZb9vaGiwkpKSrPnz57dYfvTo0VZmZqbfvLS0NOvuu++2LMuyGhsbLY/HYz355JP28oqKCisqKsr685//bFmWZe3bt8+SZO3cudMu895771kOh8P6+9//HrS+BVsoxsqyvvnHvmjRoiD2pP0Fe6y+7dChQy1+6VZUVFiRkZHWunXr7HmffvqpJcnKy8trQ2/aVyjGyrK+CSj3339/m9p+rrXnWDXZsWOHJcn64osvLMvquPsrywrNeFkW+6zvU1lZaUmyPvjgA8uy2uez1Wl+4qmrq1NhYaHS09PteWFhYUpPT1deXl6L6+Tl5fmVl6SMjAy7/KFDh+Tz+fzKuN1upaWl2WXy8vIUGxurIUOG2GXS09MVFhbW7PCXKUI1Vk0WLFig7t27a9CgQXryySd15syZYHUt6NpjrFqjsLBQ9fX1fvX069dPKSkpAdVzLoVqrJqsXr1aPXr00GWXXaacnBydOnUq4DrOlXM1VpWVlXI4HPazyzri/koK3Xg1YZ/VfBv/8R//Ibfbrcsvv9yuI9ifrQ7xNOPW+PLLL9XQ0NDsTrQJCQnav39/i+v4fL4Wy/t8Pnt507wfKhMfH++3PCIiQnFxcXYZ04RqrCTp3/7t33TllVcqLi5O27ZtU05OjkpLS/XMM8+0uV/toT3GqjV8Pp+cTmezHWWg9ZxLoRorSbrrrrvUq1cvJSUl6W9/+5tmzpyp4uJivfHGG4F14hw5F2NVU1OjmTNnasyYMfbD3jri/koK3XhJ7LO+bf369brzzjt16tQpJSYmauPGjerRo4ddR7A/W50moKBjyM7Otl8PHDhQTqdTd999t+bPn88tqHHWpk6dar8eMGCAEhMTdcMNN+jgwYP6p3/6pxC2LDTq6+s1evRoWZalF154IdTNMd4PjRf7rP83fPhwFRUV6csvv9SLL76o0aNHKz8/v1kwCZZO8xNPjx49FB4e3uwqh7KyMnk8nhbX8Xg8P1i+6b8/Vqa8vNxv+ZkzZ3T8+PHv3W6ohWqsWpKWlqYzZ87o8OHDgXbjnGiPsWoNj8ejuro6VVRUtKmecylUY9WStLQ0SdLnn3/epnraS3uOVdOX7RdffKGNGzf6HQ3oiPsrKXTj1ZJ/5H1W165dddFFF2nYsGFavny5IiIitHz5cruOYH+2Ok1AcTqdGjx4sHJzc+15jY2Nys3NldfrbXEdr9frV16SNm7caJfv06ePPB6PX5mqqirl5+fbZbxeryoqKlRYWGiX2bRpkxobG+2dpGlCNVYtKSoqUlhYWLsl8LZqj7FqjcGDBysyMtKvnuLiYpWUlARUz7kUqrFqSdOlyImJiW2qp72011g1fdkeOHBAH3zwgbp3796sjo62v5JCN14tYZ/1/xobG1VbW2vXEfTP1lmdWmuotWvXWlFRUdbKlSutffv2WVOnTrViY2Mtn89nWZZljRs3zpo1a5ZdfuvWrVZERIT11FNPWZ9++qn18MMPt3jpbGxsrPXWW29Zf/vb36xbb721xcuMBw0aZOXn51v//d//bfXt29f4y/ZCMVbbtm2zFi1aZBUVFVkHDx60/vSnP1k9e/a0xo8ff247H6D2GKuvvvrK+uSTT6x33nnHkmStXbvW+uSTT6zS0lK7zD333GOlpKRYmzZtsgoKCiyv12t5vd5z1/GzEIqx+vzzz625c+daBQUF1qFDh6y33nrLuvDCC61rrrnm3HY+QMEeq7q6OuuWW26xLrjgAquoqMjvstja2lq7no64v7Ks0IwX+6xvxqq6utrKycmx8vLyrMOHD1sFBQXWpEmTrKioKGvPnj12PcH+bHWqgGJZlrVkyRIrJSXFcjqd1tChQ63t27fby6699lprwoQJfuVfe+016+KLL7acTqd16aWXWu+8847f8sbGRuuhhx6yEhISrKioKOuGG26wiouL/cp89dVX1pgxY6yYmBjL5XJZkyZNsk6cONFufQyWcz1WhYWFVlpamuV2u60uXbpY/fv3tx5//HGrpqamXfsZDMEeq5dfftmS1Gx6+OGH7TKnT5+2/vVf/9U6//zzrejoaOu2227zCzCmOtdjVVJSYl1zzTVWXFycFRUVZV100UXWjBkzjL8PimUFd6yaLsNuafrwww/tch11f2VZ53682Gd94/Tp09Ztt91mJSUlWU6n00pMTLRuueUWa8eOHX51BPuz5bAsyzq7Yy8AAADto9OcgwIAADoPAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjPN/UuiGBeAiVdgAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "skew(p_hat_a), skew(p_hat_b)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "xV9Szdjv_mqZ",
        "outputId": "30c7b023-3342-473b-e9e5-82e088b80d54"
      },
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(1.0075073847610363, 3.132640198793857)"
            ]
          },
          "metadata": {},
          "execution_count": 6
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# 参考： 標本 CV 率の差の分布\n",
        "fig, ax = plt.subplots()\n",
        "ax.hist(p_hat_a - p_hat_b)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 517
        },
        "id": "ctlX76hhK66P",
        "outputId": "8df84864-c521-48cb-83c6-8c7a34054059"
      },
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(array([2.000e+00, 0.000e+00, 2.500e+01, 1.100e+01, 0.000e+00, 8.470e+02,\n",
              "        6.400e+01, 0.000e+00, 8.902e+03, 1.490e+02]),\n",
              " array([-0.03  , -0.0263, -0.0226, -0.0189, -0.0152, -0.0115, -0.0078,\n",
              "        -0.0041, -0.0004,  0.0033,  0.007 ]),\n",
              " <BarContainer object of 10 artists>)"
            ]
          },
          "metadata": {},
          "execution_count": 7
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGdCAYAAAAMm0nCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAl1ElEQVR4nO3df1zV9d3/8SeIHDA9IBrnQKGy1fyxXP5KosxqckmNrq3mrms6sq7N/LGgZjQNbitcbaVpZmqltXJ4u8ppXls/JmWRpi5FNNTCX8yr6aXNHWhDONkUFN7fP/ryuXkm/jgIHt7wuN9u51Z8Pu9zeL84Eo8O5xzDjDFGAAAAFgkP9QYAAACCRcAAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsE5EqDfQWhoaGnT48GF169ZNYWFhod4OAAA4D8YYffHFF0pMTFR4+JkfZ2m3AXP48GElJSWFehsAAKAZDh06pMsvv/yM59ttwHTr1k3SV18At9sd4t0AAIDz4ff7lZSU5PwcP5N2GzCNvzZyu90EDAAAljnX0z94Ei8AALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKwTEeoNAADQkvrkFoZ6C0E7MCsj1FuwDo/AAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALBOUAFTX1+vRx55RMnJyYqOjtbXv/51/epXv5IxxlljjFF+fr4SEhIUHR2ttLQ07du3L+B2qqqqlJmZKbfbrdjYWE2YMEFHjx4NWPPJJ5/ohhtuUFRUlJKSkjR79uwLGBMAALQnQQXMk08+qUWLFunZZ5/Vnj179OSTT2r27NlauHChs2b27NlasGCBFi9erJKSEl1yySVKT0/X8ePHnTWZmZnatWuXioqKtGrVKm3YsEGTJk1yzvv9fo0ePVq9e/dWaWmp5syZo1/+8pd68cUXW2BkAABguzBz6sMn53DbbbfJ4/Ho5Zdfdo6NGTNG0dHReuWVV2SMUWJioh588EH9/Oc/lyTV1NTI4/GooKBAY8eO1Z49ezRgwABt3bpVw4YNkyStXr1a3/nOd/TZZ58pMTFRixYt0i9+8Qv5fD5FRkZKknJzc/XGG29o796957VXv9+vmJgY1dTUyO12n/cXBABgtz65haHeQtAOzMoI9RbajPP9+R3UIzDXXXed1qxZoz//+c+SpI8//lgffvihbr31VknS/v375fP5lJaW5lwnJiZGKSkpKi4uliQVFxcrNjbWiRdJSktLU3h4uEpKSpw1I0eOdOJFktLT01VeXq4jR440ubfa2lr5/f6ACwAAaJ8iglmcm5srv9+vfv36qVOnTqqvr9fjjz+uzMxMSZLP55MkeTyegOt5PB7nnM/nU3x8fOAmIiIUFxcXsCY5Ofm022g8171799P2NnPmTD366KPBjAMAACwV1CMwr732ml599VUtW7ZM27Zt09KlS/XUU09p6dKlrbW/85aXl6eamhrncujQoVBvCQAAtJKgHoGZNm2acnNzNXbsWEnSwIED9X//93+aOXOm7r77bnm9XklSRUWFEhISnOtVVFRo0KBBkiSv16vKysqA2z158qSqqqqc63u9XlVUVASsafy4cc2/crlccrlcwYwDAAAsFdQjMP/85z8VHh54lU6dOqmhoUGSlJycLK/XqzVr1jjn/X6/SkpKlJqaKklKTU1VdXW1SktLnTVr165VQ0ODUlJSnDUbNmzQiRMnnDVFRUXq27dvk78+AgAAHUtQAfPv//7vevzxx1VYWKgDBw7o9ddf19NPP6077rhDkhQWFqapU6fq17/+td566y2VlZXprrvuUmJiom6//XZJUv/+/XXLLbdo4sSJ2rJlizZu3Kjs7GyNHTtWiYmJkqQf/ehHioyM1IQJE7Rr1y6tWLFC8+fPV05OTstODwAArBTUr5AWLlyoRx55RPfee68qKyuVmJioyZMnKz8/31kzffp0ffnll5o0aZKqq6s1YsQIrV69WlFRUc6aV199VdnZ2Ro1apTCw8M1ZswYLViwwDkfExOj9957T1lZWRo6dKh69uyp/Pz8gPeKAQAAHVdQ7wNjE94HBgA6Jt4Hxm6t8j4wAAAAbQEBAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6QQfMX//6V915553q0aOHoqOjNXDgQH300UfOeWOM8vPzlZCQoOjoaKWlpWnfvn0Bt1FVVaXMzEy53W7FxsZqwoQJOnr0aMCaTz75RDfccIOioqKUlJSk2bNnN3NEAADQ3gQVMEeOHNH111+vzp0765133tHu3bs1d+5cde/e3Vkze/ZsLViwQIsXL1ZJSYkuueQSpaen6/jx486azMxM7dq1S0VFRVq1apU2bNigSZMmOef9fr9Gjx6t3r17q7S0VHPmzNEvf/lLvfjiiy0wMgAAsF2YMcac7+Lc3Fxt3LhRf/rTn5o8b4xRYmKiHnzwQf385z+XJNXU1Mjj8aigoEBjx47Vnj17NGDAAG3dulXDhg2TJK1evVrf+c539NlnnykxMVGLFi3SL37xC/l8PkVGRjqf+4033tDevXvPa69+v18xMTGqqamR2+0+3xEBAJbrk1sY6i0E7cCsjFBvoc0435/fQT0C89Zbb2nYsGH6j//4D8XHx2vw4MH6zW9+45zfv3+/fD6f0tLSnGMxMTFKSUlRcXGxJKm4uFixsbFOvEhSWlqawsPDVVJS4qwZOXKkEy+SlJ6ervLych05cqTJvdXW1srv9wdcAABA+xRUwPzlL3/RokWLdOWVV+rdd9/VT3/6U91///1aunSpJMnn80mSPB5PwPU8Ho9zzufzKT4+PuB8RESE4uLiAtY0dRunfo5/NXPmTMXExDiXpKSkYEYDAAAWCSpgGhoaNGTIED3xxBMaPHiwJk2apIkTJ2rx4sWttb/zlpeXp5qaGudy6NChUG8JAAC0kqACJiEhQQMGDAg41r9/fx08eFCS5PV6JUkVFRUBayoqKpxzXq9XlZWVAedPnjypqqqqgDVN3capn+NfuVwuud3ugAsAAGifggqY66+/XuXl5QHH/vznP6t3796SpOTkZHm9Xq1Zs8Y57/f7VVJSotTUVElSamqqqqurVVpa6qxZu3atGhoalJKS4qzZsGGDTpw44awpKipS3759A17xBAAAOqagAuaBBx7Q5s2b9cQTT+h///d/tWzZMr344ovKysqSJIWFhWnq1Kn69a9/rbfeektlZWW66667lJiYqNtvv13SV4/Y3HLLLZo4caK2bNmijRs3Kjs7W2PHjlViYqIk6Uc/+pEiIyM1YcIE7dq1SytWrND8+fOVk5PTstMDAAArRQSz+JprrtHrr7+uvLw8PfbYY0pOTtYzzzyjzMxMZ8306dP15ZdfatKkSaqurtaIESO0evVqRUVFOWteffVVZWdna9SoUQoPD9eYMWO0YMEC53xMTIzee+89ZWVlaejQoerZs6fy8/MD3isGAAB0XEG9D4xNeB8YAOiYeB8Yu7XK+8AAAAC0BQQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOtcUMDMmjVLYWFhmjp1qnPs+PHjysrKUo8ePdS1a1eNGTNGFRUVAdc7ePCgMjIy1KVLF8XHx2vatGk6efJkwJp169ZpyJAhcrlcuuKKK1RQUHAhWwUAAO1IswNm69ateuGFF/Stb30r4PgDDzygP/7xj1q5cqXWr1+vw4cP6/vf/75zvr6+XhkZGaqrq9OmTZu0dOlSFRQUKD8/31mzf/9+ZWRk6Oabb9aOHTs0depU3XPPPXr33Xebu10AANCONCtgjh49qszMTP3mN79R9+7dneM1NTV6+eWX9fTTT+vb3/62hg4dqt/+9rfatGmTNm/eLEl67733tHv3br3yyisaNGiQbr31Vv3qV7/Sc889p7q6OknS4sWLlZycrLlz56p///7Kzs7WD37wA82bN68FRgYAALZrVsBkZWUpIyNDaWlpAcdLS0t14sSJgOP9+vVTr169VFxcLEkqLi7WwIED5fF4nDXp6eny+/3atWuXs+Zfbzs9Pd25jabU1tbK7/cHXAAAQPsUEewVli9frm3btmnr1q2nnfP5fIqMjFRsbGzAcY/HI5/P56w5NV4azzeeO9sav9+vY8eOKTo6+rTPPXPmTD366KPBjgMAACwU1CMwhw4d0s9+9jO9+uqrioqKaq09NUteXp5qamqcy6FDh0K9JQAA0EqCCpjS0lJVVlZqyJAhioiIUEREhNavX68FCxYoIiJCHo9HdXV1qq6uDrheRUWFvF6vJMnr9Z72qqTGj8+1xu12N/noiyS5XC653e6ACwAAaJ+CCphRo0aprKxMO3bscC7Dhg1TZmam8++dO3fWmjVrnOuUl5fr4MGDSk1NlSSlpqaqrKxMlZWVzpqioiK53W4NGDDAWXPqbTSuabwNAADQsQX1HJhu3brpqquuCjh2ySWXqEePHs7xCRMmKCcnR3FxcXK73brvvvuUmpqqa6+9VpI0evRoDRgwQOPHj9fs2bPl8/n08MMPKysrSy6XS5I0ZcoUPfvss5o+fbp+8pOfaO3atXrttddUWFjYEjMDAADLBf0k3nOZN2+ewsPDNWbMGNXW1io9PV3PP/+8c75Tp05atWqVfvrTnyo1NVWXXHKJ7r77bj322GPOmuTkZBUWFuqBBx7Q/Pnzdfnll+ull15Senp6S28XAABYKMwYY0K9idbg9/sVExOjmpoang8DAB1In1z7Hq0/MCsj1FtoM8735zd/FxIAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOkEFzMyZM3XNNdeoW7duio+P1+23367y8vKANcePH1dWVpZ69Oihrl27asyYMaqoqAhYc/DgQWVkZKhLly6Kj4/XtGnTdPLkyYA169at05AhQ+RyuXTFFVeooKCgeRMCAIB2J6iAWb9+vbKysrR582YVFRXpxIkTGj16tL788ktnzQMPPKA//vGPWrlypdavX6/Dhw/r+9//vnO+vr5eGRkZqqur06ZNm7R06VIVFBQoPz/fWbN//35lZGTo5ptv1o4dOzR16lTdc889evfdd1tgZAAAYLswY4xp7pU///xzxcfHa/369Ro5cqRqamp06aWXatmyZfrBD34gSdq7d6/69++v4uJiXXvttXrnnXd022236fDhw/J4PJKkxYsX66GHHtLnn3+uyMhIPfTQQyosLNTOnTudzzV27FhVV1dr9erV57U3v9+vmJgY1dTUyO12N3dEAIBl+uQWhnoLQTswKyPUW2gzzvfn9wU9B6ampkaSFBcXJ0kqLS3ViRMnlJaW5qzp16+fevXqpeLiYklScXGxBg4c6MSLJKWnp8vv92vXrl3OmlNvo3FN4200pba2Vn6/P+ACAADap2YHTENDg6ZOnarrr79eV111lSTJ5/MpMjJSsbGxAWs9Ho98Pp+z5tR4aTzfeO5sa/x+v44dO9bkfmbOnKmYmBjnkpSU1NzRAABAG9fsgMnKytLOnTu1fPnyltxPs+Xl5ammpsa5HDp0KNRbAgAArSSiOVfKzs7WqlWrtGHDBl1++eXOca/Xq7q6OlVXVwc8ClNRUSGv1+us2bJlS8DtNb5K6dQ1//rKpYqKCrndbkVHRze5J5fLJZfL1ZxxAACAZYJ6BMYYo+zsbL3++utau3atkpOTA84PHTpUnTt31po1a5xj5eXlOnjwoFJTUyVJqampKisrU2VlpbOmqKhIbrdbAwYMcNacehuNaxpvAwAAdGxBPQKTlZWlZcuW6c0331S3bt2c56zExMQoOjpaMTExmjBhgnJychQXFye326377rtPqampuvbaayVJo0eP1oABAzR+/HjNnj1bPp9PDz/8sLKyspxHUKZMmaJnn31W06dP109+8hOtXbtWr732mgoL7XtmOQAAaHlBPQKzaNEi1dTU6KabblJCQoJzWbFihbNm3rx5uu222zRmzBiNHDlSXq9Xf/jDH5zznTp10qpVq9SpUyelpqbqzjvv1F133aXHHnvMWZOcnKzCwkIVFRXp6quv1ty5c/XSSy8pPT29BUYGAAC2u6D3gWnLeB8YAOiYeB8Yu12U94EBAAAIBQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWiQj1BgCgOfrkFoZ6C0E7MCsj1FsA2g0egQEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUiQr0BAAA6uj65haHeQtAOzMoI6efnERgAAGCdNh0wzz33nPr06aOoqCilpKRoy5Ytod4SAABoA9pswKxYsUI5OTmaMWOGtm3bpquvvlrp6emqrKwM9dYAAECItdnnwDz99NOaOHGifvzjH0uSFi9erMLCQi1ZskS5ubkh3h0ABI/nOQAtp00GTF1dnUpLS5WXl+ccCw8PV1pamoqLi5u8Tm1trWpra52Pa2pqJEl+v791N4tWddWMd0O9haDtfDQ91FvoEBpq/xnqLXQINv43lD8bF0dr/dlovF1jzFnXtcmA+fvf/676+np5PJ6A4x6PR3v37m3yOjNnztSjjz562vGkpKRW2SNwJjHPhHoHQMvhzzPOpLX/bHzxxReKiYk54/k2GTDNkZeXp5ycHOfjhoYGVVVVqUePHgoLC2uxz+P3+5WUlKRDhw7J7Xa32O22ZczMzO0VMzNze2XzzMYYffHFF0pMTDzrujYZMD179lSnTp1UUVERcLyiokJer7fJ67hcLrlcroBjsbGxrbVFud1u6/5QXChm7hiYuWNg5o7B1pnP9shLozb5KqTIyEgNHTpUa9ascY41NDRozZo1Sk1NDeHOAABAW9AmH4GRpJycHN19990aNmyYhg8frmeeeUZffvml86okAADQcbXZgPnhD3+ozz//XPn5+fL5fBo0aJBWr1592hN7LzaXy6UZM2ac9uuq9oyZOwZm7hiYuWPoCDOHmXO9TgkAAKCNaZPPgQEAADgbAgYAAFiHgAEAANYhYAAAgHU6fMBUVVUpMzNTbrdbsbGxmjBhgo4ePXrW6xw/flxZWVnq0aOHunbtqjFjxgS86d4//vEP3XLLLUpMTJTL5VJSUpKys7NP+3sj1q1bpyFDhsjlcumKK65QQUFBa4x4mtaY+eOPP9a4ceOUlJSk6Oho9e/fX/Pnzw+4jXXr1iksLOy0i8/na5U5TxWqmaX2dT9L0v3336+hQ4fK5XJp0KBBp93GgQMHmryfN2/e3JLjNSlUM0vSJ598ohtuuEFRUVFKSkrS7NmzW2qss2qtmQ8ePKiMjAx16dJF8fHxmjZtmk6ePOmcv5jfz88995z69OmjqKgopaSkaMuWLWddv3LlSvXr109RUVEaOHCg3n777YDzxhjl5+crISFB0dHRSktL0759+wLWNOfr2pJCMXOfPn1Ouz9nzZrV4rO1GNPB3XLLLebqq682mzdvNn/605/MFVdcYcaNG3fW60yZMsUkJSWZNWvWmI8++shce+215rrrrnPOV1VVmeeff95s3brVHDhwwLz//vumb9++Abf7l7/8xXTp0sXk5OSY3bt3m4ULF5pOnTqZ1atXt9qsjVpj5pdfftncf//9Zt26debTTz81//3f/22io6PNwoULnTUffPCBkWTKy8vN3/72N+dSX1/farM2CtXM7e1+NsaY++67zzz77LNm/Pjx5uqrrz7tNvbv328kmffffz/gfq6rq2vJ8ZoUqplramqMx+MxmZmZZufOneZ3v/udiY6ONi+88EJLjtek1pj55MmT5qqrrjJpaWlm+/bt5u233zY9e/Y0eXl5zpqL9f28fPlyExkZaZYsWWJ27dplJk6caGJjY01FRUWT6zdu3Gg6depkZs+ebXbv3m0efvhh07lzZ1NWVuasmTVrlomJiTFvvPGG+fjjj813v/tdk5ycbI4dO+asac7XtaWEaubevXubxx57LOD+PHr0aKvP21wdOmB2795tJJmtW7c6x9555x0TFhZm/vrXvzZ5nerqatO5c2ezcuVK59iePXuMJFNcXHzGzzV//nxz+eWXOx9Pnz7dfPOb3wxY88Mf/tCkp6c3d5zzcjFnvvfee83NN9/sfNz4H7wjR45c+CBBCOXM7fl+njFjxlkDZvv27Rc8RzBCOfPzzz9vunfvbmpra51jDz30kOnbt+8FTHRurTXz22+/bcLDw43P53PWLFq0yLjdbmfGi/X9PHz4cJOVleV8XF9fbxITE83MmTObXP+f//mfJiMjI+BYSkqKmTx5sjHGmIaGBuP1es2cOXOc89XV1cblcpnf/e53xpjmfV1bUihmNuargJk3b14LTtK6OvSvkIqLixUbG6thw4Y5x9LS0hQeHq6SkpImr1NaWqoTJ04oLS3NOdavXz/16tVLxcXFTV7n8OHD+sMf/qAbb7wx4HOfehuSlJ6efsbbaCkXa2ZJqqmpUVxc3GnHBw0apISEBP3bv/2bNm7ceAHTnJ9QztwR7ucz+e53v6v4+HiNGDFCb731VvBDBCmUMxcXF2vkyJGKjIx0jqWnp6u8vFxHjhxpxjTn/3lbY+bi4mINHDgw4I1D09PT5ff7tWvXroDba83v57q6OpWWlgbsNTw8XGlpaWe8f871Pbd//375fL6ANTExMUpJSQmYP9iva0sJ1cyNZs2apR49emjw4MGaM2dOwK8N25oOHTA+n0/x8fEBxyIiIhQXF3fG3+P6fD5FRkae9hdFejye064zbtw4denSRZdddpncbrdeeumlgNv513cV9ng88vv9Onbs2AVMdXatPXOjTZs2acWKFZo0aZJzLCEhQYsXL9bvf/97/f73v1dSUpJuuukmbdu27cKGOodQztze7+emdO3aVXPnztXKlStVWFioESNG6Pbbb2/1iAnlzGe6nxvPtZbWmvl85rkY389///vfVV9f3+Rezjbf2dY3/vNca4L9uraUUM0sffVcr+XLl+uDDz7Q5MmT9cQTT2j69OkXPFNraZcBk5ub2+STy0697N27t9X3MW/ePG3btk1vvvmmPv30U+Xk5LTa52orM0vSzp079b3vfU8zZszQ6NGjneN9+/bV5MmTNXToUF133XVasmSJrrvuOs2bN69Zn8eGmVtaW5r5THr27KmcnBylpKTommuu0axZs3TnnXdqzpw5zbo9G2ZuaTbM3NLfzwi9nJwc3XTTTfrWt76lKVOmaO7cuVq4cKFqa2tDvbUmtdm/C+lCPPjgg/qv//qvs6752te+Jq/Xq8rKyoDjJ0+eVFVVlbxeb5PX83q9qqurU3V1dcD/wVRUVJx2Ha/XK6/Xq379+ikuLk433HCDHnnkESUkJMjr9Z72rP+Kigq53W5FR0ef/7D/X1uZeffu3Ro1apQmTZqkhx9++Jz7Hj58uD788MNzrmuKDTO31/s5WCkpKSoqKmrWdW2Y+Uz3c+O5YIV6Zq/Xe9qrXs5nngv5fm5Kz5491alTpya/tmeb72zrG/9ZUVGhhISEgDWNrzBrzte1pYRq5qakpKTo5MmTOnDggPr27duccVpXqJ+EE0qNT9T66KOPnGPvvvvueT0B7n/+53+cY3v37j3nkzvXr19vJJn9+/cbY756cudVV10VsGbcuHEX7cmdrTHzzp07TXx8vJk2bdp57yctLc3ccccdzZjk/IVy5vZ4Pzc60xNam3LPPfeYwYMHBzdEkEI5c+OTeE99pVVeXt5FexJvS8/c+CTeU1/18sILLxi3222OHz9+xv20xvfz8OHDTXZ2tvNxfX29ueyyy876hNbbbrst4FhqauppT2h96qmnnPM1NTVNPok3mK9rSwrFzE155ZVXTHh4uKmqqrqQcVpNhw4YY756qdzgwYNNSUmJ+fDDD82VV14Z8FK5zz77zPTt29eUlJQ4x6ZMmWJ69epl1q5daz766COTmppqUlNTnfOFhYVmyZIlpqyszOzfv9+sWrXK9O/f31x//fXOmsaX106bNs3s2bPHPPfccxf15bUtPXNZWZm59NJLzZ133hnwErzKykpnzbx588wbb7xh9u3bZ8rKyszPfvYzEx4ebt5///12O3N7u5+NMWbfvn1m+/btZvLkyeYb3/iG2b59u9m+fbvz6pSCggKzbNkys2fPHrNnzx7z+OOPm/DwcLNkyZJ2O3N1dbXxeDxm/PjxZufOnWb58uWmS5cuF+1l1C09c+PLqEePHm127NhhVq9ebS699NKAl1FfrO/n5cuXG5fLZQoKCszu3bvNpEmTTGxsrPMKqfHjx5vc3Fxn/caNG01ERIR56qmnzJ49e8yMGTOafElxbGysefPNN80nn3xivve97zX5MuqzfV1bUyhm3rRpk5k3b57ZsWOH+fTTT80rr7xiLr30UnPXXXddlJmbo8MHzD/+8Q8zbtw407VrV+N2u82Pf/xj88UXXzjnG18S+sEHHzjHjh07Zu69917TvXt306VLF3PHHXeYv/3tb875tWvXmtTUVBMTE2OioqLMlVdeaR566KHTXm74wQcfmEGDBpnIyEjzta99zfz2t79t5Wm/0hozz5gxw0g67dK7d29nzZNPPmm+/vWvm6ioKBMXF2duuukms3bt2osxcshmNqZ93c/GGHPjjTc2OXfjo4sFBQWmf//+pkuXLsbtdpvhw4cHvGS3NYVqZmOM+fjjj82IESOMy+Uyl112mZk1a1Zrj2uMab2ZDxw4YG699VYTHR1tevbsaR588EFz4sQJ5/zF/H5euHCh6dWrl4mMjDTDhw83mzdvds7deOON5u677w5Y/9prr5lvfOMbJjIy0nzzm980hYWFAecbGhrMI488Yjwej3G5XGbUqFGmvLw8YM25vq6t7WLPXFpaalJSUpyfW/379zdPPPHEWR9xC7UwY4xpvV9QAQAAtLx2+SokAADQvhEwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArPP/ABzvVI4uDp34AAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "skew(p_hat_a - p_hat_b)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "BbLswB5vStlC",
        "outputId": "3ed52eda-914c-4551-e135-30d785e7d47e"
      },
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "-2.667806182944919"
            ]
          },
          "metadata": {},
          "execution_count": 8
        }
      ]
    }
  ]
}