k-anonymity.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# k - Anonymity Analysis\r\n",
    "\r\n",
    "## a) Randomly generate a dataset (dataframe) with eight columns and 50,000 rows. Each column should be a categorical variable (of arbitrary name) with three levels (of arbitrary names) in approximately equal proportions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(50000, 8)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>C</th>\n",
       "      <th>I</th>\n",
       "      <th>J</th>\n",
       "      <th>K</th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>a3</td>\n",
       "      <td>b2</td>\n",
       "      <td>c3</td>\n",
       "      <td>i3</td>\n",
       "      <td>j1</td>\n",
       "      <td>k3</td>\n",
       "      <td>x1</td>\n",
       "      <td>y3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>a3</td>\n",
       "      <td>b1</td>\n",
       "      <td>c2</td>\n",
       "      <td>i3</td>\n",
       "      <td>j1</td>\n",
       "      <td>k3</td>\n",
       "      <td>x2</td>\n",
       "      <td>y1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>a3</td>\n",
       "      <td>b2</td>\n",
       "      <td>c3</td>\n",
       "      <td>i1</td>\n",
       "      <td>j3</td>\n",
       "      <td>k1</td>\n",
       "      <td>x3</td>\n",
       "      <td>y2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>a3</td>\n",
       "      <td>b2</td>\n",
       "      <td>c2</td>\n",
       "      <td>i3</td>\n",
       "      <td>j1</td>\n",
       "      <td>k2</td>\n",
       "      <td>x3</td>\n",
       "      <td>y3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>a1</td>\n",
       "      <td>b1</td>\n",
       "      <td>c3</td>\n",
       "      <td>i3</td>\n",
       "      <td>j3</td>\n",
       "      <td>k1</td>\n",
       "      <td>x2</td>\n",
       "      <td>y2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    A   B   C   I   J   K   X   Y\n",
       "0  a3  b2  c3  i3  j1  k3  x1  y3\n",
       "1  a3  b1  c2  i3  j1  k3  x2  y1\n",
       "2  a3  b2  c3  i1  j3  k1  x3  y2\n",
       "3  a3  b2  c2  i3  j1  k2  x3  y3\n",
       "4  a1  b1  c3  i3  j3  k1  x2  y2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\r\n",
    "import numpy as np\r\n",
    "\r\n",
    "dictionary = dict({\r\n",
    "    \"A\":[\"a1\",\"a2\",\"a3\"],\r\n",
    "    \"B\":[\"b1\",\"b2\",\"b3\"],\r\n",
    "    \"C\":[\"c1\",\"c2\",\"c3\"],\r\n",
    "    \"I\":[\"i1\",\"i2\",\"i3\"],\r\n",
    "    \"J\":[\"j1\",\"j2\",\"j3\"],\r\n",
    "    \"K\":[\"k1\",\"k2\",\"k3\"],\r\n",
    "    \"X\":[\"x1\",\"x2\",\"x3\"],\r\n",
    "    \"Y\":[\"y1\",\"y2\",\"y3\"]\r\n",
    "})\r\n",
    "\r\n",
    "np.random.seed(92)\r\n",
    "\r\n",
    "def categorical_df_from_dict(d: dict, size: int) -> pd.DataFrame:\r\n",
    "    \"\"\" Creates a dataframe with columns = d.keys, where each column is populated with values randomly sampled from d.values.\"\"\"\r\n",
    "\r\n",
    "    df = pd.DataFrame()\r\n",
    "\r\n",
    "    for key, val in d.items():\r\n",
    "\r\n",
    "        df[key] = np.random.choice(val, size)\r\n",
    "\r\n",
    "    return df\r\n",
    "\r\n",
    "df = categorical_df_from_dict(dictionary, size=50000)\r\n",
    "\r\n",
    "print(df.shape)\r\n",
    "display(df.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## b) Verify that the proportions of each value are similar for each of the eight columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['a1' 'a2' 'a3']\n",
      "[0.33406 0.33004 0.3359 ]\n",
      "['b1' 'b2' 'b3']\n",
      "[0.33304 0.33118 0.33578]\n",
      "['c1' 'c2' 'c3']\n",
      "[0.3289 0.3346 0.3365]\n",
      "['i1' 'i2' 'i3']\n",
      "[0.33328 0.33416 0.33256]\n",
      "['j1' 'j2' 'j3']\n",
      "[0.3339  0.33152 0.33458]\n",
      "['k1' 'k2' 'k3']\n",
      "[0.3345  0.33182 0.33368]\n",
      "['x1' 'x2' 'x3']\n",
      "[0.33586 0.33352 0.33062]\n",
      "['y1' 'y2' 'y3']\n",
      "[0.33464 0.33314 0.33222]\n"
     ]
    }
   ],
   "source": [
    "def verify_proportions(df: pd.DataFrame) -> None:\r\n",
    "    \"\"\" Prints the proportions of values in all columns in the dataframe df. \"\"\"\r\n",
    "    for col in df:\r\n",
    "        vals, counts = np.unique(df[col], return_counts=True)\r\n",
    "\r\n",
    "        print(vals)\r\n",
    "        print(counts/sum(counts))\r\n",
    "\r\n",
    "verify_proportions(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## c) How many unique rows (i.e., permutations of category levels) are possible?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total permutations: 6561\n"
     ]
    }
   ],
   "source": [
    "n = 3 # unique values in each column\r\n",
    "r = 8 # columns\r\n",
    "\r\n",
    "print(f\"Total permutations: {n**r}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['a1', 'b1', 'c1', 'i1', 'j1', 'k1', 'x1', 'y1']\n",
      "['a1', 'b1', 'c1', 'i1', 'j1', 'k1', 'x1', 'y2']\n",
      "['a1', 'b1', 'c1', 'i1', 'j1', 'k1', 'x1', 'y3']\n",
      "['a1', 'b1', 'c1', 'i1', 'j1', 'k1', 'x2', 'y1']\n",
      "['a1', 'b1', 'c1', 'i1', 'j1', 'k1', 'x2', 'y2']\n",
      "['a1', 'b1', 'c1', 'i1', 'j1', 'k1', 'x2', 'y3']\n",
      "['a1', 'b1', 'c1', 'i1', 'j1', 'k1', 'x3', 'y1']\n",
      "['a1', 'b1', 'c1', 'i1', 'j1', 'k1', 'x3', 'y2']\n",
      "['a1', 'b1', 'c1', 'i1', 'j1', 'k1', 'x3', 'y3']\n"
     ]
    }
   ],
   "source": [
    "import itertools\r\n",
    "\r\n",
    "def print_permutations(n: int) -> None:\r\n",
    "    \"\"\"\r\n",
    "    Prints the first n permutations from itertools.product and adds the letters so that the variables are labelled the same as defined in 1. a)\r\n",
    "    \"\"\"\r\n",
    "    for num, perm in enumerate(itertools.product('123', repeat=8)):\r\n",
    "        perm = list(perm)\r\n",
    "        letters = [\"a\",\"b\",\"c\",\"i\",\"j\",\"k\",\"x\",\"y\"]\r\n",
    "        for x, y in enumerate(perm):\r\n",
    "            perm[x] = letters[x]+y\r\n",
    "        print(perm)\r\n",
    "        if num >= n-1:\r\n",
    "            break\r\n",
    "\r\n",
    "print_permutations(n=9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## d) Produce a table and appropriate graph which show the frequencies (numbers of groups) by permutation group sizes up to group size of 12. That is, how many groups are unique combinations (group size = 1), how many groups are made up of a pair of matching combinations (group size = 2), how many groups are made up three the same, etc?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "size\n",
       "1      31\n",
       "2      93\n",
       "3     235\n",
       "4     454\n",
       "5     651\n",
       "6     871\n",
       "7     977\n",
       "8     949\n",
       "9     751\n",
       "10    591\n",
       "11    375\n",
       "12    272\n",
       "Name: A, dtype: int64"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "unique_perms = (\r\n",
    "    df\r\n",
    "    .groupby(df.columns.tolist(), as_index=False)\r\n",
    "    .size()\r\n",
    "    .query(\"size <= 12\")\r\n",
    ")\r\n",
    "\r\n",
    "display(unique_perms.groupby('size').count().iloc[:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAFDCAYAAACKmSIcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAc4UlEQVR4nO3de7RdZX3u8e9DotwUBQmIgAbbiAVrBSP11uIdNSpYpMWKRkvFttSjVsdp0LZYLafpaGvVU1HBS6OoNOIFLK2VYm1ttWIUvAByiBIhhUK0VdRDweDv/DHn9qxsdpKFZK357p3vZ4yMPee75lrrYZPs/ax33lJVSJIkqT27DB1AkiRJc7OoSZIkNcqiJkmS1CiLmiRJUqMsapIkSY2yqEmSJDVq8dABJmXfffetpUuXDh1DkiRpu77whS98q6qWzB5fsEVt6dKlrFu3bugYkiRJ25Xkm3ONu+tTkiSpURMrakneleSmJF8dGdsnyUVJru6/7j3y2GlJ1ie5KskxI+MPT/KV/rE3J8mkMkuSJLVkkjNqfwU8ddbYKuDiqloGXNyvk+Qw4ETg8P45ZyZZ1D/nrcApwLL+z+zXlCRJWpAmVtSq6p+B/5w1fCywpl9eAxw3Mn5uVd1aVdcA64GjkhwA7FVVn63upqTvGXmOJEnSgjbtY9T2r6obAPqv+/XjBwLXjWy3sR87sF+ePS5JkrTgtXIywVzHndU2xud+keSUJOuSrNu0adMOCydJkjSEaRe1G/vdmfRfb+rHNwIHj2x3EHB9P37QHONzqqqzqmp5VS1fsuQOlyKRJEmaV6Zd1C4AVvbLK4HzR8ZPTLJrkkPoThq4pN89+r0kj+zP9nzByHMkSZIWtIld8DbJB4DHAfsm2QicDqwG1iY5GbgWOAGgqi5Psha4AtgMnFpVt/cv9Zt0Z5DuDvxd/0eSJGnBS3cy5cKzfPny8s4EkiRpPkjyhapaPnu8lZMJJEmSNMuCvdenJE3D0lUXDh0BgA2rVwwdQdIEOKMmSZLUKIuaJElSoyxqkiRJjbKoSZIkNcqiJkmS1CiLmiRJUqO8PIck7SS8lIg0/zijJkmS1CiLmiRJUqMsapIkSY2yqEmSJDXKoiZJktQoi5okSVKjLGqSJEmNsqhJkiQ1yqImSZLUKIuaJElSoyxqkiRJjbKoSZIkNcqiJkmS1CiLmiRJUqMsapIkSY2yqEmSJDXKoiZJktQoi5okSVKjLGqSJEmNsqhJkiQ1avHQASRptqWrLhw6AgAbVq8YOoKknZwzapIkSY2yqEmSJDXKoiZJktQoi5okSVKjLGqSJEmNsqhJkiQ1yqImSZLUKIuaJElSoyxqkiRJjbKoSZIkNcqiJkmS1CiLmiRJUqMsapIkSY2yqEmSJDVqkKKW5BVJLk/y1SQfSLJbkn2SXJTk6v7r3iPbn5ZkfZKrkhwzRGZJkqRpm3pRS3Ig8D+A5VX1EGARcCKwCri4qpYBF/frJDmsf/xw4KnAmUkWTTu3JEnStA2163MxsHuSxcAewPXAscCa/vE1wHH98rHAuVV1a1VdA6wHjppuXEmSpOmbelGrqn8H/gy4FrgB+G5VfQLYv6pu6Le5Adivf8qBwHUjL7GxH5MkSVrQhtj1uTfdLNkhwP2APZOctK2nzDFWW3ntU5KsS7Ju06ZNdz2sJEnSgIbY9fkk4Jqq2lRVPwQ+DDwauDHJAQD915v67TcCB488/yC6XaV3UFVnVdXyqlq+ZMmSif0HSJIkTcMQRe1a4JFJ9kgS4InAlcAFwMp+m5XA+f3yBcCJSXZNcgiwDLhkypklSZKmbvG037CqPpfkPOCLwGbgUuAs4B7A2iQn05W5E/rtL0+yFrii3/7Uqrp92rklSZKmbepFDaCqTgdOnzV8K93s2lzbnwGcMelckiRJLfHOBJIkSY2yqEmSJDXKoiZJktQoi5okSVKjLGqSJEmNsqhJkiQ1yqImSZLUKIuaJElSowa54K0kSduydNWFQ0cAYMPqFUNH0E7OGTVJkqRGWdQkSZIaZVGTJElqlEVNkiSpURY1SZKkRlnUJEmSGmVRkyRJapRFTZIkqVEWNUmSpEZZ1CRJkhplUZMkSWqURU2SJKlRFjVJkqRGWdQkSZIaZVGTJElqlEVNkiSpURY1SZKkRlnUJEmSGmVRkyRJapRFTZIkqVEWNUmSpEZZ1CRJkhplUZMkSWqURU2SJKlRFjVJkqRGLR46gKTpWLrqwqEjALBh9YqhI0jSvOGMmiRJUqMsapIkSY2yqEmSJDXKoiZJktQoi5okSVKjLGqSJEmNsqhJkiQ1yqImSZLUKIuaJElSowYpaknuneS8JF9LcmWSRyXZJ8lFSa7uv+49sv1pSdYnuSrJMUNkliRJmrahZtTeBHy8qh4M/BxwJbAKuLiqlgEX9+skOQw4ETgceCpwZpJFg6SWJEmaoqkXtSR7Ab8IvBOgqm6rqu8AxwJr+s3WAMf1y8cC51bVrVV1DbAeOGqamSVJkoYwVlFL8pAd+J4PBDYB705yaZJ3JNkT2L+qbgDov+7Xb38gcN3I8zf2Y5IkSQvauDNqb0tySZLfSnLvu/iei4EjgbdW1RHAD+h3c25F5hirOTdMTkmyLsm6TZs23cWYkiRJwxqrqFXVY4HnAQcD65K8P8mTf8L33AhsrKrP9evn0RW3G5McANB/vWlk+4NHnn8QcP1Wcp5VVcuravmSJUt+wniSJEltWDzuhlV1dZLfA9YBbwaOSBLg1VX14TvxOv+R5Lokh1bVVcATgSv6PyuB1f3X8/unXAC8P8kbgPsBy4BLxn0/SZImaemqC4eOAMCG1SuGjqAJGKuoJXko8CJgBXAR8Myq+mKS+wGfBcYuar2XAu9LcnfgG/1r7wKsTXIycC1wAkBVXZ5kLV2R2wycWlW338n3kyRJmnfGnVH7S+BsutmzW2YGq+r6fpbtTqmqy4Dlczz0xK1sfwZwxp19H0mSpPls3KL2dOCWmZmsJLsAu1XV/62q904snSRJ0k5s3LM+/wHYfWR9j35MkiRJEzJuUdutqr4/s9Iv7zGZSJIkSYLxi9oPkhw5s5Lk4cAt29hekiRJd9G4x6i9HPhgkpnrlx0A/MpEEkmSJAkYs6hV1eeTPBg4lO5OAV+rqh9ONJkkSdJObuwL3gKPAJb2zzkiCVX1nomkkiRJ0tgXvH0v8FPAZcDMxWYLsKhJkiRNyLgzasuBw6pqzpuhS5Ikaccb96zPrwL3nWQQSZIkbWncGbV9gSuSXALcOjNYVc+aSCpJkiSNXdReO8kQkiRJuqNxL8/xT0keACyrqn9IsgewaLLRJEmSdm5jHaOW5MXAecDb+6EDgY9OKJMkSZIY/2SCU4HHADcDVNXVwH6TCiVJkqTxi9qtVXXbzEqSxXTXUZMkSdKEjFvU/inJq4HdkzwZ+CDwscnFkiRJ0rhFbRWwCfgK8BLgb4Hfm1QoSZIkjX/W54+As/s/kiRJmoJx7/V5DXMck1ZVD9zhiSRJkgTcuXt9ztgNOAHYZ8fHkSRJ0oyxjlGrqm+P/Pn3qnoj8ITJRpMkSdq5jbvr88iR1V3oZtjuOZFEkiRJAsbf9fnnI8ubgQ3AL+/wNJIkSfqxcc/6fPykg0iSJGlL4+76/J1tPV5Vb9gxcSRJkjTjzpz1+Qjggn79mcA/A9dNIpQkSZLGL2r7AkdW1fcAkrwW+GBV/fqkgkmSJO3sxr2F1P2B20bWbwOW7vA0kiRJ+rFxZ9TeC1yS5CN0dyh4NvCeiaWSJEnS2Gd9npHk74Bf6IdeVFWXTi6WJEmSxt31CbAHcHNVvQnYmOSQCWWSJEkSYxa1JKcDvwuc1g/dDThnUqEkSZI0/jFqzwaOAL4IUFXXJ/EWUtrpLV114dARfmzD6hVDR5Ak7WDj7vq8raqK7kQCkuw5uUiSJEmC8Yva2iRvB+6d5MXAPwBnTy6WJEmStrvrM0mAvwYeDNwMHAr8QVVdNOFskiRJO7XtFrWqqiQfraqHA5YzSZKkKRl31+e/JXnERJNIkiRpC+Oe9fl44DeSbAB+AIRusu2hkwomSZK0s9tmUUty/6q6FnjalPJIkiSpt70ZtY8CR1bVN5N8qKqOn0ImSZIksf1j1DKy/MBJBpEkSdKWtlfUaivLkiRJmrDtFbWfS3Jzku8BD+2Xb07yvSQ335U3TrIoyaVJ/qZf3yfJRUmu7r/uPbLtaUnWJ7kqyTF35X0lSZLmi20WtapaVFV7VdU9q2pxvzyzvtddfO+XAVeOrK8CLq6qZcDF/TpJDgNOBA4HngqcmWTRXXxvSZKk5o17HbUdKslBwArgHSPDxwJr+uU1wHEj4+dW1a1VdQ2wHjhqSlElSZIGM0hRA94I/E/gRyNj+1fVDQD91/368QOB60a229iPSZIkLWhTL2pJngHcVFVfGPcpc4zNeWJDklOSrEuybtOmTT9xRkmSpBYMMaP2GOBZ/V0OzgWekOQc4MYkBwD0X2/qt98IHDzy/IOA6+d64ao6q6qWV9XyJUuWTCq/JEnSVEy9qFXVaVV1UFUtpTtJ4JNVdRJwAbCy32wlcH6/fAFwYpJdkxwCLAMumXJsSZKkqRv3Xp/TsBpYm+Rk4FrgBICqujzJWuAKYDNwalXdPlxMSZKk6Ri0qFXVp4BP9cvfBp64le3OAM6YWjBJkqQGDHXWpyRJkrajpV2fkiRpgpauunDoCABsWL1i6AjzhjNqkiRJjbKoSZIkNcqiJkmS1CiLmiRJUqMsapIkSY2yqEmSJDXKoiZJktQoi5okSVKjLGqSJEmNsqhJkiQ1yqImSZLUKIuaJElSoyxqkiRJjbKoSZIkNcqiJkmS1CiLmiRJUqMsapIkSY2yqEmSJDXKoiZJktQoi5okSVKjLGqSJEmNsqhJkiQ1yqImSZLUKIuaJElSoyxqkiRJjbKoSZIkNcqiJkmS1KjFQweQJEmabemqC4eOAMCG1SsGfX9n1CRJkhplUZMkSWqURU2SJKlRFjVJkqRGWdQkSZIaZVGTJElqlEVNkiSpURY1SZKkRnnBWzWnlYscwvAXOpQk7dycUZMkSWqURU2SJKlRFjVJkqRGWdQkSZIaZVGTJElq1NSLWpKDk/xjkiuTXJ7kZf34PkkuSnJ1/3XvkeeclmR9kquSHDPtzJIkSUMYYkZtM/DKqvoZ4JHAqUkOA1YBF1fVMuDifp3+sROBw4GnAmcmWTRAbkmSpKmaelGrqhuq6ov98veAK4EDgWOBNf1ma4Dj+uVjgXOr6taqugZYDxw11dCSJEkDGPQYtSRLgSOAzwH7V9UN0JU5YL9+swOB60aetrEfkyRJWtAGK2pJ7gF8CHh5Vd28rU3nGKutvOYpSdYlWbdp06YdEVOSJGkwgxS1JHejK2nvq6oP98M3Jjmgf/wA4KZ+fCNw8MjTDwKun+t1q+qsqlpeVcuXLFkymfCSJElTMsRZnwHeCVxZVW8YeegCYGW/vBI4f2T8xCS7JjkEWAZcMq28kiRJQxnipuyPAZ4PfCXJZf3Yq4HVwNokJwPXAicAVNXlSdYCV9CdMXpqVd0+9dSSJElTNvWiVlX/wtzHnQE8cSvPOQM4Y2KhJEmSGuSdCSRJkhplUZMkSWqURU2SJKlRFjVJkqRGWdQkSZIaZVGTJElqlEVNkiSpURY1SZKkRlnUJEmSGmVRkyRJapRFTZIkqVEWNUmSpEZZ1CRJkhplUZMkSWqURU2SJKlRFjVJkqRGWdQkSZIaZVGTJElqlEVNkiSpURY1SZKkRlnUJEmSGmVRkyRJapRFTZIkqVGLhw6g6Vi66sKhI/zYhtUrho4gSdK84IyaJElSoyxqkiRJjbKoSZIkNcqiJkmS1CiLmiRJUqMsapIkSY2yqEmSJDXKoiZJktQoi5okSVKjLGqSJEmNsqhJkiQ1yqImSZLUKIuaJElSoyxqkiRJjbKoSZIkNWrx0AHms6WrLhw6wo9tWL1i6AiSJGkHc0ZNkiSpURY1SZKkRlnUJEmSGmVRkyRJatS8KWpJnprkqiTrk6waOo8kSdKkzYuilmQR8BbgacBhwHOTHDZsKkmSpMmaF0UNOApYX1XfqKrbgHOBYwfOJEmSNFHzpagdCFw3sr6xH5MkSVqwUlVDZ9iuJCcAx1TVr/frzweOqqqXztruFOCUfvVQ4KqpBv3J7Qt8a+gQY5ovWedLTjDrpJh1Msw6GWadjPmU9QFVtWT24Hy5M8FG4OCR9YOA62dvVFVnAWdNK9SOkmRdVS0fOsc45kvW+ZITzDopZp0Ms06GWSdjPmXdmvmy6/PzwLIkhyS5O3AicMHAmSRJkiZqXsyoVdXmJL8N/D2wCHhXVV0+cCxJkqSJmhdFDaCq/hb426FzTMh82l07X7LOl5xg1kkx62SYdTLMOhnzKeuc5sXJBJIkSTuj+XKMmiRJ0k7HoiYNLEmGzrCQJNlz6AzjSnJf//9L2haL2oD6W2M1LclPJ1meZNehs2xPksOTHJ3kPkNn2Z4kj+2vB0hVVcu/rJM8M8nLhs4xjiTHAn+SZL+hs2xPkmOAj7DlpYealOSRSZ7ff7370Hm2Jcmy/mfWovnwM1baHovaAJI8CKCqbm/5B0mSZwAfBv4U+KuZ3C1K8jTgA8ArgPckue/AkeaUZJck9wDeDpyW5Dfgx2WtuX+PSZ4CvB64Yugs25PkaOBPgPOr6qah82xL/339E+AA4JUDx9mmJM+iOyD7ScCrgAcMm2jrkhwHnAecBrwBeMl8mmEFZ9gnZT5/X5v7xbDQ9eXnsiTvh3bLWpJHA38GrKyqxwP/BawaNtXckjwOeBPw61V1HHAb8JABI21VVf2oqr4PrAHeCTw6yStmHhs03Cz934H3AqdU1UVJ7pXkAUn2GDrbVjwceEef9X5Jnpzk55Pca+hgo5I8CTgTeB6wDPiZJL84bKq59bPTpwK/WlUrgZuBhyXZL8luw6bbUp/1JcBzq+p44EvAi4BXJLnnoOG2of87enSSR0DbM+xJ9ho6w7iSHNnvuTgKuu/r0Jl+Uha1Keo/2f028HLgtiTnQLtlDVhdVZf2y6cD+zS6C/RG4CVVdUk/k/bzwG8neXuS5zT6Q28z3S6vNcBRSd6Q5I/TaeXf5beBHwIH9L8EPwq8lW52tcXv6+aR5fOAX6P79/aWJHsPE2lOi4AX9NeC3JPuVneHQ5Of+jcDuwMP7n9JPw54AfBG4Pcam63aDNwDuC9AVb0L+CawBHjGgLm2qt8TcA5daX9NkndCm2UtyS8Bn+6LZSs/o+bUT4i8k+6Wkq9K8pKBI90lTX+zF5qq+gHdL4/30+1C2G20rA2ZbQ6fo9vtOXMs3a50uzz26seaOQ6sqq6sqn/sV08Gzuxn1v4NOIHuXm+tOR/4j6q6GFgH/AawV3WamFmrqquAFcBf0M1OvJ/uF97HgeOBlsoPwCeBFyc5Fzi7qp5L9wHj+8BRgyYbUVV/X1WfSbJLVX0HuBA4PcnPtvapv6q+C7yZblfiJ4B3V9UzgXfQ3crvpweMt4U+6/uAF/XH050B/DfdbvsnDxpuDv3P1ZXA66rqFLoCfGiS86CtspZkKfA7wE10h5cc2Uq22ZIcAfwv4IVV9QLgg8CDh01111jUpqyqrq+q71fVt+im6XefKWv9VG0Tf6Gq6vaqurlfDfAd4D+ralOS5wF/lGT3wQJuRVWdUVV/1C+/G7gnbR6sfQvdD+UX05W01cD9W/vkV1Vfoitnf1xVZ/e7bt9FV9LuP2y6LVXVV+k+AP08cEg/9g26Gaw73Oh4aDOFvKo+TncM2DMam1EFoKrOozs+7dPApf3YJ+n+bbV2vNoH6D5IPAHYo6pOqqq3A/u1ttuu/3B+6cj6zVX1WGD/JG/vx1op7j8CXlNVT6Yrvn8APDzJFhfNb6S87U73Yf1L/fqlwGOSHNxIvjutqR8IO5uq+jZdWfthkq8Bf0336b8pVbW5P67quiR/TPfJ6syqumXgaFuY/Y8wyfHA/sD1wyTauqq6HrgO+H3gd6rqdXQHPzd3942quqKq3jKz3n9flwA3DJdqq/6ObhbtpCQnJzkZOAL47LCxtutLdLOXu7Qyozqqqv6Lbsby+CRP6U8wOAT48rDJtlRV362q9wEnV9UrAJK8ANgHaGKvRbY8Kevfgd9NMvqh59nAfZIcNt1kd5T/f+LbtcBl/fLr6O6/fTrdvy2S/Gz/2GDFciTrZ4AP9WOL6H7+3wh8t5+lXDZUxp+URW1g/czal4F7A79UVRuHTXRH/af8uwO/QHcsxYlV9ZWBY93BzA+JJLv2v6BfR3cyxH8Mm2yrzgaeXVUf69f/qaquGzLQtvR/D36N7izQF1TVjUNnmq3/UPEe4DnATwGHAS+qqq8Pm2zbqupDdL9QWpz9nfFZutmJ19CdYPCiqtowaKKtmCm7/d/XVXTHsP5g2FRbnEx2LkBVnUN3iZZ/nSlr/e+EzXQzloMZyfqBPtd3+98DVNXrgUvoTtRYDbwvA14SZ47v66b+0ILb6XZ/L+q3ez7w540ds7pd3kJqYP1fmLXAK6uqqU+nsyV5IfD5/iDoZiW5G90xKV/vj7NqWpI0tItjq/oZy6Ppjq372tB5For58v9/RrozKDNyaESzkjwAuFtVrW8gy550Mz0fBh4N7NofR0mS1wPPojsbeF/gJODpVXVNI1kXV9VJ/WO7VtWt/fKngAcBxwz14X07WRfRHbrzAeC7wMPoPmQ2f7mhURa1BiTZrar+e+gc2zPffqFIUkuS3I/uEie7AW8DfjhS1p5Nd8bqw4E39sdcDmaOrP89U4D6xx9Ed7jOC0eOBxvEGFk/Slconz0fPrzPZlGTJGnK+jPnzwJuq6rnJjkc+H5VfXPgaHcwkvWWqjopycPorgBwRb+rthlzZF1Gdz29c+bbTNoMi5okSQNIsi/dnV8eTXcc1eNaPE4Ztsj6KLqsR/cnRTVnJOtj+qFfaPGY2nF5MoEkSQMYOZnsXnS75ZosaTDniW9NljTYIutewPHzuaSBRU2SpEH0J5M9HXhKi2fSjzLrcNz1KUnSQObLyWRg1qFY1CRJkhrlrk9JkqRGWdQkSZIaZVGTJElqlEVN0oKQZP8k70/yjSRfSPLZ/mrv08zwyCSfS3JZkiuTvLYff1aSVdPMImlh8GQCSfNefx/SzwBrqupt/dgDgGdV1f+ete3iqto8oRxXAb9cVV/q7zN46Hy9GrqkNjijJmkheALdrXjeNjNQVd+cKWlJXpjkg0k+BnwiyT5JPprky0n+LclD++1em+RVM6+R5KtJlvZ/vpZkTf+c85LsMUeO/YAb+ve/faak9e//l/3yZSN/bklydJI9k7wryeeTXJrk2El9oyTNLxY1SQvB4cAXt7PNo4CVVfUE4A+BS6vqocCrgfeM8R6HAmf1z7kZ+K05tvkL4KokH0nykiS7zd6gqh5WVQ8Dfh9YRzcT+Brgk1X1CODxwJ8m2XOMTJIWOIuapAUnyVuSfCnJ50eGL6qq/+yXHwu8F6CqPgncJ8m9tvOy11XVv/bL5/SvsYWqeh2wHPgE8KvAx7eSbxndvQh/pap+CDwFWJXkMuBTwG7A/bf33ylp4Vs8dABJ2gEuB46fWamqU/sbM68b2eYHI8uZ4zUK2MyWH2B3m/X47O3v+CJVXwfemuRsYFOS+4w+3s+UrQVePHK/xNDdk/CquV5T0s7LGTVJC8Engd2S/ObI2FzHkM34Z+B5AEkeB3yrqm4GNgBH9uNHAoeMPOf+SR7VLz8X+JfZL5pkRX9iA8Ay4HbgO7M2ezfw7qr69MjY3wMvnXlukiO2kV3STsSiJmneq+709eOAo5Nck+QSYA3wu1t5ymuB5Um+DKwGVvbjHwL26XdB/ibwf0aecyWwsn/OPsBb53jd59Mdo3YZ3a7V51XV7TMP9meiPgf4tZETCpYDrwfuBnw5yVf7dUny8hyStD1JlgJ/U1UPGTqLpJ2LM2qSJEmNckZNkiSpUc6oSZIkNcqiJkmS1CiLmiRJUqMsapIkSY2yqEmSJDXKoiZJktSo/wcBHyIJQ3bUngAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\r\n",
    "\r\n",
    "sizes = unique_perms.loc[:, \"size\"]\r\n",
    "\r\n",
    "def plot_frequency_hist(x: pd.Series) -> tuple:\r\n",
    "    \"\"\" Creates the frequency histogram plot of the series x. \"\"\" \r\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10,5))\r\n",
    "\r\n",
    "    # define bins and binwidth\r\n",
    "    bins = np.arange(min(x.unique()), min(x.unique()) + len(x.unique()) + 1)\r\n",
    "    binwidth = 0.8\r\n",
    "\r\n",
    "    # plot hist\r\n",
    "    plt.hist(x.values, bins=bins-0.5*binwidth, width=binwidth)\r\n",
    "\r\n",
    "    # formatting\r\n",
    "    ax.set_xlabel(\"Group Size\"); ax.set_ylabel(\"Frequency\")\r\n",
    "    ax.set_xlim([min(x.unique())-1,min(x.unique())+len(x.unique())])\r\n",
    "    ax.set_xticks(bins[:-1])\r\n",
    "    ax.set_xticklabels(labels=bins[:-1],rotation=45)\r\n",
    "\r\n",
    "    return fig, ax\r\n",
    "\r\n",
    "fig, ax = plot_frequency_hist(sizes)\r\n",
    "\r\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## e) Comment upon the distribution of group sizes in d)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'poisson_gen' object has no attribute 'fit'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-48-ad0c481babd6>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     34\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     35\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 36\u001b[1;33m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0max\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mhist_fit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msizes\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdists\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'poisson'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     37\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     38\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshow\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-48-ad0c481babd6>\u001b[0m in \u001b[0;36mhist_fit\u001b[1;34m(x, dists)\u001b[0m\n\u001b[0;32m     19\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     20\u001b[0m         \u001b[0mdist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mst\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdist_name\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 21\u001b[1;33m         \u001b[0mparams\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdist\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     22\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     23\u001b[0m         \u001b[0mp\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdist\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpdf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAttributeError\u001b[0m: 'poisson_gen' object has no attribute 'fit'"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAEvCAYAAAB2a9QGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAV50lEQVR4nO3db4hd933n8fdnRxGpkxp3yezGleQdFURcEZq1GbxODXlQNyBZIcqDfSCziYNbEAY7tktCdtJ90GeLlg0hMWs8CFtpTUzM4nhZEU3jhCZhKdhG4z84kRVvB1UbTS2vp5jYIQbLwt99cG/hdjzyHEn35/tn3y8YdM7vz7nfe9BIH37n3HNTVUiSJGm4/sWoC5AkSZpGhixJkqQGDFmSJEkNGLIkSZIaMGRJkiQ1YMiSJElqYMuoC9jIRz7ykZqbmxt1GZIkSZt65pln/rGqZte3j2XImpubY3l5edRlSJIkbSrJ/9mo3cuFkiRJDXQKWUn2JHkpyUqShQ36r03yZJK3knxlXd9VSR5L8oskJ5N8cljFS5IkjatNLxcmmQHuBz4NrALHkxytqhcHhr0G3A18boNDfAv4QVX9+yRbgSsuu2pJkqQx12Ul6wZgpapOVdU54FFg/+CAqnq1qo4Dbw+2J7kS+BTwUH/cuar61TAKlyRJGmddQtY24MzA/mq/rYvfA9aAbyd5LsmDST50kTVKkiRNnC4hKxu0VcfjbwGuBx6oquuA3wDvuqcLIMnBJMtJltfW1joeXpIkaTx1CVmrwI6B/e3Ayx2PvwqsVtXT/f3H6IWud6mqw1U1X1Xzs7PvetSEJEnSROkSso4Du5Ls7N+4fgA42uXgVfUKcCbJx/pNNwMvvscUSZKkqbDppwur6nySu4AngBngSFWdSHJHv38xyUeBZeBK4J0k9wK7q+oN4EvAI/2Adgq4vc1bkSRJGh+dnvheVUvA0rq2xYHtV+hdRtxo7vPA/KWXKEmSNHl84rskSVIDY/ndhZLUwtzCsVGXwOlD+0ZdgqT3iStZkiRJDRiyJEmSGjBkSZIkNWDIkiRJasCQJUmS1IAhS5IkqQEf4SBJY2YcHjUBPm5CulyuZEmSJDVgyJIkSWrAkCVJktSAIUuSJKkBQ5YkSVIDhixJkqQGDFmSJEkNGLIkSZIaMGRJkiQ1YMiSJElqwJAlSZLUgCFLkiSpAUOWJElSA4YsSZKkBgxZkiRJDXQKWUn2JHkpyUqShQ36r03yZJK3knxlg/6ZJM8l+f4wipYkSRp3m4asJDPA/cBeYDdwa5Ld64a9BtwNfP0Ch7kHOHkZdUqSJE2ULitZNwArVXWqqs4BjwL7BwdU1atVdRx4e/3kJNuBfcCDQ6hXkiRpInQJWduAMwP7q/22rr4JfBV45yLmSJIkTbQuISsbtFWXgyf5DPBqVT3TYezBJMtJltfW1rocXpIkaWxt6TBmFdgxsL8deLnj8W8CPpvkFuCDwJVJvlNVn18/sKoOA4cB5ufnO4U4SaM3t3Bs1CUAcPrQvlGXIEn/TJeVrOPAriQ7k2wFDgBHuxy8qr5WVduraq4/78cbBSxJkqRps+lKVlWdT3IX8AQwAxypqhNJ7uj3Lyb5KLAMXAm8k+ReYHdVvdGudEmSpPHV5XIhVbUELK1rWxzYfoXeZcT3OsZPgZ9edIWSJEkTyCe+S5IkNWDIkiRJasCQJUmS1IAhS5IkqQFDliRJUgOGLEmSpAYMWZIkSQ0YsiRJkhowZEmSJDVgyJIkSWrAkCVJktSAIUuSJKkBQ5YkSVIDhixJkqQGDFmSJEkNGLIkSZIaMGRJkiQ1YMiSJElqwJAlSZLUgCFLkiSpAUOWJElSA4YsSZKkBgxZkiRJDRiyJEmSGugUspLsSfJSkpUkCxv0X5vkySRvJfnKQPuOJD9JcjLJiST3DLN4SZKkcbVlswFJZoD7gU8Dq8DxJEer6sWBYa8BdwOfWzf9PPDlqno2yW8DzyT50bq5kiRJU6fLStYNwEpVnaqqc8CjwP7BAVX1alUdB95e1362qp7tb/8aOAlsG0rlkiRJY6xLyNoGnBnYX+USglKSOeA64OmLnStJkjRpNr1cCGSDtrqYF0nyYeB7wL1V9cYFxhwEDgJcc801F3N4SdIIzC0cG3UJAJw+tG/UJUgb6rKStQrsGNjfDrzc9QWSfIBewHqkqh6/0LiqOlxV81U1Pzs72/XwkiRJY6lLyDoO7EqyM8lW4ABwtMvBkwR4CDhZVd+49DIlSZImy6aXC6vqfJK7gCeAGeBIVZ1Icke/fzHJR4Fl4ErgnST3AruBPwC+APwsyfP9Q/55VS0N/Z1IkiSNkS73ZNEPRUvr2hYHtl+hdxlxvb9l43u6JEmSpppPfJckSWrAkCVJktSAIUuSJKkBQ5YkSVIDhixJkqQGDFmSJEkNGLIkSZIaMGRJkiQ1YMiSJElqwJAlSZLUgCFLkiSpAUOWJElSA4YsSZKkBgxZkiRJDRiyJEmSGjBkSZIkNWDIkiRJasCQJUmS1IAhS5IkqQFDliRJUgOGLEmSpAa2jLoASRubWzg26hIAOH1o36hLkKSJ5EqWJElSA4YsSZKkBjqFrCR7kryUZCXJwgb91yZ5MslbSb5yMXMlSZKm0aYhK8kMcD+wF9gN3Jpk97phrwF3A1+/hLmSJElTp8tK1g3ASlWdqqpzwKPA/sEBVfVqVR0H3r7YuZIkSdOoS8jaBpwZ2F/tt3VxOXMlSZImVpeQlQ3aquPxO89NcjDJcpLltbW1joeXJEkaT11C1iqwY2B/O/Byx+N3nltVh6tqvqrmZ2dnOx5ekiRpPHUJWceBXUl2JtkKHACOdjz+5cyVJEmaWJs+8b2qzie5C3gCmAGOVNWJJHf0+xeTfBRYBq4E3klyL7C7qt7YaG6j9yJJkjQ2On2tTlUtAUvr2hYHtl+hdymw01xJkqRp5xPfJUmSGjBkSZIkNWDIkiRJasCQJUmS1ECnG98lSZpUcwvHRl0CAKcP7Rt1CXqfuZIlSZLUgCFLkiSpAUOWJElSA4YsSZKkBgxZkiRJDRiyJEmSGjBkSZIkNWDIkiRJasCQJUmS1IAhS5IkqQFDliRJUgOGLEmSpAYMWZIkSQ0YsiRJkhowZEmSJDVgyJIkSWrAkCVJktSAIUuSJKkBQ5YkSVIDhixJkqQGOoWsJHuSvJRkJcnCBv1Jcl+//4Uk1w/0/VmSE0l+nuS7ST44zDcgSZI0jjYNWUlmgPuBvcBu4NYku9cN2wvs6v8cBB7oz90G3A3MV9XHgRngwNCqlyRJGlNdVrJuAFaq6lRVnQMeBfavG7MfeLh6ngKuSnJ1v28L8FtJtgBXAC8PqXZJkqSx1SVkbQPODOyv9ts2HVNV/wB8HfglcBZ4vap+uNGLJDmYZDnJ8traWtf6JUmSxlKXkJUN2qrLmCS/Q2+Vayfwu8CHknx+oxepqsNVNV9V87Ozsx3KkiRJGl9dQtYqsGNgfzvvvuR3oTF/DPx9Va1V1dvA48AfXnq5kiRJk6FLyDoO7EqyM8lWejeuH1035ihwW/9ThjfSuyx4lt5lwhuTXJEkwM3AySHWL0mSNJa2bDagqs4nuQt4gt6nA49U1Ykkd/T7F4El4BZgBXgTuL3f93SSx4BngfPAc8DhFm9EkiRpnGwasgCqaolekBpsWxzYLuDOC8z9C+AvLqNGSZKkieMT3yVJkhowZEmSJDVgyJIkSWrAkCVJktSAIUuSJKkBQ5YkSVIDhixJkqQGOj0nS5omcwvHRl0CAKcP7Rt1CZKkhlzJkiRJasCQJUmS1IAhS5IkqQFDliRJUgOGLEmSpAYMWZIkSQ0YsiRJkhowZEmSJDVgyJIkSWrAkCVJktSAIUuSJKkBQ5YkSVIDhixJkqQGDFmSJEkNGLIkSZIaMGRJkiQ10ClkJdmT5KUkK0kWNuhPkvv6/S8kuX6g76okjyX5RZKTST45zDcgSZI0jjYNWUlmgPuBvcBu4NYku9cN2wvs6v8cBB4Y6PsW8IOquhb4BHByCHVLkiSNtS4rWTcAK1V1qqrOAY8C+9eN2Q88XD1PAVcluTrJlcCngIcAqupcVf1qeOVLkiSNpy0dxmwDzgzsrwL/rsOYbcB5YA34dpJPAM8A91TVb9a/SJKD9FbBuOaaa7rWL0nS1JhbODbqEjh9aN+oS5gaXVayskFbdRyzBbgeeKCqrgN+A7zrni6AqjpcVfNVNT87O9uhLEmSpPHVJWStAjsG9rcDL3ccswqsVtXT/fbH6IUuSZKkqdYlZB0HdiXZmWQrcAA4um7MUeC2/qcMbwRer6qzVfUKcCbJx/rjbgZeHFbxkiRJ42rTe7Kq6nySu4AngBngSFWdSHJHv38RWAJuAVaAN4HbBw7xJeCRfkA7ta5PkiRpKnW58Z2qWqIXpAbbFge2C7jzAnOfB+YvvURJkqTJ4xPfJUmSGjBkSZIkNWDIkiRJasCQJUmS1IAhS5IkqQFDliRJUgOGLEmSpAYMWZIkSQ0YsiRJkhowZEmSJDVgyJIkSWrAkCVJktSAIUuSJKkBQ5YkSVIDhixJkqQGDFmSJEkNGLIkSZIaMGRJkiQ1YMiSJElqwJAlSZLUgCFLkiSpAUOWJElSA1tGXYAkSZoscwvHRl0CAKcP7Rt1Ce/JlSxJkqQGOoWsJHuSvJRkJcnCBv1Jcl+//4Uk16/rn0nyXJLvD6twSZKkcbZpyEoyA9wP7AV2A7cm2b1u2F5gV//nIPDAuv57gJOXXa0kSdKE6LKSdQOwUlWnquoc8Ciwf92Y/cDD1fMUcFWSqwGSbAf2AQ8OsW5JkqSx1iVkbQPODOyv9tu6jvkm8FXgnUsrUZIkafJ0CVnZoK26jEnyGeDVqnpm0xdJDiZZTrK8trbWoSxJkqTx1SVkrQI7Bva3Ay93HHMT8Nkkp+ldZvyjJN/Z6EWq6nBVzVfV/OzsbMfyJUmSxlOXkHUc2JVkZ5KtwAHg6LoxR4Hb+p8yvBF4varOVtXXqmp7Vc315/24qj4/zDcgSZI0jjZ9GGlVnU9yF/AEMAMcqaoTSe7o9y8CS8AtwArwJnB7u5IlSZLGX6cnvlfVEr0gNdi2OLBdwJ2bHOOnwE8vukJNBJ/+K0nSP+cT3yVJkhowZEmSJDVgyJIkSWrAkCVJktSAIUuSJKkBQ5YkSVIDhixJkqQGDFmSJEkNGLIkSZIaMGRJkiQ1YMiSJElqwJAlSZLUgCFLkiSpAUOWJElSA4YsSZKkBgxZkiRJDRiyJEmSGjBkSZIkNWDIkiRJasCQJUmS1IAhS5IkqQFDliRJUgOGLEmSpAYMWZIkSQ10CllJ9iR5KclKkoUN+pPkvn7/C0mu77fvSPKTJCeTnEhyz7DfgCRJ0jjaNGQlmQHuB/YCu4Fbk+xeN2wvsKv/cxB4oN9+HvhyVf0+cCNw5wZzJUmSpk6XlawbgJWqOlVV54BHgf3rxuwHHq6ep4CrklxdVWer6lmAqvo1cBLYNsT6JUmSxlKXkLUNODOwv8q7g9KmY5LMAdcBT190lZIkSROmS8jKBm11MWOSfBj4HnBvVb2x4YskB5MsJ1leW1vrUJYkSdL46hKyVoEdA/vbgZe7jknyAXoB65GqevxCL1JVh6tqvqrmZ2dnu9QuSZI0trqErOPAriQ7k2wFDgBH1405CtzW/5ThjcDrVXU2SYCHgJNV9Y2hVi5JkjTGtmw2oKrOJ7kLeAKYAY5U1Ykkd/T7F4El4BZgBXgTuL0//SbgC8DPkjzfb/vzqloa6ruQJEkaM5uGLIB+KFpa17Y4sF3AnRvM+1s2vl9LkiRpqvnEd0mSpAYMWZIkSQ0YsiRJkhrodE+WRmdu4dioSwDg9KF9oy5BkqSJ4kqWJElSA4YsSZKkBgxZkiRJDRiyJEmSGjBkSZIkNWDIkiRJasCQJUmS1IAhS5IkqQFDliRJUgOGLEmSpAYMWZIkSQ0YsiRJkhowZEmSJDVgyJIkSWrAkCVJktSAIUuSJKmBLaMuYFTmFo6NugQATh/aN+oSJElSA65kSZIkNWDIkiRJasCQJUmS1ECnkJVkT5KXkqwkWdigP0nu6/e/kOT6rnMlSZKm0aYhK8kMcD+wF9gN3Jpk97phe4Fd/Z+DwAMXMVeSJGnqdFnJugFYqapTVXUOeBTYv27MfuDh6nkKuCrJ1R3nSpIkTZ0uIWsbcGZgf7Xf1mVMl7mSJElTp8tzsrJBW3Uc02Vu7wDJQXqXGgHeSvLzDrVNvPyX9+VlPgL84+Uc4H2q87K9j3V6TofL8zl8ntPh8nwO3zSd03+zUWOXkLUK7BjY3w683HHM1g5zAaiqw8BhgCTLVTXfoTZ14PkcPs/pcHk+h89zOlyez+H7/+GcdrlceBzYlWRnkq3AAeDoujFHgdv6nzK8EXi9qs52nCtJkjR1Nl3JqqrzSe4CngBmgCNVdSLJHf3+RWAJuAVYAd4Ebn+vuU3eiSRJ0hjp9N2FVbVEL0gNti0ObBdwZ9e5HRy+yPF6b57P4fOcDpfnc/g8p8Pl+Ry+qT+n6eUjSZIkDZNfqyNJktTAWIUsv4JnuJLsSPKTJCeTnEhyz6hrmgZJZpI8l+T7o65lGiS5KsljSX7R/7v6yVHXNMmS/Fn/9/3nSb6b5IOjrmnSJDmS5NXBRwkl+ZdJfpTk7/p//s4oa5w0Fzin/7X/e/9Ckv+R5KoRltjE2IQsv4KnifPAl6vq94EbgTs9p0NxD3By1EVMkW8BP6iqa4FP4Lm9ZEm2AXcD81X1cXofODow2qom0l8Ce9a1LQB/U1W7gL/p76u7v+Td5/RHwMer6g+A/w187f0uqrWxCVn4FTxDV1Vnq+rZ/vav6f3n5RP3L0OS7cA+4MFR1zINklwJfAp4CKCqzlXVr0Za1OTbAvxWki3AFVzg2YS6sKr6X8Br65r3A3/V3/4r4HPvZ02TbqNzWlU/rKrz/d2n6D1Lc6qMU8jyK3gaSjIHXAc8PeJSJt03ga8C74y4jmnxe8Aa8O3+JdgHk3xo1EVNqqr6B+DrwC+Bs/SeWfjD0VY1Nf51//mP9P/8VyOuZ9r8CfDXoy5i2MYpZHX+Ch5dnCQfBr4H3FtVb4y6nkmV5DPAq1X1zKhrmSJbgOuBB6rqOuA3eBnmkvXvE9oP7AR+F/hQks+PtirpvSX5T/Rub3lk1LUM2ziFrC5f36OLlOQD9ALWI1X1+KjrmXA3AZ9Ncpre5ew/SvKd0ZY08VaB1ar6pxXWx+iFLl2aPwb+vqrWqupt4HHgD0dc07T4v0muBuj/+eqI65kKSb4IfAb4DzWFz5Qap5DlV/AMWZLQu9flZFV9Y9T1TLqq+lpVba+qOXp/P39cVa4SXIaqegU4k+Rj/aabgRdHWNKk+yVwY5Ir+r//N+MHCYblKPDF/vYXgf85wlqmQpI9wH8EPltVb466nhbGJmT1b377p6/gOQn8d7+C57LdBHyB3orL8/2fW0ZdlLTOl4BHkrwA/FvgP4+2nMnVXxF8DHgW+Bm9f+On/qnaw5bku8CTwMeSrCb5U+AQ8Okkfwd8ur+vji5wTv8b8NvAj/r/Py2+50EmkE98lyRJamBsVrIkSZKmiSFLkiSpAUOWJElSA4YsSZKkBgxZkiRJDRiyJEmSGjBkSZIkNWDIkiRJauD/ARngzpS/9ah8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats as st\r\n",
    "\r\n",
    "def hist_fit(x: pd.Series, dists: list) -> tuple:\r\n",
    "    \"\"\" Creates density histogram plot with fitted distributions\"\"\" \r\n",
    "\r\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10,5))\r\n",
    "\r\n",
    "    # define bins and binwidth\r\n",
    "    bins = np.arange(min(x.unique()), min(x.unique()) + len(x.unique()) + 1)\r\n",
    "    binwidth = 0.8\r\n",
    "\r\n",
    "    # plot density\r\n",
    "    ax.hist(x.values, bins=bins - 0.5*binwidth,width=binwidth, density=True, label=\"Frequency Density\")\r\n",
    "\r\n",
    "    # experiment with other distributions by adding them to the list\r\n",
    "    for dist_name in dists:\r\n",
    "        xmin, xmax = ax.get_xlim()\r\n",
    "        xx = np.linspace(xmin, xmax, 100)\r\n",
    "\r\n",
    "        dist = getattr(st, dist_name)\r\n",
    "        params = dist.fit(x)\r\n",
    "\r\n",
    "        p = dist.pdf(xx, *params)\r\n",
    "        ax.plot(xx, p, linewidth=2, label=dist_name.title() + \" Distribution Fit\")\r\n",
    "\r\n",
    "    # formatting\r\n",
    "    ax.set_xlabel(\"Group Size\")\r\n",
    "    ax.set_ylabel(\"Density\")\r\n",
    "    ax.set_xlim([min(x.unique()) - 1, min(x.unique()) + 1 + len(x.unique())])\r\n",
    "    ax.set_xticks(bins[:-1])\r\n",
    "    ax.set_xticklabels(labels=bins[:-1], rotation=45)\r\n",
    "    plt.legend()\r\n",
    "\r\n",
    "    return fig, ax\r\n",
    "\r\n",
    "fig, ax = hist_fit(sizes, dists=['norm'])\r\n",
    "\r\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## f) If your random variables were, in fact, meaningful information on individuals, which group sizes are of most concern from a privacy perspective?\r\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## g) Consider the effect of missing data in the dataset you created in 1 a). How might this complicate the production of a frequency table of group sizes in 1 d)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "size\n",
       "1    16799\n",
       "2     7076\n",
       "3     2145\n",
       "4      496\n",
       "5      105\n",
       "6       22\n",
       "7        3\n",
       "Name: A, dtype: int64"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAE+CAYAAADia0YcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAY+UlEQVR4nO3dfbRldX3f8fcHhvIUIfIo5WlInKJADcKIJNoCUhWDClZIhlIdE3SQQFdMzKqDtZU0YQm1kcREiFCVAauIGAVrUBHS2DQKDBF5pkwEZTIUxkfQIjjjt3+cfeOZmTP3HmDO7Pub+36tddbd53v2b9/v2WvmzGd+e++zU1VIkiSpLVv13YAkSZKeOkOcJElSgwxxkiRJDTLESZIkNcgQJ0mS1KB5fTewue222241f/78vtuQJEma0S233PLtqtp91GtzLsTNnz+f5cuX992GJEnSjJJ8c2OveThVkiSpQYY4SZKkBhniJEmSGmSIkyRJapAhTpIkqUGGOEmSpAYZ4iRJkhpkiJMkSWqQIU6SJKlBhjhJkqQGGeIkSZIaNOfundqK+Us/13cLE/XAecf33YIkSU1zJk6SJKlBhjhJkqQGGeIkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGmSIkyRJapAhTpIkqUGGOEmSpAYZ4iRJkhpkiJMkSWqQIU6SJKlBhjhJkqQGGeIkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGjSxEJdk3yR/leTuJHcm+e2uvkuS65Lc1/189tCYs5OsSHJvklcO1Q9Pcnv32vuTpKtvm+QTXf3GJPMn9X4kSZJmk0nOxK0B3l5VzweOBM5MchCwFLi+qhYA13fP6V5bBBwMHAdcmGTrblsXAUuABd3juK5+GvC9qnoucAFw/gTfjyRJ0qwxsRBXVQ9V1d91y48BdwN7AycAy7rVlgEndssnAFdU1RNVdT+wAjgiyV7ATlX1laoq4LL1xkxt6yrg2KlZOkmSpC3ZZjknrjvM+ULgRmDPqnoIBkEP2KNbbW/gwaFhK7va3t3y+vV1xlTVGuAHwK4jfv+SJMuTLF+9evUmeleSJEn9mXiIS/JzwKeAt1XVo9OtOqJW09SnG7NuoeriqlpYVQt33333mVqWJEma9SYa4pJswyDA/feq+ouu/HB3iJTu5yNdfSWw79DwfYBVXX2fEfV1xiSZB+wMfHfTvxNJkqTZZZJXpwb4EHB3Vb1v6KVrgMXd8mLg6qH6ou6K0wMYXMBwU3fI9bEkR3bbfON6Y6a2dRJwQ3fenCRJ0hZt3gS3/RLgDcDtSW7tau8EzgOuTHIa8C3gZICqujPJlcBdDK5sPbOq1nbjzgAuBbYHru0eMAiJlydZwWAGbtEE348kSdKsMbEQV1V/w+hz1gCO3ciYc4FzR9SXA4eMqP+YLgRKkiTNJd6xQZIkqUGGOEmSpAYZ4iRJkhpkiJMkSWqQIU6SJKlBhjhJkqQGGeIkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGmSIkyRJapAhTpIkqUGGOEmSpAYZ4iRJkhpkiJMkSWqQIU6SJKlBhjhJkqQGGeIkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGmSIkyRJapAhTpIkqUGGOEmSpAYZ4iRJkhpkiJMkSWqQIU6SJKlBhjhJkqQGGeIkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGmSIkyRJapAhTpIkqUGGOEmSpAYZ4iRJkhpkiJMkSWqQIU6SJKlBhjhJkqQGGeIkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGmSIkyRJapAhTpIkqUGGOEmSpAZNLMQl+XCSR5LcMVQ7J8k/JLm1e/zq0GtnJ1mR5N4krxyqH57k9u619ydJV982ySe6+o1J5k/qvUiSJM02k5yJuxQ4bkT9gqo6tHv8JUCSg4BFwMHdmAuTbN2tfxGwBFjQPaa2eRrwvap6LnABcP6k3ogkSdJsM7EQV1VfBr475uonAFdU1RNVdT+wAjgiyV7ATlX1laoq4DLgxKExy7rlq4Bjp2bpJEmStnR9nBN3VpLbusOtz+5qewMPDq2zsqvt3S2vX19nTFWtAX4A7DrqFyZZkmR5kuWrV6/edO9EkiSpJ5s7xF0E/CJwKPAQ8EddfdQMWk1Tn27MhsWqi6tqYVUt3H333Z9Sw5IkSbPRZg1xVfVwVa2tqp8ClwBHdC+tBPYdWnUfYFVX32dEfZ0xSeYBOzP+4VtJkqSmbdYQ153jNuV1wNSVq9cAi7orTg9gcAHDTVX1EPBYkiO7893eCFw9NGZxt3wScEN33pwkSdIWb96kNpzk48DRwG5JVgLvBo5OciiDw54PAKcDVNWdSa4E7gLWAGdW1dpuU2cwuNJ1e+Da7gHwIeDyJCsYzMAtmtR7kSRJmm0mFuKq6pQR5Q9Ns/65wLkj6suBQ0bUfwyc/Ex6lCRJapV3bJAkSWqQIU6SJKlBhjhJkqQGGeIkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGmSIkyRJapAhTpIkqUGGOEmSpAYZ4iRJkho0VohLssEN6CVJktSfcWfi/jzJTUl+K8nPT7IhSZIkzWysEFdVLwVOBfYFlif5WJKXT7QzSZIkbdTY58RV1X3Au4B3AEcB709yT5J/PanmJEmSNNq458S9IMkFwN3Ay4DXVNXzu+ULJtifJEmSRpg35np/BlwCvLOqHp8qVtWqJO+aSGeSJEnaqHFD3K8Cj1fVWoAkWwHbVdX/q6rLJ9adJEmSRhr3nLgvAdsPPd+hq0mSJKkH44a47arqh1NPuuUdJtOSJEmSZjJuiPtRksOmniQ5HHh8mvUlSZI0QeOeE/c24JNJVnXP9wJ+fSIdSZIkaUZjhbiqujnJ84ADgQD3VNVPJtqZJEmSNmrcmTiAFwHzuzEvTEJVXTaRriRJkjStsUJcksuBXwRuBdZ25QIMcZIkST0YdyZuIXBQVdUkm5EkSdJ4xr069Q7gOZNsRJIkSeMbdyZuN+CuJDcBT0wVq+q1E+lKkiRJ0xo3xJ0zySYkSZL01Iz7FSN/nWR/YEFVfSnJDsDWk21NkiRJGzPWOXFJ3gJcBXywK+0NfGZCPUmSJGkG417YcCbwEuBRgKq6D9hjUk1JkiRpeuOGuCeq6smpJ0nmMfieOEmSJPVg3BD310neCWyf5OXAJ4HPTq4tSZIkTWfcELcUWA3cDpwO/CXwrkk1JUmSpOmNe3XqT4FLuockSZJ6Nu69U+9nxDlwVfULm7wjSZIkzeip3Dt1ynbAycAum74dSZIkjWOsc+Kq6jtDj3+oqj8GXjbZ1iRJkrQx4x5OPWzo6VYMZuaeNZGOJEmSNKNxD6f+0dDyGuAB4Nc2eTeSJEkay7hXpx4z6UYkSZI0vnEPp/7udK9X1fs2TTuSJEkax1O5OvVFwDXd89cAXwYenERTkiRJmt64IW434LCqegwgyTnAJ6vqzZNqTJIkSRs37m239gOeHHr+JDB/k3cjSZKksYw7E3c5cFOSTzO4c8PrgMsm1pUkSZKmNe7VqecmuRb4F13pN6rqa5NrS5IkSdMZ93AqwA7Ao1X1J8DKJAdMqCdJkiTNYKwQl+TdwDuAs7vSNsBHZxjz4SSPJLljqLZLkuuS3Nf9fPbQa2cnWZHk3iSvHKofnuT27rX3J0lX3zbJJ7r6jUnmj/2uJUmSGjfuTNzrgNcCPwKoqlXMfNutS4Hj1qstBa6vqgXA9d1zkhwELAIO7sZcmGTrbsxFwBJgQfeY2uZpwPeq6rnABcD5Y74XSZKk5o0b4p6sqmJwUQNJdpxpQFV9GfjueuUTgGXd8jLgxKH6FVX1RFXdD6wAjkiyF7BTVX2l+/2XrTdmaltXAcdOzdJJkiRt6cYNcVcm+SDw80neAnwJuORp/L49q+ohgO7nHl19b9b94uCVXW3vbnn9+jpjqmoN8ANg11G/NMmSJMuTLF+9evXTaFuSJGl2mfHq1G526xPA84BHgQOB/1RV123CPkbNoNU09enGbFisuhi4GGDhwoUj15EkSWrJjCGuqirJZ6rqcOCZBreHk+xVVQ91h0of6eorgX2H1tsHWNXV9xlRHx6zMsk8YGc2PHwrSZK0RRr3cOpXk7xoE/y+a4DF3fJi4Oqh+qLuitMDGFzAcFN3yPWxJEd2M4JvXG/M1LZOAm7ozpuTJEna4o17x4ZjgLcmeYDBFaphMEn3go0NSPJx4GhgtyQrgXcD5zE4v+404FvAyQw2dGeSK4G7gDXAmVW1ttvUGQyudN0euLZ7AHwIuDzJCgYzcIvGfC+SJEnNmzbEJdmvqr4FvOqpbriqTtnIS8duZP1zgXNH1JcDh4yo/5guBEqSJM01M83EfQY4rKq+meRTVfX6zdCTJEmSZjDTOXHDV4D+wiQbkSRJ0vhmCnG1kWVJkiT1aKbDqb+U5FEGM3Lbd8vwswsbdppod5IkSRpp2hBXVVtP97okSZL6Me73xEmSJGkWMcRJkiQ1yBAnSZLUIEOcJElSg8a97ZY0K8xf+rm+W5iYB847vu8WJEkNcSZOkiSpQYY4SZKkBhniJEmSGmSIkyRJapAhTpIkqUGGOEmSpAYZ4iRJkhpkiJMkSWqQIU6SJKlBhjhJkqQGGeIkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGmSIkyRJapAhTpIkqUGGOEmSpAYZ4iRJkhpkiJMkSWqQIU6SJKlBhjhJkqQGGeIkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGmSIkyRJapAhTpIkqUGGOEmSpAYZ4iRJkhpkiJMkSWqQIU6SJKlBhjhJkqQGGeIkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGmSIkyRJalAvIS7JA0luT3JrkuVdbZck1yW5r/v57KH1z06yIsm9SV45VD+8286KJO9Pkj7ejyRJ0ubW50zcMVV1aFUt7J4vBa6vqgXA9d1zkhwELAIOBo4DLkyydTfmImAJsKB7HLcZ+5ckSerNbDqcegKwrFteBpw4VL+iqp6oqvuBFcARSfYCdqqqr1RVAZcNjZEkSdqi9RXiCvhikluSLOlqe1bVQwDdzz26+t7Ag0NjV3a1vbvl9esbSLIkyfIky1evXr0J34YkSVI/5vX0e19SVauS7AFcl+SeadYddZ5bTVPfsFh1MXAxwMKFC0euI0mS1JJeZuKqalX38xHg08ARwMPdIVK6n490q68E9h0avg+wqqvvM6IuSZK0xdvsIS7JjkmeNbUMvAK4A7gGWNytthi4ulu+BliUZNskBzC4gOGm7pDrY0mO7K5KfePQGEmSpC1aH4dT9wQ+3X0byDzgY1X1+SQ3A1cmOQ34FnAyQFXdmeRK4C5gDXBmVa3ttnUGcCmwPXBt95AkSdribfYQV1XfAH5pRP07wLEbGXMucO6I+nLgkE3doyRJ0mw3m75iRJIkSWMyxEmSJDXIECdJktQgQ5wkSVKDDHGSJEkNMsRJkiQ1yBAnSZLUIEOcJElSgwxxkiRJDTLESZIkNcgQJ0mS1CBDnCRJUoMMcZIkSQ0yxEmSJDXIECdJktQgQ5wkSVKDDHGSJEkNMsRJkiQ1yBAnSZLUIEOcJElSgwxxkiRJDTLESZIkNcgQJ0mS1CBDnCRJUoPm9d2ApGdm/tLP9d3CRD1w3vF9tyBJs5IzcZIkSQ0yxEmSJDXIECdJktQgQ5wkSVKDDHGSJEkNMsRJkiQ1yBAnSZLUIEOcJElSgwxxkiRJDTLESZIkNcgQJ0mS1CBDnCRJUoMMcZIkSQ0yxEmSJDXIECdJktQgQ5wkSVKDDHGSJEkNMsRJkiQ1yBAnSZLUIEOcJElSgwxxkiRJDTLESZIkNWhe3w1I0iTMX/q5vluYqAfOO77vFiT1zJk4SZKkBhniJEmSGtR8iEtyXJJ7k6xIsrTvfiRJkjaHpkNckq2BDwCvAg4CTklyUL9dSZIkTV7rFzYcAayoqm8AJLkCOAG4q9euJGmW8oIPacuRquq7h6ctyUnAcVX15u75G4AXV9VZ6623BFjSPT0QuHezNtqG3YBv993ELOR+2ZD7ZDT3y2jul9HcLxtyn4y2f1XtPuqF1mfiMqK2QSqtqouBiyffTruSLK+qhX33Mdu4XzbkPhnN/TKa+2U098uG3CdPXdPnxAErgX2Hnu8DrOqpF0mSpM2m9RB3M7AgyQFJ/gmwCLim554kSZImrunDqVW1JslZwBeArYEPV9WdPbfVKg83j+Z+2ZD7ZDT3y2jul9HcLxtynzxFTV/YIEmSNFe1fjhVkiRpTjLESZIkNcgQJ0nPUJJRX3ckrSPJjn33MBsleY5/h54eQ9wc1926TJ0kz02yMMm2ffcymyQ5OMlRSXbtu5fZIslLuy8Yp6rKf4R+Jslrkvx2333MJklOAM5PskffvcwmSV4JfJp1vy5MYzLEzVFJ/hlAVa01yA0keTXwF8B7gUun9tFcl+RVwMeB3wEuS/KcnlvqVZKtkvwc8EHg7CRvhX8McnP+MzXJK4A/wNsf/qMkRwHnA1dX1SN99zNbdH9Wzgf2At7ecztNmvMfOHNRF1ZuTfIxMMgBJPkV4L8Ci6vqGOB7wNJ+u+pfkqOBPwHeXFUnAk8Ch/TYUu+q6qdV9UNgGfAh4FeS/M7Ua70217Pu79HlwJKqui7Jzkn2T7JD37317HDgv3X75J8meXmSFyfZue/G+pLkXwEXAqcCC4DnJ/mX/XbVHkPcHNOdk3EW8DbgySQfBYNc57yq+lq3/G5gFw+r8jBwelXd1M3AvRg4K8kHk5w0xw8hrmFwCGgZcESS9yV5Twbm6mfrd4CfAHt1h94/A1zEYGZ7Lv95WTO0fBXwmww+hz+Q5Nn9tNS7rYE3dt/tuiODe5ofDJ5j+lTM1Q+aOauqfsTgA+RjwO8B2w0HuT5769mNDA6lTp0nuC2wP7BTV5uT54JV1d1V9Vfd09OAC7sZua8CJzO4YfVcdTXwf6vqemA58FZgpxqYkzNyVXUvcDxwAfB1Bp8zrwY+D7wemKuB5QbgLUmuAC6pqlMY/Efxh8ARvXbWk6r6QlX9bZKtqur7wOeAdyf55+UX2I7NEDcHVdWqqvphVX0bOB3YfirIJTksyfP67XDzq6q1VfVo9zTA94HvVtXqJKcCf5hk+94anAWq6tyq+sNu+SPAs5jbJyM/DhyY5C0MAtx5wH5JTu+3rX5V1dcZBLf3VNUl3eHnDzMIcPv1210/quoOBv9pfjFwQFf7BoPZqN17bK13U//hqarPM7hjw6vn+Gz2U9L0bbf0zFXVd7p/dN6b5B4GHyrH9NxWr6pqDfDDJA8meQ/wCuBNVfV4z631JkmG/3ec5PXAnsCq/rrqV1WtSvIg8B+BM6vqs0mOAVb03Frvquouhi5s6P687A481FtT/buWwezbOUm+2dVeyCD8a+DrDC6g+i9z/MjQ2LztlgDoTsx+B/Dyqrq973761J2PsQ1wd/fz2Kq6r9+uZofuHMF/C/wu8OvdDMOclWRfYI+quqV7vtVcPZQ6Svd36TcYzEKd7L2tB0c7gJMYnLJx6Vz/vF1fkiuBf19VD/TdSwsMcaI7sfZK4O1VdVvf/cwWSd4E3Ow/PD+TZBvg5cDfd+c/iQ1nKjXQhbijGJw7eE/f/Wj28u/Q02OIEwBJtquqH/fdx2zih4okaTYzxEmSJDXIqz8kSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJG3RkuyZ5GNJvpHkliRfSfK6zdzDkUluTHJrkruTnNPVX5tk6ebsRdKWw6tTJW2xuu8p+1tgWVX9eVfbH3htVf3peuvO6+7WMYk+7gV+raq+3t2b98DurgaS9LQ5EydpS/Yy4MmpAAdQVd+cCnBJ3pTkk0k+C3wxyS5JPpPktiRfTfKCbr1zkvze1DaS3JFkfve4J8mybsxVSXYY0ccedLec6u7Te9fQ7/+zbvnWocfjSY5KsmOSDye5OcnXkpwwqR0lqT2GOElbsoOBv5thnV8GFlfVy4DfB75WVS8A3glcNsbvOBC4uBvzKPBbI9a5ALg3yaeTnJ5ku/VXqKpDq+pQBvdiXc5gBvE/ADdU1YsY3NP4vUl2HKMnSXOAIU7SnJHkA0m+nuTmofJ1VfXdbvmlwOUAVXUDsGuSnWfY7INV9b+75Y9221hHVf1nYCHwReDfAJ/fSH8LgPcyuC/tT4BXAEuT3Ar8T2A7YL+Z3qekuWFe3w1I0gTdCbx+6klVnZlkNwYzXVN+NLScEdsoYA3r/qd3u/VeX3/9DTdS9ffARUkuAVYn2XX49W6G7UrgLVW1aqif13ufWkmjOBMnaUt2A7BdkjOGaqPOWZvyZeBUgCRHA9+uqkeBB4DDuvphwAFDY/ZL8svd8inA36y/0STHdxdZACwA1gLfX2+1jwAfqar/NVT7AvDvpsYmeeE0vUuaYwxxkrZYNbj8/kTgqCT3J7kJWAa8YyNDzgEWJrkNOA9Y3NU/BezSHdY8A/g/Q2PuBhZ3Y3YBLhqx3TcwOCfuVgaHa0+tqrVTL3ZXzJ4E/ObQxQ0LgT8AtgFuS3JH91ySAL9iRJKetiTzgf9RVYf03YukuceZOEmSpAY5EydJktQgZ+IkSZIaZIiTJElqkCFOkiSpQYY4SZKkBhniJEmSGvT/AYgHVBrFOytEAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def create_missing_values(df: pd.DataFrame, pct_missing: float=0.2) -> pd.DataFrame:\r\n",
    "    \"\"\" Chooses random rows and columns and replaces the values with np.nan until the desired portion of missing data is met. \"\"\"\r\n",
    "\r\n",
    "    n = int(pct_missing*df.shape[0]*df.shape[1])\r\n",
    "\r\n",
    "    # define all possible row/column index combinations\r\n",
    "    all_rc_pairs = np.array(list(itertools.product(range(0,len(df)), range(0,8))))\r\n",
    "\r\n",
    "    # randomly select a subset of size n (no replacement)\r\n",
    "    idx = np.random.choice(len(all_rc_pairs), replace=False, size=n)\r\n",
    "    subset_rc_pairs = all_rc_pairs[idx]\r\n",
    "\r\n",
    "    # replace with np.nan\r\n",
    "    for x, y in subset_rc_pairs:\r\n",
    "        df.iat[x, y] = np.nan\r\n",
    "\r\n",
    "    return df\r\n",
    "\r\n",
    "df_missing = categorical_df_from_dict(dictionary, size=50000)\r\n",
    "\r\n",
    "df_missing = create_missing_values(df_missing, pct_missing=0.20)\r\n",
    "\r\n",
    "missing_perms = (\r\n",
    "    df_missing\r\n",
    "    .groupby(df.columns.tolist(), as_index=False, dropna=False) # experiment with dropna=True/False\r\n",
    "    .size()\r\n",
    "    .query(\"size <= 12\")\r\n",
    ")\r\n",
    "\r\n",
    "sizes_missing = missing_perms.loc[:, \"size\"]\r\n",
    "\r\n",
    "display(missing_perms.groupby('size').count().iloc[:,0])\r\n",
    "\r\n",
    "fig, ax = plot_frequency_hist(sizes_missing)\r\n",
    "\r\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## h) Imagine the code that you wrote for Part d) was to be deployed in an automated system that Mirador's customers could use independently, on potentially large volumes of their own data. Describe how you might deploy the code, and what additional considerations you might have or any changes to the code you might make. Note: it is not necessary to provide another version of the code created for d).\r\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "b3ba2566441a7c06988d0923437866b63cedc61552a5af99d1f4fb67d367b25f"
  },
  "kernelspec": {
   "display_name": "Python 3.8.5 64-bit ('base': conda)",
   "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"
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}