Wasserman_Chapter_8.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Wasserman_Chapter_8.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyMni8z2+aXQPFkFUWE2GQNI",
      "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/sajad13901/Statistics_Wasserman/blob/main/Wasserman_Chapter_8.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 144,
      "metadata": {
        "id": "5EuHOEDDemoj",
        "outputId": "8e055624-3358-47ad-9704-ee99f4bfce32",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "data size: 799\n",
            "skew_hat: 1.76\n",
            "skew_hat variance: 0.16\n"
          ]
        }
      ],
      "source": [
        "# Chapter 8, nerve data, Example 8.2\n",
        "from google.colab import files\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "from scipy.stats import norm\n",
        "\n",
        "# uploaded = files.upload()\n",
        "df = pd.read_csv('https://raw.githubusercontent.com/sajad13901/Statistics_Wasserman/main/nerve.txt', delimiter= '\\s+', index_col=False)\n",
        "data = df.waiting_time\n",
        "n, mu, sigma = data.size, data.mean(), data.std()\n",
        "\n",
        "def calc_skew(data):\n",
        "    n, mu, sigma = data.size, data.mean(), data.std()\n",
        "    tmp = 0;\n",
        "    for element in data: tmp += (element - mu)**3\n",
        "    return tmp/(n * sigma**3)\n",
        "\n",
        "# use bootstrap to calc variance\n",
        "skew_boot = []\n",
        "B = 1000\n",
        "for _ in range(B):\n",
        "    tmp_data = np.random.choice(data, size=n, replace=True)\n",
        "    skew_boot.append(calc_skew(tmp_data))\n",
        "\n",
        "print('data size:', n)\n",
        "print('skew_hat:', round(calc_skew(data),2))\n",
        "print('skew_hat variance:', round(np.std(skew_boot),2))\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Chapter 8, problem 1\n",
        "import numpy as np\n",
        "from scipy.stats import norm\n",
        "\n",
        "LSAT = np.array([576,   635,  558,  578,  666,  580,  555,  661,  651,  605,  653,  575,  545,  572,  594])\n",
        "GPA  = np.array([3.39, 3.30, 2.81, 3.03, 3.44, 3.07, 3.00, 3.43, 3.36, 3.13, 3.12, 2.74, 2.76, 2.88, 3.96])\n",
        "\n",
        "def corr(X, Y):\n",
        "    mu_x = X.mean()\n",
        "    mu_y = Y.mean()\n",
        "    return sum((X - mu_x) * (Y - mu_y)) / (sum((X - mu_x)**2) * sum((Y - mu_y)**2))**0.5\n",
        "  \n",
        "theta_hat = corr(LSAT, GPA)    \n",
        "print('estimated theta hat: %.4f' %theta_hat)\n",
        "\n",
        "# bootstrap to calc variance of theta_hat\n",
        "n = GPA.size\n",
        "B = int(1e3)\n",
        "theta_boot = []\n",
        "idx_vec = list(range(n))\n",
        "for _ in range(B):\n",
        "    sel_idx = np.random.choice(idx_vec, size=n, replace=True)\n",
        "    sel_gpa = np.take(GPA, sel_idx)\n",
        "    sel_lsat = np.take(LSAT, sel_idx)\n",
        "    theta_boot.append(corr(sel_lsat, sel_gpa))\n",
        "se_hat = np.std(theta_boot)\n",
        "print('bootstrap standard error: %.4f' %se_hat)\n",
        "\n",
        "z_95 = norm.ppf(.975)\n",
        "print('95%% interval Normal method: (%.2f, %.2f)' %(theta_hat-z_95*se_hat, theta_hat+z_95*se_hat))\n",
        "\n",
        "boot_quant_1 = np.quantile(theta_boot, .05/2)\n",
        "boot_quant_2 = np.quantile(theta_boot, 1-.05/2)\n",
        "print('95%% interval Pivotal method: (%.2f, %.2f)' %(2*theta_hat-boot_quant_2, 2*theta_hat-boot_quant_1))\n",
        "\n",
        "print('95%% interval Percentile method: (%.2f, %.2f)' %(boot_quant_1, boot_quant_2))"
      ],
      "metadata": {
        "id": "O5gs4plTjsLf",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "ddc5fd28-f963-4b84-fcec-7f830b15d361"
      },
      "execution_count": 141,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "estimated theta hat: 0.5459\n",
            "bootstrap standard error: 0.1938\n",
            "95% interval Normal method: (0.17, 0.93)\n",
            "95% interval Pivotal method: (0.16, 0.86)\n",
            "95% interval Percentile method: (0.23, 0.93)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Chapter 8, problem 2\n",
        "import numpy as np\n",
        "from scipy.stats import norm\n",
        "\n",
        "n = 50\n",
        "y = np.random.normal(0, 1, n)\n",
        "x = np.exp(y)\n",
        "\n",
        "def calc_skew(data):\n",
        "    n, mu, sigma = data.size, data.mean(), data.std()\n",
        "    tmp = 0;\n",
        "    for element in data: tmp += (element - mu)**3\n",
        "    return tmp/(n * sigma**3)\n",
        "\n",
        "t_hat = calc_skew(x)\n",
        "\n",
        "t_boot = []\n",
        "B = 10000\n",
        "for _ in range(B):\n",
        "  sel_data = np.random.choice(x, size=n, replace=True)\n",
        "  t_boot.append(calc_skew(sel_data))\n",
        "\n",
        "se_hat = np.std(t_boot)\n",
        "print('bootstrap standard error: %.4f' %se_hat)\n",
        "\n",
        "z_95 = norm.ppf(.975)\n",
        "print('95%% interval Normal method: (%.2f, %.2f)' %(t_hat-z_95*se_hat, t_hat+z_95*se_hat))\n",
        "\n",
        "boot_quant_1 = np.quantile(t_boot, .05/2)\n",
        "boot_quant_2 = np.quantile(t_boot, 1-.05/2)\n",
        "print('95%% interval Pivotal method: (%.2f, %.2f)' %(2*t_hat-boot_quant_2, 2*t_hat-boot_quant_1))\n",
        "print('95%% interval Percentile method: (%.2f, %.2f)' %(boot_quant_1, boot_quant_2))\n"
      ],
      "metadata": {
        "id": "msfhONQGjsSu",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "a6e851c4-8a7f-48c9-993c-8497d7fe35c7"
      },
      "execution_count": 142,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "bootstrap standard error: 1.1126\n",
            "95% interval Normal method: (1.95, 6.31)\n",
            "95% interval Pivotal method: (2.64, 7.13)\n",
            "95% interval Percentile method: (1.13, 5.62)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Chapter 8, problem 3\n",
        "import numpy as np\n",
        "from scipy.stats import norm\n",
        "from scipy.stats import t\n",
        "\n",
        "n = 25\n",
        "x = t.rvs(3, size=n)\n",
        "\n",
        "def T(x):\n",
        "    return (np.quantile(x, 0.75) - np.quantile(x, 0.25)) / 1.34\n",
        "\n",
        "B = 10000\n",
        "t_boot = np.empty(B)\n",
        "for i in range(B):\n",
        "  sel_data = np.random.choice(x, size=n, replace=True)\n",
        "  t_boot[i] = T(sel_data)\n",
        "\n",
        "alpha = 0.05\n",
        "t_hat = T(x)\n",
        "se_hat = np.std(t_boot)\n",
        "z_alpha = norm.ppf(1 - alpha/2)\n",
        "boot_quant_1 = np.quantile(t_boot, alpha/2)\n",
        "boot_quant_2 = np.quantile(t_boot, 1 - alpha/2)\n",
        "print('95%% interval Normal method: (%.2f, %.2f)' %(t_hat-z_alpha*se_hat, t_hat+z_alpha*se_hat))\n",
        "print('95%% interval Pivotal method: (%.2f, %.2f)' %(2*t_hat-boot_quant_2, 2*t_hat-boot_quant_1))\n",
        "print('95%% interval Percentile method: (%.2f, %.2f)' %(boot_quant_1, boot_quant_2))\n"
      ],
      "metadata": {
        "id": "Y4vJWdJMZpQe",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "b48659d4-4431-4a2a-c90a-77ad6f6ada50"
      },
      "execution_count": 156,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "95% interval Normal method: (0.52, 1.42)\n",
            "95% interval Pivotal method: (0.53, 1.46)\n",
            "95% interval Percentile method: (0.49, 1.42)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Chapter 8, problem 6\n",
        "import numpy as np\n",
        "from scipy.stats import norm\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "mu = 5\n",
        "n = 1000\n",
        "x = np.random.normal(mu, 1, n)\n",
        "\n",
        "B = 10000\n",
        "t_boot = np.empty(B)\n",
        "for i in range(B):\n",
        "  sel_data = np.random.choice(x, size=n, replace=True)\n",
        "  t_boot[i] = np.exp(sel_data.mean())\n",
        "\n",
        "t_hat = np.exp(x.mean())\n",
        "se_hat = np.std(t_boot)\n",
        "z_alpha = norm.ppf(1 - 0.05/2)\n",
        "print('95%% interval Normal method: (%.2f, %.2f)' %(t_hat-z_alpha*se_hat, t_hat+z_alpha*se_hat))\n",
        "\n",
        "# true ditribution of theta_hat = (1/y) * f_x_bar(log(y)) , f_x_bar = Norm(mu, var_x/n)  \n",
        "def norm_dist(y, mu, var):\n",
        "  return np.exp(-(y-mu)**2/(2*var))/(var*2*np.pi)**0.5\n",
        "\n",
        "y = np.linspace(100, 200, 500)\n",
        "f_y = (1/y) * norm_dist(np.log(y), mu, 1/n)\n",
        "\n",
        "plt.hist(t_boot, 100, density=True, label='bootstrap dist')\n",
        "plt.plot(y, f_y, lw=4, label='true distribution')\n",
        "plt.legend(loc='upper right')\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "nwpGL5lVV5L5",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        },
        "outputId": "d35bc51c-ce09-48a7-bd7f-164937da4327"
      },
      "execution_count": 199,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "95% interval Normal method: (137.46, 156.01)\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZwU1b338c9vFhhgZF9kHxZBQBQVIUZxjVs0GhW3a6JEvV40Jt4kmuC9N27Z9MZcjYk3PiYkLsHEaPTGRzAmUXwwMS4DdxCQRSAjDKDAwLAKzEyf54/umemq7oGB6e6qrv6+X6950efU6apf0cyP06dOnTLnHCIiEl1FQQcgIiLZpUQvIhJxSvQiIhGnRC8iEnFK9CIiEVcSdAB+vXv3dhUVFUGHISKSV+bPn7/ZOdcn3bbQJfqKigoqKyuDDkNEJK+Y2YetbdPQjYhIxCnRi4hEnBK9iEjEhW6MXkQyq76+npqaGvbs2RN0KJIBZWVlDBo0iNLS0ja/R4leJOJqamo47LDDqKiowMyCDkfawTlHbW0tNTU1DBs2rM3v09CNSMTt2bOHXr16KclHgJnRq1evg/52pkQvUgCU5KPjUD5LJXoRkYhToheRrKqrq+O///u/c3a8iooKNm/eDMCnP/3p/bb9/ve/v9/tn/3sZ6mrq6O6upqjjjrqoOJ4/fXXefPNN5vLjz76KE8++eRB7SNTlOhFgIoZs4MOIbL2l+gbGhqyeuzkRJtOa4neOUcsFmPOnDl07979kI7tT/TTp0/nmmuuOaR9tZdm3YgUiru7ZXn/29JWz5gxg1WrVjFhwgTOOusszj//fL797W/To0cPli1bxp/+9CcuuOACFi9eDMADDzzAzp07ufvuu1m1ahVf/vKX2bRpE507d+bnP/85Rx55pGf/tbW1XHXVVaxbt44TTzyR5KfmlZeXs3PnTjZs2MAVV1zB9u3baWho4Gc/+xmzZ8/mk08+YcKECYwbN47vfe97nHPOOUyePJn58+czZ84cTj311OYlWRoaGrj66qtZsGAB48aN48knn6Rz587Ny7b07t2byspKbrvtNh5//HEeffRRiouL+fWvf81PfvITXn31VcrLy7ntttuoqqpi+vTp7N69mxEjRvDLX/6SHj16cNpppzF58mTmzp1LXV0dM2fOZMqUKe3+aNSjF5Gsuu+++xgxYgRVVVX88Ic/BGDBggX8+Mc/ZsWKFft974033shPfvIT5s+fzwMPPMDNN9+c0uaee+7h5JNPZsmSJVx88cWsWbMmpc3TTz/NOeecQ1VVFQsXLmTChAncd999dOrUiaqqKmbNmgXABx98wM0338ySJUsYOnSoZx/Lly/n5ptvZunSpXTt2nW/w1EVFRVMnz6dr33ta1RVVaUk62uuuYb777+f9957j/Hjx3PPPfc0b2toaOCdd97hoYce8tS3h3r0IpJzkyZNOuA88J07d/Lmm29y2WWXNdft3bs3pd28efN4/vnnATj//PPp0aNHSpsTTjiB6667jvr6ej7/+c8zYcKEtMccOnQon/rUp9JuGzx4MCeddBIAX/jCF3j44Ye57bbb9nsO6Wzbto26ujpOPfVUAK699lrPOV5yySUAHH/88VRXVx/0/tNRj15Ecq5Lly7Nr0tKSojFYs3lpjnisViM7t27U1VV1fyzdOnSQzreKaecwrx58xg4cCDTpk1r9aJoclx+/mmNTeXk+DNx93HHjh0BKC4uztg1DPXoRQpFK2Po2XbYYYexY8eOVrf369ePjRs3UltbS3l5OS+99BLnnnsuXbt2ZdiwYTz77LNcdtllOOd47733OOaYYzzvP+WUU3j66af5j//4D15++WW2bt2acowPP/yQQYMG8c///M/s3buXBQsWcM0111BaWkp9fX2blhNYs2YNf//73znxxBN5+umnOfnkk4H4MM38+fM577zz+P3vf+857+3bt6fsp1u3bvTo0YM33niDKVOm8NRTTzX37rNFPXopaJptk329evXipJNO4qijjuL2229P2V5aWsqdd97JpEmTOOusszwXW2fNmsXMmTM55phjGDduHH/4wx9S3n/XXXcxb948xo0bx/PPP8+QIUNS2rz++uscc8wxHHvssTzzzDPceuutQPwawNFHH83VV199wPMYPXo0jzzyCGPGjGHr1q3cdNNNzce/9dZbmThxIsXFxc3tP/e5z/HCCy8wYcIE3njjDc++nnjiCW6//XaOPvpoqqqquPPOOw94/Paw5CvUYTBx4kSnB49IrlTMmE31fec3/xlFS5cuZcyYMUGHIRmU7jM1s/nOuYnp2qtHLyIScUr0IiIRp0QvkqDxeokqJXoRkYhTohcRiTglehGRiNMNUyLJlrwAc38AW6uh75Fw1r0w/LSAg8qsTF+LONC01Orqas+iZe3x0EMPceONN9K5c+dW2zz++OOcffbZDBgwoN3HO5Bp06ZxwQUXMHXqVG644Qa+/vWvM3bs2MDj8lOPXiThsuLX4dlpsHk5NO6FDQvhqUtgxZ+CDk0SHnroIXbv3r3fNo8//jjr169Pu62xsTEbYQHwi1/8otUkD/uPK9uU6EWAgWziOyW/St3gGuEPN8MnqbfVS9s1LfE7ZswYpk6d2pysX331VY499ljGjx/Pdddd17xoWbr6hx9+mPXr13P66adz+umn09jYyLRp0zjqqKMYP348Dz74IM899xyVlZVcffXVTJgwgU8++YSKigq+9a1vcdxxx/Hss8/y85//nBNOOIFjjjmGSy+9tDmWadOmMX36dCZOnMioUaN46aWXUs7DOcctt9zC6NGj+cxnPsPGjRubt5122mlUVla2Oa5cUqIXAb5a8gJlVp9+465N8NbPchtQxKRb4nfPnj1MmzaNZ555hkWLFjWvE99a/Ve/+lUGDBjA3LlzmTt3LlVVVaxbt47FixezaNEivvSlLzF16lQmTpzIrFmzqKqqolOnTkB8GYYFCxZw5ZVXcskll/Duu++ycOFCxowZw8yZM5vjrK6u5p133mH27NlMnz49ZZGyF154geXLl/P+++/z5JNPpn2wycHElStK9CK7NvP54r9668r7ectv/x+ob//KhIXKv8TvX//6V5YvX86wYcMYNWoUEF+ud968ea3W+w0fPpzVq1fzla98hT/+8Y907dq11eNfccUVza8XL17MlClTGD9+PLNmzWLJkiXN2y6//HKKioo44ogjGD58OMuWLfPsZ968eVx11VUUFxczYMAAzjjjjHbFlStK9CKLnqWjJS0H26MCbnmXOpe0ZO2eOliW+lVe2qa1JX7bo0ePHixcuJDTTjuNRx99lBtuuKHVtsnLD0+bNo2f/vSnLFq0iLvuusvTa89EnAcTV64o0Yss881COX4alHXj+UbfI9zeT105UdqmaYlfoHmJ39GjR1NdXc3KlSsBmpfrba0evEseb968mVgsxqWXXsp3v/tdFixYkNImnR07dtC/f3/q6+ubnyzV5NlnnyUWi7Fq1SpWr17N6NGjPdtPOeUUnnnmGRobG9mwYQNz585N2f+hxpVNml4pBa0bO+FD3zjr2IsA+J/Gk7iu5I8t9ateg4a9UNIxhxFmXhCrdDYt8XvdddcxduxYbrrpJsrKyvjVr37FZZddRkNDAyeccALTp0+nY8eOaeshvqzwueeey4ABA3jooYf40pe+1PzQjx/84AdAy0XVTp06Nf/nkuw73/kOkydPpk+fPkyePNmTfIcMGcKkSZPYvn07jz76KGVlZZ73Xnzxxbz22muMHTuWIUOGcOKJJ6bsf926dW2KK5fj9FqmWAqKfzniW//tDn7coeXZn8tjgxh9b3zMdtiM/8s/en8Ddn7csoMvvgAjUsdlw0zLFLdN8pz4sNMyxSIH4czi//WU/xw7vvm1owiOONv7hhWv5CIskYxqU6I3s3PNbLmZrTSzGWm2dzSzZxLb3zazikR9qZk9YWaLzGypmd2R2fBF2sE5JhV5Z1W81nist82oc73lVa9lOSgJyuOPP54XvflDccBEb2bFwCPAecBY4Coz89/+dT2w1Tk3EngQuD9RfxnQ0Tk3Hjge+Jem/wREAlf3IYdby41Qn7gOvOeGe9sMPxUs6ddk8wrYVZujADMnbEO0cugO5bNsS49+ErDSObfaObcP+C1wka/NRcATidfPAWdafF6SA7qYWQnQCdgHpD4tVyQIa97yFKtiI2nwzU+ouGse9Bvnfd/at7MdWUaVlZVRW1urZB8Bzjlqa2tTLhIfSFtm3QwE1iaVa4DJrbVxzjWY2TagF/GkfxGwAegMfM05t8V/ADO7EbgRSPtgX5Gs8CX6d138Bp2URb+GnAgfLWopr30LjvxstqPLmEGDBlFTU8OmTZuCDkUyoKysjEGDBh3Ue7I9vXIS0AgMAHoAb5jZX5xzq5MbOeceAx6D+KybLMckEudL9PNjo9O3GzwZ3nks6X351aMvLS1l2LBhQYchAWrL0M06YHBSeVCiLm2bxDBNN6AW+Cfgj865eufcRuBvQNrpPyI5tXcnsY3eC7ELYkekbzvEN1d6/QJo2JelwEQyry2J/l3gCDMbZmYdgCuBF31tXgSuTbyeCrzm4gOCa4AzAMysC/ApYBkiQft4CUXW8uVxdexwdtDKGufdBkLXgS3lxn2wSf+MJX8cMNE75xqAW4BXgKXA75xzS8zsXjO7MNFsJtDLzFYCXweapmA+ApSb2RLi/2H8yjn3XqZPQuSgfeT9Z/i+q9h/+8OP9pY3LMxsPCJZ1KYxeufcHGCOr+7OpNd7iE+l9L9vZ7p6kcD5Ev2SWMX+2/c/Bla83Or7RcJMd8ZKYdrgS/Ru6P7b9/f36JXoJX8o0UvBaJ422VgPG9/3bHv/QD16/9DNR4sglr3H0olkkhK9FJ7NH8QvqCZ87LqzmW77f0+3QdCpZ0u5fhdsWd16e5EQUaKXwrN5uae4PDa4lYZJzFKHb3zfCkTCSoleCs+mFZ7iSjewlYY+fXxL/W5anr6dSMgo0Uvh8fXo257ofXfOai695Akleik8/h59bEDb3tfnSN9+1KOX/KBELwXFiEHtB566VW7/ib55to6/R7/5A2hsSH2DSMgo0UtBGWiboWFPc3mrK6eWrm17c+ee0KVvS7lxL9R9mOEIRTJPiV4Kykjzrse30g0ArO070Di95CEleikoI229p7wy1sYLsU2U6CUPKdFLQRnhT/QHGJ9PoQuykoeU6CXyKmbMbr6gOrLIO3RzoAuxyfsA1KOXvKRELwUlZeimrXPom6T06FdALNbOqESyS4leCkY3dtLDdjaX97hS1rneB7eTLn2gU4+WcsMnsL0mQxGKZIcSvRSMofaxp7zG9cUd7K+AGfQc4a2rXdXOyESyS4leCkaFL9F/6Pod2o56+RK9VrGUkFOil4IxJFOJ3t+jV6KXkFOil4LhH7rJWI9eQzcSckr0ElnNUyIThhb5x+gPPtFXzJgNPYd7K7co0Uu4KdFLwUjt0fdtpeUB+Hv0W6v1WEEJNSV6KQid2EM/q2suN7gi1rk+h7azsm7QOWlaZuM+2La2nRGKZI8SvRSEIbbRU17velFPyaHvUOP0kkeU6KUgZOxCbJOUcXrNvJHwUqKXguDv0bc/0WuKpeQPJXopCBX2kad8KIneM4unl69Hr6EbCTEleikI/h79oUyt9Ejp0SvRS3gp0UtByPzQja9HX7dGUywltJToJfKMGANss6eu5mBXrfQr6wqderaUG/fBjg3t26dIlijRS+T1pY4O1tLbrnNd2Enn9u+4R4W3vFUPCpdwUqKXSKuYMZuBvt78Qa9B35oeQ73lrdWZ2a9IhinRS+QNylqir/CWleglpJToJfL8PfqaQ136wK+7r0dfp6EbCScleom8QbbJU1aPXgqNEr1EXmqPPltj9OrRSzgp0UvkZW3opttgsKRfoZ0fQf0nmdm3SAa1KdGb2blmttzMVprZjDTbO5rZM4ntb5tZRdK2o83s72a2xMwWmVlZ5sIXSa9luQKXvaGb4lLoNshbV7cmM/sWyaADJnozKwYeAc4DxgJXmdlYX7Prga3OuZHAg8D9ifeWAL8GpjvnxgGnAfUZi17kAHqyg062r7m8y3WkjvLMHcB/QVbj9BJCbenRTwJWOudWO+f2Ab8FLvK1uQh4IvH6OeBMMzPgbOA959xCAOdcrXNO94lLzqQftrHMHUA3TUkeaEuiHwgkPz6nJlGXto1zrgHYBvQCRgHOzF4xswVm9s10BzCzG82s0swqN23alK6JyCHJ2rBNE900JXkg2xdjS4CTgasTf15sZmf6GznnHnPOTXTOTezTJ0MXykRI7dFnPtEP85aV6CWE2pLo1wGDk8qDEnVp2yTG5bsBtcR7//Occ5udc7uBOcBx7Q1apK2yNuOmiW6akjzQlkT/LnCEmQ0zsw7AlcCLvjYvAtcmXk8FXnPOOeAVYLyZdU78B3Aq8H5mQhc5sKwtf9Ak3U1TzmX2GCLtdMBEnxhzv4V40l4K/M45t8TM7jWzCxPNZgK9zGwl8HVgRuK9W4H/Iv6fRRWwwDk3238MkWzJ+hh9l95QmrQS5r6dsHtLZo8h0k4lbWnknJtDfNglue7OpNd7gMtaee+viU+xFMmJ5Ef+ZX3oxizeq9+Y9EV1azV06ZXZ44i0g+6Mlcjqyi662u7m8l5Xyma6ZmTfnufH+odv6qozcgyRTFGil8hKnXHTC5eNf/K6aUpCToleIivrwzZNdNOUhJwSvURW1ufQN9FNUxJySvQSWVmfcdNE69JLyCnRS2TlrEfvH6PfVgONDdk5lsghUKKXyMrZGH2HztAlad+uEbb7bx4XCY4SvURWtnv0+51iqeEbCREleomkTuyht21vLje4Ij6mR0b27UnwTbTmjYSYEr1ESlMSHmC1nvoNrheNFGfvwHp+rISYEr1EUspiZmTnQmxz7149egkxJXqJpJxNrWyiMXoJMSV6iaSczbhpoqEbCTEleomk1ESf5R5910FgSdcAdm2Efbtbby+SQ0r0Ekk5H7opLoFuvkcp163J7jFF2kiJXiIp50M3oHF6CS0leomcUhroS11zOeaMDS4HDwLRzBsJKSV6iYymqY79rZYia3lu60a6U9+2h6m1jy7ISkgp0Uvk+MfnczJsA9C9wlvW0I2EhBK9RE7OVq30S3mkoHr0Eg5K9BI5KXfF5izRpxm6cS59W5EcUqKXyAls6KZLHyjt3FLetwM+2ZqbY4vshxK9RE5gQzdm0H2It27rP3JzbJH9UKKXvJZuyeCB5Piu2GT+KZaaeSMhoEQvkVJMI/19SxTnrEcPuiAroaREL5FyOFsosVhzeZPryh465i4AzaWXEFKil0gJbHy+ScrQTXVujy+ShhK9REpgM26a+Hv0GrqREFCil0gJqkff+pOm1kKsMScxiLRGiV4ixX+zVM579GVdoVPPlnKsHnZsyG0MIj5K9BIpgQ/dQJoLstW5j0EkiRK9RELT0EngF2NBc+kldJToJTKMGAPCkOg1l15CRoleIqMvdXSwlgufW105u+iU+0A0l15CRoleIiN1fD6A3jxoLr2EjhK9REbqA8EDuBALGrqR0GlTojezc81suZmtNLMZabZ3NLNnEtvfNrMK3/YhZrbTzG7LTNgiqVIfCB5Qj77bYMBayjs2QP2eYGIRoQ2J3syKgUeA84CxwFVmNtbX7Hpgq3NuJPAgcL9v+38BL7c/XJHWpfboA0r0JR2g60Bv3ba1wcQiQtt69JOAlc651c65fcBvgYt8bS4Cnki8fg4408wMwMw+D/wDWJKZkEW8mqZWBn6zVDLNpZcQaUuiHwgkd0dqEnVp2zjnGoBtQC8zKwe+BdyzvwOY2Y1mVmlmlZs2bdpfU5FWpQ7dBJjodUFWQiTbF2PvBh50zu3cXyPn3GPOuYnOuYl9+gT4yyl5y4gF96zYdHRBVkKkpA1t1gGDk8qDEnXp2tSYWQnQDagFJgNTzew/ge5AzMz2OOd+2u7IRZL0Zjsdrb65vN11ZjtdggtIc+klRNqS6N8FjjCzYcQT+pXAP/navAhcC/wdmAq85pxzwJSmBmZ2N7BTSV6yITQXYpto6EZC5ICJ3jnXYGa3AK8AxcAvnXNLzOxeoNI59yIwE3jKzFYCW4j/ZyCSNf5nxYZqfB6g5zBvecs/wLn4A8RFcqwtPXqcc3OAOb66O5Ne7wEuO8A+7j6E+ETaJDR3xTYp7welXaB+V7y8bwfs2gTlfYONSwqS7oyVSAjLzVLN3zTMoOdw78baVbkPSAQleomIINeh9w8jNevlS/RbVmc/GJE0lOglEobYRk858IuxAD1HeMtb1KOXYCjRS94rIpbSo1/j+gUUTRIN3UhIKNFL3uvHVjpaQ3N5qysPZA59yhBOL/XoJRyU6CXvDS362FP+0IVkZot/6KZ2dXyKpUiOKdFL3hti3kS/NiyJvrwvdChvKdfvgp0bW28vkiVK9JL3/BdiPwx4fN47xdJ/45SGbyT3lOgl7w01/9BNCC7ENkkZvlGil9xTope85+/Rr4kFn+ibe/W6ICshoEQvecc/uyUl0YdljB7SzKXXTVOSe0r0kte6sose1vK4g72ulI/oEWBEPv4efa0SveSeEr3ktdQZN31wYfpn7b9paoumWEruheg3QuTghW3GTYoufaDDYS3l+l2w46Pg4pGCpEQveW1omMfnIT7FMmVxM12QldxSope85h+6CV2iB+g10lvevCKYOKRgKdFLXgv90A1A79He8iYlesktJXrJa/51bkLZo+8zylvevDyYOKRgKdFLXkmeQ19KA/2p9WwPzTo3yXr7Er169JJjSvSStwbZJoqtZariBteTvXQIMKJW9BoJlvSrtr0G9u5svb1IhinRS94aZhs85VAO2wCUdIQeFd662g8CCUUKkxK95K3hvkS/KtY/oEjaQBdkJUBK9JK3htt6T3m1GxBQJG2gC7ISICV6yVsjinw9+jAn+pQevRK95I4SveQt/9DNahfmoRt/j15DN5I7SvSSlw5jN31sW3N5nyumxvUJMKID8A/dbFkNjfXBxCIFR4le8kLFjNmeOfT+8fkP3eE0UpzrsNqkYsZsKOsG5Ye3VMYaoHZlcEFJQVGil7yUD8M2/gek0G+st/zxktwFIwVNiV7y0vCi8Cf6FP3GectK9JIjSvSSl1KnVuZDoj/KW974fjBxSMFRope8lHqzVDinVnqGb/pq6EaCoUQveaeIGMPM+5SmvOjR9xkNlnTBeNta+KQuuHikYCjRS94ZbBsps5apibXuMOo4bD/vCImSjqnz6TcuDSYWKShK9JJ3jrS1nvKK2OCAIjkEKRdkFwcThxQUJXoJPf80xVG+RL/M5VOi1zi95J4SveSd0UW+Hr0bFFAkh8A/80aJXnKgTYnezM41s+VmttLMZqTZ3tHMnklsf9vMKhL1Z5nZfDNblPjzjMyGL4VotNV4ysvzeuhmCcQag4lFCsYBE72ZFQOPAOcBY4GrzMz3/ZPrga3OuZHAg8D9ifrNwOecc+OBa4GnMhW4FKYO1Kc8cCSvevRdB0Knni3l+l2wWQ8hkexqS49+ErDSObfaObcP+C1wka/NRcATidfPAWeamTnn/tc513RnyxKgk5l1zETgUphG2HpKLNZcrnG92UnnACM6SGYw4Fhv3fr/DSYWKRhtSfQDgeRB0ZpEXdo2zrkGYBvQy9fmUmCBc26v/wBmdqOZVZpZ5aZNm9oauxSg0b4LsXk1bNNk4HHeshK9ZFlOLsaa2Tjiwzn/km67c+4x59xE59zEPn1CvNSsBC6vL8Q2SenRLwgmDikYbUn064DkbtOgRF3aNmZWAnQDahPlQcALwDXOuVXtDVgKm79Hvywfe/QDfD36jxZpbXrJqrYk+neBI8xsmJl1AK4EXvS1eZH4xVaAqcBrzjlnZt2B2cAM59zfMhW0FK4xRWs85Q/ysUfftb93bfqGPbpDVrLqgIk+MeZ+C/AKsBT4nXNuiZnda2YXJprNBHqZ2Urg60DTFMxbgJHAnWZWlfjpm/GzkMhKvlmqD3X0ty3N5b2uJD8TPeiCrORUm8bonXNznHOjnHMjnHPfS9Td6Zx7MfF6j3PuMufcSOfcJOfc6kT9d51zXZxzE5J+NmbvdCRK/HfEHlX0D095qRtCPSW5DOmQpTyERBdkJYd0Z6zkjaNttae8KDY8oEgOjSfZ+3v0NZW5DUYKihK95I3xRd5E/57Lr0QPScl+4PHeDR8v1pLFkjVK9JInHEf7hm7yrUfv0bmn70EkDta+E1g4Em1K9JIX+rGVvtbS493jSvnA+e/byzNDTvSW17wZTBwSeUr0khfG+3rzS1wFjRS30jpPDP20t/zh34OJQyJPiV5CJ2WGCmnG5/N52KaJv0e/fgHU7wkmFok0JXrJC8faSk95UWxYQJFkULeB0H1IS7lxH6ybH1w8EllK9BJ6xTRyXJF3Kd8qNzKgaDJsiH/4RuP0knlK9BJ6R9oayq1lSGOz68pq1z/AiDLIP06/em4wcUikKdFL6J1QtNxTroyNBiyYYDKgYsbslusQw0/zblz7NuzZluuQJOKU6CVU0l2InehL9O/GRuUqnOzrMRR6J51PrAH+MS+4eCSSlOgl5ByTi5Z5auI9+ggZ+RlveeVfgolDIkuJXkJttK2lj7UMZex0ZSxxFcEFlA3+RP/BX8C5YGKRSFKil1A7uWiRp/xWbAwNebJiZZsNPQlKOrWUt9fApuWttxc5SEr0EkpNY/VTihZ76v8WOyqIcLKrtAyGTfHWLU+9ViFyqJToJbQ6UM8k3/j8G7HxAUWTZaPO8ZYXvxBMHBJJSvQSWicWvU9n29tc/th1Z2W+L2TWmjEXgiX9On68CDatCC4eiRQlegmts4q8D+N4rfFY8nn+/H6V94Vhp3jrljwfTCwSOUr0EkpGjM8UL/DU/Sk2MaBosiPlnoFxl3jLi3+v2TeSEUr0EhrJie9oW83htrW5vMt15M3YuCDCyp0xn4OipBlFm1fAhoXBxSORoUQvoXRhsXdt9nmxo9lLh4CiyZHOPWHEmd66+Y8HEopEixK9hE4xjVxY7F3F8aXGE1tpnd9Shm+O+ym1H8MAAAi1SURBVKK3vOhZ2LsjdwFJJCnRS+icVLTYczfsDteJv8SOCzCi7PIk+1HnwmFJK3Pu2wlVT+c+KIkUJXoJnSuLvUv1vtw4KfrDNk2KS+FYX6/+zZ9CY30w8UgkKNFLqPSnlrN90yqfj01ppXVEnXADFHdsKW9bE5+BI3KIlOglcMlDF18o+TMlFmsuL4sN5q3YmCDCyinP8M1h/eDYq70NXv8BNOzLbVASGUr0Ehpd2ckXi//sqXuy8Wwie5NUGs0J/9NfAStu2bC1GipnBhKT5D8legmN60tepqt90lze4sr5n8aTAowoQD2Hw/HTvHVzfwA7PgokHMlvSvQSCp+a8SQ3FM/x1P284QJ2UxZQRLmXMtXytBnQobylvHcbzP6G7paVg6ZEL4FJTmzfLn2KLkkLmG1x5TzReHYQYYVHeV84/d+8dctegspfBhOP5C0legnc54v+yvnF73jqftRweUH15pN5evaTp8MA3z0Ef5wB1X/NbVCS15ToJeeSE9k5d/yM75Z6e6gLY8P5TeMZuQ4rVJr/joqK4ZLHvEM4jfvgN/8Ea98NJjjJO0r0Epgh9jFPdLifctvTXLfXlfLN+huJ6Z9mi95HwIUPe+v2boMnPgfL5qR/j0iSiD18U8Kuqac69Y4f8XyHh+ht2z3b7234IsvdkCBCC53kbz7V910an2L56r0tDRo+gd9eBROvh7PugY6H5T5IyQvqNklO9WIbd5U8wTMdvpOS5B9vOJtZjWe28s7CVjFjNpz8dTjlm6kbK2fCjyfAXx+C3VtyH5yEnrmQTdWaOHGiq6ysPHBDyQsVM2bTjZ1MLFrOzInr2L3wBc/jAZv8vnEKt9f/i4Zs2qD60o9ofOkbFFua392iUhh5Jgw/DYaeBH2OhJICWSeowJnZfOdc2qfztGnoxszOBX4MFAO/cM7d59veEXgSOB6oBa5wzlUntt0BXA80Al91zr1yiOfRuh0fw/Z1SRW+X4CU34ekipT/6PzvzeT2bO77YI/te2t7j924F/Zsjy+pu3cH7N4MdWt4tcNCRhRtiLd5Dzqnucn1pw0X8UDD5RTSHbDtcsL1XP78Nn7f/ynYstq7LVYPK/4Y/4H43bU9h0H3IdClT+KnN5R2gdIyKOmU+LMs/szalB/zlrGQfUyhCib+99VeRSVw+Pj27yfJAXv0ZlYMrADOAmqAd4GrnHPvJ7W5GTjaOTfdzK4ELnbOXWFmY4HfAJOAAcBfgFHOucbWjndIPfq//Rj+fOfBvUcCV+N6c0f9DbwROzroUPJSGXv5cskf+Ernv8SXM5Zo6NIXbv/goN+2vx59W74nTwJWOudWO+f2Ab8FLvK1uQh4IvH6OeBMM7NE/W+dc3udc/8AVib2JwVsg+vJPfVf5Iy9P1KSb4c9dORHDZfDvy7ih/WXsyrW/8BvkoLUlqGbgcDapHINMLm1Ns65BjPbBvRK1L/le+9A/wHM7EbgxkRxp5ktb1P06fUGNrfj/fkmD893O/BI4ueQ5OE5t1ur52z3x/9Mc5k23xXo57x9M988pCGgoa1tCMX0SufcY8BjmdiXmVW29vUligrtfEHnXCh0zpnTlqGbdcDgpPKgRF3aNmZWAnQjflG2Le8VEZEsakuifxc4wsyGmVkH4ErgRV+bF4FrE6+nAq+5+FXeF4ErzayjmQ0DjgDeQUREcuaAQzeJMfdbgFeIT6/8pXNuiZndC1Q6514EZgJPmdlKYAvx/wxItPsd8D7QAHx5fzNuMiQjQ0B5pNDOF3TOhULnnCGhu2FKREQyS7chiohEnBK9iEjE5VWiN7NfmtlGM1ucVNfTzP5sZh8k/uyRqDcze9jMVprZe2Z2XOt7Dq9WzvmHZrYscV4vmFn3pG13JM55uZmdE0zU7ZPunJO2fcPMnJn1TpQj+zkn6r+S+KyXmNl/JtVH8nM2swlm9paZVZlZpZlNStTn/edsZoPNbK6ZvZ/4PG9N1Gc/hznn8uYHOAU4DlicVPefwIzE6xnA/YnXnwVeJr4YxqeAt4OOP4PnfDZQknh9f9I5jwUWAh2BYcAqoDjoc8jEOSfqBxOfFPAh0LsAPufTiS8b0jFR7hv1zxn4E3Be0mf7elQ+Z6A/cFzi9WHEl5YZm4scllc9eufcPOKzepIlL7/wBPD5pPonXdxbQHczy7t7xNOds3PuT865hkTxLeL3J0BElpxo5XMGeJD4DaDJMwgi+zkDNwH3Oef2JtpsTNRH+XN2QNfE627A+sTrvP+cnXMbnHMLEq93AEuJrxSQ9RyWV4m+Ff2cc4nlEfkI6Jd4nW7phpTlFyLgOuL/60OEz9nMLgLWOecW+jZF9pyBUcAUM3vbzP6fmZ2QqI/yOf8r8EMzWws8ANyRqI/UOZtZBXAs8DY5yGFRSPTNXPz7TsHMFzWzfyd+f8KsoGPJJjPrDPwbUGhLlJYAPYl/bb8d+F1iscAouwn4mnNuMPA14vfoRIqZlQO/B/7VOed5+k62clgUEv3HTV9nEn82fb2N9PILZjYNuAC4OvGPA6J7ziOIj0UvNLNq4ue1wMwOJ7rnDPEe3POJr+7vADHiC31F+ZyvBZ5PvH6WliGpSJyzmZUST/KznHNN55n1HBaFRJ+8/MK1wB+S6q9JXLn+FLAt6etRXrP4g2C+CVzonNudtCmSS0445xY55/o65yqccxXEE+BxzrmPiPDnDPwP8QuymNkooAPx1Rwj+TknrAdOTbw+A2hamD3vP+fEt7GZwFLn3H8lbcp+Dgv6SvRBXrX+DbABqCf+y3498eWQXyX+D+IvQM9EWyO+Du4qYBEwMej4M3jOK4mP3VUlfh5Nav/viXNeTmL2Qr79pDtn3/ZqWmbdRPlz7gD8GlgMLADOiPrnDJwMzCc+q+ht4PiofM6Jc3PAe0m/u5/NRQ7TEggiIhEXhaEbERHZDyV6EZGIU6IXEYk4JXoRkYhTohcRiTglehGRiFOiFxGJuP8P42kqYw6llaIAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Chapter 8, problem 7, part a\n",
        "import numpy as np\n",
        "from scipy.stats import norm\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "n = 50\n",
        "theta = 1\n",
        "x = np.random.uniform(0, theta, n)\n",
        "\n",
        "B = 10000\n",
        "t_boot = np.empty(B)\n",
        "for i in range(B):\n",
        "  sel_data = np.random.choice(x, size=n, replace=True)\n",
        "  t_boot[i] = np.max(sel_data)\n",
        "\n",
        "# true ditribution of theta_hat, find in Q2 of chapter 6 \n",
        "y = np.linspace(0.1, theta, 100)\n",
        "f_y = (n/theta**n) * y**(n-1) \n",
        "\n",
        "plt.hist(t_boot, 10, density=True, label='bootstrap dist')\n",
        "plt.plot(y, f_y, lw=2, label='true distribution')\n",
        "plt.legend(loc='upper left')\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "cMx65jwfV5pe",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "outputId": "2f081542-b370-46f1-99e0-8aaa21ce43fb"
      },
      "execution_count": 223,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAb8ElEQVR4nO3de3hV9b3n8fc3IRAUAhQjAyIGKiD3oOHSIojX4qVe0eqDByIiosfWOR4ZOfUZhLHT2jPOkbZKGSsVdNCD2Dp6RNtRhCe2XiAilChSgYkaRIGUm+WWy3f+2DubhCRks29r7/B5PQ9P9rpkrS8L+OTHb/3Wb5m7IyIimScr6AJERCQ2CnARkQylABcRyVAKcBGRDKUAFxHJUG1SebLTTjvNCwoKUnlKEZGM98EHH+xy9/xj16c0wAsKCigtLU3lKUVEMp6ZfdbUenWhiIhkKAW4iEiGUoCLiGSolPaBN6WqqoqKigoOHToUdCmSALm5ufTs2ZOcnJygSxFp9QIP8IqKCjp27EhBQQFmFnQ5Egd3p7KykoqKCnr37h10OSKtXuBdKIcOHaJr164K71bAzOjatav+NyWSIlG1wM2sHNgP1ADV7l5kZt8ClgIFQDlwk7vvjqUIhXfroT9LkdQ5kRb4he5e6O5F4eVZwAp37wusCC+LiEh9X38En74BeysSfuh4ulCuARaHPy8Gro2/nNTbs2cP8+fPT9n5CgoK2LVrFwDf/e53j7vvT3/60+Nuv+KKK9izZw/l5eUMHjz4hOpYtWoV77zzTmR5wYIFPPPMMyd0DBGJwpqFsGQibHo94YeONsAd+L9m9oGZTQ+v6+bu28OfvwK6NfWNZjbdzErNrHTnzp1xlpt4xwvw6urqpJ67foA2pbkAd3dqa2t57bXX6Ny5c0znPjbAZ8yYweTJk2M6logcR9WB0Ne2pyb80NEG+Pnufi5wOfCPZjau/kYPvdanyVf7uPuT7l7k7kX5+Y0e5Q/crFmz2LJlC4WFhcycOZNVq1YxduxYrr76agYOHNiodfvoo48yZ84cALZs2cKECRM477zzGDt2LJ988kmj41dWVnLZZZcxaNAgpk2bRv03IHXo0AGA7du3M27cOAoLCxk8eDBvv/02s2bN4uDBgxQWFjJp0iTKy8vp378/kydPZvDgwXzxxRcNWvPV1dVMmjSJAQMGMHHiRA4cCP2lqb9PaWkp48ePp7y8nAULFvDYY49RWFjI22+/zZw5c3j00UcBWLduHaNHj2bo0KFcd9117N4durUxfvx4HnjgAUaOHEm/fv14++23E/ynIdIKHfl76GtO+4QfOqqbmO6+Lfx1h5m9BIwEvjaz7u6+3cy6AzvirmZOp7gP0fRx9za76ZFHHqGsrIx169YBoZbp2rVrKSsro3fv3pSXlzf7vdOnT2fBggX07duX999/n7vvvpu33nqrwT5z587l/PPPZ/bs2SxfvpyFCxc2Os5zzz3H9773PR588EFqamo4cOAAY8eO5fHHH4/UVV5ezqeffsrixYsZPXp0o2Ns2rSJhQsXMmbMGKZOncr8+fO5//77m6y7oKCAGTNm0KFDh8g+K1asiGyfPHkyv/rVr7jggguYPXs2c+fOZd68eUDoB8Xq1at57bXXmDt3Lm+++Waz10dEgKqDoa85AbTAzexUM+tY9xm4DCgDXgGmhHebAryc8OoCMnLkyBbHMX/zzTe888473HjjjRQWFnLnnXeyffv2RvuVlJRw6623AnDllVfSpUuXRvuMGDGCp59+mjlz5rBhwwY6duzY5DnPOuusJsMb4Mwzz2TMmDEA3HrrrfzpT386bv3N2bt3L3v27OGCCy4AYMqUKZSUlES2X3/99QCcd955x/3hJiJhdV0oAbXAuwEvhYeHtQGec/c/mNka4AUzux34DLgp7mqO01JOpVNPPfqTsk2bNtTW1kaW68Y419bW0rlz50gLOR7jxo2jpKSE5cuXU1xczH333ddkf3T9uo517PC9uuX69SdifHa7du0AyM7OTvo9ApFWIdIHfkrCD91iC9zdt7r7sPCvQe7+38PrK939Ynfv6+6XuPvfEl5dCnTs2JH9+/c3u71bt27s2LGDyspKDh8+zKuvvgpAXl4evXv3ZtmyZUDoxuL69esbff+4ceN47rnnAHj99dcj/cn1ffbZZ3Tr1o077riDadOmsXbtWgBycnKoqqqK6vfx+eef8+677wKhLpnzzz8fCHWXfPDBBwD87ne/a/H33alTJ7p06RLp33722WcjrXERicGRuhZ4AAHe2nXt2pUxY8YwePBgZs6c2Wh7Tk4Os2fPZuTIkVx66aWcc845kW1Llixh4cKFDBs2jEGDBvHyy417kR566CFKSkoYNGgQv//97+nVq1ejfVatWsWwYcMYPnw4S5cu5d577wVCfexDhw5l0qRJLf4++vfvzxNPPMGAAQPYvXs3d911V+T89957L0VFRWRnZ0f2//73v89LL70UuYlZ3+LFi5k5cyZDhw5l3bp1zJ49u8Xzi0gzIn3giQ9wqz8qItmKior82Bc6bNy4kQEDBqSsBkk+/ZmK1POvfeBAJdy/GTrENhLPzD6o9xBlxEnfAhcRSaq6FngQfeAiIhKj2tqjNzHbJH4UigJcRCRZqsMjv9q0h6zEx60CXEQkWZI4BhwU4CIiyVP3GH0S5kEBBbiISPJEhhAmpwUe+CvVjlUwa3lCj1f+yJXH315ezlVXXUVZWVnc55o3bx7Tp0/nlFOav9u8aNEiLrvsMnr06BH3+VpSXFzMVVddxcSJE5k2bRr33XcfAwcODLwukZNGVd1EVokfgQJqgSfUvHnzIrMANmfRokV8+eWXTW6rqalJRlkAPPXUU82GNxy/LhGJURIf4gEFOND8VKwrVqxg+PDhDBkyhKlTp3L48OFm1//yl7/kyy+/5MILL+TCCy+kpqaG4uJiBg8ezJAhQ3jsscd48cUXKS0tZdKkSRQWFnLw4EEKCgp44IEHOPfcc1m2bBm/+c1vGDFiBMOGDeOGG26I1FJcXMyMGTMoKiqiX79+kUf663N37rnnHvr3788ll1zCjh1HJ4gcP348paWlUdclIglwJHnzoIACHAhNxXr33XezceNG8vLymD9/PocOHaK4uJilS5eyYcMGqqur+fWvf93s+h/96Ef06NGDlStXsnLlStatW8e2bdsoKytjw4YN3HbbbUycOJGioiKWLFnCunXraN8+1C/WtWtX1q5dy80338z111/PmjVrWL9+PQMGDGgw/Wx5eTmrV69m+fLlzJgxo9HkVC+99BKbNm3i448/5plnnmnyhREnUpeIxKkqefOggAIcaHoq1k2bNtG7d2/69esHHJ1Wtbn1x+rTpw9bt27lhz/8IX/4wx/Iy8tr9vw/+MEPIp/LysoYO3YsQ4YMYcmSJXz00UeRbTfddBNZWVn07duXPn36NHqBRElJCbfccgvZ2dn06NGDiy66KK66RCROCvDka24q1nh06dKF9evXM378eBYsWMC0adOa3bf+NLHFxcU8/vjjbNiwgYceeqhBKzsRdZ5IXSISpyROJQsKcKDpqVj79+9PeXk5mzdvBo5Oq9rcemg4ReuuXbuora3lhhtu4Cc/+UlkitiWpq/dv38/3bt3p6qqiiVLljTYtmzZMmpra9myZQtbt26lf//+DbaPGzeOpUuXUlNTw/bt21m5cmWj48dal4jEIIlTyUIaDiNsadhfMtRNxTp16lQGDhzIXXfdRW5uLk8//TQ33ngj1dXVjBgxghkzZtCuXbsm10No+tcJEybQo0cP5s2bx2233RZ5mcLPfvYz4OjNyPbt20d+aNT38MMPM2rUKPLz8xk1alSDUO3VqxcjR45k3759LFiwgNzc3Abfe9111/HWW28xcOBAevXqxXe+851Gx9+2bVtUdakfXCQBkjwKRdPJZoj6Y7rTnf5MRcLemA1//gVc/BCMvS/mw2g6WRGRVItMJZucR+nTrgtFmrZo0aKgSxCRE3XkJJjMKpXdOJJc+rMUqae1P0qfm5tLZWWl/uG3Au5OZWVlo5urIietJN/EDLwLpWfPnlRUVLBz586gS5EEyM3NpWfPnkGXIZIeItPJttIAz8nJoXfv3kGXISKSeJrMSkQkQ+lRehGRDKVXqomIZKjIdLJ6pZqISGZJ8ivVFOAiIsng3vrHgYuItEo1R8BrISsHsnOScgoFuIhIMiR5DDgowEVEkiMyAiU5NzBBAS4ikhxJvoEJJxDgZpZtZh+a2avh5d5m9r6ZbTazpWbWNmlViohkmjTrQrkX2Fhv+efAY+5+NrAbuD2RhYmIZLQkP0YPUQa4mfUErgSeCi8bcBHwYniXxcC1yShQRCQjJXkIIUTfAp8H/BegNrzcFdjj7tXh5QrgjKa+0cymm1mpmZVqxkEROWmkQwvczK4Cdrj7B7GcwN2fdPcidy/Kz8+P5RAiIpkn8hh98gI8mulkxwBXm9kVQC6QB/wC6GxmbcKt8J7AtqRVKSKSaZI8kRVE0QJ3939x957uXgDcDLzl7pOAlUDdK9KnAC8nrUoRkUyT5uPAHwDuM7PNhPrEFyamJBGRViAFLfATeiOPu68CVoU/bwVGJr4kEZFWIAV94HoSU0QkGdJhFIqIiMQgjcaBi4jIiTiS3PdhggJcRCQ56rpQ1AcuIpJhIl0oaTAboYiInIDITcz0HAcuIiLNOZIGT2KKiEgM6h7kaasWuIhIZkmHuVBERCQGVRpGKCKSmY6oC0VEJPPUVEFtFVgWZCfvdcEKcBGRRKs/laxZ0k6jABcRSbQUDCEEBbiISOJVJX8qWVCAi4gkXgpGoIACXEQk8VIwFzgowEVEEu9I8ieyAgW4iEjiRaaSTd4YcDjBd2KKiEhIwazlzW4rn1TXhaIWuIhIZqk/DjyJFOAiIommceAiIhlK48BFRDKUxoGLiGSoFLyRHhTgIiKJl4KXOYACXEQk8VLwOjVQgIuIJF7kSUx1oYiIZJZDe0Nfczsl9TQKcBGRRFOAi4hkqEN7Ql/bd07qaRTgIiKJFmmBBxzgZpZrZqvNbL2ZfWRmc8Pre5vZ+2a22cyWmlny3twpIpIhjFo4tC+00C4vqeeKpgV+GLjI3YcBhcAEMxsN/Bx4zN3PBnYDtyevTBGRzNCRg4BD246QndwJX1sMcA/5JryYE/7lwEXAi+H1i4Frk1KhiEgGybPwGPAk38CEKPvAzSzbzNYBO4A3gC3AHnevDu9SAZyRnBJFRDJHJ8JjwJN8AxOiDHB3r3H3QqAnMBI4J9oTmNl0Mys1s9KdO3fGWKaISGbIs3CAp0sLvI677wFWAt8BOptZXQdPT2BbM9/zpLsXuXtRfn5+XMWKiKS7PNKoC8XM8s2sc/hze+BSYCOhIJ8Y3m0K8HKyihQRyRRHW+DJ70KJ5hZpd2CxmWUTCvwX3P1VM/sY+Hcz+wnwIbAwiXWKiGSEPFLXhdJigLv7X4DhTazfSqg/XEREwiKjUNLlJqaIiESnUwpb4ApwEZEESrtx4CIiEp2jfeDqQhERySid0nUcuIiIHF9ajQMXEZHoaRSKiEiG0igUEZEMlEM1p9hhsGxo2yHp51OAi4gkSMf6/d9mST+fAlxEJEFSOQIFFOAiIgmTl8K5wEEBLiKSMKl8ChMU4CIiCZPKESigABcRSZijLXB1oYiIZBS1wEVEMlQq34cJCnARkYSJzIOiUSgiIpmlUwrfhwkKcBGRhEnlXOCgABcRSRiNAxcRyVCpfCM9KMBFRBImlXOBgwJcRCRB/GgLvF1eSs6oABcRSYD2HKat1XDIcyAnNyXnVICLiCRA3RjwvZyasnMqwEVEEqCu/3ufK8BFRDJKJ74B1AIXEck4R1vgp6TsnApwEZEEqOsD34cCXEQko9TNg7JXfeAiIpmlbgz4PvWBi4hkFvWBi4hkqLq38aTVKBQzO9PMVprZx2b2kZndG17/LTN7w8w+DX/tkvxyRUTSU7qOA68G/tndBwKjgX80s4HALGCFu/cFVoSXRUROSp0tNA48rUahuPt2d18b/rwf2AicAVwDLA7vthi4NllFioiku3z2ALDDUzMTIZxgH7iZFQDDgfeBbu6+PbzpK6BbM98z3cxKzax0586dcZQqIpK+Trc0DnAz6wD8DvjP7r6v/jZ3d8Cb+j53f9Ldi9y9KD8/P65iRUTS0SkcoqMd5JDnpN8wQjPLIRTeS9z99+HVX5tZ9/D27sCO5JQoIpLeTrfdQF3r21J23mhGoRiwENjo7v9Wb9MrwJTw5ynAy4kvT0Qk/XUL939/TWoH47WJYp8xwD8AG8xsXXjdj4FHgBfM7HbgM+Cm5JQoIpLeGrbAU6fFAHf3P9H8/wkuTmw5IiKZ52iAp7YFricxRUTidHQEigJcRCSjdAu3wL9OcReKAlxEJE6n1z3Ek+KbmApwEZE4HW2BK8BFRDJKfgBPYYICXEQkLu05RJ4d5LDnpHQqWVCAi4jEpeEcKKl7ChMU4CIicelGuP87xTcwQQEuIhKXuhZ4qocQggJcRCQu3QJ6ChMU4CIicckP6ClMUICLiMSlW0ATWYECXEQkLqcHNJUsKMBFROKiFriISIY6OgpFLXARkYyRy2Hy7EAgT2GCAlxEJGZBPoUJCnARkZjVPYW5g9T3f4MCXEQkZkH2f4MCXEQkZkGOQAEFuIhIzIJ6F2YdBbiISIwib6NXH7iISGaJTCWrFriISGbpbpVAcAHeJpCziohkiIJZy5tcn0M1vWwHtW585t1SXFWIWuAiIjHoZV/Txmr5wvM5TNtAalCAi4jE4Gz7EoAt3iOwGhTgIiIx+LYCXEQkM307SwEuIpKRIi3wWgW4iEgGcfqoC0VEJPPks4c8O8hu78Df6BhYHS0GuJn91sx2mFlZvXXfMrM3zOzT8NdgRrGLiATg7Ab936mfB7xONC3wRcCEY9bNAla4e19gRXhZROSkkA793xBFgLt7CfC3Y1ZfAywOf14MXJvgukRE0tbRIYTdA60j1j7wbu6+Pfz5K6DZ50jNbLqZlZpZ6c6dO2M8nYhI+kiHMeCQgJuY7u6AH2f7k+5e5O5F+fn58Z5ORCRwfbJC7ddMDfCvzaw7QPjrjsSVJCKSvtpziJ62iyOezRd+eqC1xBrgrwBTwp+nAC8nphwRkfTW274CoNz/EzVkB1pLNMMInwfeBfqbWYWZ3Q48AlxqZp8Cl4SXRURavXTp/4Yo5gN391ua2XRxgmsREUl76TAHSh09iSkicgLqWuBba4MdQggKcBGRE5IO84DXUYCLiETpFA7R1yqo9iw2+xlBl6MAFxGJ1nlZf6WN1VLmBfyd9kGXowAXEYnWyKxPAHi/dkDAlYQowEVEojQqayMAq2vPCbiSEAW4iEgU2nGEYbaFWjfWKMBFRDLHuVmf0s6q2ei92MepQZcDKMBFRKJS132SLv3foAAXEYnKKKu7gZke3SegABcRaVFbqhie9SmQPjcwQQEuItKiYbaFXKtiU21PdpMXdDkRCnARkRak2/jvOgpwEZEWpOMNTFCAi4gc16kcZETWJgBW1/YPuJqGFOAiIsdxZfZ7tLcjrK7tz066BF1OAwpwEZHjuDl7JQAv1IwPtpAmtPhGHhGR1qJg1vJmt5U/cmXjlTs+4dyszez39iyvGZXEymKjFriISHM+fBaA/6gZzUFyAy6mMQW4iEhTqo/A+ucBWFpzYcDFNE0BLiLSlL++Dgcq+aT2TNb7t4OupkkKcBGRpqwNdZ+Ebl5aoKU0RwEuInKsbWth85uQ3ZaXasYEXU2zFOAiIvXV1sCr/wQ4jLozreY+OZYCXESkvjULYfs6yOsJF8wKuprjUoCLiNTZ/xW89XDo8+U/h3Ydgq2nBQpwEZE6f3wQDu+DfhPgnCYe7EkzCnAREYCSR6HsRWjTHi7/V7D0HHlSnx6lF5GT3vTs/4C3ngcMrv4VdDkr6JKiogAXkbR2wvOXnKDbs1/jxznh8L7mCRh6Y9zHTBUFuIiclE5nN3NzFnF59prQiu/Pg+GTAq3pRKkPXEROKm2pYlL2m7zZbiaXZ6/hG8/l/qo74bzioEs7YWqBi0jr5w47N/Fgm//NDdklfMu+AeDNmuH816qpbKcrjwZcYiziCnAzmwD8AsgGnnL3RxJSlYicdI7X131C3OHvu6ByM+zcCJ+9C5/9GfZt445w4n1Uexbzq69hee0o0nWek2jEHOBmlg08AVwKVABrzOwVd/84UcVFbF8PXpvww4pIDNyj2anJj5EF99Dnel9H2Cdk4WRZLYaTTS3Z1NCGWtpQQw7V5FBNrlWRyxFyOQxvrIbD++HQHvhmB/x9J+zbDof3Ni6pfRee31/IczUXs8F7k8nBXSeeFvhIYLO7bwUws38HrgESH+C/nQBVBxJ+WBFJH8vaxfBNf25mfbs86Ho2nNYXziiCgjGQP4B/+fHr8ZSYdsyj+mnaxDeaTQQmuPu08PI/AKPc/Z5j9psOTA8v9gc2xV5uWjgN2BV0EWlC16IhXY+GdD2OivdanOXu+ceuTPpNTHd/Engy2edJFTMrdfeioOtIB7oWDel6NKTrcVSyrkU8wwi3AWfWW+4ZXiciIikQT4CvAfqaWW8zawvcDLySmLJERKQlMXehuHu1md0D/JHQMMLfuvtHCassfbWa7qAE0LVoSNejIV2Po5JyLWK+iSkiIsHSo/QiIhlKAS4ikqEU4M0wswlmtsnMNptZoxfjmdl9Zvaxmf3FzFaYWWZMIByDlq5Fvf1uMDM3s1Y9dCya62FmN4X/fnxkZs+lusZUiuLfSi8zW2lmH4b/vVwRRJ2pYGa/NbMdZlbWzHYzs1+Gr9VfzOzcuE7o7vp1zC9CN2W3AH2AtsB6YOAx+1wInBL+fBewNOi6g7oW4f06AiXAe0BR0HUH/HejL/Ah0CW8fHrQdQd8PZ4E7gp/HgiUB113Eq/HOOBcoKyZ7VcArxN6jn808H4851MLvGmRaQLc/QhQN01AhLuvdPe65/vfIzQOvjVq8VqEPQz8HDiUyuICEM31uAN4wt13A7j7jhTXmErRXA8H8sKfOwFfprC+lHL3EuBvx9nlGuAZD3kP6Gxm3WM9nwK8aWcAX9Rbrgiva87thH6qtkYtXovwfwPPdPcETSeX1qL5u9EP6Gdmfzaz98KzdrZW0VyPOcCtZlYBvAb8MDWlpaUTzZbj0nzgcTKzW4Ei4IKgawmCmWUB/wYUB1xKOmlDqBtlPKH/mZWY2RB33xNoVcG5BVjk7v/TzL4DPGtmg901xWi81AJvWlTTBJjZJcCDwNXufjhFtaVaS9eiIzAYWGVm5YT69V5pxTcyo/m7UQG84u5V7v7/gL8SCvTWKJrrcTvwAoC7vwvkEprc6WSU0ClIFOBNa3GaADMbDvwvQuHdmvs4j3st3H2vu5/m7gXuXkDofsDV7l4aTLlJF80UEv+HUOsbMzuNUJfK1lQWmULRXI/PgYsBzGwAoQDfmdIq08crwOTwaJTRwF533x7rwdSF0gRvZpoAM/tvQKm7vwL8D6ADsMzMAD5396sDKzpJorwWJ40or8cfgcvM7GOgBpjp7pXBVZ08UV6PfwZ+Y2b/ROiGZrGHh2S0Nmb2PKEf3qeF+/wfAnIA3H0BoXsAVwCbgQPAbXGdr5VeRxGRVk9dKCIiGUoBLiKSoRTgIiIZSgEuIpKhFOAiIhlKAS4ikqEU4CIiGer/Axv4X9I/yKyOAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}