jupyter_Approx_Monte_Carlo_confidence_interval_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d9c72de8-4fb5-486e-9794-f43dbbf67005",
   "metadata": {},
   "source": [
    "# Estimating Confidence Intervals Using Approximate Monte Carlo Simulation Iterations\n",
    "## Conditionally Bias-adjusted RMSE-based Wald Confidence Intervals with Student’s T distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d7e1ce9-6c4c-468b-9e53-dd51a1d0b6e9",
   "metadata": {},
   "source": [
    "This notebook provides an example on how the PPMF files from the 25 Approximate Monte Carlo (AMC) iterations can be used to estimate confidence intervals for a few example queries and geographies. For details on the AMC method and the results of the research, see Ashmead et al., [*An Approximate Monte Carlo Simulation Method for Estimating Uncertainty and Constructing Confidence Intervals for 2020 Census Statistics*](https://github.com/uscensusbureau/AMC_Confidence_Intervals/blob/main/Approx_Monte_Carlo_confidence_interval_paper.pdf).\n",
    "\n",
    "The method used is an RMSE-based confidence interval assuming a Student’s T distribution with 5 degrees of freedom. We use a conditional approach in\n",
    "which a bias-corrected interval was used only if some criteria was met and otherwise the non-bias-corrected versions were utilized. \n",
    "\n",
    "Bias was estimated as the \"true\" count from the PPMF0 for a given query and geography subtracted from the average count from the 25 AMC iterations for the same query and geography.\n",
    "\n",
    "For the bias-corrected interval, the estimated bias was subtracted from the point estimate of the confidence interval.\n",
    "The non-bias-corrected version did not have an adjusted interval.\n",
    "\n",
    "All three of the following criteria had to be met in order to use the bias-corrected version:\n",
    "* The PPMF0 value was greater than 5\n",
    "* The absolute value of the estimated bias divided by the estimated standard deviation of the AMC estimates was greater than or equal to 0.5\n",
    "* The estimated bias was either negative or the PPMF0 value was at least 25\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65e77e66-a718-4159-87ad-a3c6d0fcb5a9",
   "metadata": {},
   "source": [
    "### Pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "18f9f0db-c62d-4e40-ad3c-e219fea23092",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "\n",
    "def ci_wald_t_conditional(df_data, ci_level = 0.9, degrees_freedom = 5):\n",
    "    \"\"\" Calculate confidence intervals using conditionally bias-adjusted RMSE-based Wald CI with Student’s T distribution\n",
    "        Note: as our query of interest is a count, we apply a floor and ceiling function to the endpoints of the confidence intervals respectively\n",
    "        and truncate at zero to account for the outcome space being integer and non-negative\n",
    "    :df_data: a pandas dataframe aggregated to the query x geography level which contains the following columns-\n",
    "         \"PPMF0_COUNT\": the count from the PPMF0\n",
    "         \"MEAN\": the average count from all simulation iterations\n",
    "         \"STD\": the standard deviation of the counts from all the simulation iterations\n",
    "         \"MSE\": the mean squared error, where error was calculated as the count from the PPMF0 subtracted from the count from a simulation run\n",
    "    :ci_level: level for two-sided confidence interval (ex: .9)\n",
    "    :deg_free: degrees of freedom for CI using Student's T distribution\n",
    "    :return: the input pandas dataframe with a new column for the confidence interval\n",
    "    \"\"\"\n",
    "    t_value = scipy.stats.t.ppf(q=1-(1-ci_level)/2, df=degrees_freedom) \n",
    "    print(\"Calculating confidence interval using t-value\", t_value)\n",
    "    df_data[\"RMSE\"] = df_data[\"MSE\"].pow(1/2)\n",
    "    df_data[\"BIAS\"] = df_data[\"MEAN\"] - df_data[\"PPMF0_COUNT\"]\n",
    "    df_data[\"MET_CRITERIA\"] = np.where((df_data[\"PPMF0_COUNT\"] > 5) &\n",
    "                                  (df_data[\"BIAS\"].abs()/df_data[\"STD\"] >= 0.5) &\n",
    "                                  ((df_data[\"BIAS\"] < 0) | (df_data[\"PPMF0_COUNT\"] >= 25)),\n",
    "                                  True, False)\n",
    "    df_data[\"POINT_EST\"] = np.where(df_data[\"MET_CRITERIA\"], df_data[\"PPMF0_COUNT\"] - df_data[\"BIAS\"], df_data[\"PPMF0_COUNT\"])\n",
    "    df_data['CI'] = list(zip(np.floor((df_data[\"POINT_EST\"] - t_value * df_data[\"RMSE\"]).clip(lower=0)).astype(int), \n",
    "                    np.ceil(df_data[\"POINT_EST\"] + t_value * df_data[\"RMSE\"]).astype(int)))\n",
    "    return df_data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b5b75d3-c7ad-43b5-8822-4c0f5b190697",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "For this example, we have already used the 26 2010 PPMF files (the PPMF0 and the TDA outputs PPMF1 - PPMF25) to summarize counts \n",
    "for a few example queries and geographies. Though we used 2010 data, our queries were based on 2020 DHC and PL94 tables. Note the full list of queries in the 2020 DHC and PL94 tables can be found here-\n",
    "\n",
    "DHC: https://www2.census.gov/programs-surveys/decennial/2020/technical-documentation/complete-tech-docs/demographic-and-housing-characteristics-file-and-demographic-profile/2020-dhc-table-matrix.xlsx\n",
    "\n",
    "PL94: https://www2.census.gov/programs-surveys/decennial/2020/technical-documentation/complete-tech-docs/summary-file/2020Census_PL94_171Redistricting_StatesTechDoc_English.pdf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3adbb555-44e6-4715-9619-f1f3908c15a7",
   "metadata": {},
   "outputs": [
    {
     "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>PRODUCT_TABLE</th>\n",
       "      <th>QUERY</th>\n",
       "      <th>GEO</th>\n",
       "      <th>ITERATION</th>\n",
       "      <th>COUNT</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>0</td>\n",
       "      <td>1379</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1</td>\n",
       "      <td>1382</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>2</td>\n",
       "      <td>1375</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>3</td>\n",
       "      <td>1382</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>4</td>\n",
       "      <td>1373</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>5</td>\n",
       "      <td>1382</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>6</td>\n",
       "      <td>1389</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>7</td>\n",
       "      <td>1376</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>8</td>\n",
       "      <td>1380</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>9</td>\n",
       "      <td>1379</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  PRODUCT_TABLE            QUERY                      GEO  ITERATION  COUNT\n",
       "0       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200          0   1379\n",
       "1       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200          1   1382\n",
       "2       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200          2   1375\n",
       "3       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200          3   1382\n",
       "4       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200          4   1373\n",
       "5       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200          5   1382\n",
       "6       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200          6   1389\n",
       "7       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200          7   1376\n",
       "8       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200          8   1380\n",
       "9       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200          9   1379"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Our data has one row for each query, geography, and iteration (0-25)\n",
    "\n",
    "file_path = \"jupyter_data.csv\"\n",
    "df_example = pd.read_csv(file_path)\n",
    "df_example.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3beb1276-d890-4a13-8990-8d3f393feeb0",
   "metadata": {},
   "outputs": [
    {
     "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>PRODUCT_TABLE</th>\n",
       "      <th>QUERY</th>\n",
       "      <th>GEO</th>\n",
       "      <th>PPMF0_COUNT</th>\n",
       "      <th>ITERATION</th>\n",
       "      <th>COUNT</th>\n",
       "      <th>DIFF</th>\n",
       "      <th>ABS_DIFF</th>\n",
       "      <th>SQRD_ERROR</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>1</td>\n",
       "      <td>1382</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>2</td>\n",
       "      <td>1375</td>\n",
       "      <td>-4</td>\n",
       "      <td>4</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>3</td>\n",
       "      <td>1382</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>4</td>\n",
       "      <td>1373</td>\n",
       "      <td>-6</td>\n",
       "      <td>6</td>\n",
       "      <td>36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>5</td>\n",
       "      <td>1382</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>6</td>\n",
       "      <td>1389</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>7</td>\n",
       "      <td>1376</td>\n",
       "      <td>-3</td>\n",
       "      <td>3</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>8</td>\n",
       "      <td>1380</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>9</td>\n",
       "      <td>1379</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>10</td>\n",
       "      <td>1376</td>\n",
       "      <td>-3</td>\n",
       "      <td>3</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  PRODUCT_TABLE            QUERY                      GEO  PPMF0_COUNT  \\\n",
       "0       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200         1379   \n",
       "1       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200         1379   \n",
       "2       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200         1379   \n",
       "3       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200         1379   \n",
       "4       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200         1379   \n",
       "5       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200         1379   \n",
       "6       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200         1379   \n",
       "7       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200         1379   \n",
       "8       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200         1379   \n",
       "9       DHCH H4  H0040001: Total  ST-CO-TR: 24-033-807200         1379   \n",
       "\n",
       "   ITERATION  COUNT  DIFF  ABS_DIFF  SQRD_ERROR  \n",
       "0          1   1382     3         3           9  \n",
       "1          2   1375    -4         4          16  \n",
       "2          3   1382     3         3           9  \n",
       "3          4   1373    -6         6          36  \n",
       "4          5   1382     3         3           9  \n",
       "5          6   1389    10        10         100  \n",
       "6          7   1376    -3         3           9  \n",
       "7          8   1380     1         1           1  \n",
       "8          9   1379     0         0           0  \n",
       "9         10   1376    -3         3           9  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Divide the csv datafile into two - the PPMF0 data (which was input to the TDA for the AMC method) \n",
    "# and the PPMF1-PPMF25 synthetic data (which was the output from the AMC method)\n",
    "\n",
    "df_PPMF0 = df_example[df_example[\"ITERATION\"] == 0]\n",
    "df_PPMF1_25 = df_example[df_example[\"ITERATION\"] != 0]\n",
    "\n",
    "df_PPMF0 = df_PPMF0.rename(columns={\"COUNT\": \"PPMF0_COUNT\"})\n",
    "df_PPMF0 = df_PPMF0.drop(columns=[\"ITERATION\"])\n",
    "\n",
    "# Then join the files back together to calculate the differences\n",
    "\n",
    "df_calc = df_PPMF0.merge(df_PPMF1_25, on=[\"PRODUCT_TABLE\", \"QUERY\", \"GEO\"], how=\"outer\")\n",
    "df_calc[\"DIFF\"] = df_calc[\"COUNT\"] - df_calc[\"PPMF0_COUNT\"]\n",
    "df_calc[\"ABS_DIFF\"] = df_calc[\"DIFF\"].abs()\n",
    "df_calc[\"SQRD_ERROR\"] = df_calc[\"DIFF\"]**2\n",
    "\n",
    "df_calc.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fb771c2b-3776-443c-b4cf-cda0e52a9943",
   "metadata": {},
   "outputs": [
    {
     "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>PRODUCT_TABLE</th>\n",
       "      <th>QUERY</th>\n",
       "      <th>GEO</th>\n",
       "      <th>PPMF0_COUNT</th>\n",
       "      <th>MEAN</th>\n",
       "      <th>MSE</th>\n",
       "      <th>STD</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>1379.68</td>\n",
       "      <td>14.92</td>\n",
       "      <td>3.880722</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040002: Owned mortgage</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>106</td>\n",
       "      <td>106.28</td>\n",
       "      <td>5.40</td>\n",
       "      <td>2.354428</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040003: Owned no mortgage</td>\n",
       "      <td>ST-CO-TR-BG: 24-033-807200-1</td>\n",
       "      <td>17</td>\n",
       "      <td>14.20</td>\n",
       "      <td>20.08</td>\n",
       "      <td>3.570714</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040004: Rented</td>\n",
       "      <td>ST-CO-TR-BG: 24-033-807200-1</td>\n",
       "      <td>383</td>\n",
       "      <td>384.88</td>\n",
       "      <td>20.84</td>\n",
       "      <td>4.245782</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080001: Total NonVoting</td>\n",
       "      <td>STATE: 32</td>\n",
       "      <td>665025</td>\n",
       "      <td>665023.52</td>\n",
       "      <td>377.64</td>\n",
       "      <td>19.776080</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080002: NonVoting Household</td>\n",
       "      <td>STATE: 32</td>\n",
       "      <td>663773</td>\n",
       "      <td>663761.92</td>\n",
       "      <td>615.72</td>\n",
       "      <td>22.660391</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080003: Nonvoting Householder, spouse</td>\n",
       "      <td>STATE: 32</td>\n",
       "      <td>744</td>\n",
       "      <td>795.76</td>\n",
       "      <td>3226.56</td>\n",
       "      <td>23.880396</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080021: Nonvoting Institutional</td>\n",
       "      <td>STATE-COUNTY: 32-021</td>\n",
       "      <td>13</td>\n",
       "      <td>11.72</td>\n",
       "      <td>32.08</td>\n",
       "      <td>5.631163</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080029: Nonvoting Non-Institutional</td>\n",
       "      <td>STATE-COUNTY: 32-021</td>\n",
       "      <td>0</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>PL94 H1</td>\n",
       "      <td>H0010001: Total</td>\n",
       "      <td>ST-AIAN: 06-4255</td>\n",
       "      <td>1313</td>\n",
       "      <td>1313.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  PRODUCT_TABLE                                      QUERY  \\\n",
       "0       DHCH H4                            H0040001: Total   \n",
       "1       DHCH H4                   H0040002: Owned mortgage   \n",
       "2       DHCH H4                H0040003: Owned no mortgage   \n",
       "3       DHCH H4                           H0040004: Rented   \n",
       "4     DHCP PCT8                PCT0080001: Total NonVoting   \n",
       "5     DHCP PCT8            PCT0080002: NonVoting Household   \n",
       "6     DHCP PCT8  PCT0080003: Nonvoting Householder, spouse   \n",
       "7     DHCP PCT8       PCT0080021: Nonvoting Institutional    \n",
       "8     DHCP PCT8    PCT0080029: Nonvoting Non-Institutional   \n",
       "9       PL94 H1                            H0010001: Total   \n",
       "\n",
       "                            GEO  PPMF0_COUNT       MEAN      MSE        STD  \n",
       "0       ST-CO-TR: 24-033-807200         1379    1379.68    14.92   3.880722  \n",
       "1       ST-CO-TR: 24-033-807200          106     106.28     5.40   2.354428  \n",
       "2  ST-CO-TR-BG: 24-033-807200-1           17      14.20    20.08   3.570714  \n",
       "3  ST-CO-TR-BG: 24-033-807200-1          383     384.88    20.84   4.245782  \n",
       "4                     STATE: 32       665025  665023.52   377.64  19.776080  \n",
       "5                     STATE: 32       663773  663761.92   615.72  22.660391  \n",
       "6                     STATE: 32          744     795.76  3226.56  23.880396  \n",
       "7          STATE-COUNTY: 32-021           13      11.72    32.08   5.631163  \n",
       "8          STATE-COUNTY: 32-021            0       0.00     0.00   0.000000  \n",
       "9              ST-AIAN: 06-4255         1313    1313.00     0.00   0.000000  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For each query and unique geography, calculate the input statistics needed to construct the confidence intervals\n",
    "# this requires calculating the mean and standard deviation from the 25 simulation iterations, and the mean squared error of the differences\n",
    "\n",
    "# note this syntax only works for pandas >= 0.25\n",
    "df_statistics = df_calc.groupby([\"PRODUCT_TABLE\", \"QUERY\", \"GEO\", \"PPMF0_COUNT\"],  as_index=False).agg(\n",
    "    MEAN=(\"COUNT\", \"mean\"),\n",
    "    MSE=(\"SQRD_ERROR\", \"mean\"),\n",
    "    STD=(\"COUNT\", \"std\"))\n",
    "\n",
    "\n",
    "df_statistics.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "24a5fd5b-fe03-4bab-9df1-b51d4467c1bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating confidence interval using t-value 2.015048372669157\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>PRODUCT_TABLE</th>\n",
       "      <th>QUERY</th>\n",
       "      <th>GEO</th>\n",
       "      <th>PPMF0_COUNT</th>\n",
       "      <th>MEAN</th>\n",
       "      <th>MSE</th>\n",
       "      <th>STD</th>\n",
       "      <th>RMSE</th>\n",
       "      <th>BIAS</th>\n",
       "      <th>MET_CRITERIA</th>\n",
       "      <th>POINT_EST</th>\n",
       "      <th>CI</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>1379.68</td>\n",
       "      <td>14.92</td>\n",
       "      <td>3.880722</td>\n",
       "      <td>3.862642</td>\n",
       "      <td>0.68</td>\n",
       "      <td>False</td>\n",
       "      <td>1379.00</td>\n",
       "      <td>(1371, 1387)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040002: Owned mortgage</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>106</td>\n",
       "      <td>106.28</td>\n",
       "      <td>5.40</td>\n",
       "      <td>2.354428</td>\n",
       "      <td>2.323790</td>\n",
       "      <td>0.28</td>\n",
       "      <td>False</td>\n",
       "      <td>106.00</td>\n",
       "      <td>(101, 111)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040003: Owned no mortgage</td>\n",
       "      <td>ST-CO-TR-BG: 24-033-807200-1</td>\n",
       "      <td>17</td>\n",
       "      <td>14.20</td>\n",
       "      <td>20.08</td>\n",
       "      <td>3.570714</td>\n",
       "      <td>4.481071</td>\n",
       "      <td>-2.80</td>\n",
       "      <td>True</td>\n",
       "      <td>19.80</td>\n",
       "      <td>(10, 29)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040004: Rented</td>\n",
       "      <td>ST-CO-TR-BG: 24-033-807200-1</td>\n",
       "      <td>383</td>\n",
       "      <td>384.88</td>\n",
       "      <td>20.84</td>\n",
       "      <td>4.245782</td>\n",
       "      <td>4.565085</td>\n",
       "      <td>1.88</td>\n",
       "      <td>False</td>\n",
       "      <td>383.00</td>\n",
       "      <td>(373, 393)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080001: Total NonVoting</td>\n",
       "      <td>STATE: 32</td>\n",
       "      <td>665025</td>\n",
       "      <td>665023.52</td>\n",
       "      <td>377.64</td>\n",
       "      <td>19.776080</td>\n",
       "      <td>19.432962</td>\n",
       "      <td>-1.48</td>\n",
       "      <td>False</td>\n",
       "      <td>665025.00</td>\n",
       "      <td>(664985, 665065)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080002: NonVoting Household</td>\n",
       "      <td>STATE: 32</td>\n",
       "      <td>663773</td>\n",
       "      <td>663761.92</td>\n",
       "      <td>615.72</td>\n",
       "      <td>22.660391</td>\n",
       "      <td>24.813706</td>\n",
       "      <td>-11.08</td>\n",
       "      <td>False</td>\n",
       "      <td>663773.00</td>\n",
       "      <td>(663722, 663824)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080003: Nonvoting Householder, spouse</td>\n",
       "      <td>STATE: 32</td>\n",
       "      <td>744</td>\n",
       "      <td>795.76</td>\n",
       "      <td>3226.56</td>\n",
       "      <td>23.880396</td>\n",
       "      <td>56.802817</td>\n",
       "      <td>51.76</td>\n",
       "      <td>True</td>\n",
       "      <td>692.24</td>\n",
       "      <td>(577, 807)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080021: Nonvoting Institutional</td>\n",
       "      <td>STATE-COUNTY: 32-021</td>\n",
       "      <td>13</td>\n",
       "      <td>11.72</td>\n",
       "      <td>32.08</td>\n",
       "      <td>5.631163</td>\n",
       "      <td>5.663921</td>\n",
       "      <td>-1.28</td>\n",
       "      <td>False</td>\n",
       "      <td>13.00</td>\n",
       "      <td>(1, 25)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080029: Nonvoting Non-Institutional</td>\n",
       "      <td>STATE-COUNTY: 32-021</td>\n",
       "      <td>0</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.00</td>\n",
       "      <td>False</td>\n",
       "      <td>0.00</td>\n",
       "      <td>(0, 0)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>PL94 H1</td>\n",
       "      <td>H0010001: Total</td>\n",
       "      <td>ST-AIAN: 06-4255</td>\n",
       "      <td>1313</td>\n",
       "      <td>1313.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.00</td>\n",
       "      <td>False</td>\n",
       "      <td>1313.00</td>\n",
       "      <td>(1313, 1313)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>PL94 H1</td>\n",
       "      <td>H0010002: Occupied</td>\n",
       "      <td>ST-AIAN: 06-4255</td>\n",
       "      <td>1220</td>\n",
       "      <td>1219.20</td>\n",
       "      <td>29.20</td>\n",
       "      <td>5.454356</td>\n",
       "      <td>5.403702</td>\n",
       "      <td>-0.80</td>\n",
       "      <td>False</td>\n",
       "      <td>1220.00</td>\n",
       "      <td>(1209, 1231)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>PL94 H1</td>\n",
       "      <td>H0010003: Vacant</td>\n",
       "      <td>ST-AIAN: 06-4255</td>\n",
       "      <td>93</td>\n",
       "      <td>93.80</td>\n",
       "      <td>29.20</td>\n",
       "      <td>5.454356</td>\n",
       "      <td>5.403702</td>\n",
       "      <td>0.80</td>\n",
       "      <td>False</td>\n",
       "      <td>93.00</td>\n",
       "      <td>(82, 104)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>PL94 P1</td>\n",
       "      <td>P0010003: White alone</td>\n",
       "      <td>ST-ELEM: 17-42570</td>\n",
       "      <td>3536</td>\n",
       "      <td>3536.00</td>\n",
       "      <td>120.88</td>\n",
       "      <td>11.221260</td>\n",
       "      <td>10.994544</td>\n",
       "      <td>0.00</td>\n",
       "      <td>False</td>\n",
       "      <td>3536.00</td>\n",
       "      <td>(3513, 3559)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>PL94 P1</td>\n",
       "      <td>P0010004: Black alone</td>\n",
       "      <td>ST-ELEM: 17-42570</td>\n",
       "      <td>187</td>\n",
       "      <td>189.36</td>\n",
       "      <td>40.04</td>\n",
       "      <td>5.992217</td>\n",
       "      <td>6.327717</td>\n",
       "      <td>2.36</td>\n",
       "      <td>False</td>\n",
       "      <td>187.00</td>\n",
       "      <td>(174, 200)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>PL94 P1</td>\n",
       "      <td>P0010005: AIAN alone</td>\n",
       "      <td>ST-ELEM: 17-42570</td>\n",
       "      <td>8</td>\n",
       "      <td>7.44</td>\n",
       "      <td>9.68</td>\n",
       "      <td>3.123566</td>\n",
       "      <td>3.111270</td>\n",
       "      <td>-0.56</td>\n",
       "      <td>False</td>\n",
       "      <td>8.00</td>\n",
       "      <td>(1, 15)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   PRODUCT_TABLE                                      QUERY  \\\n",
       "0        DHCH H4                            H0040001: Total   \n",
       "1        DHCH H4                   H0040002: Owned mortgage   \n",
       "2        DHCH H4                H0040003: Owned no mortgage   \n",
       "3        DHCH H4                           H0040004: Rented   \n",
       "4      DHCP PCT8                PCT0080001: Total NonVoting   \n",
       "5      DHCP PCT8            PCT0080002: NonVoting Household   \n",
       "6      DHCP PCT8  PCT0080003: Nonvoting Householder, spouse   \n",
       "7      DHCP PCT8       PCT0080021: Nonvoting Institutional    \n",
       "8      DHCP PCT8    PCT0080029: Nonvoting Non-Institutional   \n",
       "9        PL94 H1                            H0010001: Total   \n",
       "10       PL94 H1                         H0010002: Occupied   \n",
       "11       PL94 H1                           H0010003: Vacant   \n",
       "12       PL94 P1                      P0010003: White alone   \n",
       "13       PL94 P1                      P0010004: Black alone   \n",
       "14       PL94 P1                       P0010005: AIAN alone   \n",
       "\n",
       "                             GEO  PPMF0_COUNT       MEAN      MSE        STD  \\\n",
       "0        ST-CO-TR: 24-033-807200         1379    1379.68    14.92   3.880722   \n",
       "1        ST-CO-TR: 24-033-807200          106     106.28     5.40   2.354428   \n",
       "2   ST-CO-TR-BG: 24-033-807200-1           17      14.20    20.08   3.570714   \n",
       "3   ST-CO-TR-BG: 24-033-807200-1          383     384.88    20.84   4.245782   \n",
       "4                      STATE: 32       665025  665023.52   377.64  19.776080   \n",
       "5                      STATE: 32       663773  663761.92   615.72  22.660391   \n",
       "6                      STATE: 32          744     795.76  3226.56  23.880396   \n",
       "7           STATE-COUNTY: 32-021           13      11.72    32.08   5.631163   \n",
       "8           STATE-COUNTY: 32-021            0       0.00     0.00   0.000000   \n",
       "9               ST-AIAN: 06-4255         1313    1313.00     0.00   0.000000   \n",
       "10              ST-AIAN: 06-4255         1220    1219.20    29.20   5.454356   \n",
       "11              ST-AIAN: 06-4255           93      93.80    29.20   5.454356   \n",
       "12             ST-ELEM: 17-42570         3536    3536.00   120.88  11.221260   \n",
       "13             ST-ELEM: 17-42570          187     189.36    40.04   5.992217   \n",
       "14             ST-ELEM: 17-42570            8       7.44     9.68   3.123566   \n",
       "\n",
       "         RMSE   BIAS  MET_CRITERIA  POINT_EST                CI  \n",
       "0    3.862642   0.68         False    1379.00      (1371, 1387)  \n",
       "1    2.323790   0.28         False     106.00        (101, 111)  \n",
       "2    4.481071  -2.80          True      19.80          (10, 29)  \n",
       "3    4.565085   1.88         False     383.00        (373, 393)  \n",
       "4   19.432962  -1.48         False  665025.00  (664985, 665065)  \n",
       "5   24.813706 -11.08         False  663773.00  (663722, 663824)  \n",
       "6   56.802817  51.76          True     692.24        (577, 807)  \n",
       "7    5.663921  -1.28         False      13.00           (1, 25)  \n",
       "8    0.000000   0.00         False       0.00            (0, 0)  \n",
       "9    0.000000   0.00         False    1313.00      (1313, 1313)  \n",
       "10   5.403702  -0.80         False    1220.00      (1209, 1231)  \n",
       "11   5.403702   0.80         False      93.00         (82, 104)  \n",
       "12  10.994544   0.00         False    3536.00      (3513, 3559)  \n",
       "13   6.327717   2.36         False     187.00        (174, 200)  \n",
       "14   3.111270  -0.56         False       8.00           (1, 15)  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now calculate the 90% confidence interval using the aggregated pandas dataframe.\n",
    "# We use 5 degrees of freedom, preferring a more conservative, larger interval\n",
    "\n",
    "df_ci = ci_wald_t_conditional(df_data=df_statistics, ci_level=.9, degrees_freedom=5)\n",
    "df_ci.head(15)\n",
    "\n",
    "# Please note the following in our output:\n",
    "# Query \"PCT0080029: Nonvoting Non-Institutional\" had a count of 0 from every AMC iteration and from the PPMF0. \n",
    "#    Therefore, the confidence interval is (0, 0)\n",
    "# Query \"H0010001: Total\" measures the total number of housing units. This is an invariant, \n",
    "#    meaning the DAS TDA does not apply noise to this statistic. \n",
    "#    Therefore, every AMC simulation contains the same count and the confidence interval has a width of 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32a5fa73-1f89-4c8e-a42b-b5b5605ab040",
   "metadata": {},
   "source": [
    "### PySpark\n",
    "Alternatively, when users are accessing all 26 of the full PPMF files, they may prefer to calculate the prior example using PySpark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a69fef03-a458-404b-b473-ba5a29a5a73a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This was run and tested on EMR 6.15.0\n",
    "\n",
    "import os, sys\n",
    "import scipy.stats\n",
    "from IPython.core.display import HTML\n",
    "\n",
    "if 'SPARK_HOME' not in os.environ:\n",
    "    os.environ['SPARK_HOME'] = '/usr/lib/spark'\n",
    "    \n",
    "sys.path.append(os.path.join(os.environ['SPARK_HOME'],'python'))\n",
    "sys.path.append(os.path.join(os.environ['SPARK_HOME'],'python','lib', 'py4j-src.zip'))\n",
    "\n",
    "from pyspark.sql import SparkSession\n",
    "import pyspark.sql.functions as F\n",
    "\n",
    "def ci_wald_t_conditional_pyspark(df_data, ci_level = 0.9, degrees_freedom = 5):\n",
    "    \"\"\" Calculate confidence intervals using conditionally bias-adjusted RMSE-based Wald CI with Student's T distribution\n",
    "        Note: as our query of interest is a count, we apply a floor and ceiling function to the endpoints of the confidence intervals respectively\n",
    "        and truncate at zero to account for the outcome space being integer and non-negative\n",
    "    :df: a pyspark dataframe aggregated to the query x geography level which contains the following columns-\n",
    "         \"PPMF0_COUNT\": the \"true\" count from the PPMF0\n",
    "         \"MEAN\": the average count from all simulation runs\n",
    "         \"STD\": the standard deviation of the counts from all the simulation runs\n",
    "         \"MSE\": the mean squared error, where error was calculated as the count from the PPMF0 subtracted from the count from a simulation run\n",
    "    :ci_level: level for two-sided confidence interval (ex: .9)\n",
    "    :deg_free: degrees of freedom for CI using Student's T distribution\n",
    "    :return: the input pyspark dataframe with a new column for the confidence interval\n",
    "    \"\"\"\n",
    "    t_value = scipy.stats.t.ppf(q=1-(1-ci_level)/2, df=degrees_freedom) \n",
    "    print(\"Calculating confidence interval using t-value\", t_value)\n",
    "    df_data = df_data.withColumn(\"RMSE\", F.sqrt(F.col(\"MSE\")))\\\n",
    "                     .withColumn(\"BIAS\", F.col(\"MEAN\") - F.col(\"PPMF0_COUNT\"))\\\n",
    "                     .withColumn(\"MET_CRITERIA\", F.when((F.col(\"PPMF0_COUNT\") > 5) &\n",
    "                                              (F.abs(F.col(\"BIAS\"))/F.col(\"STD\") >= 0.5) &\n",
    "                                              ((F.col(\"BIAS\") < 0) | (F.col(\"PPMF0_COUNT\") >= 25)),\n",
    "                                              True).otherwise(False))\\\n",
    "                     .withColumn(\"POINT_EST\", F.when(F.col(\"MET_CRITERIA\"), F.col(\"PPMF0_COUNT\") - F.col(\"BIAS\"))\n",
    "                                     .otherwise(F.col(\"PPMF0_COUNT\")))\\\n",
    "                     .withColumn(\"CI\", F.array(F.floor(F.greatest(F.lit(0), F.col(\"POINT_EST\") - t_value * F.col(\"RMSE\"))),\n",
    "                                     F.ceil(F.col(\"POINT_EST\") + t_value * F.col(\"RMSE\"))))\n",
    "    return df_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9ff58e27-a912-46e8-9600-83cdccdbe8ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Setting default log level to \"WARN\".\n",
      "To adjust logging level use sc.setLogLevel(newLevel). For SparkR, use setLogLevel(newLevel).\n",
      "24/04/11 14:28:07 WARN Client: Neither spark.yarn.jars nor spark.yarn.archive is set, falling back to uploading libraries under SPARK_HOME.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating confidence interval using t-value 2.015048372669157\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "                                                                                \r"
     ]
    },
    {
     "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>PRODUCT_TABLE</th>\n",
       "      <th>QUERY</th>\n",
       "      <th>GEO</th>\n",
       "      <th>PPMF0_COUNT</th>\n",
       "      <th>MEAN</th>\n",
       "      <th>MSE</th>\n",
       "      <th>STD</th>\n",
       "      <th>RMSE</th>\n",
       "      <th>BIAS</th>\n",
       "      <th>MET_CRITERIA</th>\n",
       "      <th>POINT_EST</th>\n",
       "      <th>CI</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040001: Total</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>1379</td>\n",
       "      <td>1379.68</td>\n",
       "      <td>14.92</td>\n",
       "      <td>3.880722</td>\n",
       "      <td>3.862642</td>\n",
       "      <td>0.68</td>\n",
       "      <td>False</td>\n",
       "      <td>1379.00</td>\n",
       "      <td>[1371, 1387]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040002: Owned mortgage</td>\n",
       "      <td>ST-CO-TR: 24-033-807200</td>\n",
       "      <td>106</td>\n",
       "      <td>106.28</td>\n",
       "      <td>5.40</td>\n",
       "      <td>2.354428</td>\n",
       "      <td>2.323790</td>\n",
       "      <td>0.28</td>\n",
       "      <td>False</td>\n",
       "      <td>106.00</td>\n",
       "      <td>[101, 111]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040003: Owned no mortgage</td>\n",
       "      <td>ST-CO-TR-BG: 24-033-807200-1</td>\n",
       "      <td>17</td>\n",
       "      <td>14.20</td>\n",
       "      <td>20.08</td>\n",
       "      <td>3.570714</td>\n",
       "      <td>4.481071</td>\n",
       "      <td>-2.80</td>\n",
       "      <td>True</td>\n",
       "      <td>19.80</td>\n",
       "      <td>[10, 29]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>DHCH H4</td>\n",
       "      <td>H0040004: Rented</td>\n",
       "      <td>ST-CO-TR-BG: 24-033-807200-1</td>\n",
       "      <td>383</td>\n",
       "      <td>384.88</td>\n",
       "      <td>20.84</td>\n",
       "      <td>4.245782</td>\n",
       "      <td>4.565085</td>\n",
       "      <td>1.88</td>\n",
       "      <td>False</td>\n",
       "      <td>383.00</td>\n",
       "      <td>[373, 393]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080001: Total NonVoting</td>\n",
       "      <td>STATE: 32</td>\n",
       "      <td>665025</td>\n",
       "      <td>665023.52</td>\n",
       "      <td>377.64</td>\n",
       "      <td>19.776080</td>\n",
       "      <td>19.432962</td>\n",
       "      <td>-1.48</td>\n",
       "      <td>False</td>\n",
       "      <td>665025.00</td>\n",
       "      <td>[664985, 665065]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080002: NonVoting Household</td>\n",
       "      <td>STATE: 32</td>\n",
       "      <td>663773</td>\n",
       "      <td>663761.92</td>\n",
       "      <td>615.72</td>\n",
       "      <td>22.660391</td>\n",
       "      <td>24.813706</td>\n",
       "      <td>-11.08</td>\n",
       "      <td>False</td>\n",
       "      <td>663773.00</td>\n",
       "      <td>[663722, 663824]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080003: Nonvoting Householder, spouse</td>\n",
       "      <td>STATE: 32</td>\n",
       "      <td>744</td>\n",
       "      <td>795.76</td>\n",
       "      <td>3226.56</td>\n",
       "      <td>23.880396</td>\n",
       "      <td>56.802817</td>\n",
       "      <td>51.76</td>\n",
       "      <td>True</td>\n",
       "      <td>692.24</td>\n",
       "      <td>[577, 807]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080021: Nonvoting Institutional</td>\n",
       "      <td>STATE-COUNTY: 32-021</td>\n",
       "      <td>13</td>\n",
       "      <td>11.72</td>\n",
       "      <td>32.08</td>\n",
       "      <td>5.631163</td>\n",
       "      <td>5.663921</td>\n",
       "      <td>-1.28</td>\n",
       "      <td>False</td>\n",
       "      <td>13.00</td>\n",
       "      <td>[1, 25]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>DHCP PCT8</td>\n",
       "      <td>PCT0080029: Nonvoting Non-Institutional</td>\n",
       "      <td>STATE-COUNTY: 32-021</td>\n",
       "      <td>0</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.00</td>\n",
       "      <td>False</td>\n",
       "      <td>0.00</td>\n",
       "      <td>[0, 0]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>PL94 H1</td>\n",
       "      <td>H0010001: Total</td>\n",
       "      <td>ST-AIAN: 06-4255</td>\n",
       "      <td>1313</td>\n",
       "      <td>1313.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.00</td>\n",
       "      <td>False</td>\n",
       "      <td>1313.00</td>\n",
       "      <td>[1313, 1313]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>PL94 H1</td>\n",
       "      <td>H0010002: Occupied</td>\n",
       "      <td>ST-AIAN: 06-4255</td>\n",
       "      <td>1220</td>\n",
       "      <td>1219.20</td>\n",
       "      <td>29.20</td>\n",
       "      <td>5.454356</td>\n",
       "      <td>5.403702</td>\n",
       "      <td>-0.80</td>\n",
       "      <td>False</td>\n",
       "      <td>1220.00</td>\n",
       "      <td>[1209, 1231]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>PL94 H1</td>\n",
       "      <td>H0010003: Vacant</td>\n",
       "      <td>ST-AIAN: 06-4255</td>\n",
       "      <td>93</td>\n",
       "      <td>93.80</td>\n",
       "      <td>29.20</td>\n",
       "      <td>5.454356</td>\n",
       "      <td>5.403702</td>\n",
       "      <td>0.80</td>\n",
       "      <td>False</td>\n",
       "      <td>93.00</td>\n",
       "      <td>[82, 104]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>PL94 P1</td>\n",
       "      <td>P0010003: White alone</td>\n",
       "      <td>ST-ELEM: 17-42570</td>\n",
       "      <td>3536</td>\n",
       "      <td>3536.00</td>\n",
       "      <td>120.88</td>\n",
       "      <td>11.221260</td>\n",
       "      <td>10.994544</td>\n",
       "      <td>0.00</td>\n",
       "      <td>False</td>\n",
       "      <td>3536.00</td>\n",
       "      <td>[3513, 3559]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>PL94 P1</td>\n",
       "      <td>P0010004: Black alone</td>\n",
       "      <td>ST-ELEM: 17-42570</td>\n",
       "      <td>187</td>\n",
       "      <td>189.36</td>\n",
       "      <td>40.04</td>\n",
       "      <td>5.992217</td>\n",
       "      <td>6.327717</td>\n",
       "      <td>2.36</td>\n",
       "      <td>False</td>\n",
       "      <td>187.00</td>\n",
       "      <td>[174, 200]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>PL94 P1</td>\n",
       "      <td>P0010005: AIAN alone</td>\n",
       "      <td>ST-ELEM: 17-42570</td>\n",
       "      <td>8</td>\n",
       "      <td>7.44</td>\n",
       "      <td>9.68</td>\n",
       "      <td>3.123566</td>\n",
       "      <td>3.111270</td>\n",
       "      <td>-0.56</td>\n",
       "      <td>False</td>\n",
       "      <td>8.00</td>\n",
       "      <td>[1, 15]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   PRODUCT_TABLE                                      QUERY  \\\n",
       "0        DHCH H4                            H0040001: Total   \n",
       "1        DHCH H4                   H0040002: Owned mortgage   \n",
       "2        DHCH H4                H0040003: Owned no mortgage   \n",
       "3        DHCH H4                           H0040004: Rented   \n",
       "4      DHCP PCT8                PCT0080001: Total NonVoting   \n",
       "5      DHCP PCT8            PCT0080002: NonVoting Household   \n",
       "6      DHCP PCT8  PCT0080003: Nonvoting Householder, spouse   \n",
       "7      DHCP PCT8       PCT0080021: Nonvoting Institutional    \n",
       "8      DHCP PCT8    PCT0080029: Nonvoting Non-Institutional   \n",
       "9        PL94 H1                            H0010001: Total   \n",
       "10       PL94 H1                         H0010002: Occupied   \n",
       "11       PL94 H1                           H0010003: Vacant   \n",
       "12       PL94 P1                      P0010003: White alone   \n",
       "13       PL94 P1                      P0010004: Black alone   \n",
       "14       PL94 P1                       P0010005: AIAN alone   \n",
       "\n",
       "                             GEO  PPMF0_COUNT       MEAN      MSE        STD  \\\n",
       "0        ST-CO-TR: 24-033-807200         1379    1379.68    14.92   3.880722   \n",
       "1        ST-CO-TR: 24-033-807200          106     106.28     5.40   2.354428   \n",
       "2   ST-CO-TR-BG: 24-033-807200-1           17      14.20    20.08   3.570714   \n",
       "3   ST-CO-TR-BG: 24-033-807200-1          383     384.88    20.84   4.245782   \n",
       "4                      STATE: 32       665025  665023.52   377.64  19.776080   \n",
       "5                      STATE: 32       663773  663761.92   615.72  22.660391   \n",
       "6                      STATE: 32          744     795.76  3226.56  23.880396   \n",
       "7           STATE-COUNTY: 32-021           13      11.72    32.08   5.631163   \n",
       "8           STATE-COUNTY: 32-021            0       0.00     0.00   0.000000   \n",
       "9               ST-AIAN: 06-4255         1313    1313.00     0.00   0.000000   \n",
       "10              ST-AIAN: 06-4255         1220    1219.20    29.20   5.454356   \n",
       "11              ST-AIAN: 06-4255           93      93.80    29.20   5.454356   \n",
       "12             ST-ELEM: 17-42570         3536    3536.00   120.88  11.221260   \n",
       "13             ST-ELEM: 17-42570          187     189.36    40.04   5.992217   \n",
       "14             ST-ELEM: 17-42570            8       7.44     9.68   3.123566   \n",
       "\n",
       "         RMSE   BIAS  MET_CRITERIA  POINT_EST                CI  \n",
       "0    3.862642   0.68         False    1379.00      [1371, 1387]  \n",
       "1    2.323790   0.28         False     106.00        [101, 111]  \n",
       "2    4.481071  -2.80          True      19.80          [10, 29]  \n",
       "3    4.565085   1.88         False     383.00        [373, 393]  \n",
       "4   19.432962  -1.48         False  665025.00  [664985, 665065]  \n",
       "5   24.813706 -11.08         False  663773.00  [663722, 663824]  \n",
       "6   56.802817  51.76          True     692.24        [577, 807]  \n",
       "7    5.663921  -1.28         False      13.00           [1, 25]  \n",
       "8    0.000000   0.00         False       0.00            [0, 0]  \n",
       "9    0.000000   0.00         False    1313.00      [1313, 1313]  \n",
       "10   5.403702  -0.80         False    1220.00      [1209, 1231]  \n",
       "11   5.403702   0.80         False      93.00         [82, 104]  \n",
       "12  10.994544   0.00         False    3536.00      [3513, 3559]  \n",
       "13   6.327717   2.36         False     187.00        [174, 200]  \n",
       "14   3.111270  -0.56         False       8.00           [1, 15]  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spark = SparkSession.builder.getOrCreate()\n",
    "\n",
    "df_spark = spark.createDataFrame(df_example)\n",
    "\n",
    "# We split the datafile into two - the PPMF0 \"true\" data and the PPMF1-PPMF25 synthetic data\n",
    "df_PPMF0_spark = df_spark.filter(\"ITERATION = 0\")\n",
    "df_PPMF1_25_spark = df_spark.filter(\"ITERATION != 0\")\n",
    "\n",
    "df_PPMF0_spark = df_PPMF0_spark.withColumnRenamed(\"COUNT\", \"PPMF0_COUNT\")\n",
    "df_PPMF0_spark = df_PPMF0_spark.drop(\"ITERATION\")\n",
    "\n",
    "# Then join the files back together & calculate the differences\n",
    "df_calc_spark = df_PPMF1_25_spark.join(df_PPMF0_spark, [\"PRODUCT_TABLE\", \"QUERY\", \"GEO\"], \"outer\")\n",
    "df_calc_spark = df_calc_spark.withColumn(\"DIFF\", F.col(\"COUNT\") - F.col(\"PPMF0_COUNT\"))\n",
    "df_calc_spark = df_calc_spark.withColumn(\"ABS_DIFF\", F.abs(F.col(\"DIFF\")))\n",
    "df_calc_spark = df_calc_spark.withColumn(\"SQRD_ERROR\", F.pow(F.col(\"DIFF\"), 2))\n",
    "\n",
    "# For each query level and unique geography, calculate the statistics needed to construct the confidence intervals\n",
    "# this includes calculating the mean, stdev, and the mean squared error of the differences\n",
    "df_statistics_spark = df_calc_spark.groupBy([\"PRODUCT_TABLE\", \"QUERY\", \"GEO\", \"PPMF0_COUNT\"]).agg(\n",
    "                    F.mean(\"COUNT\").alias(\"MEAN\"),\n",
    "                    F.mean(\"SQRD_ERROR\").alias(\"MSE\"),\n",
    "                    F.stddev(\"COUNT\").alias(\"STD\"))\n",
    "\n",
    "df_ci_spark = ci_wald_t_conditional_pyspark(df_data=df_statistics_spark, ci_level=.9, degrees_freedom=5)\n",
    "\n",
    "df_ci_spark.orderBy(\"PRODUCT_TABLE\",\"QUERY\").toPandas().head(15)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}