Confidence Intervals with Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "exceptional-vertical",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "In a series of weekly articles, I will be covering some important topics of statistics with a twist.\n",
    "\n",
    "The goal is to use Python to help us get intuition on complex concepts, empirically test theoretical proofs, or build algorithms from scratch. In this series, you will find articles covering topics such as random variables, sampling distributions, confidence intervals, significance tests, and more.\n",
    "\n",
    "At the end of each article, you can find exercises to test your knowledge. The solutions will be shared in the article of the following week.\n",
    "\n",
    "Articles published so far:\n",
    "- [Bernoulli and Binomial Random Variables with Python](https://towardsdatascience.com/bernoulli-and-binomial-random-variables-d0698288dd36)\n",
    "- [From Binomial to Geometric and Poisson Random Variables with Python](https://towardsdatascience.com/geometric-and-poisson-random-variables-with-python-e5dcb63d6b55)\n",
    "- [Sampling Distribution of a Sample Proportion with Python](https://towardsdatascience.com/sampling-distributions-with-python-f5a5f268f636)\n",
    "\n",
    "As usual the code is available on my [GitHub](https://github.com/luisroque/College_Statistics_with_Python)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "varying-dakota",
   "metadata": {},
   "source": [
    "# Confidence Intervals and Margin of Error"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "liquid-brazil",
   "metadata": {},
   "source": [
    "A group of students was trying to understand which character of *friends* was the funniest and they narrowed it down to Chandler Bing or Ross Geller. They were interested in the likelihood of Ross actually winning, since he has a very awkward sense of humor compared to the sarcastic and witty sense of humor of Chandler. They decided to setup a poll in their school. Ideally, they would ask the entire population, but shortly they understood that it would not be feasible to ask the 5,000 students their preference. Instead they decided to take random samples of 40 students and calculate the sample proportion that supported Ross. The first value that they arrived at was $\\hat p=0.61$.\n",
    "\n",
    "Notice that to continue with their study, the group of students needs to make sure that the conditions to calculate a valid confidence interval for a proportion are met. There are 3 conditions:\n",
    "* The sample has to be random\n",
    "* The sampling distribution of the sample proportions can be approximated by a normal distribution. The rule of thumb is that you need to have at least 10 successes and 10 failures\n",
    "* The samples are required to be independent. The rule of thumb is that if you are sampling without replacement your sample size should be less than 10% the population size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "animated-outdoors",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import math\n",
    "from scipy.stats import bernoulli, norm, t, skewnorm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "attractive-beatles",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_hat = 0.61\n",
    "n = 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "distinguished-character",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normal conditions:\n",
      "successes >= 10: True\n",
      "failures >= 10 : True\n",
      "--\n",
      "Independence condition:\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "print('Normal conditions:')\n",
    "print('successes >= 10: ' + str(n*p_hat >= 10))\n",
    "print('failures >= 10 : ' + str(n*(1-p_hat) >= 10))\n",
    "print('--')\n",
    "print('Independence condition:')\n",
    "print(40/5000 < 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "atomic-norwegian",
   "metadata": {},
   "source": [
    "All the conditions are met, the group of students could now focus on building the sampling distribution of the sample proportions; from this distribution they calculated the possible sample proportions they could get and their likelihoods. \n",
    "\n",
    "We already saw that the mean of the sampling distribution is the actual population proportion $p$ and the standard deviation of the sample proportions $\\sigma_{\\hat p}=\\sqrt{\\frac{p(1-p)}{n}}$.\n",
    "\n",
    "Let's start to connect these concepts with the concept of confidence intervals. What is the probability that $\\hat p = 0.61$ is within $2\\sigma_{\\hat p}$ of $p$? For a normal distribution this is approximately 95%. This is equivalent to say that there is a 95% probability that $p$ is within $2\\sigma_{\\hat p}$ of $\\hat p$. And this is the fundamental idea for a confidence interval.\n",
    "\n",
    "Now, we have a problem. We do not know $p$, so we need to use an estimate. The best estimate that we have is naturally $\\hat p$. Thus, instead of using $\\sigma_{\\hat p}$ we use the Standard Error $\\text{SE}_{\\hat p} =\\sqrt{\\frac{\\hat p(1- \\hat p)}{n}}$ to compute our confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "athletic-literature",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With 95% confidence between 0.46 and 0.76 of students prefer the awkward humor of Ross.\n"
     ]
    }
   ],
   "source": [
    "SE_hat_p = np.sqrt(p_hat*(1-p_hat)/n)\n",
    "print(f'With 95% confidence between {np.round(p_hat - 2*SE_hat_p, 2)} and {np.round(p_hat + 2*SE_hat_p, 2)} of students prefer the awkward humor of Ross.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "short-banking",
   "metadata": {},
   "source": [
    "Notice that our confidence interval calculated above can change based on what sample proportion we actually choose. If the group of students sample again 40 new students it could happen that the new sample proportion is now 0.55."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "suburban-newcastle",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With 95% confidence between 0.39 and 0.71 of students prefer the awkward humor of Ross.\n"
     ]
    }
   ],
   "source": [
    "p_hat = 0.55\n",
    "n = 40\n",
    "SE_hat_p = np.sqrt(p_hat*(1-p_hat)/n)\n",
    "print(f'With 95% confidence between {np.round(p_hat - 2*SE_hat_p, 2)} and {np.round(p_hat + 2*SE_hat_p, 2)} of students prefer the awkward humor of Ross.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daily-billy",
   "metadata": {},
   "source": [
    "We can connect this concept with margin of error. The margin of error for our first trial (knowing that we are interested in getting 95% confidence) is 2 times our $\\text{SE}_{\\hat p}$.\n",
    "\n",
    "An interesting question that often arises is \"what can you do to reduce the margin of error\"? Notice that the margin of error is dependent on $\\text{SE}_{\\hat p}$, which is inversely proportional to the sample size. So, one possible way to reduce the margin of error is to increase the sample size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "magnetic-survey",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Margin of error = 0.15732132722552272\n"
     ]
    }
   ],
   "source": [
    "print(f'Margin of error = {2*SE_hat_p}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "comfortable-demographic",
   "metadata": {},
   "source": [
    "Again, the same reasoning applies. Depending on what our sample proportion is, our margin of error could have a different value.\n",
    "\n",
    "##### The idea is that if we use this method of computing confidence intervals over and over again, it is going to produce different intervals each time (depending on the sample proportion) that include the true proportion 95% of the time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "collaborative-refund",
   "metadata": {},
   "outputs": [],
   "source": [
    "confidence_interval=0.95\n",
    "p = 0.61\n",
    "n = 50\n",
    "number_trials = 25\n",
    "\n",
    "p_hat_list = []\n",
    "SE_hat_p_list = []\n",
    "\n",
    "for i in range(number_trials):\n",
    "    s_ = bernoulli.rvs(p=p, size=n)\n",
    "    p_hat = s_[s_==1].shape[0] / s_.shape[0]\n",
    "    p_hat_list.append(p_hat)\n",
    "    SE_hat_p_list.append(2*np.sqrt(p_hat*(1-p_hat)/n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "solar-observation",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAHSCAYAAADmGtAuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWMklEQVR4nO3df6zd933X8de7PywBrWhGWmPSWi5rCassmo1L1zJAHtWgdf6oKg1lBbXVYuEh6LRI/SNV/+gmTShFgrUgICNbozJpa4dpRsMUNqLCXZnWNrvBWRM7MjWltKFXdgqFRbYsK/aHP+5xfOde+5x77znfc87nPh6Sde89P+75+Cv7mW++Pvf9qdZaAOjPy+a9AABmQ+ABOiXwAJ0SeIBOCTxApwQeoFOvGPeAqnpDkl9J8qeTXE3yUGvtn1bVzyX5u0meHz30o621x271vW6//fZ26NChXS0Y9oTvntn4eNud810HC+HJJ5/8Tmvttdt93tjAJ3kxyYdba/+1ql6d5Mmqenx03ydaa/940hc7dOhQ1tbWtrtG2Ht+/cjGx3tW57kKFkRV/c+dPG9s4Ftr60nWR5+/UFXPJrljJy8GwHC2dQ2+qg4l+cEkXxnd9KGq+mpVPVxVt93kOceraq2q1p5//vmtHgLADEwc+Kp6VZLPJbmvtfaHSR5M8v1J7srGGf4/2ep5rbWHWmsrrbWV175225eQANihiQJfVa/MRtx/tbX2SJK01s611q601q4m+aUkb5vdMgHYrrGBr6pK8qkkz7bWfmHT7Qc2Pey9SZ6Z/vIA2KlJ3kXzI0nen+TpqnpqdNtHk7yvqu5K0pJ8I8lPzWB9AOzQJO+i+d0ktcVdt3zPOwDz5SdZATol8ACdEniATgk8QKcEHqBTAg/QKYEH6JTAw5SdOnUqhw8fzqlTp+a9FPY4gYcpunDhQo4ePZrTp0/n7rvvzoULF+a9JPYwgYcpuvfee3P+/Pm01nLu3LkcO3Zs3ktiD5tkFg3MxZEjR+a9hG1ZX1/P2bNnc/Xq1STJpUuXcuLEiZw8eTIHDhwY8+w/6pPveCpJct+DR6a8yvFWV1cHf01mQ+BhSi5evJj9+/dn3759L912+fLlXLx4cY6rYi8TeBbWnj6THO3JuvrA6lyXwXJzDR6gUwIP0CmBB+iUwAN0SuABOiXwAJ0SeIBOCTxApwQeoFMCD9ApgYcFYY480ybwsABunCN/5cqVeS+JDgg8LIAb58ifOXNm3kuiA6ZJ0p0e5siff/5SXnjiidw/h3nwO7GnJ38uMIGHOdtqjvy+feu5etVlGnZH4OlOF2eT1+bBf2J1rstgubkGD9ApgQfolMADdErgATol8ACdEniATgk8QKcEHqBTAg/QKYEH6JTAA3RK4GEPsrnI3iDwsMfcuLnIhQsX5r0kZkTgYY+5cXORY8eOzXtJzIhxwTAF095k5JPveCpJct+UN/zYanOREydO5OTJkzlw4MBUX2vauhgDPTCBhz1kq81FLl++nIsXL85xVcyKwMMUTP3s8tqGHw9M+fuyp7gGD9ApgQfolMADdErgATol8ACdEniATgk8QKcEHqBTAg/QKYEH6JTAw4DMYWdIAg8DMYedoQk8DMQcdoZmmiRLadrz12dtu3PYZzUPfrfMZF8uAg8DMIedeRB4llL3Z5LmwTMFrsEDdErgATol8ACdEniATgk8QKcEHqBTAg/QKYEH6JTAA3RK4AE6JfAAnRJ4GJANPxjS2MBX1Ruq6j9X1bNVdaqqfmZ0+/dV1eNV9bXRx9tmv1xYXjb8YGiTnMG/mOTDrbUfSPL2JP+gqt6S5CNJvtBae3OSL4y+Bm7Chh8Mbey44NbaepL10ecvVNWzSe5I8p4kR0YP+9dJVpPcP5NVwg1s+DEf3Y9p7sy25sFX1aEkP5jkK0n2j+Kf1tp6Vb3uJs85nuR4khw8eHBXi4VlZcMP5qFaa5M9sOpVSX4nyT9srT1SVf+3tfaaTfd/t7V2y+vwKysrbW1tbTfrhb1htOFH7lmd5ypYEFX1ZGttZbvPm+hdNFX1yiSfS/KrrbVHRjefq6oDo/sPJDm/3RcHYHYmeRdNJflUkmdba7+w6a5Hk3xw9PkHk3x++ssDYKcmuQb/I0nen+TpqnpqdNtHk3w8yb+pqmNJvpnkb81khQDsyCTvovndJHWTu9853eUAMC1+khWgUwIP0CmBB+iUwAN0SuABOiXwAJ0SeLgF89tZZgIPN2F+O8tO4OEmzG9n2W1rXDDsVO/z26dtWvPgzW/f2wQetmB+Oz0QeAbhTHKbRvPgVx9YnesyWG6uwQN0SuABOiXwAJ0SeIBOCTxApwQeoFMCD9ApgQfolMADdErgATol8ACdEnhYYjYk4VYEHpaUDUkYR+BhSdmQhHGMC4ZNFmVjknEbfsx7Q5JpMEJ69gQelpANSZiEwMMmC3NWacMPpsA1eIBOCTxApwQeoFMCD9ApgQfolMADdErgATol8ACdEniATgk8QKcEHpidU6eSw4c3PjI4gQdm48KF5OjR5PTp5O67N75mUAIPzMa99ybnzyetJefOJebVD840SVh2CzLD/o9YX0/Onk1G8+pz6VJy4kRy8mSyiPPqF2WK6JQJPDB9Fy8m+/cnm+bV5/LljdsZjMDDsuv07JPdcw0eoFMCD9ApgQfolMADdErgATol8ACdEniATgk8QKcEHqBTAg/QKYGHBXHq1KkcPnw4p8xOZ0oEHhbAhQsXcvTo0Zw+fTp33313rly5Mu8l0QGBhwVw77335vz582mt5dy5czlz5sy8l0QHTJOkO0cWcT76Layvr+fs2bO5OpqdfunSpZx//lJeeOKJ3P/gkfkubkKrJlouJIGHObt48WL279+ffZtmp+/bt56rV12mYXcEnu50cTb560eSJKufWJ3rMlhursEDdErgATol8ACdEniATgk8QKcEHqBTAg/QKYEH6JTAA3RK4AE6JfAAnRJ4WAI2A2Enxga+qh6uqvNV9cym236uqv5XVT01+nV0tsuEvevGzUAuXLgw7yWxJCY5g/90kndtcfsnWmt3jX49Nt1lAdfcuBnIsWPH5r0klsTYccGttS9W1aEB1gKDWIYNQT75jqeSJPd87M7v2QzkxIkTOXnyZA4cODDHFU6mi9HNS2w38+A/VFUfSLKW5MOtte9u9aCqOp7keJIcPHhwFy8He89Wm4Fcvnw5Fy9enOOqWBbVWhv/oI0z+N9srR0efb0/yXeStCQ/n+RAa+3ecd9nZWWlra2t7WrBsCeMNvzIPavzXAULoqqebK2tbPd5O3oXTWvtXGvtSmvtapJfSvK2nXwfAGZnR4Gvqs0X/96b5JmbPRaA+Rh7Db6qPpPkSJLbq+q5JD+b5EhV3ZWNSzTfSPJTs1siADsxybto3rfFzZ+awVoAmCI/yQrQKYEH6JTAA3RK4AE6JfAAnRJ4gE4JPCwB8+DZCYGHBWcePDsl8LDgzINnp3YzLhiWknnwwzEPfr4EHhaYefDshsCz5yzFWeVoHvy3Hlid6zJYbq7BA3RK4AE6JfAAnRJ4gE4JPECnBB6gUwIP0CmBB+iUwAN0SuABOiXwAJ0SeOAlNhbpi8ADSWws0iOBB5LYWKRHxgXDALa7yci1DT/ue3B7z9up9fX1pd1YZCnGP8+JwAM2FumUwMMAtn2WOdrwY9WGH+yCa/AAnRJ4gE4JPECnBB6gUwIP0CmBB+iUwAN0SuABOiXwAJ0SeIBOCTxMmZnqLAqBhykyU51FIvAwRWaqs0hMk2RhbXeG+rxNc6b60PPgNzNfvR8CD1NipjqLRuBZWHv6TNI8eKbANXiATgk8QKcEHqBTAg/QKYEH6JTAA3RK4AE6JfAAnRJ4gE4JPECnBB6gUwIPTMRGJstH4IGxbGSynAQeGMtGJsvJuGCYs602Npnnhh83muZGJkPa0+OmRwQeuCUbmSwvgYc52/JM04YfTIFr8ACdEniATgk8QKcEHqBTAg/QKYEH6JTAA3RK4AE6JfAAnRJ4gE4JPEyZueksirGBr6qHq+p8VT2z6bbvq6rHq+pro4+3zXaZsBzMTWeRTHIG/+kk77rhto8k+UJr7c1JvjD6GvY8c9NZJNVaG/+gqkNJfrO1dnj09ZkkR1pr61V1IMlqa+3Ocd9nZWWlra2t7XLJ7BVbzUlfZDfOTU+Sl73sZXnTm9607bnpL82D/9JdU1zhZMxRXzxV9WRrbWW7z9vpuOD9rbX1JBlF/nW3WNjxJMeT5ODBgzt8OVh85qazaGY+D7619lCSh5KNM/hZvx792NNnkubBMwU7fRfNudGlmYw+np/ekgCYhp0G/tEkHxx9/sEkn5/OcgCYlkneJvmZJF9KcmdVPVdVx5J8PMmPVdXXkvzY6GsAFsjYa/Cttffd5K53TnktAEyRn2QF6JTAA3RK4AE6JfAAnRJ4gE4JPECnBB6gUwIPTMRGJstH4IGxbGSynAQeGMtGJstp5uOCgVvbamOTlzb8ePB77xvajRuZXLp0KSdOnMjJkye3vZHJkPb0uOkRgQduyUYmy0vgYc62PNO04QdT4Bo8QKcEHqBTAg/QKYEH6JTAA3RK4AE6JfAAnRJ4gE4JPECnBB6gUwIPS8yMdm5F4GFJmdHOOAIPS8qMdsYxTRI22Wo2+zyMmwe/rDPaNzOvffYEHpaQGe1MQuBhk4U5qzQPnilwDR6gUwIP0CmBB+iUwAN0SuABOiXwAJ0SeIBOCTxApwQeoFMCD9ApgQfolMDDgGzQwZAEHgZigw6GJvAwEBt0MDTjgllKi7Ixx6S2u0HHuA0/5mVhxikzEYGHAdigg3kQeJZS92eSNvxgClyDB+iUwAN0SuABOiXwAJ0SeIBOCTxApwQeoFMCD9ApgQfolMADdErggZeYV98XgQeSmFffI4EHkphX3yPTJGEA251fP/Q8+O3Oq18k3U8W3QWBB8yr75TAwwC2fZZpHjxT4Bo8QKcEHqBTAg/QKYEH6JTAA3RK4AE6JfAAnRJ4gE4JPECnBB6gUwIP0CmBhyVmgw5uZVeBr6pvVNXTVfVUVa1Na1HAeDboYJxpnMH/aGvtrtbayhS+FzAhG3QwjnHBsMl2N+aYlXEbfizzBh3X2Khj9nYb+JbkP1ZVS/KvWmsP3fiAqjqe5HiSHDx4cJcvByQ26GAy1Vrb+ZOr/kxr7dtV9bokjyf56dbaF2/2+JWVlba25lI9jDXa8CP3rM5zFSyIqnpyJ5fBd3UNvrX27dHH80l+I8nbdvP9AJieHQe+qv5EVb362udJ/kaSZ6a1MAB2ZzfX4Pcn+Y2quvZ9fq219ltTWRUAu7bjwLfWvp7krVNcCwBT5CdZATol8ACdEniATgk8QKcEHqBTAg/QKYEHJmL2/PIReGAss+eXk8ADY5k9v5zMg4c522oG/bh58ENa1tnz5s0LPDCG2fPLS+BhzrY80xzNg199YIv7YEKuwQN0SuABOiXwAJ0SeIBOCTxApwQeoFMCD9ApgQfolMADdErgATol8ACdEnhYAjbbYCcEHhaczTbYKYGHBWezDXbKuGD2nK022Fg01zb8uOdjdy7lZhvX2HRjvgQeFpjNNtgNgWfPWYqzytGGH9+y4Qe74Bo8QKcEHqBTAg/QKYEH6JTAA3RK4AE6JfAAnRJ4gE4JPECnBB6gUwIPAzLXnSEJPAzEXHeGJvAwEHPdGZppkiylZZjpvtn6+vq25rpfmwd/34NHBlzleEsxiZOXCDwMwFx35kHgWUrdn0mO5sGvmgfPLrgGD9ApgQfolMADdErgATol8ACdEniATgk8QKcEHqBTAg/QKYEH6JTAA3RK4GGJ2UCEWxF4WFI2EGEcgYclZQMRxjEuGDZZlI1Exm34sd0NRBZR9yOfF4DAwxKygQiTEHjYZGHOKm34wRS4Bg/QKYEH6JTAA3RK4AE6JfAAnRJ4gE4JPECnBB6gUwIP0CmBB+iUwAMTMXt++ewq8FX1rqo6U1Vnq+oj01oUsFjMnl9OOw58Vb08yb9I8u4kb0nyvqp6y7QWBiwOs+eX026mSb4tydnW2teTpKo+m+Q9SU5PY2GwV2w1g37cPPghLevs+YWZDDpHuwn8HUm+tenr55L88I0PqqrjSY4nycGDB3fxcsA8mD2/vHYT+NritvY9N7T2UJKHkmRlZeV77oe9bsszTfPgmYLd/CPrc0nesOnr1yf59u6WA8C07Cbwv5/kzVX1xqral+Qnkjw6nWUBsFs7vkTTWnuxqj6U5LeTvDzJw601b5AFWBC72pO1tfZYksemtBYApshPsgJ0SuABOiXwAJ0SeIBOCTxApwQeoFMCD9ApgQfolMADdErgATol8ACdEniATgk8QKcEHqBT1dpwu+hV1QtJzgz2govt9iTfmfciFoRjcZ1jcZ1jcd2drbVXb/dJu5oHvwNnWmsrA7/mQqqqNcdig2NxnWNxnWNxXVWt7eR5LtEAdErgATo1dOAfGvj1FpljcZ1jcZ1jcZ1jcd2OjsWg/8gKwHBcogHolMADdGomga+qd1XVmao6W1Uf2eL+qqp/Nrr/q1X1Q7NYxyKY4Fj8ndEx+GpV/V5VvXUe65y1ccdh0+P+UlVdqaofH3J9Q5rkWFTVkap6qqpOVdXvDL3GoUzw9+NPVtW/r6o/GB2Ln5zHOodQVQ9X1fmqeuYm92+/m621qf5K8vIk/z3Jn02yL8kfJHnLDY85muQ/JKkkb0/ylWmvYxF+TXgs/nKS20afv7vHYzHJcdj0uP+U5LEkPz7vdc/xz8RrkpxOcnD09evmve45HouPJvlHo89fm+T/JNk377XP6Hj8tSQ/lOSZm9y/7W7O4gz+bUnOtta+3lq7nOSzSd5zw2Pek+RX2oYvJ3lNVR2YwVrmbeyxaK39Xmvtu6Mvv5zk9QOvcQiT/JlIkp9O8rkk54dc3MAmORZ/O8kjrbVvJklrrdfjMcmxaEleXVWV5FXZCPyLwy5zGK21L2bj93cz2+7mLAJ/R5Jvbfr6udFt231MD7b7+zyWjf9C92bscaiqO5K8N8kvDriueZjkz8SfS3JbVa1W1ZNV9YHBVjesSY7FP0/yA0m+neTpJD/TWrs6zPIWzra7OYtRBbXFbTe+F3OSx/Rg4t9nVf1oNgL/V2a6ovmY5Dh8Msn9rbUrGydr3ZrkWLwiyV9M8s4kfyzJl6rqy621/zbrxQ1skmPxN5M8leSvJ/n+JI9X1X9prf3hjNe2iLbdzVkE/rkkb9j09euz8V/f7T6mBxP9PqvqLyT55STvbq3974HWNqRJjsNKks+O4n57kqNV9WJr7d8NssLhTPr34zuttQtJLlTVF5O8NUlvgZ/kWPxkko+3jYvQZ6vqfyT580meGGaJC2Xb3ZzFJZrfT/LmqnpjVe1L8hNJHr3hMY8m+cDoX4XfnuT/tdbWZ7CWeRt7LKrqYJJHkry/wzO0a8Yeh9baG1trh1prh5L82yR/v8O4J5P9/fh8kr9aVa+oqj+e5IeTPDvwOocwybH4Zjb+TyZVtT/JnUm+PugqF8e2uzn1M/jW2otV9aEkv52NfyV/uLV2qqr+3uj+X8zGuySOJjmb5GI2/ivdnQmPxceS/Kkk/3J09vpi62yC3oTHYU+Y5Fi01p6tqt9K8tUkV5P8cmtty7fOLbMJ/1z8fJJPV9XT2bhEcX9rrcsRwlX1mSRHktxeVc8l+dkkr0x23k2jCgA65SdZATol8ACdEniATgk8QKcEHqBTAg/QKYEH6NT/B5a48BVFIJ21AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "24/25=0.96\n"
     ]
    }
   ],
   "source": [
    "j=0\n",
    "_, ax = plt.subplots(1, 1, figsize=(6, 8))\n",
    "for i in range(len(p_hat_list)):\n",
    "    if (p>p_hat_list[i]-SE_hat_p_list[i]) & (p<p_hat_list[i]+SE_hat_p_list[i]):\n",
    "        # interval contains p\n",
    "        ax.errorbar(p_hat_list[i], np.arange(len(p_hat_list))[i],lolims=True, xerr=SE_hat_p_list[i], yerr=0.0, linestyle='', c='black')\n",
    "        j +=1\n",
    "    else:\n",
    "        # interval does not contain p\n",
    "        ax.errorbar(p_hat_list[i], np.arange(len(p_hat_list))[i],lolims=True, xerr=SE_hat_p_list[i], yerr=0.0, linestyle='', c='red')\n",
    "ax.axvline(0.61, color='darkorange')\n",
    "plt.xlim(0,1)\n",
    "plt.show()\n",
    "print(f'{j}/{number_trials}={np.round(j/number_trials,2)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "super-footwear",
   "metadata": {},
   "source": [
    "We are plotting the sample proportion $\\hat p$ and the interval $ \\hat p \\pm 2 \\cdot SE_{\\hat p}$. Notice that about 92% of our samples intervals contain the true proportion $p$. This number will converge to 95% as the number of samples increase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "distant-swift",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1877/2000=0.94\n"
     ]
    }
   ],
   "source": [
    "number_trials = 2000\n",
    "\n",
    "j=0\n",
    "for i in range(number_trials):\n",
    "    p_hat = bernoulli.rvs(p=p, size=n)\n",
    "    p_hat = p_hat[p_hat==1].shape[0] / p_hat.shape[0]\n",
    "    SE_hat_p = 2*np.sqrt(p_hat*(1-p_hat)/n)\n",
    "    if (p>p_hat-SE_hat_p) & (p<p_hat+SE_hat_p):\n",
    "        # interval contains p\n",
    "        j +=1\n",
    "\n",
    "print(f'{j}/{number_trials}={np.round(j/number_trials,2)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "traditional-forge",
   "metadata": {},
   "source": [
    "What if we were interested in the 99% Confidence Interval? We need to calculate the critical value, also known as $z^*$, for that specific confidence level. The critical value is nothing else than the number of standard deviations below and above the mean that we need to get to capture the desired confidence level (99%). Notice that using a z-table or when using the `norm.ppf` from `scipy`, remember that the values that you get are for a single tail confidence interval. This is not what we want, so we need to get the value for 99.5% instead (leaving 0.5% on each tail of the distribution gives the 99% of confidence)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "demonstrated-millennium",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.5758293035489004"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CI = 0.99\n",
    "critical_value = norm.ppf(CI+(1-CI)/2) # we want the critical value for a two-tail distribution\n",
    "critical_value"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "future-definition",
   "metadata": {},
   "source": [
    "There is a 99% chance that $p$ is within $2.58 \\cdot \\sigma_{\\hat p}$ of $\\hat p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "disabled-connection",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "994/1000=0.99\n"
     ]
    }
   ],
   "source": [
    "number_trials = 1000\n",
    "\n",
    "j=0\n",
    "for i in range(number_trials):\n",
    "    p_hat = bernoulli.rvs(p=p, size=n)\n",
    "    p_hat = p_hat[p_hat==1].shape[0] / p_hat.shape[0]\n",
    "    SE_hat_p = critical_value*np.sqrt(p_hat*(1-p_hat)/n)\n",
    "    if (p>p_hat-SE_hat_p) & (p<p_hat+SE_hat_p):\n",
    "        # interval contains p\n",
    "        j +=1\n",
    "\n",
    "print(f'{j}/{number_trials}={np.round(j/number_trials,2)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "subject-cannon",
   "metadata": {},
   "source": [
    "# Conditions for Valid Confidence Intervals for a Proportion\n",
    "\n",
    "Let's start by recalling the conditions for valid intervals for a proportion:\n",
    "\n",
    "* The sample has to be random\n",
    "* The sampling distribution of the sample proportions can be approximated by a normal distribution. The rule of thumb is that you need to have at least 10 successes and 10 failures\n",
    "* The samples are required to be independent. The rule of thumb is that if you are sampling without replacement your sample size should be less than 10% the population size\n",
    "\n",
    "To help us understand the implications, we are going to exemplify what happens when one of the conditions is not met. First, let's create a function to compute the confidence intervals and plot the last 50 samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "automotive-tuning",
   "metadata": {},
   "outputs": [],
   "source": [
    "def confidence_interval(p, n, number_trials, N, ci=0.95, sample='random'):\n",
    "    p_ = bernoulli.rvs(p=p, size=N)\n",
    "    p_hat_list = []\n",
    "    SE_hat_p_list = []\n",
    "    \n",
    "    if sample!='random':\n",
    "        # Inducing bias on the sampling\n",
    "        p_.sort()\n",
    "        p_ = p_[:-int(0.2*N)]\n",
    "        np.random.shuffle(p_)\n",
    "    \n",
    "    for i in range(number_trials):\n",
    "        s_ = np.random.choice(p_, n, replace=False)\n",
    "        p_hat = s_[s_==1].shape[0] / s_.shape[0]\n",
    "        p_hat_list.append(p_hat)\n",
    "        SE_hat_p_list.append(2*np.sqrt(p_hat*(1-p_hat)/n))\n",
    "\n",
    "    j=0\n",
    "    _, ax = plt.subplots(1, 1, figsize=(6, 8))\n",
    "    for i in range(len(p_hat_list)):\n",
    "        if (p>p_hat_list[i]-SE_hat_p_list[i]) & (p<p_hat_list[i]+SE_hat_p_list[i]):\n",
    "            # interval contains p\n",
    "            if i > len(p_hat_list)-50:\n",
    "                ax.errorbar(p_hat_list[i], np.arange(len(p_hat_list))[i],lolims=True, xerr=SE_hat_p_list[i], yerr=0.0, linestyle='', c='black')\n",
    "            j +=1\n",
    "        else:\n",
    "            # interval does not contain p\n",
    "            if i > len(p_hat_list)-50:    \n",
    "                ax.errorbar(p_hat_list[i], np.arange(len(p_hat_list))[i],lolims=True, xerr=SE_hat_p_list[i], yerr=0.0, linestyle='', c='red')\n",
    "    ax.axvline(0.61, color='darkorange')\n",
    "    plt.xlim(0,1)\n",
    "    plt.show()\n",
    "    print(f'{j}/{number_trials}={np.round(j/number_trials,2)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "purple-finger",
   "metadata": {},
   "source": [
    "The first example is the case where our sample is not random, i.e. there is some bias that we are introducing by the way that we are sampling from the population. Referring to our context, this could mean that our group of students is surveying people outside of the school's Comedy Club, with whom the Chandler jokes land more effectively. Our confidence suddenly drops to less than 90%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "cooked-morris",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAHSCAYAAAADqNcHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAd1UlEQVR4nO3df6zd9X3f8ee7yYxa1jSMHxYBLCgChvFKlLg0qUbG0mqkvtUYWVvcbU1WI7xlbgWTNhGmKe00OQlVq6RRlUQssLb7A4JH1lBlJCJUJpsazzUyJL6mXu5KFzwsX7KijNwbc2P7vT/O118fLtf2uef7vd/z/X7P8yFZ597v+Z5zP/7Kvq/vOd/P63wiM5EkCeCHJj0ASVJ7GAqSpJKhIEkqGQqSpJKhIEkqGQqSpNKbJz2Ac7nooovyyiuvnPQwpPZ75dDg9oLrJjsOtcIzzzzzncy8eLWPa30oXHnllezbt2/Sw5Da7/O3DG7v2D3JUaglIuJ/j/M43z6SJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlS6ZyhEBEPRcR8RBwY2vY3IuLJiPhWcXvB0H33RcRcRByKiFuHtr8zIr5Z3PepiIj6/zqSpCpGeaXw+8D7lm37MPBUZl4DPFV8T0RsBLYCNxSP+XREvKl4zGeA7cA1xZ/lzylJmrBzhkJmfg34q2WbbwP+oPj6D4B/MLT9kcx8LTNfAOaAmyLiUuAtmfn1HKzq84dDj5EktcS41xTWZ+YRgOL2kmL7ZcCLQ/sdLrZdVny9fPuKImJ7ROyLiH0vv/zymEOUJK1W3ReaV7pOkGfZvqLMfCAzN2fm5osvXvVHd0iSxjRuKBwt3hKiuJ0vth8Grhja73LgpWL75StslyS1yLih8DjwweLrDwJfHNq+NSLOi4irGFxQ3lu8xfRqRLyrmHX0gaHHSJJa4pyfkhoRDwO3ABdFxGHgN4CPA49GxJ3At4FfBMjM2Yh4FDgIHAd2ZOaJ4qk+xGAm0w8DTxR/JEktcs5QyMxfPsNdP3OG/XcCO1fYvg/YtKrRSWrM7Owsd9xxB5///Oe54YYbJj0cTYiNZkksLCywZcsWDh48yMzMDAsLC5MekibEUJDEtm3bmJ+fJzM5evQod95556SHpAlp/cprUhvdcsstkx7CG3zy3c8CcM9nblnV444cOcLc3BwnT54E4NixY+zatYv9+/dz6aWX1jzK5u3evXvSQ+gUQ0GacouLi6xfv55169aV25aWllhcXJzgqDQphoI0hlaefRZrNO/+2O6JDkPd5jUFSVLJUJAklQwFqYVmZ2fZtGkTs7Ozkx6KpoyhILWMnQFNkqEgtYydAU2Ss4/UWW3sClRVpTMwbk9hUlo5g0uGgtQmdgY0aYaCOsszzWXsKagGXlOQJJUMBUlSyVCQWsiegibFUJBaxp6CJslQkFrGnoImydlH6ix7Cq9nT0F1MBSkFrGnoEkzFNRZnmkuY09BNfCagiSpZChIkkqGgtRC9hQ0KYaC1DL2FDRJhoLUMvYUNEnOPlJn2VN4PXsKqoOhILWIPQVNmqGgzvJMcxl7CqqB1xQkSSVDQZJUMhQkSSVDQWtvdhY2bRrcSmo1Q0Fra2EBtmyBgwdhZmbwvaTWMhS0trZtg/l5yISjR8EiltRqTklto76Uso4cgbk5KIpYHDsGu3bB/v1wjiJWZzgtVj1jKGjtLC7C+vUwVMRiaWmwXVIrGQpt5NmnpAnxmoIkqWQoSJJKhoJ0Di54o2liKEhn4YI3mjaGgnQWLnijaePsI9XCBW8mr+oiO34UucBQkM7IBW80jQwF1cKzzBZwkR3VwGsKkqSSoSBJKhkKUoPsPKjtDAWpIXYe1AWGgtQQOw/qAmcfqXXsPIynak+hTs5G6y5DQWqAnQd1haGg1vEsc0z2FFQDrylIkkqGgiSpZChIPWc3QqthKEg9ZjdCq2UoSD1mN0Kr5ewjqdD1fsTynkLX1oMYl7PV6mUoSD1lN0LjMBSkQufPOO0pqAZeU5AklQwFSVLJUJAklQwFqUEWydR2hoLUEItk6gJDQWqIRTJ1gVNS1TpdL5GtxEV21BWGgtQAi2TqCkNBreNZ5pgsr6kGXlOQJJUMBUlSqVIoRMTdEXEgImYj4p5i240R8fWI+GZE/HFEvGVo//siYi4iDkXErRXHLvWWfQZNytihEBGbgLuAm4AbgZ+PiGuAzwEfzsy/BfwX4F8X+28EtgI3AO8DPh0Rb6o2fKl/7DNokqq8Urge2JOZi5l5HHgauB24Dvhasc+TwD8svr4NeCQzX8vMF4A5BoEiaYh9Bk1SldlHB4CdEXEh8H1gC7Cv2P73gS8CvwhcUex/GbBn6PGHi23SWOwzvF6begqjcJZZO40dCpn5fETcz+DVwPeA54DjwDbgUxHxEeBxYKl4SKz0NCs9d0RsB7YDbNiwYdwhSp1jn0GTVqmnkJkPAg8CRMRHgcOZ+efA3yu2XQvMFLsf5vSrBoDLgZfO8LwPAA8AbN68ecXgkDzTXMaegmpQdfbRJcXtBuD9wMND234I+LfAZ4vdHwe2RsR5EXEVcA2wt8rPlyTVq2qj+bHimsIPgB2Z+UoxTXVHcf8XgP8IkJmzEfEocJDB20w7MvNExZ8vSapRpVcKmXlzZm7MzBsz86li2+9m5rXFnw9nZg7tvzMzr87M6zLziaqDl9rEboH6wEazVAO7BeoLQ0Gqgd0C9YWfkqrG9LFXAM2slTCKJnoKzvjqP0NBqshugfrEUFBjPMtcY/YUVAOvKUiSSoaCJKlkKEg1saegPjAUpBrYU1BfGApSDewpqC+cfaTG2FNYW/YUVAdDQarInoL6xFBQYzzLXGP2FFQDrylIkkqGgiSpZChIkkqGgqSRWdDrP0NB0kgs6E0HQ0HSSCzoTQenpEprrKnS3lqW19pS0FsrTpc+zVCQdE4W9KaHoSCtscbOQi2vqQZeU5AklQwFSVLJUFA/zc7Cpk2DW0kjMxTUPwsLsGULHDwIMzOD7yWNxFBQ/2zbBvPzkAlHj4Lz6aWROftomvVx0ZsjR2BuDor59Bw7Brt2wf790IP59K/j3HqtAUNB/bK4COvXw9B8epaWBtslnZOhMM0805S0jNcUJEklQ0GSVDIUpA5zfQPVzVCQOmr5+gYnTpyY9JDUA4aC1FHL1zc4dOjQpIekHnD2kXqvqfUMmrTS+gbzLx/j1b17uXcN1lOYFNc5aJ6hIHXQSusbrFt3hJMnfQtJ1RgK6r2pOds8tZ7CJ3ZPdBjqNq8pSJJKhoIkqWQoSC1k/0CTYihILbO8f7DgehBqkKEgtczy/sGdrgehBjn7SJ01Lf2DXbt2sX//fi49x3oQn3z3swDc05GewtTMCusYQ0FqkZX6B0tLSyy6HoQaYiioszzTXOZUT+Fjuyc6DHWb1xQkSSVDQZJUMhSkBtk/UNsZClJD7B+oCwwFqSH2D9QFzj5S69g/GE+begrODOsuQ0FqgP0DdYWhoNbxLHNM9hRUA68pSJJKhoIkqWQoSJJKhoKmmmUy6fUMBU0ty2TSGxkKmlqWyaQ3ckqqzqqPRTJopkzWtDOV15ziq9UwFDSVLJNJKzMUdFaeZXaI5TXVwGsKkqSSoSBJKhkKUk3sPKgPDAWpBnYe1BeGglQDOw/qC2cfqTF2HtZWE4vsOBut/wwFqSI7D+oTQ0GN8SxzjdlTUA28piBJKhkKkqRSpVCIiLsj4kBEzEbEPcW2t0fEnoh4NiL2RcRNQ/vfFxFzEXEoIm6tOHapt+w8aFLGDoWI2ATcBdwE3Aj8fERcA/wW8O8y8+3AR4rviYiNwFbgBuB9wKcj4k2VRi/1kJ0HTVKVVwrXA3syczEzjwNPA7cDCbyl2OfHgJeKr28DHsnM1zLzBWCOQaBIGmLnQZNUZfbRAWBnRFwIfB/YAuwD7gG+EhG/zSB0frrY/zJgz9DjDxfbpLH0sfdQpfPQRE+hTs5Ga6exQyEzn4+I+4Enge8BzwHHgQ8B/zIzH4uIXwIeBH4WiJWeZqXnjojtwHaADRs2jDtEqXPsPGjSInPF38urf6KIjzI4+/8Y8NbMzIgI4LuZ+ZaIuA8gMz9W7P8V4Dcz8+tne97Nmzfnvn37ahmj1GtFT4E7dk9yFGqJiHgmMzev9nFVZx9dUtxuAN4PPMzgGsLfKXZ5L/Ct4uvHga0RcV5EXAVcA+yt8vMlSfWq2mh+rLim8ANgR2a+EhF3Ab8bEW8GjlG8DZSZsxHxKHCQwdtMOzLzRMWfL0mqUaVXCpl5c2ZuzMwbM/OpYtt/z8x3Ftt+KjOfGdp/Z2ZenZnXZeYTVQcvdY39A7WdjWapIfYP1AWGgtQQ+wfqAj8lVa1j/2A8beop2EHoLkNBaoD9A3WFoaDW8SxzTK6noBp4TUGSVDIUJEklQ0GSVDIUpA6zDKe6GQpSRy0vw5044afGqDpDQeqo5WW4Q4cOTXpI6gGnpKr3pqUMN//yMV7du5d7W1Beq4vTk5tnKEgdtFIZbt26I5w86VtIqsZQUO9NzdnmqfLaJ3ZPdBjqNq8pSJJKhoIkqWQoSDWxM6A+MBSkGriAjvrCUJBq4AI66gtnH6kxfewLQDML6IyiiUV2pmYm1xQzFKSKXEBHfWIoqDGeZa4xF9lRDbymIEkqGQqSpJKhIGki7HW0k6EgqXH2OtrLUJDUOHsd7eXsI6mjlvc+mugp1KEtvY469WlmnaEgqVH2OtrNUJA66g1np/YUVAOvKUiSSoaCJKlkKKifZmdh06bBraSRGQrqn4UF2LIFDh6EmZnB95JGYiiof7Ztg/l5yISjR8E58NLInH00zfq4vsGRIzA3B8UceI4dg127YP9+6Ogc+DPq0dx4tYehoH5ZXIT162FoDjxLS4Ptks7JUJhmnmlKWsZrCpKkkqEgSSoZCpKkkqGgqeZCL9LrGQqaWi70Ir2RoaCp5UIv0hs5JVVntXwhl77o40IvZ1pkp08LwGjtGQqaSi70Iq3MUNBZeZbZIS6yoxp4TUGSVDIUJEklQ0GqiZ0H9YGhINXAzoP6wlCQamDnQX3h7CM1xs7D2jpTT6FOzkbrP0NBqsjOg/rEUFBjPMtcY/YUVAOvKUiSSoaCJKlkKEg1saegPjAUpBrYU1BfGApSDewpqC+cfaTG2FNYW/YUVAdDQarInoL6xFBQYzzLXGP2FFQDrylIkkqGgiSpZChoqtktkF7PUNDUslsgvZGhoKllt0B6I2cf6azsFnTHmXoKzvrSahgKmkp2C6SVGQo6K88yO8SegmrgNQVJUslQkCSVDAVJUqlSKETE3RFxICJmI+KeYtvnI+LZ4s9fRsSzQ/vfFxFzEXEoIm6tNnSpGRbcNE3GvtAcEZuAu4CbgCXgyxHxpcy8Y2if3wG+W3y9EdgK3AC8DfhqRFybmScqjF9aU6cKbi+++CIzMzPMzs5y/vnnT3pY0pqp8krhemBPZi5m5nHgaeD2U3dGRAC/BDxcbLoNeCQzX8vMF4A5BoEitZYFN02bKlNSDwA7I+JC4PvAFmDf0P03A0cz81vF95cBe4buP1xsUw/0seTWtYJb1UV2nH4sqBAKmfl8RNwPPAl8D3gOOD60yy9z+lUCQKz0NCs9d0RsB7YDbNiwYdwhSpVYcNM0qlRey8wHgQcBIuKjDM7+iYg3A+8H3jm0+2HgiqHvLwdeOsPzPgA8ALB58+YVg0Pt4llmC1heUw2qzj66pLjdwCAETr0y+FngzzPz8NDujwNbI+K8iLgKuAbYW+XnS5LqVfVjLh4rrin8ANiRma8U27fy+reOyMzZiHgUOMjgbaYdzjySpHap9EohM2/OzI2ZeWNmPjW0/Z9m5mdX2H9nZl6dmddl5hNVfrbURXYe1HY2mqWGuKiPusBQkBpi50Fd4Ednq3XsPIynak+hTs5G6y5DQWqAnQd1haGg1vEsc0z2FFQDrylIkkqGgiSpZCion2ZnYdOmwa2kkRkK6p+FBdiyBQ4ehJmZwfeSRmIoqH+2bYP5eciEo0fBPoA0MmcfTbMe9gE4cgTm5qDoA3DsGOzaBfv3QwvXQKjEWVpaA4aC+mVxEdavh6E+AEtLg+2SzslQmGaeaUpaxmsKkqSSoSBJKhkKUk1cK0F9YChINXCtBPWFoSDVwLUS1BfOPlJj+rhOAjSzVsIomlhPwU+w7T9DQarItRLUJ4aCGuNZ5hpzPQXVwGsKkqSSoSBJKhkKqsZ1C6ReMRQ0PtctkHrHUND4XLdA6h1nHzWtL3P1+7xugbOkNMUMBY3HdQukXjIUmuZZqKQW85qCJKlkKEiSSoaCJKlkKEg95+I/Wg1DQeoxF//RahkKUo+5+I9WyympUqHriwAtX2SnLYv/rDU/kr1ehoLUUy7+o3EYClKh82ecLrKjGnhNQZJUMhQkSSVDQWqQnQG1naEgNcTOgLrAUJAaYmdAXeDsI7VO1/sCK2miM7C8pzBJnZ/JNcUMBakBdgbUFYaCWsezzDHZU1ANvKYgSSoZCpKkkqGgfpqdhU2bBreSRmYoqH8WFmDLFjh4EGZmBt9LGomhoP7Ztg3m5yETjh4F+wDSyJx9NM162AfgyBGYm4OiD8CxY7BrF+zfDz1aQwAAZ2lpDRgK6pfFRVi/Hob6ACwtDbZLOidDYZp5pilpGa8pSJJKhoIkqWQoSOfgGgiaJoaCdBaugaBpYyhIZ+EaCJo2zj5SLVwDYfKqrqfgp9MKDAXpjFwDQdPIUFAtPMtsAddTUA28piBJKhkKkqSSoSBJKhkK6icX2ZHGYiiof1xkRxqboaD+cZEdaWxOSZ1mPSycuciOVI2hoH5xkR2pEkNhmnmmKWkZrylIkkqGgiSpVCkUIuLuiDgQEbMRcc/Q9l+PiEPF9t8a2n5fRMwV991a5WdLbeNiPOqDsa8pRMQm4C7gJmAJ+HJEfAm4HLgN+InMfC0iLin23whsBW4A3gZ8NSKuzcwTFf8O0sSdWoznxRdfZGZmhtnZWc4///xJD0tatSqvFK4H9mTmYmYeB54Gbgc+BHw8M18DyMz5Yv/bgEcy87XMfAGYYxAoUue5GI/6osrsowPAzoi4EPg+sAXYB1wL3BwRO4FjwL/KzD8DLgP2DD3+cLFNU6KPC/FAexbjqbrIzij8iPT+GzsUMvP5iLgfeBL4HvAccLx4zguAdwE/CTwaET8OxEpPs9JzR8R2YDvAhg0bxh2i1AgX41GfVOopZOaDwIMAEfFRBmf/1wNfyMwE9kbESeCi4r4rhh5+OfDSGZ73AeABgM2bN68YHOoezzLXmIvsqAZVZx+duoi8AXg/8DDwR8B7i+3XAuuA7wCPA1sj4ryIuAq4Bthb5edLkupVtdH8WHFN4QfAjsx8JSIeAh6KiAMMZiV9sHjVMBsRjwIHGbzNtMOZR5LULpVeKWTmzZm5MTNvzMynim1LmflPMnNTZr4jM/9kaP+dmXl1Zl6XmU9UHbzUBPsHmiY2mqWzONU/OHjwIDMzMyy4NoN6zlCQzsL+gaaNn5KqWvSxg9CW/sGoqvYUnB0mMBSkM7J/oGlkKKgWnmW2gD0F1cBrCpKkkqEgSSoZClKH2aFQ3QwFqaOWdyhOnPADAlSdoSB11PIOxaFDhyY9JPWAs4/Ue9PSoZh/+Riv7t3LvWu4nkLTnNXWPENB6qCVOhTr1h3h5EnfQlI1hoJ6b2rONk/1FD6xe6LDULd5TUGSVDIUJEklQ0GSVDIUpJpYJFMfGApSDVyMR31hKEg1cDEe9YVTUtWYPpbIoD2L8VRdZGcUUzO9d4oZClJFLsajPjEU1BjPMteYi+yoBl5TkCSVDAVJUslQkBpkl0FtZyhIDbHLoC4wFKSG2GVQFzj7SK3Txz5DE12GJnoKo3KmWXcZClID7DKoKwwFtY5nmWOyp6AaeE1BklQyFCRJJUNBOge7BZomhoJ0FnYLNG0MBeks7BZo2jj7SLWwWzB5VXsKzvoSGArSGdkt0DQyFFQLzzJbwJ6CauA1BUlSyVCQJJUMBanD7FCoboaC1FHLOxQnTpyY9JDUA4aC1FHLOxSHDh2a9JDUA84+Uu9NS4di/uVjvLp3L/e2YD2FujirrXmGgtRBK3Uo1q07wsmTvoWkagwF9d7UnG2e6il8YvdEh6Fu85qCJKlkKEiSSoaCJKlkKEgdZnlNdTMUpI6yvKa1YChIHWV5TWvBKanqPctr3TU104lbxFCQOsjymtaKoaDem5qzTctrqoHXFCRJJUNBklQyFKSes8ug1TAUpB5b3mVYWFiY9JDUcoaC1GPLuwx33nnnpIeklnP2kVToep/hk+9+FoB7ip7CSl2GXbt2sX//fi699NIJjbJ+UzO7rCGGgtRTK3UZlpaWWFxcnOCo1HaGglTo/BnnqZ7Cx3ZPdBjqNq8pSJJKhoIkqWQoSB1mB0F1MxSkjnI9Ba0FQ0HqKNdT0Fpw9pF6r+v9g5W4noLWiqEgdZDrKWitGArqvak523Q9BdXAawqSpJKhIEkqGQqSpFKlUIiIuyPiQETMRsQ9xbbfjIj/ExHPFn+2DO1/X0TMRcShiLi14tglSTUb+0JzRGwC7gJuApaAL0fEl4q7P5GZv71s/43AVuAG4G3AVyPi2sx0uoQktUSVVwrXA3syczEzjwNPA7efZf/bgEcy87XMfAGYYxAokqSWqBIKB4D3RMSFEfEjwBbgiuK+X4uIb0TEQxFxQbHtMuDFoccfLra9QURsj4h9EbHv5ZdfrjBESdJqjB0Kmfk8cD/wJPBl4DngOPAZ4Grg7cAR4HeKh8RKT3OG534gMzdn5uaLL7543CFKklap0oXmzHwwM9+Rme8B/gr4VmYezcwTmXkS+A+cfovoMKdfSQBcDrxU5edLkupVdfbRJcXtBuD9wMMRMbz46+0M3mYCeBzYGhHnRcRVwDXA3io/X5JUr6ofc/FYRFwI/ADYkZmvRMR/ioi3M3hr6C+BfwaQmbMR8ShwkMHbTDuceSRJ7VIpFDLz5hW2/cpZ9t8J7KzyMyVJa8dGsySpZChIkkqGgiSpZChIkkqGgiSpZChIkkqRueInTbRGRLwKHJr0OFriIuA7kx5ES3gsTvNYnOaxOO26zPzR1T6oC2s0H8rMzZMeRBtExD6PxYDH4jSPxWkei9MiYt84j/PtI0lSyVCQJJW6EAoPTHoALeKxOM1jcZrH4jSPxWljHYvWX2iWJDWnC68UJEkNaUUoRMT7IuJQRMxFxIdXuD8i4lPF/d+IiHdMYpxNGOFY/OPiGHwjIv40Im6cxDibcK5jMbTfT0bEiYj4hSbH16RRjkVE3BIRz0bEbEQ83fQYmzLC/5Efi4g/jojnimPxq5MYZxOKJY/nI+LAGe5f/e/OzJzoH+BNwP8CfhxYx2BZz43L9tkCPMFgSc93Af9j0uOe4LH4aeCC4uufm+ZjMbTfnwD/FfiFSY97gv8u3spgrZINxfeXTHrcEzwW/wa4v/j6YgarQq6b9NjX6Hi8B3gHcOAM96/6d2cbXincBMxl5l9k5hLwCHDbsn1uA/4wB/YAb122wltfnPNYZOafZuYrxbd7GCxr2kej/LsA+HXgMWC+ycE1bJRj8Y+AL2TmtwEys6/HY5RjkcCPRkQAf51BKBxvdpjNyMyvMfj7ncmqf3e2IRQuA14c+v5wsW21+/TBav+edzI4C+ijcx6LiLiMwZKvn21wXJMwyr+La4ELImJ3RDwTER9obHTNGuVY/B5wPYM14L8J3J2DNeOn0ap/d7ah0RwrbFs+JWqUffpg5L9nRPxdBqHwt9d0RJMzyrH4JHBvZp4YnBT21ijH4s3AO4GfAX4Y+HpE7MnM/7nWg2vYKMfiVuBZ4L3A1cCTEfHfMvP/rfHY2mjVvzvbEAqHgSuGvr+cQcKvdp8+GOnvGRE/AXwO+LnM/L8Nja1poxyLzcAjRSBcBGyJiOOZ+UeNjLA5o/4f+U5mLgALEfE14Eagb6EwyrH4VeDjOXhTfS4iXgD+JrC3mSG2yqp/d7bh7aM/A66JiKsiYh2wFXh82T6PAx8orqS/C/huZh5peqANOOexiIgNwBeAX+nhWeCwcx6LzLwqM6/MzCuB/wz8ix4GAoz2f+SLwM0R8eaI+BHgp4DnGx5nE0Y5Ft9m8IqJiFgPXAf8RaOjbI9V/+6c+CuFzDweEb8GfIXBzIKHMnM2Iv55cf9nGcws2QLMAYsMzgR6Z8Rj8RHgQuDTxRny8ezhB4CNeCymwijHIjOfj4gvA98ATgKfy8wVpyl22Yj/Lv498PsR8U0Gb5/cm5m9/OTUiHgYuAW4KCIOA78B/DUY/3enjWZJUqkNbx9JklrCUJAklQwFSVLJUJAklQwFSVLJUJAklQwFSVLJUJAklf4/jV0thHVHEhsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "884/1000=0.88\n"
     ]
    }
   ],
   "source": [
    "confidence_interval(p=0.61, n=35, number_trials=1000, N=500, sample='not_random')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "stone-denmark",
   "metadata": {},
   "source": [
    "The second case is the fact that we can not assume that our sampling distribution is normal. Notice that our sampling size is 10, which does not guarantee the more than 10 successes and failures. Once again, despite the fact that we are calculating the 95% confidence interval, the percentage of times where the true proportion $p$ is inside the calculated interval is around 90%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "collectible-county",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normal conditions:\n",
      "successes >= 10: False\n",
      "failures >= 10 : False\n"
     ]
    }
   ],
   "source": [
    "n = 10\n",
    "print('Normal conditions:')\n",
    "print('successes >= 10: ' + str(n*p_hat >= 10))\n",
    "print('failures >= 10 : ' + str(n*(1-p_hat) >= 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "arbitrary-missile",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAHSCAYAAAADqNcHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfGUlEQVR4nO3df4wc533f8fc3cikkbhy7+kHQkggphuRKZCPDZlQrqFzVCSqHDKraTUKmbeyGhNS6SmAWaCGrKJwUBS07SCDZCGxDtdwk/UM/rnJjBa5tyC4ot4jZC42jbB6Vqy9RarEmeHQjuPKdqTOP3/6xw+HluvTd7s7dzjP3fgHE3c3O7j4ckPeZ3Xk++0RmIkkSwA+NewCSpPYwFCRJNUNBklQzFCRJNUNBklQzFCRJtVeNewCrufLKK/P6668f9zCk9ntppvf1dW8c7zjUCl/96le/nZlXDXq/1ofC9ddfz9GjR8c9DKn9nriz93Xv4XGOQi0REf9rmPv59pEkqWYoSJJqhoIkqWYoSJJqhoIkqWYoSJJqq4ZCRHwqIuYi4viybX8tIp6JiG9UX1+37LYHImI2ImYi4q5l298SEV+vbvtoRETzfx1J0ijW8krhd4F3rNj2fuBLmXkj8KXqZyLiFmAfsKO6z8ci4rLqPh8H7gVurP6sfExJ0pitGgqZ+WXgL1Zsvhv4ver73wP+/rLtj2fmK5n5AjAL3BYR24DXZOZXsreqz+8vu48kqSWGvaawNTNPAVRfr662XwO8uGy/k9W2a6rvV27vKyLujYijEXH0zJkzQw5RkjSopi8097tOkD9ge1+Z+Uhm7srMXVddNfBHd0iShjRsKJyu3hKi+jpXbT8JXLdsv2uBb1Xbr+2zXZLUIsOGwtPAe6rv3wN8Ztn2fRFxeUTcQO+C8mT1FtPLEfHWatbRu5fdR5LUEqt+SmpEPAbcCVwZESeBXwc+BDwZEQeAbwK/AJCZ0xHxJHACOAfcl5lL1UO9l95Mph8GPlf9kSS1yKqhkJm/dImbfvoS+x8CDvXZfhTYOdDoJGmNpqen2bt3L0888QQ7duwY93CKZaNZUvHm5+fZvXs3J06cYM+ePczPz497SMUyFCQVb//+/czNzZGZnD59mgMHDox7SMVq/cprUlfdeeedjT7ew7cfA+Dgx5t93LY7deoUs7OznD9/HoCzZ88yMTHB1NQU27ZtG/PoymMoSCrawsICW7duZcuWLfW2xcVFFhYWxjiqchkK0pgcPny42Qes1mg+/GDDj6siDfuZo15TkCTVDAVJUs1QkDSS6elpdu7cyfT09LiHogYYCpKGZj+gewwFSUOzH9A9zj5S8Zqe71+qje4p2A9ol6ZmsxkKkoZiP6CbDAUVr/H5/qWyp6AGeE1BklQzFCRJNUNBGoJz89VVhoI0IOfmq8sMBWlAzs1Xlzn7qIWcd99ebZ6bv1nXU1iNs9MGYyhIA3BuvrrOUGghz2w0FHsKaoDXFCRJNUNBklQzFKSC2ZdQ0wwFqVAr+xJLS0vjHpI6wFCQCrWyLzEzMzPuIakDnH2kVdmbaJ9+fYm5M2d5eXKS++0ptEaJMwkNBalA/foSW7ac4vx530LSaAwFrarEs51N6UJP4aHDYx2GyuY1BUlSzVCQJNUMBUlSzVCQCmZ5TU0zFKRCWV7TejAUpEJZXtN6cEqqVmV5rX0sr5WhxOnchoJUIMtrWi+GglZV4tnOpmR5TQ3wmoIkqWYoSJJqhoJUMHsKapqhIBXKnoLWg6EgFcqegtaDs4+0KnsK7WNPoQwlztwzFKQC2VPQejEUtKoSz3Y2JXsKaoDXFCRJNUNBklQzFCSpKdPTsHNn72uhDAVJasL8POzeDSdOwJ49vZ8LZChIUhP274e5OciE06fhwIFxj2gozj6S1KzN2Gs5dQpmZ6HqjXD2LExMwNQUbNs23rENyFCQpFEtLMDWrbCsN8LiYm97YQwFSc2y19IOEUPdzWsKkqSaoSBJqhkKkkbimg7dYihIGtrKNR3mC52br4sMBUlDW7mmw4FC5+brImcfqXiu99Dz8O3HADi4Qesp9FvTYWJigqmpKbYVNje/C5r6NGNDQdJQ+q3psLi4yEKBc/N1kaGg4rneQ+XCegoPHh7rMFQ2rylIkmqGgiSpZihIkmqGgjQEC1vqKkNBGpCFLXWZoSANyMKWuswpqS1kGau92lzY2ujyWimcsjwYQ0EagIUtdZ2h0EKe2WgoltfUAK8pSJJqhoIkqTZSKETE+yLieERMR8TBatutEfGViPh6RPxhRLxm2f4PRMRsRMxExF0jjl2SanZHmjF0KETETuAe4DbgVuDnIuJG4JPA+zPzbwD/GfhX1f63APuAHcA7gI9FxGWjDV+S7I40aZRXCjcDRzJzITPPAc8C7wTeCHy52ucZ4B9U398NPJ6Zr2TmC8AsvUCRpJHYHWnOKLOPjgOHIuIK4HvAbuBotf3vAZ8BfgG4rtr/GuDIsvufrLZJm1LTfZTN2lNoc3ekREOHQmY+HxEfpvdq4LvAc8A5YD/w0Yj4APA0sFjdJfo9TL/Hjoh7gXsBtm/fPuwQJW0CdkeaNVJPITMfBR4FiIgPAicz80+Av1ttuwnYU+1+kouvGgCuBb51icd9BHgEYNeuXX2DQypd430UewpaJqLfefjqRp19dHX1dTvwLuCxZdt+CPg3wCeq3Z8G9kXE5RFxA3AjMDnK80uSmjVqo/mp6prC94H7MvOlaprqfdXtnwb+A0BmTkfEk8AJem8z3ZeZSyM+vySpQaO+fXRHn20fAT5yif0PAYdGeU5Jaq3padi7F554AnbsGPdohmKjWZKaMD8Pu3fDiROwZ0/v5wIZCpLUhP37YW4OMuH0aSi0K+GnpEpq1mZcD+TUKZidhaorwdmzMDEBU1NQWFfCUJCkUS0swNatsKwrweJib3thDAVJzXI9kHYYR09BktQthoIkqWYoSEPws/vVVYaCNCA/u19dZihIA/Kz+9Vlzj5qoaY/Z1/NafNn92/W9RRW0/in0XacoSANwM/uV9cZCi3kmY2G4noKaoDXFCRJNUNBklQzFCRJNUNBKpglOjXNUJAKtbJEt7Tk6rYanaEgFWpliW5mZmbcQ1IHOCVVq7JM1z79SnRzZ87y8uQk91tea40Sp5cbClKB+pXotmw5xfnzvoWk0RgKWlWJZzub0oXy2kOHxzoMlc1rCpKkmqEgSaoZClLB7CmoaYaCVCh7CloPhoJUKHsKWg/OPtKq7Cm0jz2FMpQ4c89QkApkT0HrxVDQqko829mU7CmoAV5TkCTVDAVJUs1QkApmT0FNMxSkQtlT0HowFKRC2VPQenD2kVZlT6F97CmUocSZe4aCVCB7ClovhoJWVeLZzqZkT0EN8JqCJKlmKEiSaoaCpE6ws9EMQ0FS8VZ2Nubn58c9pGIZCpKKt7KzceDAgXEPqVjOPpLGpOn+x8O3HwPg4CbrKfTrbExMTDA1NcW2bdvGPLryGAqSitavs7G4uMjCwsIYR1UuQ0Eak8b7Hxd6Cg82/LgqUkQMdT+vKUiSaoaCJKlmKEgFc26+mmYoSIVyPQWtB0NBKpTrKWg9OPtIq3I9hfZxPYUylPgJw4aCVCDXU9B6MRS0qhLPdjYl11NQA7ymIEmqGQqSpJqhIEmqGQrSECyNqasMBWlALuiiLjMUpAG5oIu6zCmpLWRZrL3avKDLZl1kZzVOqR6MoSANwAVd1HWGQgt5ZqOhuMiOGuA1BUlSzVCQJNUMBUkjsbPRLYaCpKHZ2egeQ0HS0OxsdI+zj1Q8ex09G91TaHNnYzNqataioSBpKHY2uslQUPHsdVTsKagBXlOQJNUMBUlSbaRQiIj3RcTxiJiOiIPVtjdFxJGIOBYRRyPitmX7PxARsxExExF3jTh2adOzI6CmDR0KEbETuAe4DbgV+LmIuBH4TeDfZuabgA9UPxMRtwD7gB3AO4CPRcRlI41e2sRWdgSWlpbGPSR1wCivFG4GjmTmQmaeA54F3gkk8Jpqnx8DvlV9fzfweGa+kpkvALP0AkXSEFZ2BGZmZsY9JHXAKLOPjgOHIuIK4HvAbuAocBD4QkT8Fr3Q+alq/2uAI8vuf7LappazB9A+/ToCc2fO8vLkJPe7nkJrlDgzbuhQyMznI+LDwDPAd4HngHPAe4F/kZlPRcQvAo8CPwNEv4fp99gRcS9wL8D27duHHaLUWf06Alu2nOL8ed9C0mgis+/v5cEfKOKD9M7+HwRem5kZEQF8JzNfExEPAGTmg9X+XwB+IzO/8oMed9euXXn06NFGxih1WtVTYO/hcY5CLRERX83MXYPeb9TZR1dXX7cD7wIeo3cN4W9Xu7wd+Eb1/dPAvoi4PCJuAG4EJkd5fklSs0ZtND9VXVP4PnBfZr4UEfcAH4mIVwFnqd4GyszpiHgSOEHvbab7MtPXupLUIiOFQmbe0Wfbfwfecon9DwGHRnlOSWqt6WnYuxeeeAJ27Bj3aIZio1mSmjA/D7t3w4kTsGdP7+cCGQqS1IT9+2FuDjLh9GkodG0JPyVVUrM2Y6/l1CmYnYWqN8LZszAxAVNTUNjaEoaCJI1qYQG2boVlvREWF3vbC2MoSGpWgS3eTop+feHVeU1BklQzFCRJNUNBklQzFFQcF5aR1o+hoKKsXFhmvtCCkNRWhoKKsnJhmQOFFoSktnJK6iW4sEz79FtYZmJigqmpKbYVVhBaDw/ffgyAg2NeZKfEhWV0kaGgYvRbWGZxcZGFAgtCUlsZCpfg2Y6KUy2yc/jBw2MdhsrmNQVJUs1QkCTVDAWpYHY21DRDQSrUys7G0pKr22p0hoJUqJWdjZmZmXEPSR3g7COtys5G+/TrbMydOcvLk5PcP+aegi4qcRajoSAVqF9nY8uWU5w/71tIGo2hoFWVeLazKV3oKTx0eKzDUNm8piBJqhkKkqSaoSANwX6AuspQkAbkmg7qMkNBGpBrOqjLnH3UQvYC2qvNazq0ZT2FtnH23GAMBWkArumgrjMUWsgzGw3F9RTUAK8pSJJqhoIkqWYoSOoEuyPNMBQkFc/uSHMMBUnFszvSHGcfSWPSdB9ls/YU2twdKZGhIKlodkeaZShIY9J4H8WegpaJiKHu5zUFSVLNUJAk1QwFSVLNUJA0Ektj3WIoSBqapbHuMRQkDc3SWPc4JVXFc1Gino0ur1kaa5empjgbCpKGYmmsmwwFFc9FiSqW19QArylIkmqGgiSpZihIGok9hW4xFCQNzZ5C9xgKkoZmT6F7nH2k4tlT6LGnsLnZU5A0VvYUuslQUPHsKVTsKagBXlOQJNUMBUlSzVCQhuDcfHWVoSANyLn56jJDQRqQc/PVZc4+aiHn3bdXm+fmb3RPoRTOThuMoSANwLn56jpDoYU8s9FQ7CmoAV5TkCTVDAVJUs1QkApmX0JNMxSkQq3sSywtLY17SOoAQ0Eq1Mq+xMzMzLiHpA5w9pFWZW+iffr1JebOnOXlyUnut6fQGiXOJDQUpAL160ts2XKK8+d9C0mjMRS0qhLPdjalCz2Fhw6PdRgqm9cUJEk1Q0GSVDMUJEm1kUIhIt4XEccjYjoiDlbbnoiIY9WfP4+IY8v2fyAiZiNiJiLuGm3oknSRRb5mDH2hOSJ2AvcAtwGLwOcj4rOZuXfZPr8NfKf6/hZgH7ADeD3wxYi4KTOdLiFpJBeKfC+++CJ79uxhenqaV7/61eMeVpFGeaVwM3AkMxcy8xzwLPDOCzdGRAC/CDxWbbobeDwzX8nMF4BZeoEiSSNx4aPmjDIl9ThwKCKuAL4H7AaOLrv9DuB0Zn6j+vka4Miy209W26RNqelS4GZdZKfNCx+VaOhQyMznI+LDwDPAd4HngHPLdvklLr5KAIh+D9PvsSPiXuBegO3btw87REmbgAsfNWuk8lpmPgo8ChARH6R39k9EvAp4F/CWZbufBK5b9vO1wLcu8biPAI8A7Nq1q29wSKVrvBToIjtapvcO/uBGnX10dfV1O70QuPDK4GeAP8nMk8t2fxrYFxGXR8QNwI3A5CjPL0lq1qgfc/FUdU3h+8B9mflStX0ff/mtIzJzOiKeBE7Qe5vpPmceSVK7jPRKITPvyMxbMvPWzPzSsu3/JDM/0Wf/Q5n5hsx8Y2Z+bpTn1ublfHRp/dhoVlFWLiwzPz8/7iFJnWIoqCjOR5fWlx+dfQkuLNM+zkf/wdrSU/Cj1stmKKgYzkeX1p+hcAme7ag49hTUAK8pSJJqhoIkqWYoSAWzs6GmGQpSoVZ2NpaW/IAAjc5QkAq1srMxMzMz7iGpA5x9pFXZ2Wiffp2NuTNneXlykvs32XoKbVbiLEZDQSpQv87Gli2nOH/et5A0GkNBqyrxbGdTutBTeOjwWIehsnlNQZJUMxQkSTVDQSqYPQU1zVCQCmVPQevBUJAKZU9B68HZR1qVPYX2sadQhhJn7hkKUoHsKWi9GApaVYlnO5uSPQU1wGsKkqSaoSBJqhkKktSU6WnYubP3tVCGgiQ1YX4edu+GEydgz57ezwUyFCSpCfv3w9wcZMLp03DgwLhHNBRnH0lq1mbstZw6BbOzUPVGOHsWJiZgagq2bRvv2AZkKEjSqBYWYOtWWNYbYXGxt70whoKkZtlraYeIoe7mNQVJUs1QkCTVDAVJUs1QkArmIjtqmqEgFcpFdrQeDAWpUC6yo/XglFStykV22sdFdspQ4sfOGwpSgVxkR+vFUNCqSjzb2ZRcZEcN8JqCJKlmKEiSaoaCNAT7AeoqQ0Ea0Mp+wHyhi6lI/RgK0oBW9gMOFLqYitSPs49ayF5Ae/XrB0xMTDA1NcW2MS+m8vDtxwA4aE/hL3H23GAMBWkA/foBi4uLLBS4mIrUj6HQQp7ZaCgXegoPHh7rMFQ2rylIkmqGgiSpZihIQ7CnoK4yFKQB2VNQlxkK0oDsKajLnH3UQvYU2sueQnmczTcYQ0EagD0FdZ2h0EKe2Wgo9hTUAK8pSJJqhoIkqWYoSEOwp6CuMhSkAdlTUJcZCtKA7Cmoy5x91EL2FNrLnkJ5nM03GENBGoA9BXWdodBCntloKPYU1ACvKUiSaoaCJKlmKEiSaoaCVDBLdGqaoSAVamWJbmlpadxDUgcYClKhVpboZmZmxj0kdYBTUrUqy3Tt069EN3fmLC9PTnK/5bXWKHF6uaEgFahfiW7LllOcP+9bSBqNoaBVlXi2syldKK89dHisw1DZvKYgSaoZCpKk2kihEBHvi4jjETEdEQeXbf+1iJiptv/msu0PRMRsddtdozy3JHsKat7Q1xQiYidwD3AbsAh8PiI+C1wL3A38RGa+EhFXV/vfAuwDdgCvB74YETdlplfGpCFc6Cm8+OKL7Nmzhz/94HVcdtll4x6WCjfKK4WbgSOZuZCZ54BngXcC7wU+lJmvAGTmXLX/3cDjmflKZr4AzNILFElDsKeg9TDK7KPjwKGIuAL4HrAbOArcBNwREYeAs8C/zMw/Bq4Bjiy7/8lqm1rOnkL72FMoQ4kz94YOhcx8PiI+DDwDfBd4DjhXPebrgLcCPwk8GRE/DkS/h+n32BFxL3AvwPbt24cdotRZ9hS0XkbqKWTmo8CjABHxQXpn/zcDn87MBCYj4jxwZXXbdcvufi3wrUs87iPAIwC7du3qGxzaOCWe7WxK9hTUgFFnH124iLwdeBfwGPAHwNur7TcBW4BvA08D+yLi8oi4AbgRmBzl+SVJzRq10fxUdU3h+8B9mflSRHwK+FREHKc3K+k91auG6Yh4EjhB722m+5x5JEntMurbR3f02bYI/ONL7H8IODTKc0q6aHp6mr179/LEE0+wY9yDUSfYaJYK5XoKWg+GglQoewpaD35KqlZlT6F97CmUocSZe4aCVCB7ClovhoJWVeLZzqZkT0EN8JqCJKlmKEiSaoaCVDDXU1DTDAWpUPYUtB4MBalQ9hS0Hpx9pFXZU2gfewplKHHmnqEgFciegtaLoaBVlXi2synZU1ADvKYgSaoZCpKkmqEgSaoZCiqOhS1p/RgKKsrKwtb8/Py4hyR1iqGgoqwsbB04cGDcQ5I6xSmpl2Bhq336FbYmJiaYmppi27ZtYx7d+D18+zEADo65vOYU5rIZCipGv8LW4uIiCwsLYxyV1C2GwiV4tqPiXCivPXh4rMNQ2bymIEmqGQqSpJqhoOLYU5DWj6GgothTkNaXoaCi2FOQ1pezjy7BnkL72FP4wewpqAmGgophT0Faf4bCJXi2o+LYU1ADvKYgSaoZCpKkmqEgDcGuhLrKUJAGZFdCXWYoSAOyK6Euc/ZRC9mRaK82dyXa0lNoG2cSDsZQkAZgV0JdZyi0kGc2Goo9BTXAawqSpJqhIEmqGQrSEOwpqKsMBWlA9hTUZYaCNCB7CuoyZx+1kD2F9rKnUB5n8w3GUJAGYE9BXWcotJBnNhqKPQU1wGsKkqSaoSBJqhkKkqSaoSAVzBKdmmYoSIVaWaJbWloa95DUAYaCVKiVJbqZmZlxD0kd4JRUrcoyXfv0K9HNnTnLy5OT3G95rTVKnF5uKEgF6lei27LlFOfP+xaSRmMoaFUlnu1sShfKaw8dHuswVDavKUiSaoaCJKlmKEhSU6anYefO3tdCGQqS1IT5edi9G06cgD17ej8XyFCQpCbs3w9zc5AJp09DoYsvOftIUrM2Y6/l1CmYnYWqN8LZszAxAVNTMObFlwZlKEjSqBYWYOtWWNYbYXGxt70whoKkZtlraYeIoe7mNQVJUs1QkCTVDAVJI3FNh24xFCQNbeWaDvOFzs3XRYaCpKGtXNPhQKFz83WRs49UPNd76Hn49mMAHNyg9RT6rekwMTHB1NQU2wqbm98FTX2asaEgaSj91nRYXFxkocC5+brIUFDxXO+hcmE9hQcPj3UYKpvXFCRJNUNBklQzFCRJtZFCISLeFxHHI2I6Ig5W234jIv53RByr/uxetv8DETEbETMRcdeIY5ckNWzoC80RsRO4B7gNWAQ+HxGfrW5+KDN/a8X+twD7gB3A64EvRsRNmbk07BgkSc0a5ZXCzcCRzFzIzHPAs8A7f8D+dwOPZ+YrmfkCMEsvUCRJLTFKKBwH3hYRV0TEjwC7geuq2341Ir4WEZ+KiNdV264BXlx2/5PVtv9PRNwbEUcj4uiZM2dGGKIkaRBDh0JmPg98GHgG+DzwHHAO+DjwBuBNwCngt6u79Ptw77zEYz+Smbsyc9dVV1017BAlSQMa6UJzZj6amW/OzLcBfwF8IzNPZ+ZSZp4H/j0X3yI6ycVXEgDXAt8a5fklSc0adfbR1dXX7cC7gMciYvmHnryT3ttMAE8D+yLi8oi4AbgRmBzl+SVJzRr1Yy6eiogrgO8D92XmSxHxHyPiTfTeGvpz4J8CZOZ0RDwJnKD3NtN9zjySpHYZKRQy844+2375B+x/CDg0ynNKktaPjWZJUs1QkCTVDAVJUs1QkCTVDAVJUs1QkCTVIrPvJ020RkS8DMyMexwtcSXw7XEPoiU8Fhd5LC7yWFz0xsz80UHvVMIazTOZuWvcg2iDiDjqsejxWFzksbjIY3FRRBwd5n6+fSRJqhkKkqRaCaHwyLgH0CIei4s8Fhd5LC7yWFw01LFo/YVmSdLGKeGVgiRpg7QiFCLiHRExExGzEfH+PrdHRHy0uv1rEfHmcYxzI6zhWPyj6hh8LSL+KCJuHcc4N8Jqx2LZfj8ZEUsR8fMbOb6NtJZjERF3RsSxiJiOiGc3eowbZQ3/R34sIv4wIp6rjsWvjGOcG6Fa8nguIo5f4vbBf3dm5lj/AJcBfwr8OLCF3rKet6zYZzfwOXpLer4V+B/jHvcYj8VPAa+rvv/ZzXwslu33X4H/Avz8uMc9xn8Xr6W3Vsn26uerxz3uMR6Lfw18uPr+KnqrQm4Z99jX6Xi8DXgzcPwStw/8u7MNrxRuA2Yz888ycxF4HLh7xT53A7+fPUeA165Y4a0rVj0WmflHmflS9eMResuadtFa/l0A/BrwFDC3kYPbYGs5Fv8Q+HRmfhMgM7t6PNZyLBL40YgI4K/SC4VzGzvMjZGZX6b397uUgX93tiEUrgFeXPbzyWrboPt0waB/zwP0zgK6aNVjERHX0Fvy9RMbOK5xWMu/i5uA10XE4Yj4akS8e8NGt7HWcix+B7iZ3hrwXwfel7014zejgX93tqHRHH22rZwStZZ9umDNf8+I+Dv0QuFvreuIxmctx+Jh4P7MXOqdFHbWWo7Fq4C3AD8N/DDwlYg4kpn/c70Ht8HWcizuAo4BbwfeADwTEf8tM//vOo+tjQb+3dmGUDgJXLfs52vpJfyg+3TBmv6eEfETwCeBn83M/7NBY9toazkWu4DHq0C4EtgdEecy8w82ZIQbZ63/R76dmfPAfER8GbgV6FoorOVY/Arwoey9qT4bES8Afx2Y3JghtsrAvzvb8PbRHwM3RsQNEbEF2Ac8vWKfp4F3V1fS3wp8JzNPbfRAN8CqxyIitgOfBn65g2eBy616LDLzhsy8PjOvB/4T8M87GAiwtv8jnwHuiIhXRcSPAH8TeH6Dx7kR1nIsvknvFRMRsRV4I/BnGzrK9hj4d+fYXylk5rmI+FXgC/RmFnwqM6cj4p9Vt3+C3syS3cAssEDvTKBz1ngsPgBcAXysOkM+lx38ALA1HotNYS3HIjOfj4jPA18DzgOfzMy+0xRLtsZ/F/8O+N2I+Dq9t0/uz8xOfnJqRDwG3AlcGREngV8H/goM/7vTRrMkqdaGt48kSS1hKEiSaoaCJKlmKEiSaoaCJKlmKEiSaoaCJKlmKEiSav8PLIHZpljPL0IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "897/1000=0.9\n"
     ]
    }
   ],
   "source": [
    "confidence_interval(p=0.61, n=10, number_trials=1000, N=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "victorian-harmony",
   "metadata": {},
   "source": [
    "Finally, the last condition is the independence between samples. If the rule of thumb of the 10% is not met we can not assume independence. Once again, our confidence drops to close to 90%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "advanced-glance",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Independence condition:\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "n = 150\n",
    "N=600\n",
    "print('Independence condition:')\n",
    "print(n/N < 0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "usual-personal",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAHSCAYAAAADqNcHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdtklEQVR4nO3dfZBdd33f8fcXqExCQzB+UMC2sMPYruVt7AHFhUxNXcjURJupa5rESttAI8ZqqZKxO9MOuNMh6XQEcSYZE6YDxMU0Sf/wgzANzlBgDBmbdoIqi5Fs6cqobGOKVS+SaTwUdpFXWn37xz17dLW6q/tw7t6n837NaHbvuefe+9MZaT/n3vP77C8yE0mSAF4x6gFIksaHoSBJKhkKkqSSoSBJKhkKkqSSoSBJKr1q1APo5OKLL84rr7xy1MOQxt9LR5pfL7x2tOPQWPjGN77xvcy8pNfHjX0oXHnllezbt2/Uw5DG38O3NL/e8cQoR6ExERH/u5/H+fGRJKlkKEiSSoaCJKlkKEiSSoaCJKlkKEiSSh1DISI+ExHHI+JQy7bXR8TjEfGt4uuFLffdExFzEXEkIm5t2f7WiDhY3PfxiIjB/3UkSVV0807hj4B3r9r2IeCrmXk18NXiNhGxGdgGXF885hMR8criMZ8EdgBXF39WP6ckacQ6hkJmfg34q1WbbwP+uPj+j4F/0LL9ocx8OTOfA+aAmyLiDcBrM/Pr2VzV509aHiNJGhP9XlPYmJnzAMXXS4vtlwHPt+x3tNh2WfH96u1tRcSOiNgXEftefPHFPocoSerVoC80t7tOkOfZ3lZm3p+ZWzJzyyWX9PyrOyRJfeo3FI4VHwlRfD1ebD8KXNGy3+XAC8X2y9tslySNkX5D4THgfcX37wM+37J9W0RcEBFX0bygvLf4iOkHEfG2YtbRe1seI0kaEx1/S2pEPAjcAlwcEUeB3wJ+B3gkIt4PfAf4ZYDMbETEI8Bh4BSwMzOXi6f6AM2ZTD8GfLH4I0kaIx1DITN/dY273rXG/ruAXW227wNmehqdpIFpNBrccccdPPzww1x//fWjHo7GlI1mqQYWFhbYunUrhw8fZnZ2loWFhVEPSWPKUJBqYPv27Rw/fpzM5NixY7z//e8f9ZA0psZ+5TVp2t1yyy0DeZ6Pvf0AAHd/8uznm5+fZ25ujtOnTwNw4sQJdu/ezf79+3nDG94wkNde7YknnliX59X6MxSkKbe4uMjGjRvZsGFDuW1paYnFxcURjkrjylCQRmxgZ9XFGs1PfHRAz6da8pqCJKlkKEiSSoaCVDONRoOZmRkajcaoh6IxZChINWJfQZ0YClKN2FdQJ84+ksZIlc7CWj2FFfYV1A1DQaoJ+wrqhqEgjZFKZ9b2FDQAXlOQJJUMBUlSyVCQxoDdAY0LQ0EaMbsDGieGgjRidgc0Tpx9JJ3HoNY6WMsguwOdegqt7A9oLYaCNEJ2BzRuDAXpPCbqjNqeggbAawqSpJKhIEkqGQrSFGg0Gjz11FNOZ1VlhoI04VZ6DguLixw8eNBgUCWGgjThVnoOAEsnl+w5qBJnH0nraNg9h9Onc93XSIAJm5WlnhgK0gRr7Tm8+oLvArBx4+vtOahvhoK0joZ6Rl30FF64b4ivqanjNQVJUslQkCSVDAVJUslQkLQmF/+pH0NBUlsu/lNPhoKktlz8p56ckipNqNXFuF4W2elkkIv/dMtC3HgwFCSdw8V/6stQkCbUOWfWLrKjAfCagiSpZChIkkqGgjTl7BqoF4aCNMXsGqhXhoI0xewaqFfOPpJGZNAL8KzuKdg1UD8MBWlK2TVQPwwFaUQGflZtT0ED4DUFSVLJUJAklQwFSeew21BfhoKks9htqDdDQdJZ7DbUm7OPpAmzVr9hEOsp2G2QoSCpZLdBhoI0YdY8s7anoAHwmoIkqWQoSJJKhoKkc9hTqC9DQdJZ7CnUm6Eg6Sz2FOrN2UfShLGnoPVkKEgq2VOQoSBNGHsKWk9eU5AklQwFSVLJUJAklQwFacpZRFMvDAVpillEU68MBWmKWURTr5ySKo3IWiW0fq0ur1lEUz8MBWlKWURTPwwFaUQGflZteU0D4DUFSVLJUJAklSqFQkTcFRGHIqIREXcX226IiK9HxMGI+LOIeG3L/vdExFxEHImIWyuOXdIoNBowM9P8qqnTdyhExAxwJ3ATcAPwixFxNfBp4EOZ+TeB/wL862L/zcA24Hrg3cAnIuKV1YYvaagWFmDrVjh8GGZnm7c1Vaq8U7gO2JOZi5l5CngSuB24Fvhasc/jwD8svr8NeCgzX87M54A5moEiaVJs3w7Hj0MmHDsG9h6mTpXZR4eAXRFxEfAjYCuwr9j+94HPA78MXFHsfxmwp+XxR4ttkgblwAGosMjOec3Pw9wcFL0HTpyA3bth/35Yj96DnYeR6DsUMvPZiLiX5ruBHwJPA6eA7cDHI+LDwGPAUvGQaPc07Z47InYAOwA2bdrU7xAlDdLiImzcCC29B5aWmts1NSr1FDLzAeABgIj4CHA0M78J/L1i2zXAbLH7Uc68awC4HHhhjee9H7gfYMuWLW2DQ1IbN94I9hRUQdXZR5cWXzcB7wEebNn2CuDfAp8qdn8M2BYRF0TEVcDVwN4qry9JGqyqjeZHi2sKJ4GdmflSMU11Z3H/54D/BJCZjYh4BDhM82OmnZm5XPH1JUkDVOmdQmbenJmbM/OGzPxqse0PMvOa4s+HMjNb9t+VmW/OzGsz84tVBy+pd66voPOx0SzViOsrqBNDQaoR11dQJ/6WVGmMVFljYfV6Cqu5voK6YShINeH6CuqGoSCNkUpn1q6noAHwmoIkqWQoSJJKhoJUU/YV1I6hINWQfQWtxVCQasi+gtbi7CNpDPXTV+jUU1hhX0HnYyhINWNfQedjKEhjqK8za3sKGgCvKUiSSoaCJKlkKEiSSoaCVDOW1nQ+hoJUI5bW1ImhINWIpTV14pRUaYy4yI5GzVCQasLSmrphKEhjxEV2NGpeU5AklQwFSVLJUJBqwn6CumEoSDVgP0HdMhSkGrCfoG45+0gasSrdhFZr9RTsJ6gXhoI05ewnqBeGgjRiAzurtqegAfCagiSpZChIkkqGgqSSXQYZCpIAuwxqMhQkAXYZ1OTsI2lCdOozdFpP4XzsMmiFoSDJLoNKhoI0ITqeWdtT0AB4TUGSVDIUJEklQ0GqGbsIOh9DQaoRuwjqxFCQasQugjpx9pE0RqqsrdCpp2AXQd0wFKSasIugbhgK0hipdGZtT0ED4DUFSVLJUJAklQwFaYosLCzYQVAlhoI0JZaXl3nm4EE7CKrEUJCmxDePHOHkySU7CKrE2UfSEFTpH3Rjfn6eP3zXi+XtYXQQwB7CNDIUpCmwuLjIhg0beEUEb3rTTwF2ENQfQ0EagqGcURc9hW/fN4TX0tTymoIkqWQoSJJKhoIkqWQoSFOg0Wjw1FNP2U1QZYaCNOFWFs5ZWFzk4MGDBoMqMRSkCbeycA7A0sklS2uqxCmp0joaRmmtdeGc06fT0poqMRSkCda6cM6rL/guABs3vt7SmvpmKEjraKhn1EV57QXLa6rAawqSpJKhIEkqGQqSqmk0YGam+VUTz1CQ1L+FBdi6FQ4fhtnZ5m1NNENBUv+2b4fjxyETjh0DOxITz9lH0rRbr67E/DzMzUHRkeDECdi9G/bvh/XoSNiNGApDQVJ/Fhdh40bYsOHMtqWl5nZNLENBmnaeYasHXlOQJJUMBUlSqVIoRMRdEXEoIhoRcXex7caI2BMRByJiX0Tc1LL/PRExFxFHIuLWimOXNGCNRoOZmRkadg5qq+9QiIgZ4E7gJuAG4Bcj4mrgd4F/l5k3Ah8ubhMRm4FtwPXAu4FPRMQrK41e0sCsrMtw+PBhZmdnXZehpqq8U7gO2JOZi5l5CngSuB1I4LXFPj8JvFB8fxvwUGa+nJnPAXM0A0XSGFhZlyEzOXbsmOsy1FSV2UeHgF0RcRHwI2ArsA+4G/hyRPwezdD5uWL/y4A9LY8/WmyT1IVOazN87O0HALj7k+ffr53V6zKcOHFi3ddlcE2G8dR3KGTmsxFxL/A48EPgaeAU8AHgX2bmoxHxK8ADwM8D0e5p2j13ROwAdgBs2rSp3yFK6lLrugwrlpaWXJehhiKz7c/l3p8o4iM0z/4/CrwuMzMiAvh+Zr42Iu4ByMyPFvt/GfjtzPz6+Z53y5YtuW/fvoGMUZpqxXoK3PHEKEehMRER38jMLb0+rurso0uLr5uA9wAP0ryG8HeKXd4JfKv4/jFgW0RcEBFXAVcDe6u8viRpsKo2mh8trimcBHZm5ksRcSfwBxHxKuAExcdAmdmIiEeAwzQ/ZtqZmcsVX1+SNECVQiEzb26z7b8Db11j/13AriqvKemMRqPBHXfcwcMPP8z1ox6MpoKNZmlCre4VLC/7xlvVGQrShFrdKzhy5Mioh6Qp4G9JldZBp05BVe16BcdfPMEP9u7lg330FLphr6AeDAVpArXrFWzYMM/p036EpGoMBWkdjOSsuugpPHHfCF5bU8NrCpKkkqEgSSoZCpKkkqEgTQgXwNEwGArSBHABHA2LoSBNABfA0bA4JVWqaBRFtXYL4FRZZKcdy2r1ZChIY84FcDRMhoJU0dicUa+U1z76xEiHocnmNQVJUslQkCSVDAVJJbsQMhQkAXYh1GQoSALsQqjJ2UfShOjUh6jSU+i2CzFIYzNrS2cxFCTZhVDJUJAmRMcza3sKGgCvKUiSSoaCJKlkKEg1YxdB52MoSDViF0GdGApSjdhFUCfOPpLGSJW1GTr1FOwiqBuGglQTdhHUDUNBGiOVzqztKWgAvKYgSSoZCpKkkqEgaU12GurHUJDUlp2GejIUJLVlp6GenH0kTajVnYYq6ymsZqehvgwFSeew01BfhoI0oc45s7anoAHwmoIkqWQoSJJKhoIkqWQoSDVhEU3dMBSkGrCIpm4ZClINWERTt5ySKo1YlYV1Wq1VXrOIpl4YCtKUs4imXhgK0ogN7Kza8poGwGsKkqSSoSBJKhkKkkp2GWQoSALsMqjJUJAE2GVQk7OPpAnRqc9QZZEduwxaYShIssugkqEgTYiOZ9b2FDQAXlOQJJUMBUlSyVCQpsjCwoI9A1ViKEhTYnl5mWcOHrRnoEoMBWlKfPPIEU6eXLJnoEqcfSQNwaDWTFjL/Pw8f/iuF8vbw+gZgF2DaWQoSFNgcXGRDRs28IoI3vSmnwLsGag/hoI0BEM5oy56Ct++bwivpanlNQVJUslQkCSVDAVJ53BdhfoyFCSdxXUV6s1QkHQW11WoN2cfSRNmrc5DlfUUVriuggwFSSXXVZChIE2YNc+sXU9BA+A1BUlSyVCQJJUMBUlSqVIoRMRdEXEoIhoRcXex7eGIOFD8+XZEHGjZ/56ImIuIIxFxa7WhS1ovltfqq+8LzRExA9wJ3AQsAV+KiC9k5h0t+/w+8P3i+83ANuB64I3AVyLimsxcrjB+SQO2Ul57/vnnmZ2dpdFo8JrXvGbUw9KQVHmncB2wJzMXM/MU8CRw+8qdERHArwAPFptuAx7KzJcz8zlgjmagSBojltfqrcqU1EPAroi4CPgRsBXY13L/zcCxzPxWcfsyYE/L/UeLbZJ6YHlN66nvUMjMZyPiXuBx4IfA08Cpll1+lTPvEgCi3dO0e+6I2AHsANi0aVO/Q5TUI8trqlRey8wHgAcAIuIjNM/+iYhXAe8B3tqy+1HgipbblwMvrPG89wP3A2zZsqVtcEh1ZXlN66nq7KNLi6+baIbAyjuDnwe+mZlHW3Z/DNgWERdExFXA1cDeKq8vSRqsqr/m4tHimsJJYGdmvlRs38bZHx2RmY2IeAQ4TPNjpp3OPJKk8VL146Ob19j+T9fYvgvYVeU1JZ2r0Wiw+NRTbN68GSePqgobzdKEW+kVLCwucvDgQRfFUSWGgjThVnoFAEsnl+wVqBJ/dba0jtbqFAzK6l7B6dO57r0CsFswzQwFaYK19gpefcF3Adi48fX2CtQ3Q0FaR0M9oy56Ci/cN8TX1NTxmoIkqWQoSJJKhoKkNbmuQv0YCpLaWuk/HD58mNnZWfsPNWEoSGrLdRXqydlH0oRa3YEYxHoKK1xXob4MBUnncF2F+jIUpAl1zpm16yloALymIEkqGQqSpJKhIE05uwbqhaEgTTG7BuqVoSBNMbsG6pWzj6QRGfRaC6t7CnYN1A9DQZpSdg3UD0NBGpGBn1XbU9AAeE1BklQyFCRJJUNBqhl7CzofQ0GqEXsL6sRQkGrE3oI6cfaRNEaqdBc6radgb0HdMBSkmrC3oG4YCtIYqXRmbU9BA+A1BUlSyVCQJJUMBUlSyVCQdA4LbvVlKEg6iwW3ejMUJJ3Fglu9OSVVmjBrFdw6lde6YcFNhoKkkgU3GQrShFnzzNrymgbAawqSpJKhIEkqGQrSBLNPoEEzFKQJtbpPsLy8POohaQoYCtKEWt0nOHLkyKiHpCng7CNpHVRZLKcb7foEx188wQ/27uWDFXoK52OfoB4MBWkCtesTbNgwz+nTfoSkagwFaR2M5Kx6padw3wheW1PDawqSpJKhIEkqGQpSTdhpUDcMBakGXCNB3TIUpBpwjQR1y9lH0ogNqtOw1noKrpGgXhgK0pRzjQT1wlCQRmxgZ9Wup6AB8JqCJKlkKEiSSoaCpHPYaagvQ0HSWew01JuhIOksdhrqzdlH0oRZq9ewVk+hF3YaZChIKtlpkKEgTZg1z6ztKWgAvKYgSSoZCpKkkqEgSSoZCtKEsVim9WQoSBPEYpnWm6EgTRCLZVpvTkmVBmRQi+WspVOxbBDltXYsl9WLoSBNCItlGgZDQRqQkZ9RW17TAHhNQZJUMhQkSaVKoRARd0XEoYhoRMTdLdt/MyKOFNt/t2X7PRExV9x3a5XXljRijQbMzDS/amr0fU0hImaAO4GbgCXgSxHxBeBy4DbgZzLz5Yi4tNh/M7ANuB54I/CViLgmM5cr/h0kDdvCAmzdCs8/D7OzzWB4zWtGPSoNQJV3CtcBezJzMTNPAU8CtwMfAH4nM18GyMzjxf63AQ9l5suZ+RwwRzNQJE2a7dvh+HHIhGPHwL7E1Kgy++gQsCsiLgJ+BGwF9gHXADdHxC7gBPCvMvMp4DJgT8vjjxbbJK2H9epNzM/D3BwUfQlOnIDdu2H/fliPhXhGPaurZvoOhcx8NiLuBR4Hfgg8DZwqnvNC4G3AzwKPRMRPA9Huado9d0TsAHYAbNq0qd8hSloPi4uwcSO09CVYWmpu18Sr1FPIzAeABwAi4iM0z/6vAz6XmQnsjYjTwMXFfVe0PPxy4IU1nvd+4H6ALVu2tA0OSR14hq0+VJ19tHIReRPwHuBB4E+BdxbbrwE2AN8DHgO2RcQFEXEVcDWwt8rrS5IGq2qj+dHimsJJYGdmvhQRnwE+ExGHaM5Kel/xrqEREY8Ah2l+zLTTmUeSNF4qvVPIzJszc3Nm3pCZXy22LWXmP8nMmcx8S2b+ecv+uzLzzZl5bWZ+sergJQ2WazXIRrMkwLUa1GQoSAJcq0FN/pZUaUJ0Wq+hynoKndZqWA8j/62yastQkORaDSoZCtKE6Hhm7XoKGgCvKUiSSoaCJKlkKEhak72F+jEUJLVlb6GeDAVJbdlbqCdnH0kTanVvoUpPYTV7C/VlKEg6h72F+jIUpAl1zpm1PQUNgNcUJEklQ0GSVDIUJEklQ0HSYDQaMDPT/KqJZShIqm5hAbZuhcOHYXa2eVsTyVCQVN327XD8OGTCsWNg0W1iOSVVqosOi/T0bX4e5uagKLpx4gTs3g3798N6FN0sua0rQ0FSNYuLsHEjtBTdWFpqbtfEMRSkuvAMW13wmoIkqWQoSJJKhoKkkovqyFCQBLiojpoMBUmAi+qoydlH0oRYvajOalUW2XFRHa0wFCS5qI5KhoI0ITqeWbvIjgbAawqSpJKhIEkqGQrSFLJvoH4ZCtKUsW+gKgwFacrYN1AVzj6ShqxT36BfH3v7AV5eWuKzn315qH0DsHMwTQwFaYqcPr1s30CVGArSkK3bWXXRU3jhvnV6ftWC1xQkSSVDQZJUMhSkmrHDoPMxFKQascOgTgwFqUbsMKgTZx9JY6RKh6HTegqumaBuGApSTbhmgrphKEhjpNKZtespaAC8piBJKhkKkqSSoSBJKhkKUk1YWlM3DAWpBiytqVuGglQDltbULaekSiM2qEV31iqvWVpTLwwFacpZWlMvDAVpxAZ2Vm15TQPgNQVJUslQkCSVDAVpytlPUC8MBWmK2U9QrwwFaYrZT1CvnH0kjcig+gkrVvcU7CeoH4aCNKXsJ6gfhoI0IgM/q7anoAHwmoIkqWQoSJJKhoJUE/YV1A1DQaoB+wrqlqEg1YB9BXXL2UfSiLmegsaJoSBNOfsK6oWhII2Y6ylonHhNQZJUMhQkSSVDQZJUqhQKEXFXRByKiEZE3F1s++2I+D8RcaD4s7Vl/3siYi4ijkTErRXHLkkasL4vNEfEDHAncBOwBHwpIr5Q3H1fZv7eqv03A9uA64E3Al+JiGsyc7nfMUiSBqvKO4XrgD2ZuZiZp4AngdvPs/9twEOZ+XJmPgfM0QwUSdKYqBIKh4B3RMRFEfHjwFbgiuK+34iIZyLiMxFxYbHtMuD5lscfLbadIyJ2RMS+iNj34osvVhiiJKkXfYdCZj4L3As8DnwJeBo4BXwSeDNwIzAP/H7xkGj3NGs89/2ZuSUzt1xyySX9DlGS1KNKF5oz84HMfEtmvgP4K+BbmXksM5cz8zTwHznzEdFRzryTALgceKHK60uSBqvq7KNLi6+bgPcAD0ZE6y9TuZ3mx0wAjwHbIuKCiLgKuBrYW+X1JUmDVfXXXDwaERcBJ4GdmflSRPzniLiR5kdD3wb+GUBmNiLiEeAwzY+ZdjrzSJLGS6VQyMyb22z7tfPsvwvYVeU1JUnrx0azJKlkKEiSSoaCJKlkKEiSSoaCJKlkKEiSSpHZ9jdNjI2I+AFwZNTjGBMXA98b9SDGhMfiDI/FGR6LM67NzJ/o9UGTsEbzkczcMupBjIOI2OexaPJYnOGxOMNjcUZE7OvncX58JEkqGQqSpNIkhML9ox7AGPFYnOGxOMNjcYbH4oy+jsXYX2iWJA3PJLxTkCQNyViEQkS8OyKORMRcRHyozf0RER8v7n8mIt4yinEOQxfH4h8Xx+CZiPiLiLhhFOMchk7HomW/n42I5Yj4pWGOb5i6ORYRcUtEHIiIRkQ8OewxDksX/0d+MiL+LCKeLo7Fr49inMNQLHl8PCIOrXF/7z87M3Okf4BXAv8L+GlgA81lPTev2mcr8EWaS3q+Dfgfox73CI/FzwEXFt//Qp2PRct+fw78V+CXRj3uEf67eB3NtUo2FbcvHfW4R3gs/g1wb/H9JTRXhdww6rGv0/F4B/AW4NAa9/f8s3Mc3incBMxl5l9m5hLwEHDbqn1uA/4km/YAr1u1wtu06HgsMvMvMvOl4uYemsuaTqNu/l0A/CbwKHB8mIMbsm6OxT8CPpeZ3wHIzGk9Ht0ciwR+IiIC+Os0Q+HUcIc5HJn5NZp/v7X0/LNzHELhMuD5lttHi2297jMNev17vp/mWcA06ngsIuIymku+fmqI4xqFbv5dXANcGBFPRMQ3IuK9QxvdcHVzLP4DcB3NNeAPAndlc834Our5Z+c4NJqjzbbVU6K62WcadP33jIi/SzMU/va6jmh0ujkWHwM+mJnLzZPCqdXNsXgV8FbgXcCPAV+PiD2Z+T/Xe3BD1s2xuBU4ALwTeDPweET8t8z8f+s8tnHU88/OcQiFo8AVLbcvp5nwve4zDbr6e0bEzwCfBn4hM//vkMY2bN0ciy3AQ0UgXAxsjYhTmfmnQxnh8HT7f+R7mbkALETE14AbgGkLhW6Oxa8Dv5PND9XnIuI54G8Ae4czxLHS88/Ocfj46Cng6oi4KiI2ANuAx1bt8xjw3uJK+tuA72fm/LAHOgQdj0VEbAI+B/zaFJ4Ftup4LDLzqsy8MjOvBD4L/IspDATo7v/I54GbI+JVEfHjwN8Cnh3yOIehm2PxHZrvmIiIjcC1wF8OdZTjo+efnSN/p5CZpyLiN4Av05xZ8JnMbETEPy/u/xTNmSVbgTlgkeaZwNTp8lh8GLgI+ERxhnwqp/AXgHV5LGqhm2ORmc9GxJeAZ4DTwKczs+00xUnW5b+Lfw/8UUQcpPnxyQczcyp/c2pEPAjcAlwcEUeB3wL+GvT/s9NGsySpNA4fH0mSxoShIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkq/X/s9DXkutULDQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "904/1000=0.9\n"
     ]
    }
   ],
   "source": [
    "confidence_interval(p=0.61, n=n, number_trials=1000, N=N)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "applied-knowing",
   "metadata": {},
   "source": [
    "Now, that we saw all cases where the conditions are not met, let's create one that passes all 3 tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "independent-draft",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAHSCAYAAAADqNcHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdh0lEQVR4nO3df6zd9X3f8ee7yYxa1jSMXyKABY2Agb0SJS5LqpGxthqpXY2RtbW7rckKwlvmVjBpE2Ga0k6TQ6hakUZVErGQtd0fYDyyhipLIpLJZFPj3ToyJL6md/FKFzws7KwoI3aM4+v3/jhff31yOLbPj+853x/n+ZCse/2933Pux1/BfX3P+Xxe9xOZiSRJAD9U9wAkSc1hKEiSSoaCJKlkKEiSSoaCJKlkKEiSSm+sewDnc8kll+Q111xT9zCk5ntlpffxohvqHYca4Wtf+9q3M/PScR/X+FC45ppr2LNnT93DkJpvx229j5t31TkKNURE/O9JHufbR5KkkqEgSSoZCpKkkqEgSSoZCpKkkqEgSSqdNxQi4tMRcTgi9vUd+2sR8XREfLP4eFHf1x6IiAMRsRIRt/cdf0dEfKP42sciIqr/50iSpjHKK4XfB94zcOyDwJcz8zrgy8XfiYibgC3AuuIxH4+INxSP+QSwFbiu+DP4nJKkmp03FDLzK8BfDhy+A/iD4vM/AP5+3/HHM/O1zHwBOADcEhFXAG/KzK9mb1efP+x7jCSpISadU7g8Mw8BFB8vK45fCbzYd97B4tiVxeeDx4eKiK0RsSci9hw5cmTCIUqSxlX1RPOweYI8x/GhMvORzNyQmRsuvXTsX90hSZrQpKHwcvGWEMXHw8Xxg8DVfeddBbxUHL9qyHFJUoNMGgpPAe8vPn8/8Nm+41si4oKIuJbehPJS8RbTqxHxzmLV0fv6HiNJaojz/pbUiHgMuA24JCIOAr8BfAR4IiLuBr4F/CJAZi5HxBPAfuAksC0zV4un+gC9lUw/DHy++CNJapDzhkJm/vJZvvQzZzl/O7B9yPE9wPqxRic1wPLyMps3b2bHjh2sW7eu7uFIM2WjWTqHo0ePsnHjRvbv38+mTZs4evRo3UOSZspQkM7hrrvu4vDhw2QmL7/8MnfffXfdQ5JmqvE7r6kdbrvttrqHULlDhw5x4MABTp06BcDx48fZuXMne/fu5Yorrqh5dK/30Xc9C8B9n7htosfv2rWrsrGovQwF6SyOHTvG5Zdfzpo1a8pjJ06c4NixYzWOSpotQ0GV8C6zAYo9mnc9uKvWYajdnFOQJJUMBUlSyVCQzmN5eZn169ezvLxc91CkmTMUpHOwp6BFYyhI52BPQYvG1UeqhD2F+tlTUBUMBeks7CloERkKqoR3mQ1gT0EVcE5BklQyFCRJJUNBaiC7EaqLoSA1jN0I1clQkBrGboTq5OojtZbdiB80bU9h3lyx1kyGgtQgdiNUN0NBreWd5gB7CqqAcwqSpJKhIEkqGQpSx9l50DgMBanD7DxoXIaC1GF2HjQuVx9Jhbb3HgZ7Cm3bD2JSrkKrlqEgdZSdB03CUJAKrb/jtKegCjinIEkqGQqSpJKhIEkqGQpSx1le0zgMBanDLK9pXIaC1GGW1zQul6RKBctr7dT6pcQNYyhIHWV5TZMwFKRC6+84La+pAs4pSJJKhoIkqWQoSALsM6jHUJBkn0ElQ0GSfQaVXH0kTaCJnYbBnsKout5naP2qsjkzFKQFZ59B/QwFaQKNvPu0p6AKOKcgSSoZCpKkkqEgaWR2GbrPUJA0ErsMi8FQkDQSuwyLwdVH0ozNq9MwaU9hFHYZFoehIOm87DIsDkNBmrG53YXaU1AFnFOQJJUMBUlSyVCQzsO1+VokhoJ0Dq7N16IxFKRzcG2+Fo2rj1SJJu4vMK22rc2ftqfgWn2BoSCdlWvztYgMBVXCu8wGsKegCjinIEkqGQqSpJKhIEkqGQpaaBbTpB9kKGhhWUyTXs9Q0MKymCa9nktSdU5dLKVB+4ppozhbec3lwhqHoaCFZDFNGs5Q0Dl5l9kiltdUAecUJEklQ0GSVJoqFCLi3ojYFxHLEXFfcezmiPhqRHwjIv44It7Ud/4DEXEgIlYi4vYpxy5pBHYxNI6JQyEi1gP3ALcANwM/HxHXAZ8CPpiZfwP4z8C/Ks6/CdgCrAPeA3w8It4w3fAlnYtdDI1rmlcKNwK7M/NYZp4EngHuBG4AvlKc8zTwD4rP7wAez8zXMvMF4AC9QJE0I3YxNK5pVh/tA7ZHxMXA94CNwJ7i+N8DPgv8InB1cf6VwO6+xx8sjkmN0PZOxmBPoYtdjGFcIVetiUMhM5+PiIfovRr4LvAccBK4C/hYRHwIeAo4UTwkhj3NsOeOiK3AVoC1a9dOOkRpodnF0CSm6ilk5qPAowAR8WHgYGb+GfB3i2PXA5uK0w9y5lUDwFXAS2d53keARwA2bNgwNDikqrX+jtOegiow7eqjy4qPa4H3Ao/1Hfsh4N8AnyxOfwrYEhEXRMS1wHXA0jTfX5JUrWkbzU8WcwrfB7Zl5ivFMtVtxdc/A/wHgMxcjogngP303mbalpmrU35/SVKFpnqlkJm3ZuZNmXlzZn65OPa7mXl98eeDmZl952/PzLdm5g2Z+flpBy91ld0C1cVGs9QwdgtUJ0NBahi7BaqTvyVVrdX2XsEw03QLzrafQlO1frVXRxkKUoPYLVDdDAW1lneaA+wpqALOKUiSSoaCJKlkKEgamf2J7jMUJI3E/sRiMBQkjcT+xGJw9ZE0Y/PqU8yyp9D1vRlcyXaGoSDpvOxPLA5DQZqxud2F2lNQBZxTkCSVDAVJUslQkCSVDAVJI7O81n2GgqSRWF5bDIaCpJFYXlsMLkmVZszyWvNZXjvDUJB0XpbXFoehIM2Y5TW1iXMKkqSSoSBJKhkKWmiuu5d+kKGgheW6e+n1DAUtLNfdS6/n6iOd07zW2M9bF9fdn62n4Bp8jcNQ0EJy3b00nKGgc/Ius0XsKagCzilIkkqGgiSpZCiom5aXYf363kdJIzMU1D1Hj8LGjbB/P2za1Pu7pJEYCuqeu+6Cw4chE15+GewfSCNz9dEi62IH4dAhOHAAiv4Bx4/Dzp2wdy+0tH9wVq4M0wwYCuqWY8fg8suhr3/AiRO945LOy1BYZN5pShrgnIIkqWQoSJJKhoLUYu4HoaoZClJLDe4Hsbq6WveQ1AGGgtRSg/tBrKys1D0kdYCrj9R5XdwTYth+EIePHOfVpSXuH9hPoc38Lb3zZyhILTRsP4g1aw5x6pRvIWk6hoI6b2HuNk/vp/DwrlqHoXZzTkGSVDIUJEklQ0E6D7sAWiSGgnQOg12Ao+7NoI4zFKRzGOwC3O3eDOo4Vx+pEovSBdi5cyd79+7ligbuzfDRdz0LwH0T9hQWZpWWzslQkM5iWBfgxIkTHHNvBnWYoaBKeJfZAKd7Cg/uqnUYajfnFCRJJUNBklQyFCRJJUNBOg/La1okhoJ0DpbXtGgMBekcLK9p0bgkVZWwvFY/y2uqgqEgnYXlNS0iQ0GV8C6zASyvqQLOKUiSSoaCJKlkKEgC7GOox1CQZB9DJUNBkn0MlVx9JE2gib2MSXsKbetjjMuVceMxFKQFZx9D/QwFaQKNvPu0p6AKOKcgSSoZCpKk0lShEBH3RsS+iFiOiPuKY2+LiN0R8WxE7ImIW/rOfyAiDkTESkTcPuXYpc6yM6C6TBwKEbEeuAe4BbgZ+PmIuA74LeDfZubbgA8VfycibgK2AOuA9wAfj4g3TDV6qYPsDKhO07xSuBHYnZnHMvMk8AxwJ5DAm4pzfgx4qfj8DuDxzHwtM18ADtALFEl97AyoTtOsPtoHbI+Ii4HvARuBPcB9wBcj4rfphc5PFedfCezue/zB4pg0kSZ2BaY1TWdg2v0U5q2RK7g0eShk5vMR8RDwNPBd4DngJPAB4F9k5pMR8UvAo8DPAjHsaYY9d0RsBbYCrF27dtIhSq1jZ0B1i8yhP5fHf6KID9O7+38QeHNmZkQE8J3MfFNEPACQmQ8W538R+M3M/Oq5nnfDhg25Z8+eSsYodVrRU2DzrjpHoYaIiK9l5oZxHzft6qPLio9rgfcCj9GbQ/jbxSk/DXyz+PwpYEtEXBAR1wLXAUvTfH9JUrWmbTQ/WcwpfB/YlpmvRMQ9wO9GxBuB4xRvA2XmckQ8Aeyn9zbTtsxcnfL7S5IqNFUoZOatQ479d+AdZzl/O7B9mu8p6Yzl5WU2b97Mjh07WFf3YNQJNpqllhrsM6yu+sJb0zMUpJYa7DOsrKzUPSR1gL8lVZ23KH2Gw0eO8+rSEve3pKcwCrsM82coSC00rM+wZs0hTp3yLSRNx1BQ5y3M3ebp/RQe3lXrMNRuzilIkkqGgiSpZChIkkqGgiTAjX3UYyhIcmMflQwFSW7so5JLUqUJNLEQN+kmO9Ns7NMGC7MkuSKGgrTg3NhH/QwFaQKNvPs8XV57cFetw1C7OacgSSoZCpKkkqEgVcR1/uoCQ0GqgOv81RWGglQB1/mrK1x9pLlp4tr+KjRlnf+kPYVxNHLVlSplKEhTcp2/usRQ0Nx4lzlj9hRUAecUJEklQ0GSVDIUpBazG6GqGQpSSw12I1ZXV+sekjrAUJBaarAbsbKyUveQ1AGuPlLndbEfMawbcfjIcV5dWuL+GfYU5s0Va/NnKEgtNKwbsWbNIU6d8i0kTcdQUOctzN3m6Z7Cw7tqHYbazTkFSVLJUJAklQwFSYCdB/UYCpLcD0IlQ0GS+0Go5OojaQJN7D5Mup9CU/aDmJWFWX1WEUNBWnDuB6F+hoI0gUbefbqfgirgnIIkqWQoSJJKhoIkqWQoSBqZBbfuMxQkjcSC22IwFCSNxILbYnBJqjRj8yq6TVpeG4UFt8VhKEg6Lwtui8NQkGZsbnehltdUAecUJEklQ0GSVDIUpBazN6CqGQpSSw32BlZXV+sekjrAUJBaarA3sLKyUveQ1AGuPlLnNXFDnGkN6w0cPnKcV5eWuH8GPYW62B+YP0NBaqFhvYE1aw5x6pRvIWk6hoI6b2HuNk/3FB7eVesw1G7OKUiSSoaCJKlkKEgNZP9AdTEUpIZx3wLVyVCQGsZ9C1QnVx+ptRalfzDqvgWz3E9hFhZmVVjLGApSg7hvgepmKKi1vNMc4H4KqoBzCpKkkqEgSSoZClLH2XnQOAwFqcPsPGhchoLUYXYeNC5XH0mFtvceBnsK03Qe2sRVaNUyFKSOsvOgSRgKUqH1d5z2FFQB5xQkSSVDQZJUMhQkSaWpQiEi7o2IfRGxHBH3Fcd2RMSzxZ+/iIhn+85/ICIORMRKRNw+3dAlzZtFuO6beKI5ItYD9wC3ACeAL0TE5zJzc985vwN8p/j8JmALsA54C/CliLg+M1enGL+kOTldhHvxxRfZtGkTy8vLXHjhhXUPSxWb5pXCjcDuzDyWmSeBZ4A7T38xIgL4JeCx4tAdwOOZ+VpmvgAcoBcoklrAItximGZJ6j5ge0RcDHwP2Ajs6fv6rcDLmfnN4u9XArv7vn6wOCZ12rxKcbPcZKfrRbjWL0eu0MShkJnPR8RDwNPAd4HngJN9p/wyZ14lAMSwpxn23BGxFdgKsHbt2kmHKKkiFuEWx1Tltcx8FHgUICI+TO/un4h4I/Be4B19px8Eru77+1XAS2d53keARwA2bNgwNDiktpjbXajlNVVg2tVHlxUf19ILgdOvDH4W+LPMPNh3+lPAloi4ICKuBa4Dlqb5/pKkak37ay6eLOYUvg9sy8xXiuNb+MG3jsjM5Yh4AthP722mba48kqRmmeqVQmbempk3ZebNmfnlvuP/JDM/OeT87Zn51sy8ITM/P833ltrIdf5qOhvN0py44Y3awFCQ5sR1/moDf3W2Gqftm90MM491/rPsKYzLdf/tZShIc+A6f7WFoaDG8S5zQvYUVAHnFCRJJUNBklQyFKSK2EFQFxgKUgXsIKgrDAWpAnYQ1BWuPtLcdLF/AM3Za2AePQVXhnWfoSBNyQ6CusRQ0Nx4lzlj9hRUAecUJEklQ0GSVDIUpBazG6GqGQpSSw12I1ZX3chQ0zMUpJYa7EasrKzUPSR1gKuP1Hld7EcM60YcPnKcV5eWuL8B+ylUxRVr82coSC00rBuxZs0hTp3yLSRNx1BQ5y3M3ebpnsLDu2odhtrNOQVJUslQkCSVDAWpxewpqGqGgtRS9hQ0C4aC1FL2FDQLrj5S59lTaK+FWTnWIIaC1EL2FDQrhoI6b2HuNu0pqALOKUiSSoaCJKlkKEiSSoaC1GKW11Q1Q0FqKctrmgVDQWopy2uaBZekqvMsr7XXwiwnbhBDQWohy2uaFUNBnbcwd5uW11QB5xQkSSVDQZJUMhQkjcxeRPcZCpJGMtiLOHr0aN1D0gwYCpJGMtiLuPvuu+sekmbA1UfSjM2rJ/HRdz0LwH0z6CkM60Xs3LmTvXv3csUVV1T+/eZtYVaojcBQkHRew3oRJ06c4NixYzWOSrNgKEgzNre70NM9hQfn9P3USc4pSJJKhoIkqWQoSA1kH0B1MRSkhrEPoDoZClLD2AdQnVx9pNZalH0SRu0DzLKnMAt2A5rJUJAaxD6A6mYoqLW80xxgT0EVcE5BklQyFCRJJUNBUuXsWbSXoSCpUvYs2s1QkFQpexbt5uojqUGm6V40oafQ1X0XFmmlm6EgqTL2LNrPUJAaZKo7UnsKqoBzCpKkkqEgSSoZCpKkkqEgtZglMVXNUJBaarAktrq6WveQ1AGGgtRSgyWxlZWVuoekDnBJqjpvUTbjOXzkOK8uLXF/SzbZGcUilcaawlCQWmhYSWzNmkOcOuVbSJqOoaDOW5i7zdPltYd31ToMtZtzCpKkkqEgSSpNFQoRcW9E7IuI5Yi4r+/4r0fESnH8t/qOPxARB4qv3T7N95Y0GrsMGsfEcwoRsR64B7gFOAF8ISI+B1wF3AH8RGa+FhGXFeffBGwB1gFvAb4UEddnpjNj0oyc7jK8+OKLbNq0ieXlZS688MK6h6UGm+aVwo3A7sw8lpkngWeAO4EPAB/JzNcAMvNwcf4dwOOZ+VpmvgAcoBcokmbEDW80rmlWH+0DtkfExcD3gI3AHuB64NaI2A4cB/5lZv4pcCWwu+/xB4tjUiO0vc8wuMlOVze8GbQwq8vmZOJQyMznI+Ih4Gngu8BzwMniOS8C3gn8JPBERPw4EMOeZthzR8RWYCvA2rVrJx2itNDc8EaTmKqnkJmPAo8CRMSH6d393wh8JjMTWIqIU8Alxdeu7nv4VcBLZ3neR4BHADZs2DA0OKSqtf6O0012VIFpVx+dnkReC7wXeAz4I+Cni+PXA2uAbwNPAVsi4oKIuBa4Dlia5vtLkqo1baP5yWJO4fvAtsx8JSI+DXw6IvbRW5X0/uJVw3JEPAHsp/c20zZXHklSs0z79tGtQ46dAP7xWc7fDmyf5ntKOmN5eZnNmzezY8cO1tU9GHWCjWappdxPQbNgKEgt5X4KmgV/S6o6r+39g2HcT0GzYihILeR+CpoVQ0GdtzB3m+6noAo4pyBJKhkKkqSSoSA1kHsgqC6GgtQwg/2Do0eP1j0kLRBDQWoY90BQnVx9pNZalP7BqHsgDO6n0HQLsyqsZQwFqUHcA0F1MxTUWt5pDnA/BVXAOQVJUslQkCSVDAVJUslQkDrOIpzGYShIHWYRTuMyFKQOswincbkkVSq0vQw3WF6bpgjXJi5NrpahIHWURThNwlCQCq2/47S8pgo4pyBJKhkKkqSSoSA1kN0C1cVQkBrGboHqZChIDWO3QHVy9ZFaq+29gmHcZEd1MxSkBrFboLoZCmot7zQH2FNQBZxTkCSVDAVJUslQkDrOzoPGYShIHWbnQeMyFKQOs/Ogcbn6SCq0vffgfgqqgqEgdZSdB03CUJAKrb/jtKegCjinIEkqGQqSpJKhILWYHQRVzVCQWmqwg7C6ulr3kNQBhoLUUoMdhJWVlbqHpA5w9ZE6r+39g2GGdRAOHznOq0tL3N+S/RRG0foVYS1kKEgtNKyDsGbNIU6d8i0kTcdQUOctzN3m6Z7Cw7tqHYbazTkFSVLJUJAklQwFSVLJUJDmyLKZms5QkObEDW/UBoaCNCdueKM2cEmqGmdRymZVb3gzuMlOnRZmGXAHGQrSHLjhjdrCUFDjeJc5ITfZUQWcU5AklQwFSVLJUJAE2KFQj6EgyQ6FSoaCJDsUKrn6SJpAE7sUk/YU5tGhqJOr2cZjKEgLzg6F+hkK0gQaefdpT0EVcE5BklQyFCRJJUNB6jj7BxqHoSB1mP0DjctQkDrM/oHG5eojqdDE7sE4BnsKXe8fnNbIlWAtZihIHWX/QJMwFKRC6+847SmoAs4pSJJKhoIkqWQoSJJKU4VCRNwbEfsiYjki7iuO/WZE/J+IeLb4s7Hv/Aci4kBErETE7VOOXZJUsYknmiNiPXAPcAtwAvhCRHyu+PLDmfnbA+ffBGwB1gFvAb4UEddn5uqkY5AkVWuaVwo3Arsz81hmngSeAe48x/l3AI9n5muZ+QJwgF6gSJIaYppQ2Ae8OyIujogfATYCVxdf+7WI+HpEfDoiLiqOXQm82Pf4g8Wx14mIrRGxJyL2HDlyZIohSpLGMXEoZObzwEPA08AXgOeAk8AngLcCbwMOAb9TPCSGPc1ZnvuRzNyQmRsuvfTSSYcoSRrTVBPNmfloZr49M98N/CXwzcx8OTNXM/MU8O858xbRQc68kgC4Cnhpmu8vSarWtKuPLis+rgXeCzwWEf2/VOVOem8zATwFbImICyLiWuA6YGma7y9Jqta0v+biyYi4GPg+sC0zX4mI/xgRb6P31tBfAP8UIDOXI+IJYD+9t5m2ufJIkpplqlDIzFuHHPuVc5y/Hdg+zfeUJM2OjWZJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVInPob5pojIh4FVipexwNcQnw7boH0RBeizO8Fmd4Lc64ITN/dNwHtWGP5pXM3FD3IJogIvZ4LXq8Fmd4Lc7wWpwREXsmeZxvH0mSSoaCJKnUhlB4pO4BNIjX4gyvxRleizO8FmdMdC0aP9EsSZqfNrxSkCTNSSNCISLeExErEXEgIj445OsRER8rvv71iHh7HeOchxGuxT8qrsHXI+JPIuLmOsY5D+e7Fn3n/WRErEbEL8xzfPM0yrWIiNsi4tmIWI6IZ+Y9xnkZ4f+RH4uIP46I54pr8at1jHMeii2PD0fEvrN8ffyfnZlZ6x/gDcD/An4cWENvW8+bBs7ZCHye3pae7wT+R93jrvFa/BRwUfH5zy3yteg7778C/wX4hbrHXeN/F2+mt1fJ2uLvl9U97hqvxb8GHio+v5TerpBr6h77jK7Hu4G3A/vO8vWxf3Y24ZXCLcCBzPzzzDwBPA7cMXDOHcAfZs9u4M0DO7x1xXmvRWb+SWa+Uvx1N71tTbtolP8uAH4deBI4PM/Bzdko1+IfAp/JzG8BZGZXr8co1yKBH42IAP4qvVA4Od9hzkdmfoXev+9sxv7Z2YRQuBJ4se/vB4tj457TBeP+O++mdxfQRee9FhFxJb0tXz85x3HVYZT/Lq4HLoqIXRHxtYh439xGN1+jXIvfA26ktwf8N4B7s7dn/CIa+2dnExrNMeTY4JKoUc7pgpH/nRHxd+iFwt+a6YjqM8q1+Chwf2au9m4KO2uUa/FG4B3AzwA/DHw1InZn5v+c9eDmbJRrcTvwLPDTwFuBpyPiv2Xm/5vx2Jpo7J+dTQiFg8DVfX+/il7Cj3tOF4z074yInwA+BfxcZv7fOY1t3ka5FhuAx4tAuATYGBEnM/OP5jLC+Rn1/5FvZ+ZR4GhEfAW4GehaKIxyLX4V+Ej23lQ/EBEvAH8dWJrPEBtl7J+dTXj76E+B6yLi2ohYA2wBnho45yngfcVM+juB72TmoXkPdA7Oey0iYi3wGeBXOngX2O+81yIzr83MazLzGuA/Af+8g4EAo/0/8lng1oh4Y0T8CPA3gefnPM55GOVafIveKyYi4nLgBuDP5zrK5hj7Z2ftrxQy82RE/BrwRXorCz6dmcsR8c+Kr3+S3sqSjcAB4Bi9O4HOGfFafAi4GPh4cYd8Mjv4C8BGvBYLYZRrkZnPR8QXgK8Dp4BPZebQZYptNuJ/F/8O+P2I+Aa9t0/uz8xO/ubUiHgMuA24JCIOAr8B/BWY/GenjWZJUqkJbx9JkhrCUJAklQwFSVLJUJAklQwFSVLJUJAklQwFSVLJUJAklf4/9/17A/6EVUoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "947/1000=0.95\n"
     ]
    }
   ],
   "source": [
    "confidence_interval(p=0.61, n=35, number_trials=1000, N=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "original-liverpool",
   "metadata": {},
   "source": [
    " Our confidence converges effectively to 95%. We can feel confident about our confidence."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "guilty-twelve",
   "metadata": {},
   "source": [
    "# Estimating the Population Mean\n",
    "\n",
    "We have been solving the problem of estimating the population proportion, which percentage of the population prefers Ross's humor over Chandler's. A different problem is the estimation of the mean of a population. Let's see the main differences. When we estimated the confidence interval for the population proportion we defined:\n",
    "\n",
    "$$\\hat p \\pm z^* \\cdot \\sqrt{\\frac{\\hat p (1-\\hat p)}{n}}$$\n",
    "\n",
    "Following the same principle, we would define the population mean as,\n",
    "\n",
    "$$\\bar x \\pm z^* \\cdot \\frac{s}{\\sqrt{n}}$$\n",
    "\n",
    "Notice that we do not know our population standard deviation, so we use our best estimate, the sample standard deviation. Unfortunately, if we were to use this approach to calculate our confidence interval, we would be underestimating the actual interval. To achieve the 95% confidence we need to use a different critical value based on a t-distribution.\n",
    "\n",
    "\n",
    "$$\\bar x \\pm t^* \\cdot \\frac{s}{\\sqrt{n}}$$\n",
    "\n",
    "Let's do a couple of experiments to prove our point. First, we define the confidence interval for the mean based on the sample standard deviation but calculating the critical value from a normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "finnish-approval",
   "metadata": {},
   "outputs": [],
   "source": [
    "def confidence_interval_mean(μ, σ, n, number_trials, N, ci=0.95, sample='random'):\n",
    "    x_ = norm.rvs(loc=μ, scale=σ, size=N)\n",
    "    x_hat_list = []\n",
    "    SE_hat_x_list = []\n",
    "    \n",
    "    if sample!='random':\n",
    "        # Inducing bias on the sampling\n",
    "        x_.sort()\n",
    "        x_ = x_[:-int(0.2*N)]\n",
    "        np.random.shuffle(x_)\n",
    "    \n",
    "    for i in range(number_trials):\n",
    "        s_ = np.random.choice(x_, n, replace=False)\n",
    "        x_hat = np.mean(s_)\n",
    "        x_hat_list.append(x_hat)\n",
    "        SE_hat_x_list.append(norm.ppf(ci+(1-ci)/2)*np.std(s_)/np.sqrt(n))\n",
    "\n",
    "    j=0\n",
    "    _, ax = plt.subplots(1, 1, figsize=(6, 8))\n",
    "    for i in range(len(x_hat_list)):\n",
    "        if (μ>x_hat_list[i]-SE_hat_x_list[i]) & (μ<x_hat_list[i]+SE_hat_x_list[i]):\n",
    "            # interval contains p\n",
    "            if i > len(x_hat_list)-50:\n",
    "                ax.errorbar(x_hat_list[i], np.arange(len(x_hat_list))[i],lolims=True, xerr=SE_hat_x_list[i], yerr=0.0, linestyle='', c='black')\n",
    "            j +=1\n",
    "        else:\n",
    "            # interval does not contain p\n",
    "            if i > len(x_hat_list)-50:    \n",
    "                ax.errorbar(x_hat_list[i], np.arange(len(x_hat_list))[i],lolims=True, xerr=SE_hat_x_list[i], yerr=0.0, linestyle='', c='red')\n",
    "    ax.axvline(μ, color='darkorange')\n",
    "    #plt.xlim(0,1)\n",
    "    plt.show()\n",
    "    print(f'{j}/{number_trials}={np.round(j/number_trials,2)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "jewish-heart",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAHSCAYAAAADqNcHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeh0lEQVR4nO3df6zdd33f8ec7QCgQvLIlsTw7lumUeGBDgLitWX/Ma1Wa2hJJRyOTVqSqrXlkaIJuQqQIBWlT1W7S2mBtTWSlSYjKgmsIS6KRoozVmA0DdbDBvk6gF1rgLne+oWEjs2cuid/743z93cnN/XHOueec74/zfEhX99zv+ZHvV8nN63vO9/26n8hMJEkCuKTqHZAk1YehIEkqGQqSpJKhIEkqGQqSpJKhIEkqvbTqHVjJ5Zdfnps2bap6N6TR+v7XO99fs7na/VBrPP7449/LzCv6fV7tQ2HTpk0cO3as6t2QRuvgjs733Yer3Au1SER8e5Dn+fGRJKlkKEiSSoaCJKlkKEiSSoaCJKlkKEiSSiuGQkRcFRF/HhFPRMRURLy32P63I+KxiPjL4vtrup7zOxExHRFfj4hf7tp+XUScLO7bHxExmsOSJA2il3cKzwH/MjNfB2wH3hMRrwduAz6bmVcDny1+prjvncAW4HrgjyLiJcVr3QnsA64uvq4f4rFIklZpxVDIzNnM/Epx+1ngCWA9cAPw0eJhHwVuLG7fAHw8M3+YmX8FTAM/FRHrgDWZeTQ7K/vc3/UcSVIN9HVNISI2AW8GvgSszcxZ6AQHcGXxsPXAd7ueNlNsW1/cXrh9sX/Ovog4FhHHnn766X52UZK0Cj2HQkRcBnwSeF9m/mC5hy6yLZfZ/uKNmQcyc1tmbrviir7/dIckaUA9hUJEvIxOIHwsMx8sNp8pPhKi+D5XbJ8Brup6+gbgqWL7hkW2S5JqopfpowD+GHgiM/+g666Hgd8sbv8m8FDX9ndGxMsj4rV0Lih/ufiI6dmI2F685i1dz5Ek1UAvfyX1Z4B3AScj4kSx7YPA7wN/GhF7ge8ANwFk5lRE/Clwms7k0nsy8/niebcC9wGvAB4tviRJNbFiKGTmf2Px6wEAv7jEc34X+N1Fth8Dtvazg9JqTU1NsXv3bg4ePMiWLVuq3h2p1mw0q9XOnj3Lzp07OX36NLt27eLs2bNV75JUa4aCWm3Pnj3Mzc2RmZw5c4a9e/dWvUtSrdV+5TUtb8eOHVXvQm3Nzs4yPT3NhQsXADh//jyHDh3i+PHjrFu3ruK9e6E73noCgPfduaPS/Th8+HCl/3xVz1BQa507d461a9dy6aWXltvm5+c5d+5chXsl1Zuh0HCe2bVEsUbz4d87XOluSF5TkCSVDAVJUslQUHWmpmDr1s53SbVgKKgaZ8/Czp1w+jTs2tX5WVLlDAVVY88emJuDTDhzBuwPSLXg9FHTNbGnMDsL09NQ9Ac4fx4OHYLjx6Fm/YGeOAGmFjEUNH7nzsHatdDVH2B+vrNdUqUMhabzLFXSEHlNQZJUMhQkSSVDQRqTqakptm7dypS9DNWYoSCNges6qCkMBWkMXNdBTeH0kWqjrWtD9LKuQ13WUxgG/3JvsxkK0oi5roOaxFBQbUz0GabrKagmvKYgSSoZCpKkkqEgDcjegdrIUJAGYO9AbWUoSAOwd6C2cvpIq9bWfsFSeukd9KuJPYWJnhZrMUNB6pO9A7WZoaBV84xxCOwpqCa8piBJKhkKkqSSoSBJKhkK0ohYblMTGQrSCFhuU1MZCtIIWG5TUzmSqoFMWmGtH4OU25pYXhsmx5rrw1CQhsxym5rMUNBAPLMbMstrqgmvKUiSSoaCJKlkKKhZpqZg69bOd0lDZyioOc6ehZ074fRp2LWr87OkoTIU1Bx79sDcHGTCmTPg7L80dE4fTYqm9wpmZ2F6GorZf86fh0OH4PhxGHBhm1pwiks1YyioGc6dg7VroWv2n/n5znZJQ2MoTArPSCX1wGsKkqSSoSBJKhkKklbFdSPaxVCQNDDXjWgfQ0HSwFw3on2cPpJWaRhrSzRxPYVB1o1os7b85WBDQdJAXDeinQwFaZWGcoboegqqCa8pSJJKhoIkqWQoSGPgLL+awlCQRsxZfjWJoSCNmLP8ahKnj9STYcziT6JeZ/mb2FMYt7b0AOrOUJBGyFl+NY2hoJ54ljZi9hRUE15TkCSVDAVJUslQkCSVDAVpAJbR1FaGgtQny2hqM0NB6pNlNLWZI6kamTYW3ka1sEwdy2uOIU8mQ0Hqg2U0tZ2hoJHxTLMPltdUE15TkCSVDAVJUmnFUIiIeyJiLiJOdW27NiKORsTJiHgkItYU2y+NiHuL7V+NiB1dz7mu2D4dEfsjIkZxQFKb2IfQuPXyTuE+4PoF2+4GbsvMNwCfAt5fbP8nAMX2XwL+XURc/GfcCewDri6+Fr6mpC72IVSFFUMhM48AzyzYvBk4Utx+DHhHcfv1wGeL580B/wvYFhHrgDWZeTQzE7gfuHG1Oy+1mX0IVWHQ6aNTwNuBh4CbgKuK7V8FboiIjxfbriu+XwBmup4/A6wf8J+tCdfG/sPCnsKo+hB158Ra9QYNhT3A/oi4HXgYmC+23wO8DjgGfBv4AvAcsNj1g1zqxSNiH52Pmti4ceOAuyg1l30IVWWgUMjMJ4G3AUTENcCuYvtzwG9ffFxEfAH4S+D7wIaul9gAPLXM6x8ADgBs27ZtyfDQZGrl2aQ9BdXEQCOpEXFl8f0S4EPAXcXPr4yIVxW3fwl4LjNPZ+Ys8GxEbC+mjm6h89GTJKlGVnynEBEPADuAyyNiBvgwcFlEvKd4yIPAvcXtK4HPRMQF4H8A7+p6qVvpTDK9Ani0+JIk1ciKoZCZNy9x10cWeexf05lMWux1jgFb+9k5SYOZmppi9+7dHDx4kC1btlS9O2oQG81Sy9hv0GoYClLL2G/QavhXUjV2bewZrNaw1lOY1H5D1do0EWcoSC1iv0GrZSho7Np0VjU09hRUE15TkCSVDAVJUslQkBrM9RY0bIaC1FD2ETQKhoLUUPYRNApOH6nR2tJ56LenYB/hhZxoGx5DQWog+wgaFUNBjdaaM0R7CqoJrylIkkqGgiSpZChIkkqGgtRSFts0CENBaiGLbRqUoSC1kMU2DcqRVLVeEwpuw1pkByy2XdSaceUxMxSklrHYptUwFNR6jThjtLymmvCagiSpZChIkkqGgjQm9gbUBIaCNAb2BtQUhoI0BvYG1BROH2nomtALGKdeegPD7Ck0VSOmxCaAoSCNmL0BNYmhoKHzjG8A9hRUE15TkCSVDAVJUslQkIbADoLawlCQVskOgtrEUJBWyQ6C2sTpI72IPYPeDWvtgknvKTixVh+GgrQKdhDUNoaCXsSztgrYU1BNeE1BklQyFCRJJUNBY+Ecv9QMhoJGzjl+qTkMBY2cc/xSc7R6+sh5++oNa46/7aroKThlpsW0OhRUPef4pWZpdSh4JqTGsKegmvCagiSpZChIkkqGgqRF2S2ZTIaCpBexWzK5DAVJL2K3ZHK1evpI6keVvZY6radgt2R5bZ9qNBQkvYDdkslmKEiFSs8A7SmoJrymIEkqGQqSpJKhIEkqGQrSKlnyUpsYCtIqWPJS2xgK0ipY8lLbOJKqyjV1MaRhlrzqUF5reylLvTEUpAFZ8lIbGQqqnGeoWF5TbXhNQZJUMhQkSSVDQdJY2euoN0NB0tjY66g/Q0HS2NjrqD+nj6QhGrRzUYeewqi5eE9/qprKMxQkjYW9jmYwFKQhGvjszp6CasJrCpKkkqEgSSqtGAoRcU9EzEXEqa5t10bE0Yg4GRGPRMSaYvvLIuKjxfYnIuJ3up5zXbF9OiL2R0SM5pA0SZx5l4arl3cK9wHXL9h2N3BbZr4B+BTw/mL7TcDLi+3XAf80IjYV990J7AOuLr4WvqbUF2fepeFbMRQy8wjwzILNm4Ejxe3HgHdcfDjwqoh4KfAKYB74QUSsA9Zk5tHMTOB+4MbV774mmTPv0vANOn10Cng78BCddwdXFds/AdwAzAKvBH47M5+JiG3ATNfzZ4D1A/6z1YOmrlHQq7bNvPfSU/CvyWocBg2FPcD+iLgdeJjOOwKAnwKeB/4u8Brg8xHxX4DFrh/kUi8eEfvofNTExo0bB9xFtZkz79JoDBQKmfkk8DaAiLgG2FXc9evAn2Xmj4C5iPjvwDbg88CGrpfYADy1zOsfAA4AbNu2bcnw0NI8q2wYewqqiYFGUiPiyuL7JcCHgLuKu74D/EJ0vArYDjyZmbPAsxGxvZg6uoXOR0+SpBrpZST1AeAosDkiZiJiL3BzRHwDeJLOGf+9xcP/A3AZnWsOfwHcm5lfK+67lc7U0jTwTeDRYR6IJGn1Vvz4KDNvXuKujyzy2P9D58LzYq9zDNja195JYmpqit27d3Pw4EG2bNlS9e6o5Ww0SzVmF0PjZihINWYXQ+PmX0lVI7S9d7FYT6FtXYzVcJpufAwFqabsYqgKhoIaofVnivYUVBNeU5AklQwFSVLJUJAklQwFSc0zNQVbt3a+a6gMBUnNcvYs7NwJp0/Drl2dnzU0hoKkZtmzB+bmIBPOnAELfUPlSKo0yZpWCpydhelpKAp9nD8Phw7B8ePQpEJfjUesDQVJzXHuHKxdC12FPubnO9s1FIaCNMlqfMaqanhNQZJUMhQkSSVDQaqhqakptm7dypRz+BozQ0GqGRfWUZUMBalmXFhHVXL6SI0wKYvs7L5980QvrNP6P5HeAIaCVCMurKOqGQpqhNafQRaL7HzXRXZUMa8pSJJKhoIkqWQoSCNgz0BNZShIQ2bPQE1mKEhDZs9ATeb0kUai7b2CpczOzg7UM7jYU3jfnTvGsJfj0fqJsZYyFKQhsmegpjMUNBKeJfap6CkctqeginlNQZJUMhQkSSVDQWoIuw8aB0NBagC7DxoXQ0FqALsPGhenj1R7k9B5WK6nMGj3oe2ccBsNQ0GqObsPGidDQbU3EWeE9hRUE15TkCSVDAVJUslQkCSVDAVJ4zc1BVu3dr6rVgwFSeN19izs3AmnT8OuXZ2fVRuGgqTx2rMH5uYgE86cAYt4teJIqlR3bSrvzc7C9DQURTzOn4dDh+D4cWhDEa8F49OGgqTxOXcO1q6FriIe8/Od7aoFQ0Gquxacfao5vKYgSSoZCpKkkqEgjZEL5ajuDAVpTFwoR01gKEhj4kI5agKnj1QrbV1QZ6WFcpZbZKeJJuLPnbeUoSCNgQvlqCkMBdXKxJ5husiOasJrCpKkkqEgSSoZCpJK9ihkKEgC7FGow1CQBNijUIfTR2qktvUZqu4prNSj0HA0YbrOUJBkj0IlQ0GN1IQzrr7YU1BNeE1BklQyFCRJJUNBEmBHQR2GgiQ7CioZCpLsKKjk9JFaq0ldhip7CnYUBte6KTgMBWni2VFQN0NBrdWoszh7CqoJrylIkkqGgiSpZChIkkorhkJE3BMRcxFxqmvbtRFxNCJORsQjEbGm2P4bEXGi6+tCRLypuO+64vHTEbE/ImJkRyXVmCUx1Vkv7xTuA65fsO1u4LbMfAPwKeD9AJn5scx8U2a+CXgX8NeZeaJ4zp3APuDq4mvha0qtZ0lMdbdiKGTmEeCZBZs3A0eK248B71jkqTcDDwBExDpgTWYezcwE7gduHHCfpcayJKa6G3Qk9RTwduAh4CbgqkUesxu4obi9Hpjpum+m2KYJ16SC2WotVxI7+I/PANUtsjMujRoTnlCDhsIeYH9E3A48DMx33xkRPw2cy8yL1yEWu36QS714ROyj81ETGzduHHAXpXqxJKYmGCgUMvNJ4G0AEXENsGvBQ95J8dFRYQbY0PXzBuCpZV7/AHAAYNu2bUuGh5rPM8eC5TXVxEAjqRFxZfH9EuBDwF1d911C5yOlj1/clpmzwLMRsb2YOrqFzkdPkqQa6WUk9QHgKLA5ImYiYi9wc0R8A3iSzhn/vV1P+XlgJjO/teClbqUztTQNfBN4dAj7L0kaohU/PsrMm5e46yNLPP4wsH2R7ceArf3snNQmU1NT7N69m4MHD7Jly5aqd0dalI1maQzsJ6gpDAVpDOwnqCn809kayCT1C1arl0VsqlxkpymcVBsPQ0EaMfsJahJDQQPxrG3I7CmoJrymIEkqGQqSpJKhIKlvrgnRXoaCpL7YuWg3Q0FSX+xctJvTR9IqDaOz0ZSeQi+dizabhKk7Q0FSz+xctJ+hIK3SUM4e7SmoJrymIEkqGQqSpJKhoInhbL20MkNBE8HZeqk3hoImgrP1Um+cPtKKmr52QhNm6/vtKUzCvLyqYSio9Zytl3pnKGhFnpWOgT0F1YTXFCRJJUNBklQyFKQJY19DyzEUpAliX0MrMRSkCWJfQytx+kit1pSOxTjWU2hCX6NqTtoZCtLEsK+hXhgKarXGnPnZU1BNeE1BklQyFCRJJUNBklQyFKRVsgymNjEUpFWwDKa2MRSkVbAMprZxJFVD1ZSy2DAMsww2jvLaKDRm5Fc9MxSkAVkGUxsZChoqzxwHZHlNNeE1BUlSyVCQJJUMBWmI7Cyo6QwFaUjsLKgNDAVpSOwsqA2cPtJQTVJPodtqOwtN7Sksx0m0ZjIUpCGws6C2MBQ0VJ4dDsiegmrCawqSpJKhIEkqGQqSFmXnYjIZCpJexM7F5DIUJL2InYvJ5fSRtIxx9S7q1FMY5joRTTeJ03SGgqQXsHMx2QwFaRljO1O0p6Ca8JqCJKlkKEiSSoaCtErO86tNDAVpFZznV9sYCtIqOM+vtnH6qEEmda2CuhrmPH+degptNIl9g0EZCtKAnOdXGxkKDeLZTovZU1BNeE1BklQyFCRJJUNBklQyFKQJZ/lO3QwFaYJZvtNChoI0wSzfaSFHUqUhGrRgWEV5rQ2L6TimPXyGgjShLN9pMYaCNEQDn7laXlNNeE1BklQyFCRJpRVDISLuiYi5iDjVte3aiDgaEScj4pGIWNN13xuL+6aK+3+s2H5d8fN0ROyPiBjNIUnNcbEj4Cio6qKXdwr3Adcv2HY3cFtmvgH4FPB+gIh4KfAnwLszcwuwA/hR8Zw7gX3A1cXXwteUJkp3R+DkyZM8//zzVe+StHIoZOYR4JkFmzcDR4rbjwHvKG6/DfhaZn61eO7fZObzEbEOWJOZRzMzgfuBG4ew/1JjdXcE5n80z9e//vWqd0kaeProFPB24CHgJuCqYvs1QEbEZ4ArgI9n5r8F1gMzXc+fKbappVwQaHkLOwIXLiRzTz/Nr27e3JiOQFXsJozWoKGwB9gfEbcDDwPzXa/3s8BPAueAz0bE48APFnmNXOrFI2IfnY+a2Lhx44C7KNXXwo7Aj738f3Ih046AKjdQKGTmk3Q+KiIirgF2FXfNAJ/LzO8V930aeAud6wwbul5iA/DUMq9/ADgAsG3btiXDQ/Xl2Vyfip7Cd+84XOluSAONpEbElcX3S4APAXcVd30GeGNEvLK46PwPgdOZOQs8GxHbi6mjW+h89CRJqpFeRlIfAI4CmyNiJiL2AjdHxDeAJ+mc8d8LkJnfB/4A+AvgBPCVzPzPxUvdSmdqaRr4JvDocA9FkrRaK358lJk3L3HXR5Z4/J/Q+bho4fZjwNa+9k5qsampKXbv3s3BgwfZUvXOSAUbzVIFFq5jYEdBdWEoSBVYuI6BHQXVhX8lVWNhb+H/W2wdg7mnz/Psl7/MB8a4nkLdOcFWDUNBGrPF1jG49NJZLlzwIyRVz1DQWHjWt4KL6yn84eFKd0PymoIkqWQoSJJKhoKkcl2HqampqndFFTMUpAm3sDPhgj+TzVCQJtzCzsTevXur3iVVyOkjaYgG7WPc8dYTALxvzD2FxToThw4d4vjx441a18HptuExFKQJtlhnYn5+3nUdJpihIA3RwGesF3sKvzfg86Uh8ZqCJKlkKEiSSoaCJKlkKEhjYkFMTWAoSGNgQUxNYShIY2BBTE3hSKp64iI5g+ulIFZVea3JLKyNhqEgjZgFMTWJoaCeeFY2YpbXVBNeU5AklQwFSVLJUNDEsjcgvZihoIlkb0BanKGgiWRvQFqc00c1ZCdgtOq4sEzbewpOrzWHoaCJY29AWpqhUEOeVU0gewqqCa8pSJJKhoIkqWQoSC1lD0ODMBSkFrKHoUEZClIL2cPQoJw+krpU1REZZk+hjj2MYXAqbzwMBall7GFoNQwFqUtlZ6P2FFQTXlOQJJUMBUlSyVDQxHKOX3oxQ0ETyTl+aXGGgiaSc/zS4pw+mgCuz/BCdZzjr+N6CvYCJpOhoInjHL+0NENhAnjG1wD2FFQTXlOQJJUMBUlSyVCQJJUMBalmLNWpSoaCVCOW6lQ1Q0GqEUt1qpojqRoaS3KDu+OtJ/jh/Dyf+MQPa1WqawNHsvtjKEg1ceHC85bqVDlDQUPjGdkqFOW1p/7wcKW7IXlNQZJUMhQkSSVDQVoFOwVqG0NBGpCdArWRoSANyE6B2sjpowaxB1Afw16op46L7DSZk3CDMxSkAbhQj9rKUGgQz35azEV2VBNeU5AklQwFSVLJUJAazq6EhslQkBrMroSGzVCQGsyuhIbN6SPV3iT0MwbpKQy7K9F2Tu/1xlCQGsquhEbBUFDtTcQZnj0F1YTXFCRJJUNBklQyFCRJpRVDISLuiYi5iDjVte3aiDgaEScj4pGIWFNs3xQR/zciThRfd3U957ri8dMRsT8iYjSHJEkaVC/vFO4Drl+w7W7gtsx8A/Ap4P1d930zM99UfL27a/udwD7g6uJr4WtKkiq2Yihk5hHgmQWbNwNHituPAe9Y7jUiYh2wJjOPZmYC9wM39r23kqSRGvSawing7cXtm4Cruu57bUQcj4jPRcTPFdvWAzNdj5kpti0qIvZFxLGIOPb0008PuIuSpH4NGgp7gPdExOPAq4H5YvsssDEz3wz8C+A/FtcbFrt+kEu9eGYeyMxtmbntiiuuGHAXJUn9Gqi8lplPAm8DiIhrgF3F9h8CPyxuPx4R3wSuofPOYEPXS2wAnhp8tyVJozDQO4WIuLL4fgnwIeCu4ucrIuIlxe2foHNB+VuZOQs8GxHbi6mjW4CHhrD/kqQhWvGdQkQ8AOwALo+IGeDDwGUR8Z7iIQ8C9xa3fx74VxHxHPA88O7MvHiR+lY6k0yvAB4tviRJNbJiKGTmzUvc9ZFFHvtJ4JNLvM4xYGtfeydJGisbzZKkkqEgSSoZCpKkkqEgSSoZCpKkkqEgSSpF5+/T1VdEPA18u+r9WMblwPeq3olV8hjqwWOoh7Ycw6sys++/E1T7UKi7iDiWmduq3o/V8BjqwWOoh0k/Bj8+kiSVDAVJUslQWL0DVe/AEHgM9eAx1MNEH4PXFCRJJd8pSJJKhkIPIuKeiJiLiFMrPO4nI+L5iPi1ce1br3o5hojYEREnImIqIj43zv3rxUrHEBF/KyIeiYivFsfwW+Pex5VExFUR8ecR8USxj+9d5DEREfsjYjoivhYRb6liX5fS4zH8RrHvX4uIL0TEtVXs61J6OYaux9by97rXY+j79zoz/Vrhi846EW8BTi3zmJcA/xX4NPBrVe9zv8cA/Dhwms5yqgBXVr3PAxzDB4F/U9y+AngGuLTq/V6wj+uAtxS3Xw18A3j9gsfspLPeSADbgS9Vvd8DHMM/AF5T3P6VJh5DcV9tf697/PfQ9++17xR6kJlH6PwPZjn/nM5aEnOj36P+9XAMvw48mJnfKR5fu+Po4RgSeHWxut9lxWOfG8e+9SozZzPzK8XtZ4EngPULHnYDcH92fBH48YhYN+ZdXVIvx5CZX8jM7xc/fpEXLsdbuR7/PUCNf697PIa+f68NhSGIiPXAr1IsS9pQ1wCviYjDEfF4RNxS9Q4N4N8Dr6Oz/vdJ4L2ZeaHaXVpaRGwC3gx8acFd64Hvdv08w+L/w6rcMsfQbS81XmlxqWNo0u/1Mv8e+v69XnHlNfXkDuADmfl85yS1kV4KXAf8Ip0lU49GxBcz8xvV7lZffhk4AfwC8PeAxyLi85n5g0r3ahERcRmdM9D3LbJ/i/1HVLsxwRWO4eJj/hGdUPjZce5br1Y4hjtowO/1CsfQ9++1oTAc24CPF//hXA7sjIjnMvM/VbpX/ZkBvpeZZ4GzEXEEuJbO55RN8VvA72fnw9PpiPgr4O8DX652t14oIl5G55f4Y5n54CIPmQGu6vp5A513P7XRwzEQEW8E7gZ+JTP/Zpz714sejqH2v9c9/rfU1++1Hx8NQWa+NjM3ZeYm4BPAP6vTfzg9egj4uYh4aUS8EvhpOp9RNsl36JwRERFrgc3AtyrdowWK6x1/DDyRmX+wxMMeBm4pppC2A/87M2fHtpMr6OUYImIj8CDwrjq+2+zlGOr+e93jf0t9/177TqEHEfEAsAO4PCJmgA8DLwPIzNp/3ggrH0NmPhERfwZ8DbgA3J2Zy47gjlsP/x7+NXBfRJyk8xHMBzKzbn/t8meAdwEnI+JEse2DwEYoj+PTdCaQpoFzdN4B1Ukvx3A78HeAPyrOtJ/Lev2RuV6Ooe5WPIZBfq9tNEuSSn58JEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpNL/A3hjXSXSjs/qAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1837/2000=0.92\n"
     ]
    }
   ],
   "source": [
    "confidence_interval_mean(2, 0.5, 12, 2000, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pacific-occasion",
   "metadata": {},
   "source": [
    "What we were expecting happened, as the true mean is only contained in our confidence intervals 92% of the time. This is not the intended, as we are calculating the 95% confidence intervals for the mean. \n",
    "\n",
    "In the second experiment, the critical value used is now calculated from a t-distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "integrated-montgomery",
   "metadata": {},
   "outputs": [],
   "source": [
    "def confidence_interval_mean_t(μ, σ, n, number_trials, N, ci=0.95, sample='random'):\n",
    "    x_ = norm.rvs(loc=μ, scale=σ, size=N)\n",
    "    x_hat_list = []\n",
    "    SE_hat_x_list = []\n",
    "    \n",
    "    if sample!='random':\n",
    "        # Inducing bias on the sampling\n",
    "        x_.sort()\n",
    "        x_ = x_[:-int(0.2*N)]\n",
    "        np.random.shuffle(x_)\n",
    "    \n",
    "    for i in range(number_trials):\n",
    "        s_ = np.random.choice(x_, n, replace=False)\n",
    "        x_hat = np.mean(s_)\n",
    "        x_hat_list.append(x_hat)\n",
    "        SE_hat_x_list.append(t.ppf(ci+(1-ci)/2, df=n-1)*np.std(s_)/np.sqrt(n))\n",
    "\n",
    "    j=0\n",
    "    _, ax = plt.subplots(1, 1, figsize=(6, 8))\n",
    "    for i in range(len(x_hat_list)):\n",
    "        if (μ>x_hat_list[i]-SE_hat_x_list[i]) & (μ<x_hat_list[i]+SE_hat_x_list[i]):\n",
    "            # interval contains p\n",
    "            if i > len(x_hat_list)-50:\n",
    "                ax.errorbar(x_hat_list[i], np.arange(len(x_hat_list))[i],lolims=True, xerr=SE_hat_x_list[i], yerr=0.0, linestyle='', c='black')\n",
    "            j +=1\n",
    "        else:\n",
    "            # interval does not contain p\n",
    "            if i > len(x_hat_list)-50:    \n",
    "                ax.errorbar(x_hat_list[i], np.arange(len(x_hat_list))[i],lolims=True, xerr=SE_hat_x_list[i], yerr=0.0, linestyle='', c='red')\n",
    "    ax.axvline(μ, color='darkorange')\n",
    "    #plt.xlim(0,1)\n",
    "    plt.show()\n",
    "    print(f'{j}/{number_trials}={np.round(j/number_trials,2)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "lucky-malpractice",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAHSCAYAAAAAMxtqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdjklEQVR4nO3df6zdd33f8ec7gCkQvLI5QZ4dY1YlHjgQwG5r1h/zWpWmRCKpaGTSilRNNI8MTdBNiBQhkDZVLZPGgrU1kZUmISoLLhCaRCONMlZjNgypAYN9E5saSuEuHg4NG1l+YHL93h/n6/Tk5tx7zj33nPP98Xk+pKt77vece/L5KtHL33z8ep9vZCaSpDKcU/cCJEmzY+hLUkEMfUkqiKEvSQUx9CWpIIa+JBXk+XUvYJh169bl5s2b616GtHo/ON77/rIt9a5Dnbdu3Truu++++zLz0sXPNT70N2/ezKFDh+pehrR6+3b2vu/aX+cqVIiIWDfouNs7klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSBDQz8iLoiIv4iIhyJiLiLeVR3/+xFxf0T8VfX9ZX2/83sRcSIijkfEr/Yd3xYRR6rn9kRETOe0JEmDjHKl/zTwbzLzVcAO4J0R8WrgeuCzmXkh8NnqZ6rn3gZsBS4F/iginle9143AbuDC6us5gwOSpOkZGvqZeTIzv1I9fgx4CNgAXA58tHrZR4ErqseXAx/PzB9l5l8DJ4CfiYj1wNrMPJi9O7fc3vc7kqQZWNGefkRsBl4PfAl4eWaehN4fDMD51cs2AN/t+7X56tiG6vHi44P+Obsj4lBEHHrkkUdWskRJ0jJGDv2IOBf4FPDuzPzhci8dcCyXOf7cg5l7M3N7Zm4/77zzRl2iJGmIkUI/Il5AL/A/lpl3Voe/V23ZUH0/VR2fBy7o+/WNwMPV8Y0DjkuSZmSU9k4Afww8lJkf7nvqbuC3q8e/DdzVd/xtEfHCiHglvb+wfaDaAnosInZU73l13+9IkmZglE/Z/Dng7cCRiDhcHXsf8IfAn0bEtcB3gCsBMnMuIv4UeJBe8+edmblQ/d51wG3Ai4B7qy9J0owMDf3M/B8M3o8H+OUlfuf3gd8fcPwQcPFKFijVZW5ujl27drFv3z62bt1a93KkiXAiVxrg8ccf581vfjMPPvggl112GY8//njdS5ImwtCXBrjmmms4deoUmcn3vvc9rr322rqXJE1E4++c1WY7d+6sewkaw8mTJzlx4gRnzpwB4KmnnuITn/gEX/3qV1m/fv3Y73vDGw8D8O4bd05glRrF/v37615C4xj60iJPPPEEL3/5y1mzZs0zx06fPs0TTzxR46qkyTD0p8irDD1LdY/c/X+wv9ZlqGzu6UtSQQx9SSqIoS8tMjc3x8UXX8zc3FzdS5EmztCX+tjPV9cZ+lIf+/nqOts7DWXHf/am1c8/y57+aGy9TZehL1Xs56sEhn5DebXTQfb01QDu6UtSQQx9SSqIoS+tgB1+tZ2hL43IDr+6wNCXRmSHX11ge6dD7PZPzyQ6/Pb0R2d7bXoMfWkEdvjVFYZ+h3h11HD29NUA7ulLUkEMfUkqiKEvrZBdfbWZoS+tgF19tZ2hL62AXX21ne2dgtnrX5nVdvXt6ffYMquXoS+NyK6+usDQL5hXXDNmT18N4J6+JBXE0Jekghj6klQQQ18S4NBZKQx9SQ6dFcTQl+TQWUGsbEpjGGewranDWZO4QUyTWU1+NkNfKpxDZ2Ux9KUxjHX16HCWGsA9fUkqiKEvSQUx9KUWs1uvlTL0pZayW69xGPpSS9mt1zhs76gWJd7AZZI9/a536+tSQqff0JdayG69xmXoqxYlXFE9hz19NYB7+pJUEENfkgpi6EvLsAevrjH0pSXYg1cXGfrSEuzBq4ts7zRUiT32JplGD76pn6ffJkW2vibM0JcGsAevrjL0G8ormg6yp68GcE9fkgpi6EtSQQx9tY7deWl8hr5axe68tDqGvlrF7ry0OrZ3WqjUDn/bP0N+Gj19W15aKUNfrWF3Xlo9Q7+FvLprKXv6agD39CWpIIa+JBXE0Jekghj6UsM4fKZpMvSlBnH4TNNm6EsN4vCZps3KpkZW6lDYpAwbzmr78FkblVh/NvSlhnD4TLNg6GtkJV4VTZTDWWoA9/QlqSCGviQVZGjoR8QtEXEqIo72HbskIg5GxJGIuCci1lbH10TErdXxr0XEzr7f2VYdPxEReyIipnFCmj575FJ7jXKlfxtw6aJjNwPXZ+ZrgE8D76mO/3OA6vivAP8hIs7+M24EdgMXVl+L31MtYI9carehoZ+ZB4BHFx3eAhyoHt8PvLV6/Grgs9XvnQL+D7A9ItYDazPzYGYmcDtwxWoXr9mzRy6127jtnaPAW4C7gCuBC6rjXwMuj4iPV8e2Vd/PAPN9vz8PbBjznz0ye+WTZY98dUa9iYotKU3TuKF/DbAnIj4A3A2cro7fArwKOAT8DfAF4Glg0P59LvXmEbGb3lYQmzZtGnOJmjR75FL7jRX6mXkMeBNARFwEXFYdfxr43bOvi4gvAH8F/ADY2PcWG4GHl3n/vcBegO3bty/5h8MwXjGpUezpqwHGqmxGxPnV93OA9wM3VT+/OCJeUj3+FeDpzHwwM08Cj0XEjqq1czW9rSFJ0gwNvdKPiDuAncC6iJgHPgicGxHvrF5yJ3Br9fh84L6IOAP8L+DtfW91Hb0m0IuAe6svSdIMDQ39zLxqiac+MuC136bX7Bn0PoeAi1eyOKmL5ubm2LVrF/v27WPr1q11L0eFcSJXmqGFhQXnHFQrQ1+aoWPHjzvnoFr5KZsaizMQK/ehbQ/w5JNP8tRTvZ+dcxifzbzxGfrSjJw5s8CaNWt4xSv+LuCdc9CsGfoai1daY6h6+t++YX+ty1DZ3NOXpIIY+pJUEENfRfPeACqNoa9ieW8AlcjQV7G8N4BKZHunY+zPj6aOewOM+nn6bWB7q70MfRXJewOoVIZ+x3gF1mB+nr4awD19SSqIoS9JBTH0Jakghr40ZWcHwJwDUBMY+tIU9Q+AHTlyhIWFhbqXpMIZ+tIU9Q+Anf7xaY4fP173klQ4K5uaiRKHxhYPgJ05k5x65BF+fcuWom6aYo24WQx9aUoWD4D9xAv/N2cyHQBTrQx9zYRXezwznPVdb6KiGrmnL0kFMfQlqSCGvlQobyBTJkNfKpA3kCmXoS8VyBvIlMv2jlqrbd3/ptxEpY4byHRdm9pphr5UGG8gUzZDX63VpqsrwJuoqBHc05ekghj6klQQQ1/qY3ddXWfoSxW76yqBoS9V7K6rBLZ3Oqht/fUmmEV3vSk9/WlrXauqMIa+hN11lcPQ7yCvtBrKnr4awD19SSqIoS9JBTH0pZZxlkCrYehLLeIsgVbL0JdaxFkCrZbtHdWqpJmC1fb0/Rz81bHV1mPoSy3hLIEmwdBXrYq6+rKnrwZwT1+SCmLoS1JBDH2pjx14dZ2hL1XswKsEhr5UsQOvEtje6bCSOvCr5efpr05RLayWM/Ql7MCrHIZ+h3n11TD29NUA7ulLUkEMfUkqiKEvSQUx9KWaORCmWTL0pRo5EKZZM/SlGjkQplmzsqllOeA1OYuHs7wpyuxZYzb0pdo4EKY6GPpalldGE+RwlhrAPX1JKoihL0kFMfQ1FrvlUjsZ+loxu+VSexn6WjG75VJ72d7BLvpK2C0f39me/ut21bsOlc3Q14rYLZfazdDHLrpmpOrpS3VyT1+SCmLoS1JBhoZ+RNwSEaci4mjfsUsi4mBEHImIeyJibXX8BRHx0er4QxHxe32/s606fiIi9kRETOeU1HR2/KX6jHKlfxtw6aJjNwPXZ+ZrgE8D76mOXwm8sDq+DfgXEbG5eu5GYDdwYfW1+D1VADv+Ur2Ghn5mHgAeXXR4C3Cgenw/8NazLwdeEhHPB14EnAZ+GBHrgbWZeTAzE7gduGL1y1fb2PGX6jVue+co8BbgLnpX9xdUxz8JXA6cBF4M/G5mPhoR24H5vt+fBzaM+c+unb3+8ZTe8V/8efqrZetM4xg39K8B9kTEB4C76V3RA/wMsAD8Q+BlwOcj4r8Bg/bvc6k3j4jd9LaC2LRp05hLVNPY8ZfqN1boZ+Yx4E0AEXERcFn11G8Cf56ZPwZORcT/BLYDnwc29r3FRuDhZd5/L7AXYPv27Uv+4VAXr7A0Fj9PXw0wVmUzIs6vvp8DvB+4qXrqO8AvRc9LgB3Ascw8CTwWETuq1s7V9LaGJEkzNEpl8w7gILAlIuYj4lrgqoj4BnCM3hX7rdXL/zNwLr09/78Ebs3Mr1fPXUev9XMC+CZw7yRPRJI03NDtncy8aomnPjLgtf+P3l/sDnqfQ8DFK1qdpFaYm5tj165d7Nu3j61bt9a9HC3DiVxJq+LsRbsY+pJWxdmLdvFTNqUhJjWXMemefhOUPnsxjrrbf4a+pLE5e9E+hr40xMSuzOzpqwHc05ekghj6klQQQ1+SCmLoSzXwRjKqi6EvzZjDTKqToS/NmMNMqpOVzYJ5M5jZuuGNh/nR6dN88pM/cpipIeoelKqDoS/N0JkzCw4zqVaGfsFKvMqpVTWc9fB/3F/rMlQ29/QlqSCGviQVxNCXpsg+vprG0JemZHEff2Fhoe4lSYa+NC2L+/jHjx+ve0mS7R0tzy7/eAbdXOTUI0/x2AMP8N4O3URllmybTYahL03BoJuLrFlzkjNn3OJRvQx9Lcurqwk6exMVe/qqkXv6klQQQ1+SCmLoSy3nLIBWwtCXWszP5tdKGfpSi/nZ/Fop2zvqpCbOF9zwxsMAvHtCPf1BswClfDa/rbLxGfpSSw2aBfCz+TWMoa9OauSV4Nme/h/sr3UZKpt7+pJUEENfkgpi6EurYEdebWPoS2OyI682MvSlMdmRVxvZ3tHENbEjP2njdOQn3dOvQyNbUVoRQ18agx15tZWhr4nzanAJ9vTVAO7pS1JBDH1JKoihL0kFMfRVC4eapHoY+po5h5qk+hj6mjmHmqT6WNkcooRBo1kq+cYfsxzOsjarpRj6mimHmqR6GfpDeMWkiXE4Sw3gnr4kFcTQl6SCGPpSYZyRKJuhLxXEGQkZ+lJBnJGQ7R21WpvmKOq+iUrJMxKT0oU2n6EvFcIZCYGhr5Zr1ZWXPX01gHv6klQQQ1+SCmLoSw1jj17TZOhLDWKPXtNm6EsNYo9e02Z7R7VoU79+Uob19O3Rj65Vra2GMfSlhrBHr1kw9FWLIq/U7OmrAdzTl6SCGPqSVBBDXxqDXXq1laEvrZBderWZoS+tkF16tZntncKV2JdfjdV06ev+PP0mKbK91RCGvrQCdunVdoZ+4bzimiF7+moA9/QlqSCGviQVxNCXpIIMDf2IuCUiTkXE0b5jl0TEwYg4EhH3RMTa6vhvRcThvq8zEfG66rlt1etPRMSeiIipnZXUEQ6BadJGudK/Dbh00bGbgesz8zXAp4H3AGTmxzLzdZn5OuDtwLcz83D1OzcCu4ELq6/F7ympj0NgmoahoZ+ZB4BHFx3eAhyoHt8PvHXAr14F3AEQEeuBtZl5MDMTuB24Ysw1S0VwCEzTMG5l8yjwFuAu4ErgggGv2QVcXj3eAMz3PTdfHVPLOdw1upUMZ3lDldkrpb48buhfA+yJiA8AdwOn+5+MiJ8FnsjMs38PMGj/Ppd684jYTW8riE2bNo25RKm9HALTtIwV+pl5DHgTQERcBFy26CVvo9raqcwDG/t+3gg8vMz77wX2Amzfvn3JPxxUv1KujibC4Sw1wFiVzYg4v/p+DvB+4Ka+586ht+Xz8bPHMvMk8FhE7KhaO1fT2xqSJM3QKJXNO4CDwJaImI+Ia4GrIuIbwDF6V+y39v3KLwLzmfmtRW91Hb3Wzwngm8C9E1i/JGkFhm7vZOZVSzz1kSVevx/YMeD4IeDilSxOmoa5uTl27drFvn372Lp1a93LkWbKiVwVxe67Smfoqyh231U6P1q5YKV17OvuvjfhJiq2rWToqxh23yVDv2he9c2YPX01gHv6klQQQ1+SCmLoSx3nZ/Krn6EvdZhzCVrM0Jc6zLkELWZ7R6tSWtd/NWbd0697LkErN4tGnaEvdZRzCRrE0Neq2PVfAXv6agD39CWpIIa+JBXE0FeR7K6rVIa+imN3XSUz9FUcu+sqme2dQpXar6+zu96Ez9NfzPZVeQx9FcXuukpn6BfKK7wa2NNXA7inL0kFMfQlqSCGvlSxu68SGPoSdvdVDkNfwu6+ymF7R89RWod/Vt39Jvb0+9noKoOhr+LZ3VdJDH09h1d8U2JPXw3gnr4kFcTQl6SCGPqSVBBDX9JIHF7rBkNf0lAOr3WHoS9pKIfXusPKplqvLcNkTR/OWkqdN54p2bSq04a+pGU5vNYthr5arzXDZA5nqQHc05ekghj6klQQQ19agr10dZGhLw1gL11dZehLA9hLV1fZ3lFreu6zMq1eelt7+ktpTWtKz2LoS4vYS1eXGfryim1W7OmrAdzTl6SCGPqSVBBDX2oZ5we0Goa+1CLOD2i1DH2pRZwf0GrZ3tHU2P9/ttX29P1c++kroclm6Est4fyAJsHQ19SUcNW0Ivb01QDu6UtSQQx9SSqIoS/pOZwF6C5DX9KzOAvQbYa+pGdxFqDbbO9IlWnPFbTh8/SdBVhaV9pohr6kZzgL0H2GvlSZ+pWcPX01gHv6klQQQ1+SCmLoS1JBDH0Vz0EklcTQV9EcRFJpDH0VzUEklcbKZsN5I5LpmfUgUhuGs2ahK0NObWXoq1gOIqlEhn7DeVXUIQ5nqQHc05ekghj6klSQoaEfEbdExKmIONp37JKIOBgRRyLinohY2/fca6vn5qrnf6I6vq36+URE7ImImM4pSQLnDzTYKFf6twGXLjp2M3B9Zr4G+DTwHoCIeD7wJ8A7MnMrsBP4cfU7NwK7gQurr8XvKWlCnD/QUoaGfmYeAB5ddHgLcKB6fD/w1urxm4CvZ+bXqt/928xciIj1wNrMPJiZCdwOXDGB9UsawPkDLWXc9s5R4C3AXcCVwAXV8YuAjIj7gPOAj2fmvwc2APN9vz9fHZNWpU1zDLPq6XsjlJUprSE3buhfA+yJiA8AdwOn+97v54GfBp4APhsRXwZ+OOA9cqk3j4jd9LaC2LRp05hLlMrk/IGWM1boZ+Yxels5RMRFwGXVU/PA5zLz+9VznwHeQG+ff2PfW2wEHl7m/fcCewG2b9++5B8OUquu0uzpqwHGqmxGxPnV93OA9wM3VU/dB7w2Il5c/aXuPwUezMyTwGMRsaNq7VxNb2tIkjRDo1Q27wAOAlsiYj4irgWuiohvAMfoXbHfCpCZPwA+DPwlcBj4Smb+1+qtrqPX+jkBfBO4d7KnIkkaZuj2TmZetcRTH1ni9X9Cbztn8fFDwMUrWp1UgLm5OXbt2sW+ffvYunVr3ctRxzmRK9XIPr1mzdCXamSfXrPmp2yq8drUxV/O4p5+aX36VjWtOszQl2pin151MPTVeJ25QrSnrwZwT1+SCmLoS1JBDH2pMH7OftkMfakgzgXI0JcK4lyAbO+oMbrSx1/KrD5PfymlzQU0RdPaZ4a+VAjnAgSGvhqkaVdEE2dPXw3gnr4kFcTQl6SCGPqSVBBDXyqYg1rlMfSlQjmoVSZDXyqUg1plsrKpos1yIKzu4ax+Dmo9W+frwn0MfalADmqVy9BX0WZ6hedwlhrAPX1JKoihL0kFMfTVWXbQpecy9NVJdtClwQx9dZIddGkw2zt6jrbfzKSpHfTV9PRL6pFrugx9dY4ddGlphr6ew6vKKbGnrwZwT1+SCmLoS1JBDH1pxpwfUJ0MfWmGFhYWnB9QrQx9aYaOHT/u/IBqZXtHU9P2vv+kfWjbAzz55JM89VTv56bMDzSN7bHpMvSlGTlzZoE1a9bwilf8XcA7P6BZM/Q1NV6xLVL19L99w/5al6GyuacvSQUx9CWpIIa+GsHuujQbhr5q52ffS7Nj6Kt2fva9NDudbu/YE2++pn72/TSs5vP0m8ZmVnt1OvTVfH72vTRbnQ59r0bUKH6evhrAPX1JKoihL0kFMfQlqSCGvqRncVCu2wx9Sc9wUK77DH1Jz3BQrvs6XdmU+tU9rNf04aySBuUmqW3VcENfEuCgXCkMfRWj9isyh7PUAO7pS1JBDH1JKoihL/Wxo66uM/Slih11lcDQlyp21FUC2zsNV3e3vBSz6Kg3vadfl9pbVYUx9CXsqKschn7DeRXUIfb01QDu6UtSQQx9SSqIoa9Osm8vDWboq3Ps20tLM/TVOfbtpaXZ3mkZe/vLa/Jnwje9p29TrAyGvjrFvr20PEO/ZbwaazF7+moA9/QlqSCGviQVxNCXpIIMDf2IuCUiTkXE0b5jl0TEwYg4EhH3RMTa6vjmiHgyIg5XXzf1/c626vUnImJPRMR0TkmStJRRrvRvAy5ddOxm4PrMfA3waeA9fc99MzNfV329o+/4jcBu4MLqa/F7SpKmbGjoZ+YB4NFFh7cAB6rH9wNvXe49ImI9sDYzD2ZmArcDV6x4tZKkVRl3T/8o8Jbq8ZXABX3PvTIivhoRn4uIX6iObQDm+14zXx0bKCJ2R8ShiDj0yCOPjLlESdJi44b+NcA7I+LLwEuB09Xxk8CmzHw98K+B/1Lt9w/av8+l3jwz92bm9szcft555425REnSYmMNZ2XmMeBNABFxEXBZdfxHwI+qx1+OiG8CF9G7st/Y9xYbgYfHX7YkaRxjXelHxPnV93OA9wM3VT+fFxHPqx7/I3p/YfutzDwJPBYRO6rWztXAXRNYvyRpBYZe6UfEHcBOYF1EzAMfBM6NiHdWL7kTuLV6/IvAv42Ip4EF4B2ZefYvga+j1wR6EXBv9SVJmqGhoZ+ZVy3x1EcGvPZTwKeWeJ9DwMUrWp0kaaKcyJWkghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFid7nnzVXRDwC/M2E3m4d8P0JvVeTdPG8unhO4Hm1TVvP6/sAmfmcTzNufOhPUkQcysztda9j0rp4Xl08J/C82qaL5+X2jiQVxNCXpIKUFvp7617AlHTxvLp4TuB5tU3nzquoPX1JKl1pV/qSVLTOhX5E3BIRpyLi6JDX/XRELETEb8xqbasxynlFxM6IOBwRcxHxuVmub1zDzisi/l5E3BMRX6vO63dmvcaViogLIuIvIuKhas3vGvCaiIg9EXEiIr4eEW+oY60rMeJ5/VZ1Pl+PiC9ExCV1rHUlRjmvvte2KjcGysxOfdH7TP83AEeXec3zgP8OfAb4jbrXPInzAn4SeJDe7SoBzq97zRM6r/cBH6oenwc8Cqype91Dzmk98Ibq8UuBbwCvXvSaN9O7p0QAO4Av1b3uCZ3XPwFeVj3+ta6cV/Vc63Jj0FfnrvQz8wC9YFjOv6L3uf+npr+iyRjhvH4TuDMzv1O9vhXnNsJ5JfDS6o5r51avfXoWaxtXZp7MzK9Ujx8DHgI2LHrZ5cDt2fNF4CcjYv2Ml7oio5xXZn4hM39Q/fhFnn2b1EYa8d8XtDA3Bulc6A8TERuAX6e6xWOHXAS8LCL2R8SXI+Lquhc0If8JeBW9eyofAd6VmWfqXdLoImIz8HrgS4ue2gB8t+/neQYHTSMtc179rqVld8hb6ry6lBtj3Ri95W4A3puZC72Lx854PrAN+GV6t6Q8GBFfzMxv1LusVftV4DDwS8BPAfdHxOcz84e1rmoEEXEuvSvDdw9Y76D/+FpRpRtyXmdf88/ohf7Pz3JtqzHkvG6gI7lRYuhvBz5e/YtbB7w5Ip7OzD+rdVWrNw98PzMfBx6PiAPAJfT2J9vsd4A/zN6m6omI+GvgHwMP1Lus5UXEC+gFyMcy884BL5kHLuj7eSO9/5tptBHOi4h4LXAz8GuZ+bezXN+4RjivzuRGcds7mfnKzNycmZuBTwL/so3/4ga4C/iFiHh+RLwY+Fl6e5Nt9x16//dCRLwc2AJ8q9YVDVH9/cMfAw9l5oeXeNndwNVVi2cH8H8z8+TMFjmGUc4rIjYBdwJvb8v/ZY5yXl3Kjc5d6UfEHcBOYF1EzAMfBF4AkJmt3Y8bdl6Z+VBE/DnwdeAMcHNmLltbbYIR/n39O+C2iDhCb0vkvZnZ9E89/Dng7cCRiDhcHXsfsAmeOa/P0GvwnACeoPd/NE03ynl9APgHwB9VV8VPZ/M/sGyU8+oMJ3IlqSDFbe9IUskMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCvL/AYn2ZMcc2jdiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1892/2000=0.95\n"
     ]
    }
   ],
   "source": [
    "confidence_interval_mean_t(2, 0.5, 12, 2000, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "stone-captain",
   "metadata": {},
   "source": [
    "Our confidence converges effectively to 95%. Once again, we can feel confident about our confidence."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "stainless-sellers",
   "metadata": {},
   "source": [
    "# Conditions for valid t intervals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "advised-installation",
   "metadata": {},
   "source": [
    "We already saw the conditions for valid intervals for a proportion. In the case of t intervals, the same rules apply. The difference is on how to validate if our distribution can be considered normal. :\n",
    "\n",
    "* The sample has to be random\n",
    "* The sampling distribution of the sample mean can be approximated by a normal distribution. There are three ways to achieve it: sample size is bigger than 30 (central limit theorem applies), original distribution is normal or the original distribution is symmetric.\n",
    "* The samples are required to be independent. The rule of thumb is that if you are sampling without replacement your sample size should be less than 10% the population size\n",
    "\n",
    "We leave the experiments to prove that the bullets presented as an exercise for this week."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "innovative-infection",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "You will get the solutions on next week's article."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "authorized-pattern",
   "metadata": {},
   "source": [
    "1. Change the function `confidence_interval_mean_t` and build 5 different experiments to calculate a 95% interval for a sample mean, 3 of them where the conditions for a t interval are not met and 2 where they are. For the 3 where the conditions are not met define the following scenarios: sampling not random, original distribution not approximately normal and independence not met. For the 2 cases where conditions are met define: one case where the original distribution is normal and another where the original distribution is skewed but the sampling mean is normally distributed. For these last two cases, the true population mean should be contained in the calculated confidence intervals 95% of the time. \n",
    "\n",
    "*Hint: you might find useful to use the function `skewnorm` from `scipy`. Below, you have a modified version of a normal distribution skewed by a skewness parameter, mean and standard deviation.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "outside-transaction",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "target mean=2.0000 actual mean=2.0000\n",
      "target stdev=0.5000 actual stdev=0.5000\n"
     ]
    }
   ],
   "source": [
    "# code adapted from https://stackoverflow.com/questions/49367436/scipy-skewnorm-mean-not-matching-theory\n",
    "\n",
    "skew = 4.0\n",
    "mean = 2\n",
    "stdev = 0.5\n",
    "\n",
    "delta = skew / math.sqrt(1. + math.pow(skew, 2.))\n",
    "adjStdev = math.sqrt(math.pow(stdev, 2.) / (1. - 2. * math.pow(delta, 2.) / math.pi))\n",
    "adjMean = mean - adjStdev * math.sqrt(2. / math.pi) * delta\n",
    "\n",
    "print('target mean={:.4f} actual mean={:.4f}'.format(mean, float(skewnorm.stats(skew, loc=adjMean, scale=adjStdev, moments='mvsk')[0])))\n",
    "print('target stdev={:.4f} actual stdev={:.4f}'.format(stdev, math.sqrt(float(skewnorm.stats(skew, loc=adjMean, scale=adjStdev, moments='mvsk')[1]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "plain-revolution",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPJUlEQVR4nO3dbYylZX3H8e9PoD5Rg2QHXFlwaLJtiiatZEKxJIYUmxIhLi/ErIl2Y0k2NrZq20QXm5T0BcmaNkb7nI1Q10jRjdqy8aFK1xrTF6IDUgVWK9EtrGzZ8Vlro13898XcNONwZufMuc/MOefa7yeZnPvhOnP/z33Db6+57oeTqkKS1JanTboASdL4Ge6S1CDDXZIaZLhLUoMMd0lq0NmTLgBg27ZtNT8/P+kyJGmm3Hvvvd+sqrlB66Yi3Ofn51lcXJx0GZI0U5L851rrHJaRpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGTcUdqpqc+X0fHbj82P7rtrgSSeNkz12SGrRuuCe5PcnJJA+sWPZnSb6c5ItJ/jHJeSvW3Zzk4SRfSfJbm1S3JOk0hhmWeQ/wV8B7Vyy7G7i5qk4leTtwM/DWJJcBu4EXAs8H/iXJL1bVE+Mt+8y21lAKOJwiadm6Pfeq+gzw7VXLPllVp7rZzwI7uuldwPur6sdV9XXgYeCKMdYrSRrCOMbcfwf4eDd9EfDoinXHu2VPkWRvksUki0tLS2MoQ5L0pF7hnuSPgVPAHU8uGtCsBr23qg5U1UJVLczNDXzWvCRpRCNfCplkD3A9cE1VPRngx4GLVzTbATw2enmSpFGM1HNPci3wVuAVVfWjFasOA7uTPD3JpcBO4HP9y5QkbcS6PfckdwJXA9uSHAduYfnqmKcDdycB+GxVvb6qHkxyCHiI5eGaN3iljCRtvXXDvapePWDxbadpfytwa5+iJEn9eIeqJDXIcJekBvngsMb4IDBJYM9dkppkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIR/5OsbUe3ytJ67HnLkkNMtwlqUGGuyQ1yHCXpAZ5QlUb4ne0SrPBnrskNWjdnnuS24HrgZNV9aJu2fnAB4B54Bjwqqr6TrfuZuAm4AngjVX1iU2pXBviZZXSmWWYnvt7gGtXLdsHHKmqncCRbp4klwG7gRd27/mbJGeNrVpJ0lDWDfeq+gzw7VWLdwEHu+mDwA0rlr+/qn5cVV8HHgauGE+pkqRhjTrmfmFVnQDoXi/oll8EPLqi3fFu2VMk2ZtkMcni0tLSiGVIkgYZ9wnVDFhWgxpW1YGqWqiqhbm5uTGXIUlntlHD/fEk2wG615Pd8uPAxSva7QAeG708SdIoRg33w8CebnoPcNeK5buTPD3JpcBO4HP9SpQkbdQwl0LeCVwNbEtyHLgF2A8cSnIT8AhwI0BVPZjkEPAQcAp4Q1U9sUm1S5LWsG64V9Wr11h1zRrtbwVu7VOUJKkf71CVpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDfLLOjSQjwiWZps9d0lqkD33KdBCL9mv35Omiz13SWqQ4S5JDTLcJalBhrskNcgTqpoIT8BKm8ueuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KBe4Z7kD5I8mOSBJHcmeUaS85PcneSr3etzx1WsJGk4I4d7kouANwILVfUi4CxgN7APOFJVO4Ej3bwkaQv1HZY5G3hmkrOBZwGPAbuAg936g8ANPbchSdqgkcO9qr4B/DnwCHAC+F5VfRK4sKpOdG1OABcMen+SvUkWkywuLS2NWoYkaYA+wzLPZbmXfinwfODZSV4z7Pur6kBVLVTVwtzc3KhlSJIG6DMs8zLg61W1VFX/C3wY+HXg8STbAbrXk/3LlCRtRJ9wfwS4MsmzkgS4BjgKHAb2dG32AHf1K1GStFEjf81eVd2T5IPAfcAp4AvAAeBc4FCSm1j+B+DGcRQqSRper+9QrapbgFtWLf4xy714SdKEeIeqJDXIcJekBhnuktSgXmPu0nrm93100iVIZyR77pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDvUNVUWeuO1mP7rxtLe+lMYc9dkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGeZ27ZoLf6CRtjD13SWqQPXc16XQ9fe9e1ZnAnrskNahXuCc5L8kHk3w5ydEkL0lyfpK7k3y1e33uuIqVJA2n77DMu4B/rqpXJvk54FnA24AjVbU/yT5gH/DWntuRxsaHjelMMHLPPclzgJcCtwFU1U+q6rvALuBg1+wgcEO/EiVJG9VnWOYXgCXg75N8Icm7kzwbuLCqTgB0rxcMenOSvUkWkywuLS31KEOStFqfcD8buBz426p6MfDfLA/BDKWqDlTVQlUtzM3N9ShDkrRan3A/Dhyvqnu6+Q+yHPaPJ9kO0L2e7FeiJGmjRg73qvov4NEkv9QtugZ4CDgM7OmW7QHu6lWhJGnD+l4t8/vAHd2VMl8DXsfyPxiHktwEPALc2HMbkqQN6hXuVXU/sDBg1TV9fq8kqR/vUJWkBvlsmS3kkw0lbRV77pLUIMNdkhpkuEtSgwx3SWqQJ1Q3gSdOJU2aPXdJapDhLkkNMtwlqUGGuyQ1yHCXpAZ5tYzU8btV1RJ77pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBXgoprcNLJDWL7LlLUoMMd0lqkOEuSQ0y3CWpQb1PqCY5C1gEvlFV1yc5H/gAMA8cA15VVd/pu51p5DcuSZpW4+i5vwk4umJ+H3CkqnYCR7p5SdIW6hXuSXYA1wHvXrF4F3Cwmz4I3NBnG5Kkjevbc38n8BbgpyuWXVhVJwC61wsGvTHJ3iSLSRaXlpZ6liFJWmnkcE9yPXCyqu4d5f1VdaCqFqpqYW5ubtQyJEkD9DmhehXwiiQvB54BPCfJ+4DHk2yvqhNJtgMnx1GoJGl4I/fcq+rmqtpRVfPAbuBTVfUa4DCwp2u2B7ird5WSpA3ZjGfL7AcOJbkJeAS4cRO2IU0tn0WjaTCWcK+qTwOf7qa/BVwzjt8rSRqNd6hKUoN85K80Iu9Q1jSz5y5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3yDtUheCeipFljz12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIG9ikrbIWjfDHdt/3RZXojOB4S5NmKGvzeCwjCQ1yHCXpAYZ7pLUoJHH3JNcDLwXeB7wU+BAVb0ryfnAB4B54Bjwqqr6Tv9SJYFj9BpOn577KeCPquqXgSuBNyS5DNgHHKmqncCRbl6StIVGDveqOlFV93XTPwCOAhcBu4CDXbODwA09a5QkbdBYLoVMMg+8GLgHuLCqTsDyPwBJLljjPXuBvQCXXHLJOMrozS/l0DTxv0f10fuEapJzgQ8Bb66q7w/7vqo6UFULVbUwNzfXtwxJ0gq9wj3JOSwH+x1V9eFu8eNJtnfrtwMn+5UoSdqokcM9SYDbgKNV9Y4Vqw4De7rpPcBdo5cnSRpFnzH3q4DXAl9Kcn+37G3AfuBQkpuAR4Abe1UoSdqwkcO9qv4NyBqrrxn190qS+vMOVUlqkOEuSQ0y3CWpQT7PXWrc6W6G8nk07bLnLkkNMtwlqUEOy0hnMB8f3C577pLUIHvuUiN8iqRWsucuSQ0y3CWpQYa7JDXIcJekBp2RJ1Q98SSpdWdkuEs6Pa9/n30Oy0hSg+y5S9o0/gUwOfbcJalB9twlDc2LEWaH4S5pyzlcs/kclpGkBjXRc7cXIEk/y567JDWoiZ67pDaM64Stf7Ub7pIatNGh2nF+ifi0DBNv2rBMkmuTfCXJw0n2bdZ2JElPtSk99yRnAX8N/CZwHPh8ksNV9dBmbG8tXpMraaVJZsJW9+g3q+d+BfBwVX2tqn4CvB/YtUnbkiStsllj7hcBj66YPw782soGSfYCe7vZHyb5ypi2vQ345ph+16TM+mew/smy/jHK2zf8lg3VP8LvX+kFa63YrHDPgGX1MzNVB4ADY99wslhVC+P+vVtp1j+D9U+W9U/WtNS/WcMyx4GLV8zvAB7bpG1JklbZrHD/PLAzyaVJfg7YDRzepG1JklbZlGGZqjqV5PeATwBnAbdX1YObsa0Bxj7UMwGz/hmsf7Ksf7Kmov5U1fqtJEkzxWfLSFKDDHdJatDMhnuS25OcTPLAGuuT5C+6xx98McnlW13j6QxR/9VJvpfk/u7nT7a6xrUkuTjJvyY5muTBJG8a0Gba9/8wn2Gaj8Ezknwuyb939f/pgDZTewyGrH9q9/+TkpyV5AtJPjJg3WT3f1XN5A/wUuBy4IE11r8c+DjL19xfCdwz6Zo3WP/VwEcmXecatW0HLu+mfx74D+CyGdv/w3yGaT4GAc7tps8B7gGunJVjMGT9U7v/V9T4h8A/DKpz0vt/ZnvuVfUZ4NunabILeG8t+yxwXpLtW1Pd+oaof2pV1Ymquq+b/gFwlOW7klea9v0/zGeYWt1+/WE3e073s/rqiKk9BkPWP9WS7ACuA969RpOJ7v+ZDfchDHoEwsz8z9t5Sfdn68eTvHDSxQySZB54Mcs9r5VmZv+f5jPAFB+DbkjgfuAkcHdVzdQxGKJ+mOL9D7wTeAvw0zXWT3T/txzu6z4CYcrdB7ygqn4F+EvgnyZbzlMlORf4EPDmqvr+6tUD3jJ1+3+dzzDVx6CqnqiqX2X5DvArkrxoVZOpPgZD1D+1+z/J9cDJqrr3dM0GLNuy/d9yuM/0IxCq6vtP/tlaVR8DzkmybcJl/b8k57AcindU1YcHNJn6/b/eZ5j2Y/Ckqvou8Gng2lWrpv4YwNr1T/n+vwp4RZJjLD/19jeSvG9Vm4nu/5bD/TDw290Z6yuB71XViUkXNawkz0uSbvoKlo/VtyZb1bKurtuAo1X1jjWaTfX+H+YzTPkxmEtyXjf9TOBlwJdXNZvaYzBM/dO8/6vq5qraUVXzLD9e5VNV9ZpVzSa6/2f2a/aS3Mny2fRtSY4Dt7B8Uoaq+jvgYyyfrX4Y+BHwuslUOtgQ9b8S+N0kp4D/AXZXdwp+ClwFvBb4UjdmCvA24BKYjf3PcJ9hmo/BduBglr8Y52nAoar6SJLXw0wcg2Hqn+b9P9A07X8fPyBJDWp5WEaSzliGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQ/wE5cWG1rAe2ygAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Original skewed distribution\n",
    "\n",
    "plt.hist(skewnorm.rvs(a = skew, loc=adjMean, scale=adjStdev, size=2000), bins=50);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "surrounded-frontier",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPdklEQVR4nO3dfYylZXnH8e+vILaKCIaBrrx0absRkWglU0q1MTbbRhTj0jYYTNUNYjZNrMWmti72D/5oSNbYmLaxttkAdY1EpYphW2Ir2WqosYCDIAILSgRhy5YdxdfaoKtX/zgPdTrMsGfOc86cmXu+n2RzzvNy9lx3Zvd3rrmfl5OqQpLUlp+ZdgGSpPEz3CWpQYa7JDXIcJekBhnuktSgo6ddAMCJJ55YmzdvnnYZkrSu3H777d+oqpmltq2JcN+8eTNzc3PTLkOS1pUkX19um9MyktQgw12SGmS4S1KDDHdJatARwz3JNUkOJbl7wbr3JrkvyV1JPpnk+AXbLk/yQJL7k7xqQnVLkp7GMJ37B4HzF627CTi7ql4MfAW4HCDJWcDFwIu613wgyVFjq1aSNJQjhntV3Qw8vmjdp6vqcLd4C3Bq93wb8NGqeqKqHgQeAM4dY72SpCGMY879LcCnuuenAI8s2HagW/cUSXYkmUsyNz8/P4YyJElP6hXuSf4cOAxc++SqJXZb8obxVbW7qmaranZmZskLrCRJIxr5CtUk24HXAlvrp9/4cQA4bcFupwKPjl6e9PQ277xxyfUP7bpglSuR1paROvck5wPvAl5XVT9YsGkvcHGSZyY5A9gC3Na/TEnSShyxc0/yEeCVwIlJDgBXMDg75pnATUkAbqmqP6iqe5JcB9zLYLrmbVX140kVL0la2hHDvaresMTqq59m/yuBK/sUJUnqxytUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatDIX7MnrVd+NZ82Ajt3SWqQ4S5JDTLcJalBhrskNcgDqloXljsIKmlpdu6S1CA7dzVplE7fUyTVEjt3SWqQnbsmqoVuuIUxaOMx3LWmeOBUGo8jTsskuSbJoSR3L1j3vCQ3Jflq93jCgm2XJ3kgyf1JXjWpwiVJyxtmzv2DwPmL1u0E9lXVFmBft0ySs4CLgRd1r/lAkqPGVq0kaShHDPequhl4fNHqbcCe7vke4MIF6z9aVU9U1YPAA8C54ylVkjSsUc+WObmqDgJ0jyd1608BHlmw34Fu3VMk2ZFkLsnc/Pz8iGVIkpYy7lMhs8S6WmrHqtpdVbNVNTszMzPmMiRpYxs13B9LsgmgezzUrT8AnLZgv1OBR0cvT5I0ilHDfS+wvXu+HbhhwfqLkzwzyRnAFuC2fiVKklbqiOe5J/kI8ErgxCQHgCuAXcB1SS4FHgYuAqiqe5JcB9wLHAbeVlU/nlDtkqRlHDHcq+oNy2zausz+VwJX9ilKktSP95aRpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNOnraBWhj2rzzxmmXIDXNzl2SGmS4S1KDek3LJPlj4K1AAV8GLgGeBXwM2Aw8BLy+qr7Vq0qteU6zSGvLyJ17klOAPwJmq+ps4CjgYmAnsK+qtgD7umVJ0irqOy1zNPBzSY5m0LE/CmwD9nTb9wAX9nwPSdIKjRzuVfWfwF8CDwMHge9U1aeBk6vqYLfPQeCkcRQqSRreyHPuSU5g0KWfAXwb+Mckb1zB63cAOwBOP/30UcuQpma54wwP7bpglSuRnqrPtMxvAQ9W1XxV/Qi4HngZ8FiSTQDd46GlXlxVu6tqtqpmZ2ZmepQhSVqsT7g/DJyX5FlJAmwF9gN7ge3dPtuBG/qVKElaqZGnZarq1iQfB74IHAbuAHYDxwLXJbmUwQfAReMoVJI0vF7nuVfVFcAVi1Y/waCLlyRNifeW0ZI8WCitb95+QJIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapD3ltng/GJrqU127pLUIMNdkhpkuEtSgwx3SWqQB1SlMfOLTrQW2LlLUoPs3KVVYkev1WTnLkkNMtwlqUGGuyQ1yHCXpAZ5QFUr4r1opPXBzl2SGtSrc09yPHAVcDZQwFuA+4GPAZuBh4DXV9W3+ryPtBF56qT66Nu5/zXwL1V1JvASYD+wE9hXVVuAfd2yJGkVjRzuSY4DXgFcDVBVP6yqbwPbgD3dbnuAC/uVKElaqT6d+y8C88A/JLkjyVVJng2cXFUHAbrHk5Z6cZIdSeaSzM3Pz/coQ5K0WJ9wPxo4B/i7qnop8N+sYAqmqnZX1WxVzc7MzPQoQ5K0WJ9wPwAcqKpbu+WPMwj7x5JsAugeD/UrUZK0UiOHe1X9F/BIkhd0q7YC9wJ7ge3duu3ADb0qlCStWN+LmN4OXJvkGOBrwCUMPjCuS3Ip8DBwUc/3kCStUK9wr6o7gdklNm3t8/dKkvrxClVJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhrkNzFtEH6DkrSxGO7SlPnBq0lwWkaSGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIa1DvckxyV5I4k/9wtPy/JTUm+2j2e0L9MSdJKjKNzvwzYv2B5J7CvqrYA+7plSdIq6hXuSU4FLgCuWrB6G7Cne74HuLDPe0iSVq5v5/5XwJ8BP1mw7uSqOgjQPZ601AuT7Egyl2Rufn6+ZxmSpIVGDvckrwUOVdXto7y+qnZX1WxVzc7MzIxahiRpCUf3eO3LgdcleQ3ws8BxST4MPJZkU1UdTLIJODSOQiU9vc07b1xy/UO7LljlSrQWjNy5V9XlVXVqVW0GLgb+rareCOwFtne7bQdu6F2lJGlF+nTuy9kFXJfkUuBh4KIJvIe0YS3XoY+yv119u8YS7lX1WeCz3fNvAlvH8fdKkkbjFaqS1CDDXZIaZLhLUoMmcUBVE7bSA2qSNh47d0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapBXqK5hXokqaVR27pLUIMNdkhpkuEtSg5xzX0V+gbGk1WLnLkkNsnNfAzwrRtK42blLUoMMd0lqkNMy0gbmQf522blLUoMMd0lqkOEuSQ0aOdyTnJbkM0n2J7knyWXd+ucluSnJV7vHE8ZXriRpGH0698PAn1TVC4HzgLclOQvYCeyrqi3Avm5ZkrSKRg73qjpYVV/snn8P2A+cAmwD9nS77QEu7FmjJGmFxjLnnmQz8FLgVuDkqjoIgw8A4KRlXrMjyVySufn5+XGUIUnq9A73JMcCnwDeUVXfHfZ1VbW7qmaranZmZqZvGZKkBXqFe5JnMAj2a6vq+m71Y0k2dds3AYf6lShJWqk+Z8sEuBrYX1XvW7BpL7C9e74duGH08iRJo+hz+4GXA28Cvpzkzm7du4FdwHVJLgUeBi7qVaEkacVGDveq+hyQZTZvHfXvlST15xWqktQgw12SGmS4S1KDvJ+7pKfwPu/rn+EuaWiG/vrhtIwkNcjOfQKW626kVtnRrz127pLUIMNdkhpkuEtSgwx3SWqQB1QlTYwHWqfHzl2SGmS4S1KDDHdJapDhLkkNMtwlqUGeLSNp1XkWzeTZuUtSgwx3SWqQ0zI9ePdHSWuVnbskNWhDdu4r7bg9yCNNlwdgV87OXZIa1ETnPulPdefWpdXh/7XxaSLcl+M/FGljchpngtMySc5Pcn+SB5LsnNT7SJKeaiKde5KjgL8Ffhs4AHwhyd6quncS7ydpY1rpb+fT7OhX+70n1bmfCzxQVV+rqh8CHwW2Tei9JEmLTGrO/RTgkQXLB4BfW7hDkh3Ajm7x+0nun1AtfZ0IfGPaRUxAq+OCdsfmuCYk75nYX33EsfV8719YbsOkwj1LrKv/t1C1G9g9ofcfmyRzVTU77TrGrdVxQbtjc1zrzzTHNqlpmQPAaQuWTwUendB7SZIWmVS4fwHYkuSMJMcAFwN7J/RekqRFJjItU1WHk/wh8K/AUcA1VXXPJN5rFaz5qaMRtTouaHdsjmv9mdrYUlVH3kuStK54bxlJapDhLkkNMtyBJNckOZTk7mW2PzfJPyX5UpJ7klyy2jWOYohxnZDkk0nuSnJbkrNXu8ZRJDktyWeS7O9+HpctsU+S/E13+4u7kpwzjVpXasixnZnkP5I8keSd06hzpYYc1+93P6u7knw+yUumUetKDTm2bd247kwyl+Q3Jl5YVW34P8ArgHOAu5fZ/m7gPd3zGeBx4Jhp1z2Gcb0XuKJ7fiawb9o1DzmuTcA53fPnAF8Bzlq0z2uATzG45uI84NZp1z3GsZ0E/CpwJfDOadc8xnG9DDihe/7qxn5mx/LTY5wvBu6bdF127kBV3cwgsJfdBXhOkjD4IT0OHF6N2voYYlxnAfu6fe8DNic5eTVq66OqDlbVF7vn3wP2M7gqeqFtwIdq4Bbg+CSbVrnUFRtmbFV1qKq+APxoCiWOZMhxfb6qvtUt3sLg+pg1b8ixfb+6ZAeezaKLOifBcB/O+4EXMrgQ68vAZVX1k+mWNBZfAn4XIMm5DC5lXhf/oZ6UZDPwUuDWRZuWugXG4g+ANe1pxrauDTmuSxn85rWuPN3YkvxOkvuAG4G3TLoWw304rwLuBJ4P/Arw/iTHTbOgMdkFnJDkTuDtwB2sg99InpTkWOATwDuq6ruLNy/xknVz3u8RxrZuDTOuJL/JINzftZq19XWksVXVJ6vqTOBC4C8mXU/TX9YxRpcAu7pfqx5I8iCDOerbpltWP90/wEtgcAASeLD7s+YleQaD/0jXVtX1S+yybm+BMcTY1qVhxpXkxcBVwKur6purWV8fK/mZVdXNSX4pyYlVNbEbptm5D+dhYCtANyf9AuBrU61oDJIc390eAuCtwM3roUvsPoiuBvZX1fuW2W0v8OburJnzgO9U1cFVK3JEQ45t3RlmXElOB64H3lRVX1nN+voYcmy/3O1Hd+bWMcBEP7y8QhVI8hHglQxuz/kYcAXwDICq+vskzwc+yOCoeBh08R+eSrErMMS4fh34EPBj4F7g0gUHtNas7jSyf2dw/OPJYx/vBk6H/xtbGBwrOR/4AXBJVc1NodwVGXJsPw/MAcd1+3yfwdkZa/aDechxXQX8HvD1bvvhWgd3ixxybO8C3szgIPj/AH9aVZ+baF2GuyS1x2kZSWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIa9L++yBaFli2CbAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Approximately normal distribution of the sample mean because sample \n",
    "# size is bigger than 30 (CTL applies)\n",
    "\n",
    "plt.hist(np.mean([skewnorm.rvs(a = skew, loc=adjMean, scale=adjStdev, size=35) for _ in range(2000)], axis=1), bins=50);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "floating-malawi",
   "metadata": {},
   "source": [
    "### Answers from last week"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sized-singapore",
   "metadata": {},
   "source": [
    "1. Rick takes an SRS of 75 citizens in a distant planet to see what proportion of citizens sampled are satisfied with their standard of living. Suppose that 60% of the 1,000,000,000 citizens that live in the planet are satisfied with their standard of living.\n",
    "What are the mean and standard deviation of the sampling distribution of the proportion of citizens who are satisfied with their standard of living?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "extra-stroke",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6\n",
      "0.0565685424949238\n"
     ]
    }
   ],
   "source": [
    "mu_hat_p = 0.6\n",
    "print(mu_hat_p)\n",
    "\n",
    "sigma_p_hat = np.sqrt(0.6*(1-0.6)/75)\n",
    "print(sigma_p_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "impaired-disposition",
   "metadata": {},
   "source": [
    "2. A certain planet with over 1,000,000 households has a mean household income of \\\\$1,000,000 with a standard deviation of \\\\$150,000. Rick plans to take random samples of 700 households and calculate the sample mean income. Calculate the mean and standard deviation of the sampling distribution of $\\bar x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "yellow-hacker",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "μ_x_bar = 1000000\n",
      "σ_x_bar = 5669.467095138409\n"
     ]
    }
   ],
   "source": [
    "μ = 1000000\n",
    "σ = 150000\n",
    "n = 700\n",
    "\n",
    "print(f'μ_x_bar = {μ}')\n",
    "print(f'σ_x_bar = {σ/n**(1/2)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "basic-composer",
   "metadata": {},
   "source": [
    "3. Rick is performing quality control tests on different portal guns since there is some variability in the manufactoring process. A certain gun has a target thickness of 5mm. The distribution of thicknesses is skewed to the right with a mean of 5mm and a standard deviation of 1mm. A quality control check on this part involves taking a random sample of 35 points and calculating the mean thickness of those points. What is the shape of the sampling distribution of the sample mean thickness? What is the probability that the mean thickness in the sample is within 0.2mm of the target value?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "diverse-dependence",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Since n = 35 >= 30, the central limit theorem applies.\n",
    "# Even though the population is skewed to the right, the sample means\n",
    "# are normally distributed due to the sample size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "normal-clarity",
   "metadata": {},
   "outputs": [],
   "source": [
    "μ = 5\n",
    "σ = 1\n",
    "n = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "strong-romania",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9544997361036418"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm.cdf(5.2, μ, σ/n**(1/2)) - norm.cdf(4.8, μ, σ/n**(1/2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}