examples/census_income_covshift.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "237cb9c0-9b68-42ac-aa97-3d3f5e9b4bdc",
   "metadata": {},
   "source": [
    "# Relationship between age and income in a covariate-shifted population\n",
    "\n",
    "The goal is to investigate the relationship between age and income using US census data. The target of inference is the linear regression coefficient when regressing yearly income in dollars on age, while controlling for sex. The data from California in the year 2019 is downloaded through the Folktables interface (1). A shift in the distribution of the covariates is introduced between the gold-standard and unlabeled data sets by randomly sampling the unlabeled data set with sampling weights 0.8 for females and 0.2 for males. Predictions of income are made by training a gradient boosting tree via XGBoost (2) on the previous year’s data.\n",
    "\n",
    "1. F. Ding, M. Hardt, J. Miller, L. Schmidt, “Retiring adult: New datasets for fair machine learning” in Advances in Neural Information Processing Systems 34 (2021), pp. 6478–6490.\n",
    "2. T. Chen, C. Guestrin, “XGBoost: A scalable tree boosting system” in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (2016), pp. 785–794."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "598bdc10-1a86-4258-84d8-8d2d4f0a3662",
   "metadata": {},
   "source": [
    "### Import necessary packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3c78cb6c-791d-499f-b8e9-57a17dd7ec11",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "\n",
    "sys.path.append(os.path.abspath(os.path.join(os.getcwd(), os.pardir)))\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from ppi_py.datasets import load_dataset\n",
    "from ppi_py import ppi_ols_ci, classical_ols_ci\n",
    "from statsmodels.regression.linear_model import OLS\n",
    "from scipy.optimize import brentq\n",
    "from tqdm import tqdm\n",
    "from utils import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cf90ae6",
   "metadata": {},
   "source": [
    "### Import the census income data set\n",
    "\n",
    "Load the data. The data set contains reported income (```Y```), predicted income (```Yhat```), and age and sex (```X```)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a6da3138",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_folder = \"./data/\"\n",
    "data = load_dataset(dataset_folder, \"census_income\")\n",
    "Y_total = data[\"Y\"]\n",
    "Yhat_total = data[\"Yhat\"]\n",
    "X_total = data[\"X\"][:, 0:1]  # Age\n",
    "Z_total = data[\"X\"][:, 1:2] - 1  # Sex"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8969f9db",
   "metadata": {},
   "source": [
    "### Problem setup\n",
    "\n",
    "Specify the error level (```alpha```), range of values for the labeled data set size (```ns```), and number of trials (```num_trials```)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5b3c8f29",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "coordinate = 0\n",
    "weights = np.array([0.8, 0.2])  # Women, men in the new population\n",
    "inD_weights = np.array(\n",
    "    [1 - Z_total.mean(), Z_total.mean()]\n",
    ")  # Women, men in the old population\n",
    "n_total = Y_total.shape[0]  # Total number of labeled examples\n",
    "n_inD = 10000\n",
    "ns = np.linspace(100, 2000, 10).astype(\n",
    "    int\n",
    ")  # Test for different numbers of labeled ballots\n",
    "num_trials = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e9580e6-0ced-499f-92b0-f5a3f717f865",
   "metadata": {},
   "source": [
    "Generate \"in-distribution\" and \"out-of-distribution\" data sets via weighted sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c1960205",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate indexes for inD and ooD datasets.\n",
    "randperm = np.random.permutation(n_total)\n",
    "inD_idxs, temp_idxs = randperm[:n_inD], randperm[n_inD:]\n",
    "probability_vector = np.array(\n",
    "    [weights[0] if z == 0 else weights[1] for z in Z_total[temp_idxs]]\n",
    ")\n",
    "probability_vector /= probability_vector.sum()\n",
    "size = int(\n",
    "    (n_total - n_inD) * min(weights[0] / weights[1], weights[1] / weights[0])\n",
    ")\n",
    "shifted_idxs = np.random.choice(\n",
    "    temp_idxs, size=size, replace=False, p=probability_vector\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "72590739",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate inD and ooD datasets\n",
    "X_inD, X_ooD = X_total[inD_idxs], X_total[shifted_idxs]\n",
    "Y_inD, Y_ooD = Y_total[inD_idxs], Y_total[shifted_idxs]\n",
    "Yhat_inD, Yhat_ooD = Yhat_total[inD_idxs], Yhat_total[shifted_idxs]\n",
    "Z_inD, Z_ooD = Z_total[inD_idxs], Z_total[shifted_idxs]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c7c3ba4-b8f2-4094-a81a-48b5b794e508",
   "metadata": {},
   "source": [
    "Compute the ground-truth value of the estimand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cfe10ab3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute ground truth\n",
    "true_theta = OLS(Y_ooD, exog=X_ooD).fit().params[coordinate]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83ce18be",
   "metadata": {},
   "source": [
    "### Construct intervals\n",
    "\n",
    "Form confidence intervals for all methods and problem parameters. A dataframe with the following columns is formed:\n",
    "1. ```method``` (one of ```PPI```, ```Classical```, and ```Imputation```)\n",
    "2. ```n``` (labeled data set size, takes values in ```ns```)\n",
    "3. ```lower``` (lower endpoint of the confidence interval)\n",
    "4. ```upper``` (upper endpoint of the confidence interval)\n",
    "5. ```trial``` (index of trial, goes from ```0``` to ```num_trials-1```)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "812f8fd5",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████| 10/10 [05:45<00:00, 34.58s/it]\n"
     ]
    }
   ],
   "source": [
    "# Run prediction-powered inference and classical inference for many values of n\n",
    "results = []\n",
    "for i in tqdm(range(ns.shape[0])):\n",
    "    for j in range(num_trials):\n",
    "        # Prediction-Powered Inference\n",
    "        n = ns[i]\n",
    "        rand_idx = np.random.permutation(n_inD)\n",
    "        _X = X_inD[rand_idx[:n]]\n",
    "        _Y = Y_inD[rand_idx[:n]]\n",
    "        _Z = Z_inD[rand_idx[:n]]\n",
    "        _Yhat = Yhat_inD[rand_idx[:n]]\n",
    "        importance_weights = np.array(\n",
    "            [\n",
    "                (\n",
    "                    weights[0] / inD_weights[0]\n",
    "                    if z == 0\n",
    "                    else weights[1] / inD_weights[1]\n",
    "                )\n",
    "                for z in _Z\n",
    "            ]\n",
    "        )\n",
    "        ppi_ci = ppi_ols_ci(\n",
    "            _X, _Y, _Yhat, X_ooD, Yhat_ooD, w=importance_weights, alpha=alpha\n",
    "        )\n",
    "\n",
    "        # Classical interval\n",
    "        classical_ci = classical_ols_ci(\n",
    "            _X, _Y, w=importance_weights, alpha=alpha\n",
    "        )\n",
    "\n",
    "        # Append results\n",
    "        results += [\n",
    "            pd.DataFrame(\n",
    "                [\n",
    "                    {\n",
    "                        \"method\": \"PPI\",\n",
    "                        \"n\": n,\n",
    "                        \"lower\": ppi_ci[0][coordinate],\n",
    "                        \"upper\": ppi_ci[1][coordinate],\n",
    "                        \"trial\": j,\n",
    "                    }\n",
    "                ]\n",
    "            )\n",
    "        ]\n",
    "        results += [\n",
    "            pd.DataFrame(\n",
    "                [\n",
    "                    {\n",
    "                        \"method\": \"Classical\",\n",
    "                        \"n\": n,\n",
    "                        \"lower\": classical_ci[0][coordinate],\n",
    "                        \"upper\": classical_ci[1][coordinate],\n",
    "                        \"trial\": j,\n",
    "                    }\n",
    "                ]\n",
    "            )\n",
    "        ]\n",
    "\n",
    "# Imputed CI\n",
    "imputed_ci = classical_ols_ci(X_ooD, Yhat_ooD, alpha=alpha)\n",
    "results += [\n",
    "    pd.DataFrame(\n",
    "        [\n",
    "            {\n",
    "                \"method\": \"Imputation\",\n",
    "                \"n\": np.nan,\n",
    "                \"lower\": imputed_ci[0][coordinate],\n",
    "                \"upper\": imputed_ci[1][coordinate],\n",
    "                \"trial\": 0,\n",
    "            }\n",
    "        ]\n",
    "    )\n",
    "]\n",
    "\n",
    "df = pd.concat(results, axis=0, ignore_index=True)\n",
    "df[\"width\"] = df[\"upper\"] - df[\"lower\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d15ba288",
   "metadata": {},
   "source": [
    "### Plot results\n",
    "\n",
    "Plot:\n",
    "1. Five randomly chosen intervals from the dataframe for PPI and the classical method, and the imputed interval;\n",
    "2. The average interval width for PPI and the classical method, together with a scatterplot of the widths from the five random draws."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6077b2c4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnsAAADMCAYAAAALZGVzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyVUlEQVR4nO3deXSc13nn+e9TexVQhR0kwc2iqIVaKEuiJFoiFcqxLG+xY59JnDiZdOLYSjyx4+7O6aS9TTqTuHs603Fixx3nSHYmdrrtsdtLd+wTLbZFSiRlWaIWU5tFkaLAnVgIoAqo9a2680cViwAJEiCBQhUKv885dXDr1lsvngviRT289733mnMOEREREWlOvnoHICIiIiK1o2RPREREpIkp2RMRERFpYkr2RERERJqYkj0RERGRJqZkT0RERKSJBWZ4XeuySKOzegdQY7oGpZE19fX3tre9zT344IP1DkPkQmZ1Dc6U7ImIyCUqPP9AtRy8/u11jEQuxdDQUL1DEJkXSvZERGqkuP/xalnJnojUi+7ZExEREWli6tkTEamRwHVvq3cIUgOpZJKhgQFy2SzhSITu3l7iiUS9wxI5LyV7IiI1ErjijnqHIPMslUxypL8ffyBAKBymUChwpL+fVWvXKuGThqVhXBERkVkaGhjAHwgQDAYxM4LBIP5AgKGBgXqHJnJeSvZERERmKZfNEghMHRQLBALkstk6RSQyMyV7IiIisxSORPA8b0qd53mEI5E6RSQyM92zJyJSI4Vn/1e1HLzxPXWMROZLd28vR/r7gXKPnud5FD2PFStX1jkykfNTz56ISI0UX99TfUhziCcSrFq7lmAwSD6XIxgManKGNDz17ImIiFyEVDDFK/FXGYmM0BHsIByMEkfJnjQuJXsiIjUSeOO76x2CzLNjmaPsGNpOxB+lPdBOuphmx9B2tnXfRV9UQ7nSmJTsiYjUSOCyW+odgsyzvcm9tORbiI1HsTwEQi34W33sTe5VsicNS/fsiYiIzFIqmaRlOIoVgSBYEVqGo6SSyXqHJnJeSvZERERmqT3dhucrgB8wwA+er0B7uq3eoYmcl5I9ERGRWeqik4J5FEoFnINCqUDBPLrorHdoIuele/ZERGok/9S3quXQLb9ax0iag5ndBvxn59w2M3sj8LdAEcgBv+WcO2lmXwDuAFKVt70HyAP/Deit1P8r59zgpcTQ0dKJZX0MF4fIFLNE/RGWB1fQHmmfU9tEakk9eyIiNVI68nz1IXNjZn8MfBk4vVXF54GPOee2Ad8F/qRSfxNwj3NuW+UxBnwEeN45txX4GvDpS42ju7eXiAuzJrSWa+PXsia0logL093be6mnFKk5JXsiIrIYHADeN+n5rznnnquUA0DWzHzAFcB9ZrbbzD5YeX0L8GCl/ADwlksNQosqy2KkYVwREeDJUz+d/5Pe+fYz5Vqcv0Zu7byt3iGcwzn3HTN7w6TnxwHM7Hbgo8CdQAvlod3PUZ5Csd3M9gAJYKzy1hRw3tkUZnYvcC/AmjVrpj0mnkgouZNFRcmeiAjQGeqkJ6yhOIBsMUvEH5n5wDozs/cDnwLe6ZwbNDM/8HnnXLry+iPADUASiFfeFgdGz3dO59x9wH0AmzZtcrWLXmThaBhXREQWHTP7Tco9etucc69Vqq8EdpmZ38yClIdvnwF2A++oHPN2YOdCxytST+rZExGRRaXSg/cF4BDwXTMDeNQ596dm9t+BJ4AC8DXn3ItmdhD4qpntojwz9wN1Cl2kLpTsiYjIouCcex3YXHk67cJ2zrm/BP7yrLo08Cs1DU6kgWkYV0QEMP05FJEmpZ49ERHg8tbL5/2c2e99plqOvPfP5/38IiKzof/KioiIiDQx9eyJiNRIcPNv1DsEEREleyKyuD0z+jQpLzXzgfUQnFQeOnHBQ3+he1tNQxGRpUvJnogsahFflMsS6+odxpxlvAzRQLTeYYhIE9I9eyIiIiJNTMmeiIiISBNTsiciIiLSxHTPnogsahF/pN4hyBJzLHOUvcm9jBRG6Ah2sDGxkb7oynqHJXJeSvZEZFFb19K4kzO0qHLzOZY5yo6h7UT8UdoD7aSLaXYMbWdb911K+KRhKdkTEamR4JYP1jsEmWd7k3uJ+KPE/DGA6te9yb1K9qRhKdkTEakRf89l9Q5B5tlIYYT2QPuUuqgvykhhpD4BicyCJmiIiIjMUkewg0wpM6UuU8rQEeyoU0QiM1PPnojIWXYO72SsMFrvMBZcW7CdrV1b6x1GQ9uY2MiOoe1AuUcvU8qQLWbY3LG5zpGJnJ+SPRGRs8T8Ma5uvbreYSy4Q5lD9Q6h4fVFV7Kt+64ps3E3d2zW/XrS0JTsiYiIXIS+6Eold7Ko6J49ERERkSamZE9ERESkiSnZExE5y+m105aapdpukWane/ZERM6yIb6h3iHURU+4p94hiEgNqGdPREQWBTO7zcx2VMrrzWyXme00sy+Zma9S/2Ez22NmT5jZuyp1UTP7TuXYfzEzZbWypCjZExGRhmdmfwx8GYhUqj4HfNo5txUw4D1mthz4Q+AO4B7gP5lZGPgI8Hzl2K8Bn17o+EXqScO4IrJkfOvINxnIDdQ7jLr56OUfq3cIc3EAeB/wT5XnNwOPVsoPAG8FisBu51wOyJnZfmAjsAX4y0nHfmahghZpBEr2RGTJSAQTbO58U73DqJvR/CjtofZ6h3FJnHPfMbM3TKoy55yrlFNAG5AAxiYdM1396bppmdm9wL0Aa9asmZfYRepNw7giIrIYlSaV48AokKyUL1R/um5azrn7nHObnHObenp0a580ByV7IiKyGD1rZtsq5bcDO4Enga1mFjGzNmAD8AKwG3jHWceKLBkaxhURkcXoj4D7zSwEvAx82zlXNLMvUE7mfMCnnHNZM/sS8FUz2wXkgQ/ULWqROlCyJyJLRiKQqHcIMgfOudeBzZXyPuAXpjnmfuD+s+rSwK8sQIgiDUnJnogsGbd33VHvEKQJFItFxkZGiCcSBEOheocjMiMleyIiIhfhpb17GR4cpDUe56bbbsPn0+3v0tj0GyoiInIRMuk0AOOpFKOnTtU5GpGZqWdPRJaET7747xf+m3q5M+VAeF5PvTK6kj9Yt6gXSV60Oru7SU9MAHDy+HE6u7vrHJHIhSnZE5El4brE9dzScWu9w5hXw7lhusJd9Q5jyVm2YgVH+vsBGBoYoFgs4vf76xyVyPlpGFdEROQitMbjxFpagPJkjaGBpbsFnywOSvZEREQugpnRu2JF9fnJ48frGI3IzJTsiYiIXKRly5dXy6eGh8nn83WMRuTClOyJyJIQMN2iLPMnGouRaG8vP3GOwRMn6hqPyIXor5+ILAm/uur99Q5BmsyyFStIjo4C5aHclWvW1DcgkfNQz56IiMgl6Fm2DDMDIDk2Vl1/T6TRqGdPRKRGst//82o58kufqWMkUguhUIjO7m6GBweBcu/eGy6/vM5RiZxLPXsiIrXi5c88pCktO2tWrnOujtGITE89eyIiIhchlUwyNDBALpslGArh8/spFYtk0mlSySSJtrZ6hygyhZI9EZEaCb/r0/UOoSGZ2YeBfw1EAQOcc25dXYOapVQyyZH+fvyBAKFwGM/zCIVCZDMZoNy7p2RPGo2GcUVEasSC4epDpvh94B3ABuDqytdFYWhgAH8gQDAYxMwIBoO0tLZWXx84cYJSqVTHCEXOpWRPREQW2pBzrt85lzv9qHdAs5XLZgkEpg6KxVpa8PnKH6eFfJ7RU6fqEZrIeWkYV0REFoSZ/cdKMWRmDwHPAA7AOffJugV2EcKRCPl8Hi/oESGCz3wUi0Va4nFSY2NAeSi3s7u7zpGKnKFkT0REFsorZ309bdFMYe3u7eWb/d/gpeLLrPOt4257C6VikVVr1vDy888DMDgwwJWehz+gj1hpDPpNFBGpkez3zqytF3nvn1/gyKXBOfdVADP7onPuo6frzexrwNfqFthFiCcSDAQGwYPXSq9xMnSSG1fdTGs8Tv9rr5GemKBULDI0ODhlWRaRetI9eyIisiDM7A/M7DjwYTM7VnkcB1bWO7aLcVXi6mp5f+QA8UQCMztnzT2RRqFkT0REFoRz7r8651YAf+ac66s8VjjnfvFSzmdmv21mOyqPJ8wsa2Y3mdnRSfXvrxz7YTPbUznuXXNpx5aurdXyi8kXOZUfBqB3UrJ3aniYfG7RzDuRJmczrPa9aO6jkCXL6h1AjekalEZ2Udefmf3W+V5zzs1pGNfM/ivwM6AEtDnn/mrSa8uBHwKbgAiwC9g00yzgTZs2uT179kz72t8f/BL7xsu3Ht7V/WZ+acW7AXj2yScZGx0FYP3VV7NqzZq5NEtkJrO6BtWzJyIiC2VD5fHbwPuB1cD7KuVLZmabgGudc/cBNwPvNLPHzOwrZhYHbgV2V5Z5GQP2Axvn8j23Turde2LkCfKl8pZ4vbMcyk0lkxzcv5+fv/ACB/fvJ5VMziUckQtSsiciIgvCOfcJ59wngLxz7p3Ouc86534ZCM7x1J8E/qxSfhL4d865O4HXgD8FEsDYpONTwLTbXJjZvZXh3j2Dg4Pn/YYb4tfQFeoCIFNM8/To0wD0LluGWbmzJTU2RjqdPue9p3fhKBQKhMJhCoUCR/r7lfBJzSjZExGRhdZrZu0AZtYFdF3qiSrnudo5t71S9T3n3NOny8CNQBKIT3pbHBid7nzOufucc5ucc5t6enrO+3195uOOzi3V57uGH8M5RzAUmrLG3sA0vXvT7cLhDwQYGhiYucEil0BLr4hI0xoaGKBQKNQ7jIYRDAbp7u2tdxgAnwX2mFmScq/bB+dwrjuBH016/pCZfcw59yTwi8DTlHv7PmtmESBMeSj5hTl8TwBu67yNB08+QN7lOZ49zoGJA6xvXc+yFSsYrvQKnjx+nLXr1lV7+6C8C0coPHULvUAgQC6bnWtIItNSsiciTcsfCBDXpvRTeJ53znZfC8059x3gO2bWCww754pzON1VlIdrT/sI8EUzywMngHudc0kz+wKwk/KI1qecc3POrKL+GDd3bOInpx4HYOfwY6xvXU9XTw9+v59isUgmnSaVTJKY9HsYjkQoFAoEg2dGrz3PIxyJzDUkkWkp2RMRkQVxejFlM3sayE2qxzl3+6Wc0zn3/5z1/BngnHM55+4H7r+U73EhW7u2VpO9F5LPM5IfoSPUQc+yZZw4dgwo9+5NTva6e3s50t8PlHv0PM+j6HmsWLmolhuURUT37ImIyEI5vY3IZUA/8F3gw8Cv1y2iOVoeWcEVLVcA4HDsPrULYMoCywPHj1MqlarP44kEq9auJRgMks/lCAaDrFq7lngisbDBy5Khnj0REVkQzrmTla+dZrYBeDdwHzAAvLeesc3Flq6tvDrxKgA/PfUEb+29h/bOTkLhMPlcjkKhwMjwMF2TJnzEEwkld7Jg1LMnIk3L7/fXOwSZhpndALwTeHOl6uU6hjNn1yauoyPYAcBEcYLnRp/FzOhdvrx6zMkTJ+oVnoh69kSkeSU0OaNRPUZ5UsWnnHP/Uu9g5spnPu7o2sIPTnwfKE/UuKXjVpatWFG9N2/o5Em8DRuqk2OOZY6yN7mXkcIIHcEONiY20hfVPXtSG+rZExGRhdYF/Btgq5n92My+Ue+ALsaxzFEePPkA3zjydR48+QDHMke5rWMzQSvPrj2aPcrB9EFa43FiLS0AlEql6nIsxzJH2TG0nXQxTXugnXQxzY6h7RzLHK1bm6S5KdkTEakRl0lWHzJFO7ASWAvEKE/WWBTOl6iNFUa5qf3m6nG7hh/DzKZM1DhZmZ27N7mXiD9KzB/DzIj5Y0T8UfYm9y54e2Rp0DCuiEiN5B48sypI5L1/foEjl5yHKO9u8Vnn3Iv1DuZiTE7UgOrXvcm9bO26k5+OPFF+PraX0cIoy1as4OD+/QCcOnWKfC7HSGGE9kD7lPNGfVFGCiML1xBZUpTsiUjTOdzfIB1F1//mmfICxhQKhab0KDUa59zNMx/VmC6UqPVF+1jXcjmvTRygRInHh3fzjuXvpK29nbHRUXCOgRMn6Ih0kC6mq4kiQKaUqU7yEJlvSvZEpOnEWlpoaW2tdxh1dfYODTI/OoIXTtS2dt3JaxMHAHji1E+4u/etLOvrKyd7lBdY3njDRnYdfozSRJCQFyQfKOC15Nm4evOCt0eWBt2zJyIiMksbExvJFjOki2mcc6SLabLFDBsTGwG4LnEd7cF2AMaL4zw39iw9y5ZV98ZNJZMEkwEuT68j5EJkfDlCLsTl6XXEC/F6NUuanJI9ERGRWeqLrmRb913E/DFGvVFi/hjbuu+qLpviNz+3d95RPX7n8E4CgcCUBZWP9PcTDye4PLGe69qu4/LEeuLhBEMDAwveHlkaNIwrIiJyEfqiKy+4Jt7mzjfx8MBDeM7jSOYw/enX6V2+vJrMjadSJNrbp7wnEAiQy2ZrGbYsYerZExERmUetgVZuar+p+nzn8E66enrwVxZULhaLZNLpKe/xPI9wJLKgccrSoZ49EWk6Xd3d9Q4BgFJqqFr2xRsjJlkYW7ru5MmRJwH42dhzvGfFe+hZtowTR8sLJ0+MjxMKhwkEAnieR9HzWLFSO2hIbSjZExGpkfyPPl8ta529pWVVdBWXxS7jYPpgeRmWU49z24rbqsleLp+lP9dPIZkjGAlzed964olEnaOWZqVhXBERkRrY0rW1Wv7JqcdpbYsTCocBKHklstEstjZIqmeCJ7JPaLs0qRn17ImIVBzp7yeXy83fCa/61TPlffvm77wLJBwOs2rt2nqHsWhtbLuBxPE2kt4YKS/F3uTPWLZiBYdffx2A2FiEQrubsgvHhSZ+iFwqJXsiIhXZTIZ4WxvxtrZ6h9Iw8vk8oVCo3mEsSn7zc0fXHTxw8l+A8kSND634UDXZ8yUNig782i5NakvDuCIiIjWyufNN+M0PwKFMP8O+U5QiDgArWTnhQ9ulSW0p2RMREamReCDOjW03Vp/vPrWT3uXLqs/9pyDtTd2FQ2S+KdkTERGpoS1dd1bLz449y4oVq6rP/eM+2g61siWxRffrSc0o2RMRqYhEo1DZw1RkvqyJrWFttDzRpeiKPJ1+mp7ly6uv+1PGoWcOMnDiRL1ClCanCRoiIhXzPfM0+73PVMtaZ682zOxZYKzy9CDwWeAfAQe8APyBc65kZh8Gfg/wgL9wzv1gIePc0n0n/Yf/CYDHTz3OJ6/5FKFQiKOHDgHgFQq8tHcvQwMDXLFhA8FgcMZzppJJhgYGyGWzhCMRunt7tVafTEs9eyIisiiZWQTAObet8vgd4HPAp51zWwED3mNmy4E/BO4A7gH+k5mFFzLWGxI3EA+UE7GkN8aL4y9wxdVXc8PNN0/ZJm3gxAmeevxxhgcHL3i+VDLJwVdfZXhoiFQyyfDQEAdffZVUMlnTdsjipGRPRKRGrL2v+pCauAGImdnDZvaImW0GbgYerbz+APAW4FZgt3Mu55wbA/YDCzobIuAL8KbON1Wf7xzeCUBHVxe3vOlNLO878zuSz+V4/tlneeWll/A8b9rzHT18mEwmgzmH3+/HnCOTyXD08OHaNkQWJQ3jiojUSPiuj9Q7hGaXBv4L8GXgCsrJnTnnXOX1FNAGJDgz1Du5/hxmdi9wL8CaNWvmNdjbO2/nx4M/ouiKvJ4+yJHMYVZFVxMIBrn6uuvo7u3llZdeopDPA3D8yBFGhoe5+rrraO+YuixLamwMv9+Pz19e1sX8fvyVepGzKdkTkbp59qmn6h3CgorFYlx17bX1DqOZ7AP2V5K7fWY2TLln77Q4MAokK+Wz68/hnLsPuA9g06ZNbrpjLlUi2MYNiTfyzNjTQLl379dXfaD6endvL4n2dva99BJDAwNAeaHv5556itVr1/KG9evxV5I7KN+UOCX2+QxWmoqSPRGpm/bOTrp6euodxoLKZrNEJt2jJXPyQeB64P8wsz7KPXgPm9k259wO4O3AduBJ4LOVe/zCwAbKkzcW3JburdVk79nRZ/il5e+mNdBafT0UCnHtDTdw8vhxXv35zylWhnEP9/czPDTEhuuvJ55IkGhrY3RkhKIZPp+PUqlEqVg8pwdQBHTPnoiILF5fAdrNbBfwTcrJ38eBPzOznwAh4NvOuRPAF4CdwCPAp5xz2XoEvDa6ltXR1QB4zuOJUz855xgzY3lfH7fcfjsdXV3V+vTEBM/89Ke8fuAAK1atwsI+xgpjDKcHGSuMYWEffatXL1hbZPGwM7c2TEu9wtLomn1RtKa+Bg8eOLDkevZCoVAz9ew19fW3adMmt2fPnnk/71MjT/KNI18HIEiQWwO3sG3Fm+lq7z7nWOccxw4f5sC+fZRKpWp9uDXKwfjrWAYCnh8vUKSQ8HjL2rdqcealZVbXoHr2REREFtB633qiRAEoUGC39zifO/xXPHzkQfKl/JRjzYyVa9aw6U1vItF2Zk5JbjzD8hM9FP0lTvYOcaJzkKN2jN3Duxe0LbI4KNkTkbox7VYhS9Do4Cnu8b+VuJ2ZM5Ihw4MjD/IfX/ksu4Z34pWmLrkSa2nhxltv5bIrrqheNz7nY9XwCtYdXUv3eCdhi9CfeX0hmyKLhIZxZbFr9mxB1+Ailtv5lWo5vPV36xhJzTT19VerYdyfv/ACoXCYEiVeLr7MU94ext34lGPag+3c3ftWbmm/lYBv6lzK8VSKx55+hFg+OqU+7y8wkBjk/Td8gGAoNO9xS0Oa1TWoZE8Wu6b+sEHXoDS2pr7+apXsHdy/n0KhUN0SzXMez+ef59nSs0yQnnJsZ7CTt/bew80dm/DbmWVXvnXomxSOZVmRXIbvrEE6n8/Hsr4+Vq1ZQ0trK7WmbdvqSsmeLAlN/WGDrkFpbE19/dUq2Uslkxzp78cfCBAIBPA8j6LnsWx1H897z/PjwR8x7k3t6esOdfPW3nu4qf1mfObjWOYoP3r9YcJjIVozLUS8MH7nP+d7dXR1sWrNGjq7u2ty28T52rJq7VolfAtDyZ4sCU39YYOuwRntfOSReofQ0FrjcW685ZZanb6pr79aJXtw4d6wXCnH7uFdbB98hInixJT39YZ7uaf3bayzdew/+AojpREyLkuUCK35FqxoZNLpc75fNBZj1dq1LF+xAn9g/pbYPbh/P6PZUYaLQ2SKWaL+CF3+btoj7Vy2fv28fR85LyV7siQ09YcNugZn9PILL0zZV1TOFYlEiMZitTh1U19/tUz2ZiNbzLJreCfbh7aTKZ41vGudbLKbuTJ0ZbXHrlAoEAgE6Ozq4sihQ9VdOCYLBAKsWLWKlatXE4lGz3n9Yj3zsyc56h0l4AsQsCCeK+CVPFYGVnLTDbfO+fwyo1ldg9pBQ0REpAFF/BHe0ns3d3Rt4bGhR3l0aAfZUnkt6FPuFA+7H/JM7lluDtzEOv86AoEA+VyO9s5O2js7yaTTHD18mONHjlAsFgHwPI/Dr7/O4f5+enp7WbV2LYm2tkse4h3mFEEXIOAr338YtCBWLNdL41CyJyIi0sCi/ij3LHsbW7vv5NHBHTw2/Ci5Ug6AITfEQ4WHCRVCXG7ruCZ4Dc45zIxoLMb6q67iDZdfzomjRzly6BDZTKZ8UucYPHmSwZMniScSrFi1iq6eHsLh8EXFNhobo2OkDYqUF3MrQaAUZCQ2Nr8/BJkTDePKYtfUw0joGpyRhnFnpmHcS1PLYdxjmaPsTe5lpDBCR7CDjYmNs975Ytwb5+GjD/JE8gk8vHNe7wp1san9FjZ13EJX6Mx2a845hgcHOXLoEKOnpu95a2ltpXf5crp6emhpbZ2xx+/Bkw+QG88SG49ieXAhSLdmCLdGeNuyt8+qPTInumdPloSm/rBB1+CMXnnpJXqXL693GA1Nyd6lqVWydyxzlB1D24n4o0R9UTKlDNlihm3dd13UVmfHR46x88SjvOS9TJLktMesi61jU8et3NB2A1H/mXv0xlMpDu7fz/Dg4HnPH4lG6erpobunh7aODny+c/dhmK+2yCVTsidLQlN/2KBrUBpbU19/tUr2Hjz5AOlimpj/TAJ++vml9IY55+hPv85To0/x3OizZEqZc44JWJDrE9ezqeMWrmq9Cp/5OLh/P9lslnw2SzqdJpfNnvd7BAIBOru76e7tpbOri0BljUCYWy+lzJmSPVkSmvrDBl2D0tia+vqrVbL3jSNfpz3QPmWI1DnHqDfKr6/6wJzOXSgVeDH1IntGnuTnqZ9TonTOMYlAgpvab6Z7qIMO6ySbTuN5Hj6fD/P5yGezFAqF6qSOs5kZ7R0ddPX20tXTQ3QeZvXKJVOyJ0tCU3/YoGtQGltTX3+N3rM3U49aqpDimbGn2TPyFEezR6c9R7vXztrcalbmVhBzMTCjs7ubKzdsYHRkhOGBAYYGBy/Y69fS2kp3by/tnZ3E4/EpvX5Sc0r2ZElo6g8bdA1KY2vq66+R79m72HMcyxzjqdEneWb0aVJeatpzxr04Pflu1oUv483Xv5WIPwKUex0nxscZqiR+48np7w88LRKNEk8kaE0kiMfjtCYShBZgr94lum2bkj1ZEpr6wwZdg4ta8fDPqmX/6hvqGEnNNPX116izceHSeweLrsgrqVfYM/oUe0d/RsnOHeYF8OFjbWwtV7ReyVWtV7Emtra6N282m2V4cJDhwUFGhoeZIY8AIByJ0BqPT0kCQ+HwvG3htoS3bdOiyiIi9VTY8+1quUmTPblEfdGVc5rEMFIYoT3QPqUu6osyUhi54Pv85ueaxDVck7iGHYd/yAH/a5wIn2TEP0LJziRtJUocTB/kYPogDw88RNgX5vKW9VzZeiVXtl5FX2UXjtGREfa//DLZXI5SsUipNH3ymMtmyVWSxNOCoVA5+asM/eYyGTzPIxKNXnSv3NDAAP5AgGBlCPn016GBgWZP9mZFyZ6IiMgi0xHsOKdnL1PK0BHsmPU5vHCR1alVrMguJ295UoFxRgOjDEdOkfRNHerNlXK8lHqRl1IvApAItHFl65W0jrWQcC3EI62Yz0epWKTgeeWevNZWUskk4+PjuGmSwEI+z6mhIU4NDU2p9/l8HDt8mI7ubtra2ojGYkRiMSKRyLTLvwDVewpHUyk8zyMQCBCNxcidZ5LJ+TTrULCSPRGRGvGtur7eIUiT2pjYyI6h7QBT7tnb3LF51uc41nqCeCZGsBgkXAwTcAFaSy20tXXwvvX/G6+O72Pf+D72jb9yTo9h0htjz+hT5ScJaCu10VvqoaPUTps/QZvXxk3XlPfGLZVKpCcmGE8mSaVSjCeTjKdS553tWyqVKJVKDBw/zsDx42deMCuvGRmNEonFiEajRGMxorEYzjmGhgcplAo4V8LMx3h6nN6eZbP+eUweCg6FwxQKBY709zfFULCSPRGRGgnd8qv1DkGaVF90Jdu675py39/mjs0XNTScDKY4lRimPdNOyAuSDxQYjY4SCIZoDbRyY/tN3Nh+UzmRyg+xb/wV9o3vY//4q+es5TfmG2PMV9kirTIXY/srO1kZWUlfpI++6EpW9qxkWV8fZoZzjkw6TSqZ5JWfv4TnFTBn2IVuQXOObCZT3vJtmh1AHA4DzHw4HHkvz8jESHX7uJkMDQxQKpXITuodDIXDTTEUrGRPRERkEZrrfX8GnPQPcCoxQsACeM6jUCqwilVTjzOjJ9xDT7iHO7q2UHIlDmcO8+r4Pp47+Qwn3MlpJ3oM54cYzg+xN3lmolLUHysnf5G+ciLYvpKsP4fPB76AvzzZowTOK2EYfb2ryGQyZNJp8rncDO2pJHTuTDk7mmbX9u20tLaWHy0t1XLorH2Ax1MpMuk0Pr8ffyBAsVRiYnz8vPch1tJ8Dycr2RMRkUXJzILAPwBvAMLAXwBHgO8Dr1YO+5Jz7ptm9mHg9wAP+Avn3A8WPuLG4oB4II7niniuQMCCRAKRGZcA8Fl5pu7a2Fra8228duBVMoEsaV+aZDBFyp8iE8jipjlTppjmwMR+Dkzsr9ZZu9FabCVRitNSaiFaihDxR7CQjzuuuYuwr5yUFYtFspXEL1Pp4cuk02TSadLpifP2ChY9j+ToKMnR0Sn1wWCwmvjFWlvJZjKkvTQTxQk85xGwAK20Eixc3LqBc03UajGcrGRPREQWq98Ehp1z/7uZdQHPAv8X8Dnn3F+dPsjMlgN/CGwCIsAuM/uhc+7CXUVNLuKPkC6miVrsTM+eK1TX15uNJ3I/hU5HT6abeDFOZ6mLwfAQpajjfSvfx7HMMY5mj3Ise5SjmaNkS+cuzuzMkQqkSHHu+n9PvPhTov4Y7cE22oMdtAXb6Ah20JZoo72rg+XBlbQF2/nBE/+TWDpCyRzmIEAAf8mPYfiYflJHoVBgdGSE0ZGp9yNGieDMUcKRJ0eyNEb/a68RDAYJBINTvgaDQXx+f3WYOJVMcvDVV8kXCrhSifHxcVJjY1x2xRWzTtRqMbNYyZ6ISI14B5+qlgOX3VLHSJrW/wC+Pem5B9wMXGVm76Hcu/evgVuB3ZXkLmdm+4GNwFMsYauiqwhZiKSXJFvKEPFF6Qp00RvpnfU5kt4Y8Wicwdhwtc45x3gxxaroalZFV0+pHymMcCx7jKOZIxzLHuNY9hjD+aHpTl2VKabJFNMczx4/7zG+uI9INEykFCHqRYmUIgRdgHy8yO+u/10yE2kmxsfPPCYmKJ1ngogPHzjwn4477zi4f/+0x0J5mDsYChEIBKrbzPkrW88ZkM/leP2117j62msJBAIz3j+Yy2bPGWIOBAIX3MVkJkr2RERqxHvun6tlJXvzzzk3DmBmccpJ36cpD+d+2Tn3tJl9CvhT4DlgbNJbU0DbdOc0s3uBewHWrFlTs9gbwcbERnbkt9MX7Zsyo3djYuOsz5EItJEpZoj6z+yPmy3lSATO/fGaGZ2hTjpDnVyXuK5a/8qrL/HKxCu8XnqdNGmKvhLmM/K+PBNMUHQzL59SokQ6kCFNBkJTe+r+zwOfoTvUXb7vsKuXnr4e+kJraSslIAvpSvJ3aOAg4WL4wpNEpuGcI5/LTbmn8OxEcnhggN0DA/h8PsKRCOFIhEjlazgSIRwOT6mbmJigkMtVJ4oEw2FaWlouKq7JlOyJLFbP/eWZcmw5XPlb9YtFpE7MbDXwPeDvnHNfN7N259xo5eXvAX8LPAbEJ70tDowyDefcfcB9UN5Bo0ZhN4T5mNG7pXMLPzj5fQAivjDZUo5cKcvdPXfP+hwT6XGsCFf41hOwIJ4r4BU8VgZW8saNm5jwJhj1RhnNjzBaGGO0MMJoYZSxwmj5qzd2wYSw6IqczJ3kZO7kOa+FfWF6wr30tPUwUDhBLB+jM99BrBjF5/yUfCXwwTXLr8UrFChUHpPL060heD6lUql6n+H5+Hw+Sq6Ew1H0l/D5fIQyIXp6Z9/jejYleyKLVesaaLvizPOJY9DSV7945Bz+N2yqdwhNzcyWAQ8DH3XO/bhS/ZCZfcw59yTwi8DTwJPAZ80sQrnnbwPwQj1ibjRzndH7xo4bAdh1ahdJb4xEoI27e+6u1s/GMKdwRUeqlKpOjIi4CMOBU/jMRzwYJx6Ms3rSkPBkJVfi+8f/mRPZ44wWRsmUMjjnysuvuDz5Uv683ztXynEkc5gjmcMQpPyY1IFmzggR4kVeJhqOEo1FifqjRP2x8ldfhIi1EHFhQqUwB/btI5j343N+/PjxO1+5p9DA7wucd+h4SnsqyaNhBDw/xWCRZDjJybETLFuxYtY/18mU7ImI1EjwxvfUO4Rm90mgA/iMmX2mUvdvgb8xszxwArjXOZc0sy8AOwEf8Cnn3KXfACVTvLHjxotK7s52PHyC6HgYfOWEiGI5CRtNJGf1fp/5uKXjFnYMbacn0jtlSHpb9110hroYzA8ymBtgMDdYfuTL5ekmjEzmzJEjx0BuYHaNmbyBiYNwKUSoFCZEiKt6ryLha6O11EqsFCHkhfAX/JTyRXK5XHVLuXOWemkxLOBjIHVuz+RsKdkTEZFFyTn3ceDj07x0+zTH3g/cX/Og5KIlgymSiSQ9mW6CxQAFv8dA6yAEZ3/v3ExD0qujq8/pGXTOMe6NM5gfYCA3yP7RfRweO0Tal6FgBQoUcDaHkXyDnD9Pzl/uWXx85PFpD/Obn0RrgraONhL+NkZPjZDItxIpRIh4YULpMP6wn0w4M+37Z0PJnoiIiNRN1B9lODTEsegJghag4DwKxTxd/u6LOs/FDkmbWXWIeF3L5YwWRoi7OLHxKJaDUsgxHpvAF/OzuXMzmWLmzKOUmfq88jg8cQhXLOFZkaJvdvvyFl2RkcLImS3pTg8nA+bgjckb6El3UYxf+uLOSvZEFivzz3yMiEiDWxldSchCjHlj1SVguiPd9ER6FjSOkcII7S3tFFvP9OZFXYxRb5TlkdndK/fFA3+LN56nPdNOwPOTDmYYiYxAxLi9awtj3ijJQpKkN8ZYIUmyMHbO1nOTOYPxwDhhF6Ylr9m4IkvP5b9S7whkBt6ru6vlwBV31DESkcY1H0vAzIeOYAfpYpqYP1aty5QydAQ7LvCuqSL+CKeiaUZiY9WFqs356Qx1srV767TvyZVylQSwnPw9cWA3OX952ZlgKUi7v52WYAtWuPS2KdkTEakR74UHq2UleyLTm48lYObDxsRGdgxtB5iSdG7u2Dzrc1zKQtXl5V/Kew8DDAUGyOQyZC2D5zz85qdYLBGNRM97jpko2RMREZG6musSMPMVw1yTzvnopezo6uTU66/i/OD3+SkVHZlimr6+S//5KNkTEakR//pzJoWKSAOba9I5HwnjUHCYUg/4xgyfZ5QCjlKnYyg4PPObz0PJnohIjQSvf3u9QxCRBTbXhHGkMEJ7WwfWfmbpmdP7Cl8q3yW/U0RERETmVUew45wZuhc7UeRsSvZEREREGsTGxEayxQzpYhrnHOlies6zk5XsiYiIiDSI0/f9xfzlNf5i/hjbuu+a09Cw7tkTEamRwsuPVMvBDW+uYyQispjM9+xkJXsiIjVS/Pn2alnJnojUi4ZxRURERJqYOedmPkpE6sLMHgQm7wbeDQzVKZx6WGrthcXV5iHn3NvqHUStmNkg0F/vOC6g0X9XGj0+aPwYZ4pvVtegkj2RRcTM9jjnNtU7joWy1NoLS7PNcmka/Xel0eODxo9xvuLTMK6IiIhIE1OyJyIiItLElOyJLC731TuABbbU2gtLs81yaRr9d6XR44PGj3Fe4tM9eyIiIiJNTD17IiIiIk1MyZ5IAzGzoJl93cweN7OdZna1ma03s12V518yM1/l2A+b2R4ze8LM3lXv2C+Wmd1mZjsq5Vm30cyiZvadyrH/YmY9dWzGRTmrzW+stGGHmT1kZssq9U3VZpkbM3u28juyw8z+30b5e7AYrt+zYrzJzI5O+lm+v14xVv7O/1Pl/E+a2btr/jN0zumhhx4N8gDeA3yrUr4b+A7wz8C2St3fA+8FlgPPA2Gg7XS53vFfRDv/uBLzE5Xns24j8G+B/1A59teAz9e7PZfY5keBN1bKvwd8rtnarMecf2ciwLNn1dX978FiuH6nifFDwB+ddUxdYgR+B/ibSrkLOFTrn6F69kQayz4gUPlfXQIoADdTTgwAHgDeAtwK7HbO5ZxzY8B+YGMd4r1UB4D3TXp+MW3cAjx41rGLwdlt/jXn3HOVcgDI0nxtlrm5AYiZ2cNm9oiZbaYx/h4shut3uhjfaWaPmdlXzCxexxj/B/CZSc89avwzVLIn0ljGgTcAPwfuB75AeSLV6ZlUKcr/w0sAY5Ped7p+UXDOfYdyInvaxbRxcv2iaffZbXbOHQcws9uBjwJ/TZO1WeYsDfwX4B7g94H/TgP8PVgM1+80MT4J/Dvn3J3Aa8Cf1itG59y4cy5VSTi/DXyaGv8MleyJNJZ/AzzknLuS8v/qvwqEJr0eB0aBZKV8dv1iVZpUnqmNk+sXdbsr9w39PfBO59wgS6DNclH2Af/Nle0DhoFlk15vlL8Hi+H6/Z5z7unTZeBG6hijma0GtgP/5Jz7OjX+GSrZE2ksI5z5H9spIAg8a2bbKnVvB3ZS/l/qVjOLmFkbsAF4YWFDnVcX08bdwDvOOnbRMbPfpNyjt80591qluqnbLBftg8BfAZhZH+UenYcb8O/BYrh+HzKzWyvlXwSerleMlclYDwN/4pz7h0p1TX+GgXmLXkTmw18D/2BmOyn36H0S2APcb2Yh4GXg2865opl9gfJF7gM+5ZzL1ivoefBHzLKNZvYl4KtmtgvIAx+oW9SXyMz8lIfoDwHfNTOAR51zf9qsbZZL8hXgHyv/7o5y8jdE4/09WAzX70eAL5pZHjgB3OucS9Ypxk8CHcBnzOz0vXsfB75Qq5+hFlUWERERaWIaxhURERFpYkr2RERERJqYkj0RERGRJqZkT0RERKSJKdkTERERaWJK9kRERESamJI9ERERkSamZE9ERESkiSnZExEREWliSvZEREREmpiSPREREZEmpmRPRBqSmV1mZt8xsx1mttvM/s7M4pXX/oOZ/f4073m7mf3YzH5oZo+Y2W/UKK7nzOyrZvbLZvaCmf3hfH8fEZH5Eqh3ACIiZzOzKPDPwIeccz+t1P0r4BvAuy7w1r8HbnDOjVYSw5+Z2Q+dcwPzGN4dwI+dc39kZl8GPuGc+/48nl9EZF4p2RORRvRO4NHTiR6Ac+6rZvYRM1t3gfedBD5uZt8GXgI2OOdykw8ws9uAzwMGHAV+A7ga+FugCGSBDzvnDpnZx4APAA74/4D/CXwaiJnZIcqJ561mNuSc+8k8tFuk6ZjZbwPvAGLA5cB/ds79Yz1jWmo0jCsijWgdcGCa+oPAmgu8792UP1C+ARwHPmFmdtYx9wG/45y7DfgRsAG4H/ioc+4XgL8DPmdm1wDvB7ZUHr8MRIH/G/i6c+7zwIPAHyvRE5lRm3PuXZSv0X9f72CWGiV7ItKIjgJvmKb+CuDQdG8wsw5grXPuT5xzG4Gbgbdx7rDvMufcywDOub9zzj0D9Dnnnqu8/hhwLXAdsBb4MfAI0AWsn0ObRJay5ypfDwOROsaxJCnZE5FG9L+Au83s1tMVZvYhYNA599p53hMGvmVmqyvPjwMngNxZxx0zsysq5/wTM3tvpW5j5fVfAPYBrwAvAnc557YB/wg8P9eGiSxRrt4BLGW6Z09EGo5zbtzMfgn4azProvy3ai/w65MO+0QlAQRIOefuqtxj910z8wA/8APn3MNnnf73gH8wsxLlhPBvgNeBL1aGfD3gd51zr5nZj4FdZhYGnqTc4ygisqiYc0q2RURERJqVhnFFREREmpiSPREREZEmpmRPREREpIkp2RMRERFpYkr2RERERJqYkj0RERGRJqZkT0RERKSJKdkTERERaWL/P3Ig4P0uZig0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "make_plots(\n",
    "    df,\n",
    "    \"./plots/census_income_covshift.pdf\",\n",
    "    intervals_xlabel=\"OLS coeff\",\n",
    "    true_theta=true_theta,\n",
    "    empty_panel=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb197baf-efcc-49cd-98b0-359dbbc9e517",
   "metadata": {},
   "source": [
    "### Power experiment\n",
    "\n",
    "For PPI and the classical approach, find the smallest value of ```n``` such that the method has power 80% against the null $H_0: \\theta^* < 800$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8858298f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find n such that we reject H0: OLS coeff < 800 with probability 80% using a test at level alpha\n",
    "num_experiments = 100\n",
    "list_rand_idx = [np.random.permutation(n_inD) for i in range(num_experiments)]\n",
    "\n",
    "\n",
    "def _to_invert_ppi(n):\n",
    "    n = int(n)\n",
    "    nulls_rejected = 0\n",
    "    # Data setup\n",
    "    for i in range(num_experiments):\n",
    "        rand_idx = list_rand_idx[i]\n",
    "        _X = X_inD[rand_idx[:n]]\n",
    "        _Y = Y_inD[rand_idx[:n]]\n",
    "        _Z = Z_inD[rand_idx[:n]]\n",
    "        _Yhat = Yhat_inD[rand_idx[:n]]\n",
    "        importance_weights = np.array(\n",
    "            [\n",
    "                (\n",
    "                    weights[0] / inD_weights[0]\n",
    "                    if z == 0\n",
    "                    else weights[1] / inD_weights[1]\n",
    "                )\n",
    "                for z in _Z\n",
    "            ]\n",
    "        )\n",
    "\n",
    "        ppi_ci = ppi_ols_ci(\n",
    "            _X, _Y, _Yhat, X_ooD, Yhat_ooD, w=importance_weights, alpha=alpha\n",
    "        )\n",
    "\n",
    "        if ppi_ci[0][coordinate] > 800:\n",
    "            nulls_rejected += 1\n",
    "    return nulls_rejected / num_experiments - 0.8\n",
    "\n",
    "\n",
    "def _to_invert_classical(n):\n",
    "    n = int(n)\n",
    "    nulls_rejected = 0\n",
    "    # Data setup\n",
    "    for i in range(num_experiments):\n",
    "        rand_idx = list_rand_idx[i]\n",
    "        _X = X_inD[rand_idx[:n]]\n",
    "        _Y = Y_inD[rand_idx[:n]]\n",
    "        _Z = Z_inD[rand_idx[:n]]\n",
    "        _Yhat = Yhat_inD[rand_idx[:n]]\n",
    "        importance_weights = np.array(\n",
    "            [\n",
    "                (\n",
    "                    weights[0] / inD_weights[0]\n",
    "                    if z == 0\n",
    "                    else weights[1] / inD_weights[1]\n",
    "                )\n",
    "                for z in _Z\n",
    "            ]\n",
    "        )\n",
    "\n",
    "        classical_ci = ppi_ols_ci(\n",
    "            _X,\n",
    "            _Y,\n",
    "            _Yhat,\n",
    "            X_ooD,\n",
    "            Yhat_ooD,\n",
    "            w=importance_weights,\n",
    "            alpha=alpha,\n",
    "            lam=0,\n",
    "        )\n",
    "\n",
    "        if classical_ci[0][coordinate] > 800:\n",
    "            nulls_rejected += 1\n",
    "    return nulls_rejected / num_experiments - 0.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f3ad00cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The PPI test requires n=225 labeled data points to reject the null.\n",
      "The classical test requires n=318 labeled data points to reject the null.\n"
     ]
    }
   ],
   "source": [
    "n_ppi = int(brentq(_to_invert_ppi, 100, 2000, xtol=1))\n",
    "n_classical = int(brentq(_to_invert_classical, 100, 10000, xtol=1))\n",
    "print(\n",
    "    f\"The PPI test requires n={n_ppi} labeled data points to reject the null.\"\n",
    ")\n",
    "print(\n",
    "    f\"The classical test requires n={n_classical} labeled data points to reject the null.\"\n",
    ")"
   ]
  }
 ],
 "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}