notebooks/diffpriv.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Local differential privacy for privacy-preserving data analysis.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import numpy.random as npr\n",
    "import itertools as it\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "npr.seed(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1: Randomised response.\n",
    "\n",
    "Let us consider the so-called \"randomised response\" mechanism for collecting responses to a potentially embarrassing or private question. In this mechanism, the respondents are instructed to as follows:\n",
    "\n",
    "1. Flip a coin.\n",
    "2. If **tails**, then respond truthfully.\n",
    "3. If **heads**, then flip a second coin and respond *Yes* if heads and *No* if tails.\n",
    "\n",
    "i. Calculate the probability of the randomised response being \"yes\" and the probability of the randomised response being \"no\", conditional to true response being \"yes\"/\"no\".\n",
    "\n",
    "ii. Verify that this mechanism satisfies local $\\epsilon$-DP with $\\epsilon = \\ln 3$.\n",
    "\n",
    "iii. Assuming a population of independent respondents who have the probability $\\pi$ to answer yes. Design a method to estimate $\\pi$ from the noisy responses. Simulate the mechanism with $n$ participants for $n = 10, 100, 1000, 10000$ with $\\pi = 0.01, 0.1$. How accurately can you estimate $\\pi$ in the different cases?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I.  \n",
    "$Pr(yes|truth=yes) = Pr(yes,random|truth=yes) + Pr(yes,non-random|truth=yes) = $  \n",
    "$Pr(yes|random)Pr(coin = random) + Pr(yes|truth=yes,coin=non-random)Pr(coin=non-random) =  $  \n",
    "$0.5\\cdot0.5 + 1 \\cdot0.5 = 0.75$\n",
    "\n",
    "Similarly, $Pr(no|truth=yes) = 0.25$, $Pr(no|truth=no) = 0.75$, $Pr(yes|truth=no) = 0.25$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "II.  \n",
    "\n",
    "$\\frac{Pr(yes|truth=yes)}{Pr(yes|truth=no)}= 3 \\leq e^{\\epsilon} \\Leftrightarrow \\ln(3) \\leq \\epsilon$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(yes | truth = yes) = 0.75, P(no | truth = yes) = 0.25\n",
      "P(no | truth = no) = 0.75, P(yes | truth = no) = 0.25\n"
     ]
    }
   ],
   "source": [
    "#i) Monte Carlo simulation\n",
    "n = 10000\n",
    "\n",
    "truth = 'yes'\n",
    "random_choice = np.where(npr.uniform(size=n)>0.5, 'yes', 'no')\n",
    "choice = np.where(npr.uniform(size=n)>0.5, truth, random_choice)\n",
    "ctr = (choice==truth).sum()\n",
    "print(f'P({truth} | truth = {truth}) = {ctr/n:.2f}, P({\"no\" if truth == \"yes\" else \"yes\"} | truth = {truth}) = {(n-ctr)/n:.2f}')\n",
    "\n",
    "truth = 'no'\n",
    "random_choice = np.where(npr.uniform(size=n)>0.5, 'yes', 'no')\n",
    "choice = np.where(npr.uniform(size=n)>0.5, truth, random_choice)\n",
    "ctr = (choice==truth).sum()\n",
    "print(f'P({truth} | truth = {truth}) = {ctr/n:.2f}, P({\"no\" if truth == \"yes\" else \"yes\"} | truth = {truth}) = {(n-ctr)/n:.2f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7511"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "truth = 'yes'\n",
    "random_choice = np.where(npr.uniform(size=n)>0.5, 'yes', 'no')\n",
    "choice = np.where(npr.uniform(size=n)>0.5, truth, random_choice)\n",
    "(choice==truth).sum()/n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ii)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\frac{\\mathcal{P}(\\text{yes}|\\text{truth=yes})}{\\mathcal{P}(\\text{yes}|\\text{truth=no})} = \\frac{0.5\\cdot (1_{truth} + 0.5\\cdot(0+1_{truth}))}{0.5\\cdot(0_{truth}+0.5\\cdot(1+0_{truth}))} = e^{\\epsilon} = 3 \\Leftrightarrow \\epsilon = \\ln(3) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n: 10, pi: 0.01, estimate: 1.30\n",
      "n: 100, pi: 0.01, estimate: 0.90\n",
      "n: 1000, pi: 0.01, estimate: 0.99\n",
      "n: 10000, pi: 0.01, estimate: 0.98\n",
      "n: 10, pi: 0.1, estimate: 0.90\n",
      "n: 100, pi: 0.1, estimate: 1.00\n",
      "n: 1000, pi: 0.1, estimate: 0.88\n",
      "n: 10000, pi: 0.1, estimate: 0.91\n"
     ]
    }
   ],
   "source": [
    "#iii)\n",
    "n_, pi_ = [10, 100, 1000, 10000], [0.01, 0.1]\n",
    "for pi, n in it.product(pi_, n_):\n",
    "    truth = np.where(npr.uniform(size=n)>pi, 'yes', 'no')\n",
    "    random = np.where(npr.uniform(size=n)>0.5, 'yes', 'no')\n",
    "    choice = np.where(npr.uniform(size=n)>0.5, truth, random)\n",
    "    i = (choice=='yes').sum()\n",
    "    print(f'n: {n}, pi: {pi}, estimate: {max(round((i/n-.25)*2,4),0):.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the $\\pi = 0.1$ even low numbers of samples are sufficient to be rather close to the true value; whereas even 10000 samples do not suffice to be within a $+- 10\\%$ bound."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2: Generalising randomised response\n",
    "\n",
    "Consider a generalised randomised response mechanism, where each person flips $k$ coins in step 1 and responds truthfully if either  \n",
    "i. all of them display tails  \n",
    "ii. any one of them displays tails.\n",
    "\n",
    "Verify that both of these produce $\\epsilon$-DP mechanisms for each $k$. Evaluate and plot $\\epsilon$ as a function of $k$ for each case and $k = 1, \\dots, 10$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verification.  \n",
    "Case 1:  \n",
    "$Pr(yes|truth=yes) = 0.5^{k} + (1-0.5^{k})\\cdot 0.5 = 0.5(0.5^{k} + 1)$, using symmetry properties:\n",
    "$Pr(yes|truth=no) = 1 - Pr(yes|truth=yes)$  \n",
    "$\\frac{Pr(yes|truth=yes)}{Pr(yes|truth=no)} = \n",
    "\\frac{0.5(0.5^{k} + 1)}{1-(0.5(0.5^{k} + 1))} = \n",
    "\\frac{0.5(0.5^{k} + 1)}{1-(0.5^{k}0.5 - 0.5)} = \n",
    "\\frac{0.5(0.5^{k} + 1)}{0.5(1-0.5^{k})} = \n",
    "\\frac{0.5^{k} + 1}{1-0.5^{k}} \n",
    "\\leq \n",
    "e^{\\epsilon}$  \n",
    "Obviously, this produces an $\\epsilon$-DP mechanisms for each $k$, where privacy increases ($\\epsilon$ decreases) with increasing $k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Case 2: Inversion of case 1 given symmetric probabilities of tails and heads.  \n",
    "$Pr(yes|truth=yes) = (1-0.5^{k}) + 0.5^{k}\\cdot 0.5 = 1+0.5^{k}(0.5-1) =1-0.5\\cdot 0.5^{k}$\n",
    ", using symmetry properties:\n",
    "$Pr(yes|truth=no) = 1 - Pr(yes|truth=yes)$  \n",
    "$\\frac{Pr(yes|truth=yes)}{Pr(yes|truth=no)} = \n",
    "\\frac{1-0.5\\cdot 0.5^{k}}{1-(1-0.5\\cdot 0.5^{k})} = \n",
    "\\frac{1-0.5\\cdot 0.5^{k}}{0.5\\cdot 0.5^{k}} =\n",
    "\\frac{2-0.5^{k}}{0.5^{k}} \\leq e^{\\epsilon}$  \n",
    "Obviously, this produces an $\\epsilon$-DP mechanisms for each $k$, where privacy decreases ($\\epsilon$ increases) with increasing $k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_random(k, mode='i'):\n",
    "    if mode=='i':\n",
    "        return 1-.5**k\n",
    "    if mode=='ii':\n",
    "        return 1-(1-.5**k) #ofc could be just .5**k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAm4klEQVR4nO3deXxV5Z3H8c/vZiGQBJAEElaDimICCIqIe9wgghVrbcEZUWutXcZabZ0ZZzrV1tra1mmttrZTam3tVMGlahlFFtEUt1pFEdkUZJEgOwiELdtv/rgXSGKAkNyTk3vv9/165eW55zz33N99XpJvznOWx9wdERFJXZGwCxARkXApCEREUpyCQEQkxSkIRERSnIJARCTFKQhERFKcgkBSgpn9s5nNrPfazey4Nq7hbDN7vy0/M/a5K83swrb+XEkcCgJpd2K/uHabWWW9n1+1Zp/u/oi7j4pXjS2s4WV3PyHMGkSakh52ASIH8Rl3fyHsIkRSgY4IJKGY2bVm9qqZ/crMtpnZEjO7oNH25Wa2w8xWmNk/11v/ykH22cXM/mRmG81slZn9l5lF6r/PzP7bzLbG9nnxIeo72czeiX3+E2b2mJndFdtWamYVseV/N7MnG733PjO7v15NvzeztWa2xszuMrO0ltTU6DNOjLW/sjntJTUoCCQRnQZ8COQDdwBPmVk3M8sG7gcudvdc4AxgXjP290ugC3AMcC5wNfDFRp/3fuzzfgr83sys8U7MLBN4Gvgj0A2YDHz2IJ85BRhjZrmx96YBXwAejW3/I1ADHAcMA0YB1x9pTY3qOxmYAXzD3Scfqq2kFgWBtFfPmNkn9X6+XG/bBuAX7l7t7o8R/YU4NratDhhkZh3dfa27LzzUh8R+AU8A/sPdd7j7SuBnwMR6zVa5++/cvRZ4GOgJFDSxu5FEh1vvj9X2FPCPpj7X3VcBb3MgKM4Hdrn7382sABgD3OzuO919A3BvrM4jrWmfs4GpwNXu/uwh2kkKUhBIe3WZu3et9/O7etvWeMOnJa4Cern7TmA88FVgrZk9Z2YDD/M5+UBGbB/199e73ut1+xbcfVdsMaeJffVqorbVh/jsR4F9QzT/xIGjgaNjNa3dF4TAb4EeLahpn68Cr7l7+SHaSIpSEEgi6t1oGKQf8DGAu89w94uI/oW8BPhdE++vbxNQTfSXb/39rWlBXWubqK3vIdo/AZSaWR+iRwb7gmA1sBfIrxeEnd29pAU17fNVoJ+Z3duKfUiSUhBIIuoB3GRmGWb2eeBEYJqZFZjZuNi5gr1AJdGhooOKDa08DvzQzHLN7GjgW8CfW1DX60AtcKOZpZvZOGDEIT57I1AO/AFY4e6LY+vXAjOBn5lZZzOLmNmxZnZuC2raZwdQBpxjZj9uxX4kCSkIpL36v0b3ETxdb9sbwACif83/ELjC3TcT/f/5W0SPDrYQPfH7tWZ81jeAncBy4BWif5k/dKQFu3sVcDnwJeAT4CrgWaKhdDCPAhdy4Ghgn6uBTGARsBV4kuhRTou5+yfARcDFZvYDADN73sz+szX7lcRnmphGEomZXQtc7+5nhV1Lc5jZG8D/uPsfwq5F5GB0RCASR2Z2rpkVxoaGrgGGANPDrkvkUHRnsUh8nUD0nEM20aGmK2Jj/iLtloaGRERSnIaGRERSXMINDeXn53tRUVHYZbTKzp07yc7ODruMdkP9cYD6oiH1R0Ot6Y+5c+ducvfuTW1LuCAoKirirbfeCruMVikvL6e0tDTsMtoN9ccB6ouG1B8NtaY/zGzVwbZpaEhEJMUpCEREUpyCQEQkxSXcOQIRkfqqq6upqKhgz549YZcSuC5durB48eJDtsnKyqJPnz5kZGQ0e78KAhFJaBUVFeTm5lJUVMRh5uZJeDt27CA3N/eg292dzZs3U1FRQf/+/Zu9Xw0NiUhC27NnD3l5eUkfAs1hZuTl5R3x0ZGCQEQSnkLggJb0RcoEwdsfbeXHzy9Bj9QQEWkoZYJg4Zpt/M/fPmTZhsqwSxGRJGNmXHXVVftf19TU0L17dy655JIj2k9RURGbNm1qdvvrr7+eRYsWHdFnNCVlgmBUSSEAMxauO0xLEZEjk52dzYIFC9i9ezcAs2bNonfv3od5V+s9+OCDFBcXt3o/KRMEBZ2zGNavK9MVBCISgDFjxvDcc88BMHnyZK688sr927Zs2cJll13GkCFDGDlyJPPnzwdg8+bNjBo1ipKSEq6//voGQ9d//vOfGTFiBEOHDuUrX/kKtbW1n/rM0tLSuDxyJ6UuHy0rKeTu55ewessu+nbrFHY5IhJn3/+/hSz6eHtc91ncqzN3fKbksO0mTJjAnXfeySWXXML8+fO57rrrePnllwG44447GDZsGM888wwvvvgiV199NfPmzeP73/8+Z511FrfffjvPPfccv//97wFYvHgxjz32GK+++ioZGRl8/etf55FHHuGzn/1sXL/bPilzRAAwOjY8NHPR+pArEZFkM2TIEFauXMnkyZMZM2ZMg22vvPIKEydOBOD8889n8+bNbN++nTlz5uw/tzB27FiOOuooAGbPns3cuXM59dRTGTp0KLNnz2b58uWB1Z5SRwRF+dkMLMxlxoJ1fOms5t9sISKJoTl/uQfp0ksv5dZbb6W8vJzNmze3eD/uzjXXXMPdd9/dYP2OHTtaW2KTUuqIAKJHBW+u2sLGHXvDLkVEksx1113HHXfcweDBgxusP/vss3nkkUeA6KOk8/Pz6dy5M+eccw6PPvooAM8//zxbt24F4IILLuDJJ59kw4YNQPQcw6pVB32KdKulXBCUDSrEHWZpeEhE4qxPnz7cdNNNn1r/ve99j7lz5zJkyBBuu+02Hn74YSB67mDOnDmUlJTw1FNP0a9fPwCKi4u56667GDVqFEOGDOGiiy5i7drgpr5OqaEhgIGFuRyd14kZC9fxT6f1C7scEUkClZWfvj+ptLR0/yQy3bp145lnnvlUm7y8PGbOnNnkPsePH8/48eMbrGs8NFReXt6iehtLuSMCM2N0SSGvfbiJbburwy5HRCR0KRcEED1PUF3rvLRkQ9iliIiELiWDYFjfrvTI7aC7jEWShJ4hdkBL+iIlgyASiQ4Plb+/kd1Vn75bT0QSR1ZWFps3b1YYcGA+gqysrCN6X8qdLN5ndEkh//v3VcxZunH/jWYiknj69OlDRUUFGzduDLuUwO3Zs+ewv+T3zVB2JFI2CE47phtdOmYwY8E6BYFIAsvIyDii2bgSWXl5OcOGDYv7fgMbGjKzh8xsg5ktOMh2M7P7zWyZmc03s5ODqqUpGWkRLjyxgBcWr6e6tq4tP1pEpF0J8hzBH4GyQ2y/GBgQ+7kB+E2AtTRpdEkB2/fU8PflLb8VXEQk0QUWBO4+B9hyiCbjgD951N+BrmbWM6h6mnLO8d3pmJHG9AW6ekhEUleY5wh6A6vrva6IrfvUfdRmdgPRowYKCgridjcdQEk3eHbeai7ouolIG817WllZGdfvkOjUHweoLxpSfzQUVH8kxMlid58ETAIYPny477ttOx62dV3DN6fMo8sxJ3HK0d3itt9DKS8vJ57fIdGpPw5QXzSk/mgoqP4I8z6CNUDfeq/7xNa1qfMG9iAjzTQ8JCIpK8wgmApcHbt6aCSwzd2De7zeQXTOyuDM4/KZsXC9bkgRkZQU5OWjk4HXgRPMrMLMvmRmXzWzr8aaTAOWA8uA3wFfD6qWwykrKeSjLbtYvDaYSR9ERNqzwM4RuPuVh9nuwL8E9flH4sLiAiJPv8f0heso7tU57HJERNpUSj5rqLH8nA4ML+rGDJ0nEJEUpCCIKSsp5P31O1ixaWfYpYiItCkFQczoQdHnDenR1CKSahQEMb27dmRw7y66jFREUo6CoJ6yQYXMW/0J67btCbsUEZE2oyCoZ9/jqGcu0lGBiKQOBUE9x/XI4dju2RoeEpGUoiBopGxQIW+s2MLWnVVhlyIi0iYUBI2UlfSkts55YfH6sEsREWkTCoJGBvXuTO+uHXUZqYikDAVBI2bGqJIC5izdROXemrDLEREJnIKgCWUlhVTV1PG39zeGXYqISOAUBE0YXtSNvOxMpmt4SERSgIKgCWkR46LiAl5cvJ491bVhlyMiEigFwUGMHlTIzqpaXvtwU9iliIgESkFwEGccm0duh3RmLNBlpCKS3BQEB9EhPY3zT+zBrMXrqamtC7scEZHAKAgOYXRJIVt2VvHmyq1hlyIiEhgFwSGce3x3OqRHdHOZiCQ1BcEhZHdI55zjuzNj4TqiUyyLiCQfBcFhjC4pZO22Pcyv2BZ2KSIigVAQHMaFJ/YgLWIaHhKRpKUgOIyunTI5/Zg8pi/Q8JCIJCcFQTOMHlTI8k07WbahMuxSRETiTkHQDKOKCwA0c5mIJCUFQTMUdM7i5H5dmaG5jEUkCSkImqlsUCEL1mxn9ZZdYZciIhJXgQaBmZWZ2ftmtszMbmtiez8ze8nM3jGz+WY2Jsh6WmN0SSGArh4SkaQTWBCYWRrwAHAxUAxcaWbFjZr9F/C4uw8DJgC/Dqqe1jo6L5uBhbnMXKiH0IlIcgnyiGAEsMzdl7t7FTAFGNeojQOdY8tdgI8DrKfVygYV8uaqLWzcsTfsUkRE4saCujbezK4Aytz9+tjricBp7n5jvTY9gZnAUUA2cKG7z21iXzcANwAUFBScMmXKlEBqPpzVO+r47qu7ubYkk9K+GS3eT2VlJTk5OXGsLLGpPw5QXzSk/mioNf1x3nnnzXX34U1tS29VVa13JfBHd/+ZmZ0O/K+ZDXL3Bs99dvdJwCSA4cOHe2lpadtXGq2DB5eUs7w6m++VjmjxfsrLywnrO7RH6o8D1BcNqT8aCqo/ghwaWgP0rfe6T2xdfV8CHgdw99eBLCA/wJpaxcwoKynk9Q83sW13ddjliIjERZBB8CYwwMz6m1km0ZPBUxu1+Qi4AMDMTiQaBBsDrKnVRg8qpLrWeWnJhrBLERGJi8CCwN1rgBuBGcBiolcHLTSzO83s0lizbwNfNrN3gcnAtd7OH+gztE9XeuR20F3GIpI0Aj1H4O7TgGmN1t1eb3kRcGaQNcRbJGKMLinkybkV7K6qpWNmWtgliYi0iu4sboGyQYXsrq5lztJ2PYolItIsCoIWGNG/G107ZTBDw0MikgQUBC2QkRbhgoEFvLB4PdW1dYd/g4hIO6YgaKGyQYVs31PD35dvDrsUEZFWURC00NkD8umUmaarh0Qk4SkIWigrI43SE7ozc9F66ura9RWvIiKHpCBohdElhWzcsZd3Vm8NuxQRkRZTELTC+QN7kJkW0fCQiCQ0BUEr5GZlcOZxeUxfuI52fkO0iMhBKQhaaXRJIau37Gbx2h1hlyIi0iIKgla6sLiAiMF0TWEpIglKQdBK+TkdOLWom+4yFpGEpSCIg9Elhby/fgcrNu0MuxQRkSOmIIiD0YMKAZih4SERSUAKgjjo3bUjQ/p00WWkIpKQFARxMrqkkHmrP2Hdtj1hlyIickQUBHEyuiQ6PDRzkY4KRCSxKAji5LgeORzXI0fDQyKScBQEcVRWUsgbK7awZWdV2KWIiDSbgiCORpcUUlvnvLB4fdiliIg0m4Igjgb17kzvrh2ZqctIRSSBKAjiyMwYXVLInKWbqNxbE3Y5IiLNoiCIs7JBhVTV1FH+/oawSxERaRYFQZydcvRR5GVnMmOhzhOISGJQEMRZWsQYVVLAi4vXs6e6NuxyREQOS0EQgNElheysquW1DzeFXYqIyGEpCAJwxrH55HZIZ8YCDQ+JSPsXaBCYWZmZvW9my8zstoO0+YKZLTKzhWb2aJD1tJXM9Ajnn9iDWYvXU1NbF3Y5IiKHFFgQmFka8ABwMVAMXGlmxY3aDAD+AzjT3UuAm4Oqp62VlRSyZWcVb67cGnYpIiKHFOQRwQhgmbsvd/cqYAowrlGbLwMPuPtWAHdPmmsuzz2hOx3SI5qjQETavfQA990bWF3vdQVwWqM2xwOY2atAGvA9d5/eeEdmdgNwA0BBQQHl5eVB1Bt3xd2MqW+v4tzcDZjZ/vWVlZUJ8x3agvrjAPVFQ+qPhoLqjyCDoLmfPwAoBfoAc8xssLt/Ur+Ru08CJgEMHz7cS0tL27bKFtqcW8G3n3iXbscN46S+XfevLy8vJ1G+Q1tQfxygvmhI/dFQUP0R5NDQGqBvvdd9YuvqqwCmunu1u68APiAaDEnhghN7kB4xpmt4SETasSMOAjOLmFnnZjR9ExhgZv3NLBOYAExt1OYZokcDmFk+0aGi5UdaU3vVtVMmI4/JY8aCdbh72OWIiDSpWUFgZo+aWWczywYWAIvM7F8P9R53rwFuBGYAi4HH3X2hmd1pZpfGms0ANpvZIuAl4F/dfXNLv0x7NHpQIcs37WTZhsqwSxERaVJzjwiK3X07cBnwPNAfmHi4N7n7NHc/3t2Pdfcfxtbd7u5TY8vu7t9y92J3H+zuU1r2Ndqv0cUFmKGZy0Sk3WpuEGSYWQbRIJjq7tWAxjqaoUfnLE7ud5TOE4hIu9XcIPgtsBLIJnplz9HA9qCKSjajSwpY+PF2Vm/ZFXYpIiKf0qwgcPf73b23u4+JDeesAs4LuLakMbqkEEA3l4lIu9Tck8V5Zna/mb1tZnPN7D6gS8C1JY2j87I5sWdnBYGItEvNHRqaAmwEPgdcEVt+LKiiktHokgLeWrWVjTv2hl2KiEgDzQ2Cnu7+A3dfEfu5CygIsrBkUzaoEHeYtUiPphaR9qW5QTDTzCbEbiaLmNkXiN4DIM10QkEuRXmddPWQiLQ7zQ2CLwOPAHtjP1OAr5jZDjPT1UPNYGaMHlTIa8s2sbNaV96KSPvR3CDoAlwL/MDdM4Ai4EJ3z3X35jxuQohePVRT57y7UXMZi0j70dwgeAAYCVwZe70D+FUgFSWxoX26UtC5A2+tqwm7FBGR/ZobBKe5+78AewBiE8lkBlZVkopEjMuG9ebtDbW8slQT24tI+9DcIKiOTT3pAGbWHdBkvC1w8wXH0yvbuOXxeWyu1KWkIhK+5gbB/cDTQA8z+yHwCvCjwKpKYh0z0/ja0Cy27a7mX5+cr8dTi0jomvuIiUeAfwPuBtYCl7n7E0EWlsz65kb4z4sH8uKSDTz82sqwyxGRFNfsqSrdfQmwJMBaUso1ZxTx8tJN/Oj5JZx2TB4n9tTFVyISjiCnqpRDMDN+esUQunbM4BuT32F3lS4pFZFwKAhClJfTgZ9/YSgfbqzkzmcXhV2OiKQoBUHIzhqQzw3nHMPkf3zE9AVrwy5HRFKQgqAd+PZFJzCkTxf+/S/v8fEnu8MuR0RSjIKgHchMj3D/hGHU1NZx82PzqK3TJaUi0nYUBO1EUX42d44bxD9WbOGBl5aFXY6IpBAFQTty+cm9GTe0F/fNXsrcVVvCLkdEUoSCoB0xM+66bBC9umZx0+R5bNtdHXZJIpICFATtTG5WBvdNGMa67Xv4ztPv6REUIhI4BUE7dHK/o/jWRcfz7Py1PDG3IuxyRCTJKQjaqa+eeywjj+nG96Yu5MONlWGXIyJJTEHQTqVFjF+MH0ZmeoSbJr/D3ho9gkJEghFoEJhZmZm9b2bLzOy2Q7T7nJm5mQ0Psp5EU9gli59+bggLP97OPdPfD7scEUlSgQVBbCKbB4CLgWLgSjMrbqJdLvBN4I2gaklko0oKmTjyaB58ZQXl728IuxwRSUJBHhGMAJa5+3J3rwKmAOOaaPcD4CfEpsGUT/vO2BM5oSCXW594l407NKuZiMSXBXV5opldAZS5+/Wx1xOJzn18Y702JwPfcffPmVk5cKu7v9XEvm4AbgAoKCg4ZcqUKYHU3FYqKyvJyck5ovdU7Kjj+6/vZmC3NG45pQMRs4Cqa3st6Y9kpb5oSP3RUGv647zzzpvr7k0Ovzd7Ypp4M7MI8HPg2sO1dfdJwCSA4cOHe2lpaaC1Ba28vJyWfAfPX8l3/7qQ5elHc/3Zx8S/sJC0tD+SkfqiIfVHQ0H1R5BDQ2uAvvVe94mt2ycXGASUm9lKYCQwVSeMD+6qkUdzUXEBP5m+hAVrtoVdjogkiSCD4E1ggJn1N7NMYAIwdd9Gd9/m7vnuXuTuRcDfgUubGhqSKDPjJ58bQrfsTG6a8g67qmrCLklEkkBgQeDuNcCNwAxgMfC4uy80szvN7NKgPjfZdcvO5N7xQ1mxaSffn6pZzUSk9QI9R+Du04BpjdbdfpC2pUHWkkzOODafr517LL8u/5Bzju/O2CE9wy5JRBKY7ixOULdcdDxD+3bltqfmU7F1V9jliEgCUxAkqIy06Kxm7vDNKfOoqa0LuyQRSVAKggTWL68Td102iLmrtnL/i5rVTERaRkGQ4C4b1pvLh/XmVy8u5R8rNKuZiBw5BUESuPOyQfTr1ombp7zDtl2a1UxEjoyCIAnkdEjnvgnD2LBjL7c9NV+zmonIEVEQJImT+nbl1tEn8PyCdUx5c3XY5YhIAlEQJJEbzj6Gs47L5/v/t5BlG3aEXY6IJAgFQRKJRIyff+EkOmWm843J89hTrVnNROTwFARJpkfnLO65YgiL127nJ9OXhF2OiCQABUESuuDEAq49o4g/vLqSF5esD7scEWnnFARJ6raLBzKwMJdbn5jPhu2a/E1EDk5BkKSyMtL45ZXD2FVVw7cef5e6Ol1SKiJNUxAksQEFuXz3kmJeWbaJ3728POxyRKSdUhAkuX8a0Y+ykkLumfE+8ys+CbscEWmHFARJzsz48ecG0z23AzdNfofKvZrVTEQaUhCkgK6dorOardqyizv+ujDsckSknVEQpIiRx+Rx43nH8Ze3K/jrvDVhlyMi7YiCIIV884IBnNyvK//19AJWb9GsZiISpSBIIelpEe6bMAyAm6a8Q7VmNRMRFAQpp2+3Tvzo8sG889EnfOfp9/Q8IhFREKSiz5zUi6+VHsvjb1Uw5r6XefujrWGXJCIhUhCkqH8vG8gj15/G3po6rvjNa9z9/GIdHYikKAVBCjvzuHym33w240/tx2//tpxLfvkK767+JOyyRKSNKQhSXG5WBndfPpiHrxvBzr01XP6b17hnxhL21ujoQCRVKAgEgHOP7870m8/h8mG9eeClD7n0l6+yYM22sMsSkTagIJD9unTM4J7Pn8RD1w5n664qLnvgVe6d9QFVNbrMVCSZBRoEZlZmZu+b2TIzu62J7d8ys0VmNt/MZpvZ0UHWI81z/sACZt1yLpee1Iv7Zi/lsgdeZdHH28MuS0QCElgQmFka8ABwMVAMXGlmxY2avQMMd/chwJPAT4OqR45Ml04Z/Hz8UCZNPIUNO/Yw7oFX+OXspboJTSQJBXlEMAJY5u7L3b0KmAKMq9/A3V9y933POvg70CfAeqQFRpUUMvOWcykb1JOfzfqAy3/9Gh+s3xF2WSISR+YezMxVZnYFUObu18deTwROc/cbD9L+V8A6d7+riW03ADcAFBQUnDJlypRAam4rlZWV5OTkhF3GEXtzXQ1/WriX3TXw2QEZlBVlkBaxVu83UfsjCOqLhtQfDbWmP84777y57j68qW3praoqTszsKmA4cG5T2919EjAJYPjw4V5aWtp2xQWgvLycRPwOpcAXK/fy3WcW8MSCdSzdnc1/f/4kjuvRun+oidofQVBfNKT+aCio/ghyaGgN0Lfe6z6xdQ2Y2YXAd4BL3X1vgPVIHOTndODX/3wy9185jJWbdzLm/pf53Zzl1GpOZJGEFWQQvAkMMLP+ZpYJTACm1m9gZsOA3xINgQ0B1iJxZGZcelIvZt5yDuce350fTlvMF377Ois27Qy7NBFpgcCCwN1rgBuBGcBi4HF3X2hmd5rZpbFm9wA5wBNmNs/Mph5kd9IO9cjNYtLEU7h3/EksXb+Di++bw0OvrKBORwciCSXQcwTuPg2Y1mjd7fWWLwzy8yV4ZsZnh/XhjGPzue0v87nz2UVMX7iOe64YwtF52WGXJyLNoDuLJS4KOmfx0LWncs8VQ1j88XbKfvEyf3p9pY4ORBKAgkDixsz4/PC+zLjlHIYXHcXtf13IVb9/Q9NiirRzCgKJu15dO/Kn60Zw9+WDeXf1J5T9Yg6PvvERQd2zIiKtoyCQQJgZV47ox4xbzuGkvl35z6ff4+qH/sHHn+wOuzQRaURBIIHqc1Qn/vyl0/jBuBLeWrmV0ffO4fE3V+voQKQdURBI4CIRY+LpRcy4+RxO7NWZf/vLfK7745us27Yn7NJEBAWBtKF+eZ2Y8uWR3PGZYl5fvplR9/6Np96u0NGBSMjaxbOGJHVEIsYXz+xP6Qk9uPWJd/nW4+/SK9v4fM0HjB3Sk+MLcsMuUSTlKAgkFP3zs3n8K6fz5NzV/OGlhdz/4lLum72UAT1yGDO4p0JBpA0pCCQ0aRFj/Kn9KNi5nOJTRjJjwTqenb+2QSiMHdKTsYN7MkChIBIYBYG0Cz1ys5h4ehETTy9iw/Y9TF+4jufmr+W+2Uv5xQtLOb4gdqSgUBCJOwWBtDs9Omdx9elFXF0vFJ5tFApjB/di7JBCjuuhUBBpLQWBtGuNQ+H5Bet47r21/GL2B9z7wgecUJAbO6egUBBpKQWBJIwenbO45owirjmjiPXb9zB9QXT46NOh0LPVs6aJpBIFgSSkgkah8Px7a5n23roGoTB2SE/GDFYoiByOgkASXkHnLK49sz/Xntl/fyg8995a7n3hA34+6wMGFh44Uji2u0JBpDEFgSSV+qGwbtsenl+wlmnvreXnsw6EwtjBPRmjUBDZT0EgSauwSxZfPLM/X6wXCs/NX8vPZn3Az+qFwtghPTlGoSApTEEgKaF+KKzdtpvn31vHtPcahsJJfboyoCCH4wtyGVCQQ2HnLMws7NJFAqcgkJTTs0tHrjurP9eddSAUZi9ZzwuL1/PYW6v3t8vNSmdAjxwG9MhVQEhSUxBISqsfCgCbK/eydEMlS9fv4IP1lSzdsEMBIUlPQSBST15OB/JyOjDymLwG6xUQkswUBCLNoICQZKYgEGmFeAfEpt117K2ppUN6Wlt/FUlhCgKRALQmIG7923RyOqSTl5NJt+xM8rIzycvuQLec2HJOJt2yO9RbzlRwSKsoCETa0OECYvqrb5Pfu4jNO6vYsrOKzZVVrPlkD++t2caWnVVU1zY9rWdOh/RoaMTCIrrcocnlbtmZZGUoOOQABYFIO7AvIPZ8lEFp6YAm27g72/fUsGVnFVt27mVTZVVsuYpNlXv3L8cjOLp2yqBjRjqdMtPolJlGVkZabDmdrIyIzm0kmUCDwMzKgPuANOBBd/9xo+0dgD8BpwCbgfHuvjLImkQSlZnRpWMGXTpm0D8/+7Dt4xkcjXXMaBwQjcOiqfX7ltMbrO+YmbZ/fx0z08hKTyMSUdC0pcCCwMzSgAeAi4AK4E0zm+rui+o1+xKw1d2PM7MJwE+A8UHVJJJKWhMcn+yqYnd1LburatlVVdvEck2T6zdVVrGrahd7quvYFWuzt6buiGvvmBENBWqryfnHS6SnGZlpEdLTjPRI5MByWoTM2Lr0NCMjLUJGbH1GJPo6PbYuI/aejEi9Nmn12kTqtUmLNGifFjHSIkbEIBIx0syImBGJEFsf/UmLbbMIpMVem9Vfbp8BF+QRwQhgmbsvBzCzKcA4oH4QjAO+F1t+EviVmZm7N+/PEhGJm/rBAYcPjuaqq3N2V0fDYk/sv7uqapoVNKsqPia/e1eq65ya2jqqa53q2jpqap2qmjp2VtXG1kfXVdfVUV3j1NQ1alt75GEUlP1hYU0sNwqVSIToshmRiHFRzxpKA6gpyCDoDayu97oCOO1gbdy9xsy2AXnApvqNzOwG4AaAgoICysvLAyq5bVRWVib8d4gn9ccBqdYXaUBu7Ge/jNhPNlRmVJOTs+0I9hipt+cD3J06h9rYT00d1Nb5geV92+q8weuaWJvaOnCgzqM/0eXoPj22ri62/cBrP7Dsjd6/v329feDUeVPv8f3LkZqaQP7/SIiTxe4+CZgEMHz4cC8tLQ23oFYqLy8n0b9DPKk/DlBfNKT+aCio/ogcvkmLrQH61nvdJ7auyTZmlg50IXrSWERE2kiQQfAmMMDM+ptZJjABmNqozVTgmtjyFcCLOj8gItK2Ahsaio353wjMIDpg95C7LzSzO4G33H0q8Hvgf81sGbCFaFiIiEgbCvQcgbtPA6Y1Wnd7veU9wOeDrEFERA4tyKEhERFJAAoCEZEUpyAQEUlxCgIRkRRniXa1ppltBFaFXUcr5dPo7ukUp/44QH3RkPqjodb0x9Hu3r2pDQkXBMnAzN5y9+Fh19FeqD8OUF80pP5oKKj+0NCQiEiKUxCIiKQ4BUE4JoVdQDuj/jhAfdGQ+qOhQPpD5whERFKcjghERFKcgkBEJMUpCNqQmfU1s5fMbJGZLTSzb4ZdU9jMLM3M3jGzZ8OuJWxm1tXMnjSzJWa22MxOD7umMJnZLbF/JwvMbLKZZYVdU1sxs4fMbIOZLai3rpuZzTKzpbH/HhWvz1MQtK0a4NvuXgyMBP7FzIpDrils3wQWh11EO3EfMN3dBwInkcL9Yma9gZuA4e4+iOij7FPpMfV/BMoarbsNmO3uA4DZsddxoSBoQ+6+1t3fji3vIPoPvXe4VYXHzPoAY4EHw64lbGbWBTiH6BwduHuVu38SalHhSwc6xmYv7AR8HHI9bcbd5xCdo6W+ccDDseWHgcvi9XkKgpCYWREwDHgj5FLC9Avg34jO453q+gMbgT/EhsoeNLPssIsKi7uvAf4b+AhYC2xz95nhVhW6AndfG1teBxTEa8cKghCYWQ7wF+Bmd98edj1hMLNLgA3uPjfsWtqJdOBk4DfuPgzYSRwP/RNNbPx7HNGA7AVkm9lV4VbVfsSm9I3btf8KgjZmZhlEQ+ARd38q7HpCdCZwqZmtBKYA55vZn8MtKVQVQIW77ztCfJJoMKSqC4EV7r7R3auBp4AzQq4pbOvNrCdA7L8b4rVjBUEbMjMjOga82N1/HnY9YXL3/3D3Pu5eRPQk4IvunrJ/8bn7OmC1mZ0QW3UBsCjEksL2ETDSzDrF/t1cQAqfPI+ZClwTW74G+Gu8dqwgaFtnAhOJ/vU7L/YzJuyipN34BvCImc0HhgI/Crec8MSOjJ4E3gbeI/q7KmUeN2Fmk4HXgRPMrMLMvgT8GLjIzJYSPWL6cdw+T4+YEBFJbToiEBFJcQoCEZEUpyAQEUlxCgIRkRSnIBARSXEKApFWMrOi+k+JFEk0CgIRkRSnIBCJIzM7JvbQuFPDrkWkudLDLkAkWcQeDzEFuNbd3w27HpHmUhCIxEd3os9+udzdU/kZQZKANDQkEh/biD4o7aywCxE5UjoiEImPKuCzwAwzq3T3R8MuSKS5FAQiceLuO2MT7syKhcHUsGsSaQ49fVREJMXpHIGISIpTEIiIpDgFgYhIilMQiIikOAWBiEiKUxCIiKQ4BYGISIr7f7IW6h/g0Cr9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "vec_random = np.vectorize(p_random)\n",
    "x = np.arange(1,11)\n",
    "y = (1-vec_random(x)) + vec_random(x)*0.5\n",
    "eps = y/(1-y)\n",
    "eps = np.log(eps)\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x,eps, label = 'Mode i')\n",
    "ax.set(xlabel='k', ylabel='eps',\n",
    "       title='Epsilon given k.')\n",
    "ax.grid()\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAo8klEQVR4nO3dd3xUVd7H8c+PEGqoASK9SO8lYEUCVsCCuq5ir+jjrrrPrj4qoliw7apr37WXpakIiCDFlkUUdSlCEgLSIUAooSWB9PP8kUHKUgLJnTsz+b5fr7yYyZ3c8+O8km9uzrlzjjnnEBGRyFPB7wJERMQbCngRkQilgBcRiVAKeBGRCKWAFxGJUAp4EZEIpYCXsGdm15jZrAOeOzNrHeQa+prZsmC2GWh3jZmdE+x2JTwo4CWoAoG018yyDvh4tTTndM6Ncc6dV1Y1nmAN3znn2vlZg8ihKvpdgJRLFznnvvK7CJFIpyt4CRlmdqOZfW9mr5rZLjNbamZnH3J8lZllmtlqM7vmgM/POcI5a5nZh2a21czWmtkIM6tw4NeZ2XNmtiNwzoFHqa+nmS0MtP+JmX1kZqMCxxLMLC3w+H4zm3DI175kZi8fUNM7ZrbJzDaY2SgzizqRmg5po0Pg9UNL8nqJfAp4CTWnACuBesBIYKKZ1TWz6sDLwEDnXA3gdOCXEpzvFaAW0AroB1wP3HRIe8sC7f0VeMfM7NCTmFklYBLwPlAXGAdceoQ2xwODzKxG4GujgN8DYwPH3wcKgNZAD+A84NbjremQ+noCM4G7nHPjjvZaKT8U8OKHyWa284CP2w44tgV40TmX75z7iOKgGxw4VgR0NrOqzrlNzrmUozUSCNargAedc5nOuTXA88B1B7xsrXPuLedcIfAB0BCIO8zpTqV4SPPlQG0TgZ8P165zbi2wgP2/AAYAe5xzP5pZHDAI+JNzLts5twX4e6DO461pn77AFOB659zUo7xOyhkFvPhhiHOu9gEfbx1wbIM7eAW8tUAj51w2cCVwB7DJzKaZWftjtFMPiA6c48DzNT7gefq+B865PYGHMYc5V6PD1Lb+KG2PBfYNlVzN/qv35oGaNu37BQe8ATQ4gZr2uQP4wTmXeJTXSDmkgJdQ0/iQ4YhmwEYA59xM59y5FF/RLgXeOszXH2gbkE9xqB54vg0nUNemw9TW9Civ/wRIMLMmFF/J7wv49UAuUO+AX3A1nXOdTqCmfe4AmpnZ30txDolACngJNQ2Au80s2syuADoAX5hZnJldEhiLzwWyKB6yOaLAEMfHwJNmVsPMmgN/BkafQF1zgULgj2ZW0cwuAfocpe2tQCLwHrDaOZca+PwmYBbwvJnVNLMKZnaymfU7gZr2yQQuAM4ys2dKcR6JMAp48cPnh9wHP+mAYz8BbSi++n4S+J1zLoPi79U/U3w1v53iCdP/KUFbdwHZwCpgDsVX0u8eb8HOuTzgMuAWYCdwLTCV4l82RzIWOIf9V+/7XA9UApYAO4AJFP9VcsKcczuBc4GBZvYEgJlNN7PhpTmvhDfThh8SKszsRuBW59yZftdSEmb2E/BP59x7ftcicji6ghcpITPrZ2YnBYZobgC6AjP8rkvkSPROVpGSa0fxmH51iod8fhcYUxcJSRqiERGJUBqiERGJUCE1RFOvXj3XokULv8solezsbKpXr+53GSFBfXEw9cfB1B/7laYv5s+fv805V/9wx0Iq4Fu0aMG8efP8LqNUEhMTSUhI8LuMkKC+OJj642Dqj/1K0xdmtvZIxzREIyISoRTwIiIRSgEvIhKhQmoM/nDy8/NJS0sjJyfH71JKpFatWqSmpgalrSpVqtCkSROio6OD0p6IhJeQD/i0tDRq1KhBixYtOMaeByEhMzOTGjVqeN6Oc46MjAzS0tJo2bKl5+2JSPgJ+SGanJwcYmNjwyLcg8nMiI2NDZu/bEQk+EI+4AGF+xGoX0TkaMIi4EVEItXPq7fzxao8T86tgC8BM+Paa6/97XlBQQH169fnwgsvPK7ztGjRgm3btpX49bfeeitLliwBYNCgQezcufO42hOR0LUjO4/7PlnE79+Yy7frC9iTV1DmbYT8JGsoqF69OsnJyezdu5eqVavy5Zdf0rhx42N/YSm9/fbbvz3+4osvPG9PRLznnOPTBRt46otUdu/N5/Z+regZnU61SmUfx7qCL6FBgwYxbdo0AMaNG8fQoUN/O7Z9+3aGDBlC165dGTBgAIsXLwYgIyOD8847j06dOnHrrbdy4Mqdo0ePpk+fPnTv3p3bb7+dwsLC/2ozISHht6UbjvfqX0RCz4otWQx960fu/WQRLWKrMfXuM3lwYAcqV/RmPi2sruAf+zyFJRt3l+k5OzaqyciLjr3f8VVXXcXjjz/OhRdeyOLFi7n55pv57rvvABg5ciQ9evRg8uTJTJ06leuvv55ffvmFxx57jDPPPJNHHnmEadOm8c477wCQmprKRx99xPfff090dDR33nknY8aM4frrry/T/5uIhIac/EJe/3YF//j3SqpGR/HUpV24qndTKlTw9kaJsAp4P3Xt2pU1a9Ywbtw4Bg0adNCxOXPm8OmnnwLQr18/MjIy2L17N7Nnz2bixIkADB48mDp16gDw9ddfM3/+fHr37g3A3r17adCgQRD/NyISLHOWb2PE5CTWZOxhSPdGPDS4I/VrVA5K22EV8CW50vbSxRdfzL333ktiYiIZGRknfB7nHDfccANPP/10GVYnIqFka2Yuo6Yt4bNfNtIithqjbzmFM9vUC2oNGoM/DjfffDMjR46kS5cuB32+b9++jBkzBoDvvvuOevXqUbNmTc466yzGjh0LwPTp09mxYwcAZ599NhMmTGDLli1A8Rj+2rVHXPFTRMJIUZFjzE9rOfv5RKYnpXP32W2Y8aezgh7uEGZX8H5r0qQJd9999399/tFHH+Xmm2+ma9euVK5cmQ8++AAoHpsfOnQonTp14vTTT6dZs2YAdOzYkVGjRnHeeedRVFREdHQ0r732Gs2bNw/q/0dEytbS9N0Mn5jEgnU7ObVVXUYN6ULrBjG+1aOAL4GsrKz/+lxCQsJvC/TXrVuXyZMnAwevRRMbG8usWbMOe84rr7ySK6+88qjtJiYm/vZ4zZo1x123iATHnrwCXvpqOW/PWU2tqtE8f0U3LuvZ2Pd3myvgRURK4evUzTzyWQobdu7lyvimPDCwPXWqV/K7LEABLyJyQtJ35fDY5ylMT06nTYMYPr79NPq0rOt3WQcJi4B3zvn+p04oOvCNUyISHIVFjg9+WMPzs5ZRUOS47/x23Na3FZUqht49KyEf8FWqVCEjI0NLBh9i33rwVapU8bsUkXJjcdpOhk9KInnDbvq1rc8Tl3SmWWw1v8s6opAP+CZNmpCWlsbWrVv9LqVEcnJygha6+3Z0EhFvZebk8/ysX/lw7hpiYyrz6tU9GNylYchfdIZ8wEdHR4fVjkWJiYn06NHD7zJEpAw455ienM5jn6ewJTOX605tzr3nt6NmlfDYJjPkA15ExA/rt+/hkc+S+XbZVjo2rMkb18XTvWltv8s6Lp4FvJm1Az464FOtgEeccy961aaISGnlFxbxzpzVvPjVr1QwY8TgDtx4egsqRoXeJOqxeBbwzrllQHcAM4sCNgCTvGpPRKS05q/dzkOTklmansm5HeN49OJONK5d1e+yTliwhmjOBlY657TgioiEnF178nlmxlLG/byORrWq8OZ1vTiv00l+l1VqFox7qc3sXWCBc+7VwxwbBgwDiIuL6zV+/HjP6/FSVlYWMTH+rT0RStQXB1N/HCwU+sM5x9xNhYxfmktWPpzbvCKXtq5EFY824DiS0vRF//795zvn4g93zPOAN7NKwEagk3Nu89FeGx8f7/btYBSuEhMTf1ujprxTXxxM/XEwv/tj9bZsHp6czJwV2+jWtDZPXdqZTo1q+VJLafrCzI4Y8MEYohlI8dX7UcNdRCQYcgsK+WfiKl5LXEHlqAo8cUknrj6lOVEe767kh2AE/FBgXBDaERE5qrkrM3hochKrtmZzYdeGPHJhRxrUjNx3g3sa8GZWHTgXuN3LdkREjmZ7dh5PTkvl0wVpNK1blfdv6k1Cu8jfJtPTgHfOZQOxXrYhInIkzjk+mZfGU9NTycop4M6Ek7lrQBuqVoryu7Sg0DtZRSQiLd+cyUOTkvl5zXZ6t6jDk5d2oW1cDb/LCioFvIhElJz8Ql79ZgVvzF5JtUoVefbyLlzRqykVInAS9VgU8CISMWb/upWHP0tmbcYeLuvRmOGDO1AvprLfZflGAS8iYW9LZg5PTE3l80UbaVWvOmNvPYXTW9fzuyzfKeBFJGwVFTnG/ryOZ2csJTe/iD+d04Y7+p1MlejyMYl6LAp4EQlLqZt2M3xSEgvX7eT0k2MZNaQzreprKYgDKeBFJKzsySvgxa+W886c1dSuGs0Lv+/GpT0ah/zuSn5QwItI2PhqyWZGTklhw869XNW7KQ8MbE/tapX8LitkKeBFJORt2rWXx6YsYUZKOm3jYvjkjtPo3aKu32WFPAW8iISswiLHBz+s4flZyygoctx3fjtu69uKShXDb3clPyjgRSQkLU7byfBJSSRv2E2/tvV54pLONIut5ndZYUUBLyIhJTMnn+dn/cqHc9cQG1OZV6/uweAuDTWJegIU8CISEpxzzEhO59HPU9iSmct1pzbn3vPbUbNKtN+lhS0FvIj4bv32PYycksI3S7fQoWFN3rgunu5Na/tdVthTwIuIb/ILi3h3zmpe/Go5ZjBicAduPL0FFaM0iVoWFPAi4osVOwp55pU5LE3P5NyOcTx6cSca167qd1kRRQEvIkG1a08+z85cyrifcjipFrxxXS/O73SS32VFJAW8iASFc44pizbyxNQlbM/O47zmFXn+5n7EVFYMeUU9KyKeW5uRzYjJyXy3fBtdm9Ti/Zv6sG35QoW7x9S7IuKZvIIi3py9kle+WUF0VAUeu7gT157anKgKRuJyv6uLfAp4EfHEz6u3M3xSEiu2ZDG4S0MeuagjcTWr+F1WueJpwJtZbeBtoDPggJudc3O9bFNE/LUjO49npi/lo3nraVy7Ku/d2Jv+7Rv4XVa55PUV/EvADOfc78ysEqCFJEQilHOOiQs28OQXqezam8/t/Vpxz9ltqFZJAwV+8aznzawWcBZwI4BzLg/I86o9EfHPyq1ZPDw5mR9WZtCjWW2eurQLHRrW9Luscs+cc96c2Kw78CawBOgGzAfucc5lH/K6YcAwgLi4uF7jx4/3pJ5gycrKIiZG24aB+uJQkdgfeYWOaavymbYqn+go+H3bSvRrWpEKJVgYLBL740SVpi/69+8/3zkXf7hjXgZ8PPAjcIZz7iczewnY7Zx7+EhfEx8f7+bNm+dJPcGSmJhIQkKC32WEBPXFwSKtP35YuY0Rk5JZtS2bi7s1YsSFHWhQo+STqJHWH6VRmr4wsyMGvJeDY2lAmnPup8DzCcADHrYnIkGQkZXLk1+kMnHBBprVrcaHN/fhrLb1/S5LDsOzgHfOpZvZejNr55xbBpxN8XCNiIShoiLHJ/PX8/T0pWTnFvDH/q3544DWVImO8rs0OQKvp7fvAsYE7qBZBdzkcXsi4oHlmzN5aFIyP6/ZTp8WdXny0s60iavhd1lyDJ4GvHPuF+CwY0MiEvpy8gt59ZsVvDF7JdUrV+Svl3fld72aUKGCdlcKB7pBVUQOa/avW3n4s2TWZuzhsp6NeWhQB2JjKvtdlhwHBbyIHGRLZg6jpqYyZdFGWtWrzthbT+H01vX8LktOgAJeRIDiSdSxP6/j2RlLyc0v4k/ntOGOfidrEjWMKeBFhKXpuxk+MYkF63ZyWqtYRl3amZPr601I4U4BL1KO7ckr4KWvl/P2d6upVTWa56/oxmU9G2MleCeqhD4FvEg59e3SLYyYnMyGnXu5Mr4pDwxsT53qlfwuS8qQAl6knNm8O4fHPk/hi6R0WjeI4ePbT6NPy7p+lyUeUMCLlBOFRY7RP67lbzOXkV9YxH3nt+O2vq2oVLGC36WJRxTwIuVAysZdDJ+YxKK0XfRtU49RQzrTPLa632WJxxTwIhEsO7eAF7/6lXe/X0OdatG8dFV3Lu7WSJOo5YQCXiRCfbVkM498lszGXTkM7dOMBy5oT61q0X6XJUGkgBeJMOm7cnh0SgozUtJpGxfDhKGnEd9Ck6jlkQJeJEIUFjn+NXcNz836lfzCIv7vgnbceqYmUcszBbxIBEjesIvhk5JYnLaLs9rWZ9QlnWkWqz3uyzsFvEgYy84t4IUvf+W971dTt3plXh7ag4u6NtQkqgAKeJGw9eWSzYwMTKJec0oz/u+C9tSqqklU2U8BLxJmNu3ay6NTUpiZspl2cTX49Ooe9GquSVT5bwp4kTBRWOT4cO4anpu5jELnuP+C9tzatyXRUZpElcNTwIuEgaS04knUpA276Ne2Pk9oElVKQAEvEsKycgt4YdavvP/DamJjKvPq1T0Y3EWTqFIyCniREDUrJZ2RU1JI3108iXrf+ZpEleOjgBcJMRt3Fk+izlqymfYn1eDVq3vSq3kdv8uSMORpwJvZGiATKAQKnHPxXrYnEs4KCov4YO5aXphVPIn6wMD23HKmJlHlxAXjCr6/c25bENoRCVtJabt4cNJikjfspn+7+jx+SWea1tUkqpSOhmhEfJSVW8Dzs5bxwQ9rqBdTmdev6cnAzidpElXKhDnnvDu52WpgB+CAN5xzbx7mNcOAYQBxcXG9xo8f71k9wZCVlUVMjHajB/XFoQ7tj/mbCxi9JI+duY4BzSpyeZtKVIsuP8Gu74/9StMX/fv3n3+k4W+vA76xc26DmTUAvgTucs7NPtLr4+Pj3bx58zyrJxgSExNJSEjwu4yQoL442L7+2LBzLyM/S+Gr1OJJ1Kcv60KPZuVvElXfH/uVpi/M7IgB7+kQjXNuQ+DfLWY2CegDHDHgRSJZYZHj7e9W8cKXv+IcDB/UnpvO0CSqeMezgDez6kAF51xm4PF5wONetScSyhan7eTxH3NYuzuVAe0b8PglnWhSR5Oo4i0vr+DjgEmByaKKwFjn3AwP2xMJOVm5BTw3cxkfzl1DzUrGP67pyQWaRJUg8SzgnXOrgG5enV8k1B34TtTrTm3OqdW2MrBLQ7/LknJEt0mKlLFNu4onUfe9E/W1a3rSs1kdEhMT/S5NyhkFvEgZOXBP1IKiIr0TVXyngBcpAykbdzF8UjKL1u/UnqgSMhTwIqWwJ6+AF79azjtzVlOnWjQvXdWdi7s10iSqhAQFvMgJ+nbZFkZMSmbDzr0M7dOU+y9oT+1qlfwuS+Q3CniR47QlM4fHP1/C1MWbaN0gho9vP40+LbUnqoQeBbxICRUVOcb9Zx3PTF9KbkERfz63Lbf3a0XlilF+lyZyWMcd8GZWAYhxzu32oB6RkPTr5kwenJjE/LU7OK1VLE9e2plW9bVQloS2EgW8mY0F7qB4447/ADXN7CXn3N+8LE7Ebzn5hbzyzXLe+PcqalSpyHNXdOPyno01iSphoaRX8B2dc7vN7BpgOvAAMB9QwEvEmrN8GyMmJ7EmYw+X92zCQ4M7ULe6JlElfJQ04KPNLBoYArzqnMs3M+/WGRbxUUZWLqOmpTJp4QZa1qvO2FtP4fTW9fwuS+S4lTTg3wDWAIuA2WbWHNAYvEQU5xyfzE/jqS9Syc4t4O4Brbmzf2uqRGsSVcJTiQLeOfcy8PIBn1prZv29KUkk+FZuzWL4xCR+Wr2d3i3q8NSlXWgTV8PvskRKpaSTrLHASOBMirffm0Px2u4Z3pUm4r3cgkL+kbiS179dSZXoCjx9WReujG9KhQqaRJXwV9IhmvEU78R0eeD5NcBHwDleFCUSDD+tymD4pCRWbs3m4m6NePjCjtSvUdnvskTKTEkDvqFz7okDno8ysyu9KEjEazv35PH0F0v5aN56mtatyvs39SahXQO/yxIpcyUN+FlmdhXwceD574CZ3pQk4g3nHJ/9spEnpi5h59587uh3Mvec3YaqlTSJKpGppAF/G3AP8K/A8ygg28xuB5xzrqYXxYmUlbUZ2YyYnMx3y7fRvWltRl/WhQ4N9W0rka2kAV+L4nH3ls65x82sGcXDNj95V5pI6eUXFvHm7FW8/PVyKkVV4IlLOnH1Kc2J0iSqlAMlDfjXgCJgAMV3z2QCnwK9PapLpNQWrtvBgxOTWJqeycDOJ/HoxZ2Iq1nF77JEgqakAX+Kc66nmS0EcM7tMDO9Z1tCUlZuAc/NXMYHc9dwUs0qvHV9POd2jPO7LJGgK2nA55tZFMX3wGNm9Sm+oj+mwNfNAzY45y48oSpFSuirJZt5+LNk0nfncMNpLbj3/HbEVNaq2FI+lfQ7/2VgEtDAzJ6k+C6aESX82nuAVEAzWuKZLbtzeOzzJUxL2kS7uBq8dk1Pejar43dZIr4q6VIFY8xsPnA2YMAQ51zqsb7OzJoAg4EngT+XplCRwykqcnw0bz1PfZFKbkER953fjmFntSI6qoLfpYn4zpzzblFIM5sAPA3UAO493BCNmQ0DhgHExcX1Gj9+vGf1BENWVhYxMdoIArzvi41ZRXyQksuyHUV0qFuBGzpV5qTqoRvs+t44mPpjv9L0Rf/+/ec75+IPd8yzwUkzuxDY4pybb2YJR3qdc+5N4E2A+Ph4l5BwxJeGhcTERML9/1BWvOqLvIIi/vnvlbw6dwVVK0Xx18s7c0V8k5DfhEPfGwdTf+znVV94Oft0BnCxmQ0CqlC8C9Ro59y1HrYpEW7+2u088GkSy7dkcVG3Rjyi9WNEjsizgHfOPQg8CBC4gr9X4S4nandOPn+bsYzRP62lUa2qvHdjb/q31/oxIkej+8ck5M1MSeeRz5LZmpnLTae35C/ntaW6bn0UOaag/JQ45xKBxGC0JZFj8+4cHvksmZkpm+nQsCZvXhdPt6a1/S5LJGzoMkhCTlGRY+zP63h2+lLyCou4/4L23Nq3pW59FDlOCngJKcs3Z/LgxCTmrd3BGa1jeerSLjSPre53WSJhSQEvISG3oJDXv13J64krqF65Is9d0Y3LezYO+VsfRUKZAl589/Pq7Tw4cTErt2YzpHsjRlzYkXoxuvVRpLQU8OKbXXvzeWb6Usb9vI4mdbR1nkhZU8BL0DnnmJGczsgpKWzLyuW2vi3533PbUq2Svh1FypJ+oiSoNu3ay8OTU/gqdTOdGtXknRt606VJLb/LEolICngJisIix+gf1/K3mcsoKCpi+KD23HxGSyrq1kcRzyjgxXPL0jN5YOJiFq7bSd829XhySBeaxVbzuyyRiKeAF8/kFTqem7mMf/57JTWrRvP3K7sxpLtufRQJFgW8eOKnVRk88v1e0ves4LKejRkxuCN1q2sbX5FgUsBLmcrMyefZGUsZ/eM66lc1/nVLH/q2qe93WSLlkgJeysy3S7fw0KQk0nfncOuZLeldZbPCXcRHCngpte3ZeTwxdQmTFm6gTYMYPv2f0+nRrA6JiVv8Lk2kXFPAywlzzjEtaRMjP0th19587jm7DXf2P5nKFaP8Lk1EUMDLCdq8O4eHJycza8lmujapxZjbTqH9STX9LktEDqCAl+PinOPjeesZNS2VvAK9YUkklCngpcTWZezhwUmL+X5FBqe0rMuzl3elRT2t1S4SqhTwckyFRY73f1jDczOXEVXBePLSzgzt3YwKFfSGJZFQpoCXo1q+OZP7P13MgnU76d+uPk9e2oVGtav6XZaIlIACXg4rv7CIfyau5JVvVlC9chQvXtmdS7o30jIDImHEs4A3syrAbKByoJ0JzrmRXrUnZScpbRf3TVjE0vRMLurWiJEXaYclkXDk5RV8LjDAOZdlZtHAHDOb7pz70cM2pRRy8gv5+1e/8tbsVdSvUZm3ro/n3I5xfpclIifIs4B3zjkgK/A0OvDhvGpPSuenVRk8MDGJ1duyGdqnKQ8M7ECtqtF+lyUipWDFOezRyc2igPlAa+A159z9h3nNMGAYQFxcXK/x48d7Vk8wZGVlERMT43cZJba3wPHJr3l8s66A+lWNmzpXpmNs2bwTNdz6wmvqj4OpP/YrTV/0799/vnMu/nDHPA343xoxqw1MAu5yziUf6XXx8fFu3rx5ntfjpcTERBISEvwuo0S+XbaFhyYmsWl3Djef0ZK/nFe2+6KGU18Eg/rjYOqP/UrTF2Z2xIAPyl00zrmdZvYtcAFwxICX4NgRWBxs4gGLg/VsVsfvskSkjHl5F019ID8Q7lWBc4FnvWpPjs05xxdJ6YyckszOPfncfXYb/qDFwUQilpdX8A2BDwLj8BWAj51zUz1sT45iy+4cRgQWB+vSuBb/uuUUOjTU4mAikczLu2gWAz28Or+UjHOOT+al8cS0JeQVFPHgwPbccqYWBxMpD/RO1gi2fvseHpyYxJwV2+gTWByspRYHEyk3FPARqLDI8eHcNfx1RvHiYKOGdObqPlocTKS8UcBHmJVbs7jvk0VaHExEFPCRoiiwpO+zM5ZStZIWBxMRBXxESNuxh/s+WczcVRkMaN+AZy7rQoOaVfwuS0R8poAPY/u2z3tiaioAf728K1fEN9FVu4gACviwtWV3Dg9OTOLrpVs4tVVd/va7bjStW83vskQkhCjgw9DUxRsZMTmZvXmFjLyoIzec1kJ3yIjIf1HAh5Ed2Xk8MiWFzxdtpFvT2rzw+26cXF+r8YnI4Sngw8S3S7fwf58uZueePO49ry139DtZ70YVkaNSwIe4rNwCRk1dwvj/rKddXA3ev6k3nRrV8rssEQkDCvgQNndlBvdNWMTGnXu5o9/J/O+5bbTyo4iUmAI+BOXkF/LXGct49/vVtIitxid3nEav5nX9LktEwowCPsQsWr+TP3/8Cyu3ZnP9ac15YGD7Mt1lSUTKDyVHiMgrKOKVb5bzeuJKGtSozOhbTuHMNvX8LktEwpgCPgQsS8/kzx//QsrG3VzeswmPXNSRWlWj/S5LRMKcAt5HhUWOt75bxQuzfqVm1Yq8cV0vzu90kt9liUiEUMD7ZM22bP7yySLmr93BBZ1O4slLOxMbU9nvskQkgijgg8w5x+gf1/LUF0uJjjIt6ysinlHAB9HGnXu5/9PFfLd8G2e1rc+zl3ehYS1txiEi3lDAB4FzjokLNvDo5ykUFjlGDenMNac001W7iHjKs4A3s6bAh0Ac4IA3nXMvedVeqNqWlcvwiUnMWrKZ3i3q8NwV3Wgeq42vRcR7Xl7BFwB/cc4tMLMawHwz+9I5t8TDNkPKjOR0HpqURGZOAcMHteeWM1sRpWV9RSRIPAt459wmYFPgcaaZpQKNgYgP+F1783l0SgqTFm6gc+OajPt9d9rG1fC7LBEpZ4IyBm9mLYAewE/BaM9PydsKeODvs9malcvdZ7fhrgGtidayviLiA3POeduAWQzwb+BJ59zEwxwfBgwDiIuL6zV+/HhP6/FKXqFj/NI8vllfQMPqxm1dK9OqVvle+TErK4uYGG1Iso/642Dqj/1K0xf9+/ef75yLP9wxTwPezKKBqcBM59wLx3p9fHy8mzdvnmf1eGXV1izuHLOApemZnN+iIi/dcg5Vost3uAMkJiaSkJDgdxkhQ/1xMPXHfqXpCzM7YsB7eReNAe8AqSUJ93A1dfFG7p+wmEoVK/DeTb2xTUsU7iISErwcHD4DuA4YYGa/BD4GedheUOUWFDLys2T+OHYh7U6qwbS7+9K/XQO/yxIR+Y2Xd9HMASLynsD12/fwh7ELWJy2i1vPbMn9A9trIlVEQo7eyXqcvlyymb98/AsO+Oe1vbigs1Z/FJHQpIAvofzCIp6buYw3Zq+iU6OavH5NT70jVURCmgK+BDbt2stdYxcyb+0Orj21GSMGd9REqoiEPAX8Mcz+dSt/+ugXcvILeemq7lzSvbHfJYmIlIgC/ggKixwvfb2cV75ZTpsGMbx+TS9aN9CbMkQkfCjgD2NrZi5/+mgh36/I4PKeTXhiSCeqVVJXiUh4UWod4qdVGdw1biG79ubz18u78vveTf0uSUTkhCjgA4qKHP+cvZLnZi6jeWx1Pri5Dx0a1vS7LBGRE6aAB3Zk5/GXTxbxzdItDO7akGcu60KNKtF+lyUiUirlPuAXrtvBH8cuZEtmDo9f0onrTm2urfREJCKU24B3zvHe92t4enoqcTWrMOGO0+nWtLbfZYmIlJlyGfC7c/K5f8Jipienc06HBjx/RXdqVdOQjIhElnIX8Ckbd/GHMQtYv2Mvwwe157a+rTQkIyIRqdwEvHOO8f9Zz8gpKdStVomPhp1KfIu6fpclIuKZchHw2bkFjJiczKSFG+jbph4vXtmd2JjKfpclIuKpiA/45Zsz+Z8xC1i5NYs/n9uWP/RvTVQFDcmISOSL6ICftDCN4ROTqV45itG3nMIZrev5XZKISNBEZMDn5Bfy2OcpjPt5PX1a1uWVoT2Iq1nF77JERIIq4gJ+zbZs7hyzgCWbdnNnwsn8+dy2VNR2eiJSDkVUwE9P2sR9ExYTVcF498Z4BrSP87skERHfRETA5xUU8fT0VN77fg3dm9bm1at70KRONb/LEhHxVdgH/K49+Vz/3s8sWr+Tm89oyQMD21OpooZkREQ8C3gzexe4ENjinOvsVTs1qlSked1q3HFWKwZ2aehVMyIiYcfLK/j3gVeBDz1sgwoVjJeH9vCyCRGRsOTZWIZzbjaw3avzi4jI0ZlzzruTm7UAph5tiMbMhgHDAOLi4nqNHz/es3qCISsri5gYbc4N6otDqT8Opv7YrzR90b9///nOufjDHfN9ktU59ybwJkB8fLxLSEjwt6BSSkxMJNz/D2VFfXEw9cfB1B/7edUXut1ERCRCKeBFRCKUZwFvZuOAuUA7M0szs1u8aktERP6bZ2PwzrmhXp1bRESOTUM0IiIRytPbJI+XmW0F1vpdRynVA7b5XUSIUF8cTP1xMPXHfqXpi+bOufqHOxBSAR8JzGzeke5JLW/UFwdTfxxM/bGfV32hIRoRkQilgBcRiVAK+LL3pt8FhBD1xcHUHwdTf+znSV9oDF5EJELpCl5EJEIp4EVEIpQCvgyYWVMz+9bMlphZipnd43dNocDMosxsoZlN9bsWP5lZbTObYGZLzSzVzE7zuyY/mdn/Bn5Oks1snJlV8bumYDKzd81si5klH/C5umb2pZktD/xbpyzaUsCXjQLgL865jsCpwB/MrKPPNYWCe4BUv4sIAS8BM5xz7YFulOM+MbPGwN1AfGCfiCjgKn+rCrr3gQsO+dwDwNfOuTbA14HnpaaALwPOuU3OuQWBx5kU/wA39rcqf5lZE2Aw8LbftfjJzGoBZwHvADjn8pxzO30tyn8VgapmVhGoBmz0uZ6gOsJud5cAHwQefwAMKYu2FPBlLLCLVQ/gJ59L8duLwP8BRT7X4beWwFbgvcBw1dtmVt3vovzinNsAPAesAzYBu5xzs/ytKiTEOec2BR6nA3FlcVIFfBkysxjgU+BPzrndftfjFzO7ENjinJvvdy0hoCLQE/iHc64HkE0Z/fkdjgJjy5dQ/IuvEVDdzK71t6rQ4orvXS+T+9cV8GXEzKIpDvcxzrmJftfjszOAi81sDTAeGGBmo/0tyTdpQJpzbt9fdBMoDvzy6hxgtXNuq3MuH5gInO5zTaFgs5k1BAj8u6UsTqqALwNmZhSPsaY6517wux6/OecedM41cc61oHgC7RvnXLm8SnPOpQPrzaxd4FNnA0t8LMlv64BTzaxa4OfmbMrxpPMBpgA3BB7fAHxWFidVwJeNM4DrKL5S/SXwMcjvoiRk3AWMMbPFQHfgKX/L8U/gL5kJwAIgieIMKldLFhxht7tngHPNbDnFf+U8UyZtaakCEZHIpCt4EZEIpYAXEYlQCngRkQilgBcRiVAKeBGRCKWAFzkKM2tx4Kp/IuFEAS8iEqEU8CIlZGatAguG9fa7FpGSqOh3ASLhILDUwHjgRufcIr/rESkJBbzIsdWneG2Qy5xz5XkdGQkzGqIRObZdFC+SdabfhYgcD13BixxbHnApMNPMspxzY/0uSKQkFPAiJeCcyw5sZPJlIOSn+F2TyLFoNUkRkQilMXgRkQilgBcRiVAKeBGRCKWAFxGJUAp4EZEIpYAXEYlQCngRkQj1/0+69f0bfudzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "vec_random = np.vectorize(p_random)\n",
    "x = np.arange(1,11)\n",
    "y = (1-vec_random(x,mode='ii')) + vec_random(x,mode='ii')*0.5\n",
    "eps = y/(1-y)\n",
    "eps = np.log(eps)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x,eps, label = 'Mode ii')\n",
    "ax.set(xlabel='k', ylabel='eps',\n",
    "       title='Epsilon given k.')\n",
    "ax.grid()\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3: How much money do you have on your bank account?\n",
    "\n",
    "To answer a question about a scalar value, we need a different mechanism. Here we use the Laplace mechanism, which works as follows:\n",
    "\n",
    "* Evaluate the *sensitivity* of the function $f$ of interest:\n",
    "$$ \\Delta f = \\max\\limits_{\\mathcal{D} \\sim \\mathcal{D}'} | f(\\mathcal{D}) - f(\\mathcal{D}') |, $$\n",
    "where $\\mathcal{D} \\sim \\mathcal{D}'$ denotes pairs of data sets that differ by a single sample.\n",
    "* Compute $\\mathcal{M}(\\mathcal{D}) = f(\\mathcal{D}) + \\frac{\\Delta f}{\\epsilon} \\eta$\n",
    "where $\\eta \\sim \\operatorname{Laplace}(0, 1)$\n",
    "\n",
    "In this exercise we apply the Laplace mechanism in locally differentially private mode, i.e. each individual is assumed to form a different single-observation data set that is protected using this mechanism.\n",
    "\n",
    "In order to estimate the mean of the distribution of data held by the users, we can use the following moment identities with $x$ denoting the user data and $z$ denoting the Laplace noise added for DP:\n",
    "$$ \\mathrm{E}[x + z] = \\mathrm{E}[x] + \\mathrm{E}[z], $$\n",
    "$$ \\mathrm{Var}[x + z] = \\mathrm{Var}[x] + \\mathrm{Var}[z]. $$\n",
    "\n",
    "i) What is the probability $Pr(\\mathcal{M}(\\mathcal{D})=s)$ for some outcome $s$? Also, find the ratio $\\frac{Pr(\\mathcal{M}(\\mathcal{D})=s)}{Pr(\\mathcal{M}(\\mathcal{D}')=s)}$.\n",
    "\n",
    "ii) (Optional) Verify that this mechanism satisfies $\\epsilon$-DP.\n",
    "\n",
    "iii) Assuming a population of independent respondents whose bank account balances follow the exponential distribution $\\mathrm{Exp}(\\lambda)$ with $\\lambda = 1/5000$. Simulate the mechanism for $n$ participants with $n = 10, 100, 1000, 10000$. How accurately can you estimate the mean and the standard deviation of the bank account balance distribution?\n",
    "Repeat each simulation 100 times to get an idea of the accuracy.\n",
    "\n",
    "Hint: you will need to set an upper bound on the allowed balances to get a bounded sensitivity. Larger values can be handled e.g. by setting them to the maximum. In order to guarantee DP, this bound should be selected independently of the actual data at hand. (Mechanisms that do not require such an a priori bound exist, but are more complicated. If you are interested, see e.g. https://dl.acm.org/citation.cfm?id=1250803.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>n</th>\n",
       "      <th>True Mean</th>\n",
       "      <th>True Std.</th>\n",
       "      <th>Noisy Mean</th>\n",
       "      <th>Noisy Std.</th>\n",
       "      <th>Estimated Std.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>10.0</td>\n",
       "      <td>5052.18</td>\n",
       "      <td>4303.96</td>\n",
       "      <td>5416.01</td>\n",
       "      <td>13773.56</td>\n",
       "      <td>3773.56</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>100.0</td>\n",
       "      <td>4996.74</td>\n",
       "      <td>4898.50</td>\n",
       "      <td>4966.94</td>\n",
       "      <td>14579.81</td>\n",
       "      <td>4579.81</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1000.0</td>\n",
       "      <td>4983.62</td>\n",
       "      <td>4988.29</td>\n",
       "      <td>5032.42</td>\n",
       "      <td>14932.69</td>\n",
       "      <td>4932.69</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>10000.0</td>\n",
       "      <td>4997.40</td>\n",
       "      <td>4991.35</td>\n",
       "      <td>4998.94</td>\n",
       "      <td>15011.66</td>\n",
       "      <td>5011.66</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(6)\n",
    "\n",
    "simulations = 100\n",
    "participants = [10, 100, 1000, 10000]\n",
    "inv_scale = 1/5000\n",
    "eps = 10\n",
    "upper_bound = 100000\n",
    "lower_bound = 0\n",
    "mu = sensitivity = upper_bound-lower_bound\n",
    "report = pd.DataFrame()\n",
    "\n",
    "for n in participants:\n",
    "    D = npr.exponential(scale = 1/inv_scale, size=(n,simulations))\n",
    "    np.clip(D,lower_bound,upper_bound)\n",
    "    avg_, std_ = D.mean(axis=0), D.std(axis=0)\n",
    "\n",
    "    eta = npr.laplace(loc=0.0, scale=1.0, size=(n,simulations))\n",
    "    D = npr.exponential(scale = 1/inv_scale, size=(n,simulations))\n",
    "    D += (mu/eps)*eta\n",
    "    np.clip(D,lower_bound,upper_bound)\n",
    "    avg, std = D.mean(axis=0), D.std(axis=0)\n",
    "    report = report.append(pd.Series([n,avg_.mean(),std_.mean(),avg.mean(),std.mean(),std.mean()-mu/eps]),ignore_index=True)\n",
    "    \n",
    "report.columns = ['n','True Mean', 'True Std.', 'Noisy Mean', 'Noisy Std.', 'Estimated Std.']\n",
    "display(HTML(report.round(2).to_html()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first row reports the difference between the mean of the actual data and the mean of 100 simulated added-noise data set means as well as its deviations.  \n",
    "The actual mean and the mean of 100 means is reported below.  \n",
    "The third row reports the two drifferent standard deviations in the same manner.  \n",
    "The fourth row shows the estimated standard deviation of the noisy dataset, where sensitivity and $\\epsilon$ information\n",
    "is used to gauge the noise variance. The mean of 100 data sets is reported. Furthermore the variation of the estimated standard deviation is reported as well.\n",
    "\n",
    "The mean approaches the true mean rather quick over multiple simulations, the bias standard deviation shows that only with n=10000, the value becomes somewhat stable and small. (Standard deviation reflects less than 10\\% of the actual value) As the noise just adds a fixed variance, the estimate is pretty accurate over multiple simulations right away, but for stability reasons needs a bit to approach acceptable levels (n=1000). Obviously, if we would not have any information regarding the sensitivity/epsilon used, we could not converge towards the true variation in the data (obviously?).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "i)  \n",
    "$$Pr(\\mathcal{M}(\\mathcal{D})=s) = \n",
    "Pr(f(\\mathcal{D}) + \\frac{\\Delta f}{\\epsilon} \\eta=s) = \n",
    "Pr(\\eta = \\frac{\\epsilon}{\\Delta f}(s - f(\\mathcal{D})), \\eta \\sim \\operatorname{Laplace}(0, 1)$$  \n",
    "$$Pr\\left(\\eta = \\frac{\\epsilon}{\\Delta f}(s - f(\\mathcal{D}))\\right) = \n",
    "{\\displaystyle {\\frac {1}{2}}\\exp \\left(-{\\frac {|\\epsilon}{\\Delta f}(s - f(\\mathcal{D}) |}\\right)}$$\n",
    "$$\\frac{Pr(\\eta = \\frac{\\epsilon}{\\Delta f}(s - f(\\mathcal{D}))}{Pr(\\eta = \\frac{\\epsilon}{\\Delta f}(s - f(\\mathcal{D'}))} =   \n",
    "\\frac{{\\displaystyle \\exp \\left(-|\\frac{\\epsilon}{\\Delta f}(s - f(\\mathcal{D})) |\\right)}}\n",
    "{{\\displaystyle \\exp \\left(-|\\frac{\\epsilon}{\\Delta f}(s - f(\\mathcal{D'})) |\\right)}} $$\n",
    "$$= {\\displaystyle \\exp \\left(-|\\frac{\\epsilon}{\\Delta f}(s - f(\\mathcal{D})) | + |\\frac{\\epsilon}{\\Delta f}(s - f(\\mathcal{D'})) |\\right)} $$\n",
    "$$= {\\displaystyle \\exp \\left(\\frac{\\epsilon}{\\Delta f}\\left(|s - f(\\mathcal{D'}) | - |s - f(\\mathcal{D}) |\\right)\\right)} $$\n",
    "$$ = {\\displaystyle \\exp \\left(\\frac{\\epsilon}{\\Delta f}\\left(f(\\mathcal{D}) - f(\\mathcal{D'}) \\right)\\right)}$$\n",
    "\n",
    "ii)\n",
    "\n",
    "$${\\displaystyle \\exp \\left(\\frac{\\epsilon}{\\Delta f}\\left(f(\\mathcal{D}) - f(\\mathcal{D'}) \\right)\\right)} \\leq \\exp(\\epsilon) $$  \n",
    "$$ \\Leftrightarrow {\\displaystyle           \\frac{\\epsilon}{\\Delta f}\\left(f(\\mathcal{D}) - f(\\mathcal{D'}) \\right)} \\leq \\epsilon$$  \n",
    "$$ \\Leftrightarrow {\\displaystyle           \\frac{\\left(f(\\mathcal{D}) - f(\\mathcal{D'}) \\right)}{\\Delta f}} \\leq 1$$  \n",
    "$$ \\Leftrightarrow {\\displaystyle           \\frac{\\left(f(\\mathcal{D}) - f(\\mathcal{D'}) \\right)}{\\max\\limits_{\\mathcal{D} \\sim \\mathcal{D}'} | f(\\mathcal{D}) - f(\\mathcal{D}') |}} \\leq 1$$  \n",
    "\n",
    "By definition the following is true:\n",
    "$$ {\\displaystyle f(\\mathcal{D}) - f(\\mathcal{D'}) \\leq \\max\\limits_{\\mathcal{D} \\sim \\mathcal{D}'} | f(\\mathcal{D}) - f(\\mathcal{D}') |}$$\n",
    "<div align='right'>$\\square$</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4: Does your age affect your account balance?\n",
    "\n",
    "We can apply the Laplace mechanism to fit a linear regression model to model the bank account balance as a function of age. We know that the respondents are of age 20-70 years and their account balances are at most 150000.\n",
    "\n",
    "Linear regression fit to scalar $x_i, y_i$ involves fitting the model\n",
    "$$ y_i = \\alpha + \\beta x_i + \\epsilon_i, $$\n",
    "where $\\beta$ is the regression coefficient and $\\alpha$ is the intercept. Assuming regression errors $\\epsilon_i$ are normally distributed, the maximum likelihood estimates of the parameters are\n",
    "$$ \\hat{\\beta}= \\frac{\\sum_{i = 1}^n (x_i - \\bar{x})(y_i - \\bar{y}) }{ \\sum_{i = 1}^n (x_i - \\bar{x})^2} \\\\\n",
    "   \\hat{\\alpha} = \\bar{y} - \\hat{\\beta} \\bar{x},$$\n",
    "where $\\bar{x} = \\frac{1}{n} \\sum_{i = 1}^n x_i$ and $\\bar{y} = \\frac{1}{n} \\sum_{i = 1}^n y_i$.\n",
    "\n",
    "i. Use the data in \"balances.csv\" loaded below, and release $x$ (the age) and $y$ (the bank account balance) using Laplace mechanism with total privacy loss $\\epsilon=1$ and fit the regression model to the noisy data. Repeat the process 100 times and evaluate how the model compares with one learned from clean data. Because you are releasing the composition of two queries, you will need to use $\\epsilon/2$ as the $\\epsilon$ for each query to reach total privacy loss of $\\epsilon$.\n",
    "\n",
    "ii. Design a two-step protocol that first estimates the means $\\bar{x}$ and $\\bar{y}$, and then uses these to release the terms $(x_i - \\bar{x}) (y_i - \\bar{y})$, $(x_i - \\bar{x})^2$ using the Laplace mechanism. Fit the regression model using the noisy information. Compare with the above. In order to make the results comparable, the total privacy loss should be $\\epsilon=1$. Because you are releasing the composition of four queries, you will need to use $\\epsilon/4$ as the $\\epsilon$ for each query to reach total privacy loss of $\\epsilon$.\n",
    "\n",
    "iii. Repeat the above cases with total privacy loss $\\epsilon = 10$ and compare the results.\n",
    "\n",
    "Hint: you can use the known range of values of $x$ and $y$ to bound the sensitivites. For the two-step protocol, you can also use the published private means $\\bar{x}$ and $\\bar{y}$ to evaluate the sensitivity for the second step, as long as you do not use any other features of the specific data set.\n",
    "\n",
    "Comment: a better solution would take into account the knowledge of the noise added. Developing algorithms like this is in many cases still an open research problem. Furthermore, the privacy budget can be split between different queries arbitrarily, not necessarily evenly like here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('../data/balances.csv', header=0, index_col=0)\n",
    "\n",
    "x = data.values[:,0]\n",
    "y = data.values[:,1]\n",
    "N = len(data)\n",
    "eps = 1\n",
    "#dF_x = x.max() - x.min()\n",
    "#dF_y = y.max() - y.min()\n",
    "dF_x = 50\n",
    "dF_y = 100000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "i. Use the data in \"balances.csv\" loaded below, and release $x$ (the age) and $y$ (the bank account balance) using Laplace mechanism with total privacy loss $\\epsilon=1$ and fit the regression model to the noisy data. Repeat the process 100 times and evaluate how the model compares with one learned from clean data. Because you are releasing the composition of two queries, you will need to use $\\epsilon/2$ as the $\\epsilon$ for each query to reach total privacy loss of $\\epsilon$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "b, a = [], []\n",
    "np.random.seed(1)\n",
    "for _ in range(100):\n",
    "    x_prv = x + npr.laplace(0,1,size=N)*(dF_x/(eps/2))\n",
    "    y_prv = y + npr.laplace(0,1,size=N)*(dF_y/(eps/2))\n",
    "    #y_prv = np.clip(y_prv, -150000, 150000)\n",
    "    b_hat = sum((x_prv-x_prv.mean())*(y_prv-y_prv.mean()))/sum((x_prv-x_prv.mean())**2)\n",
    "    a_hat = y_prv.mean() - b_hat*x_prv.mean()\n",
    "    b.append(b_hat)\n",
    "    a.append(a_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated parameters on privatized data (100 runs):\n",
      "beta: 15.14 with sd: 62.21139294190902\n",
      "alpha: 49514.46 with sd: 10198.189293584526\n"
     ]
    }
   ],
   "source": [
    "b_hat = np.mean(b)\n",
    "a_hat = np.mean(a)\n",
    "b_hat, a_hat \n",
    "print(f'Estimated parameters on privatized data (100 runs):')\n",
    "print(f'beta: {round(b_hat,2)} with sd: {np.std(b)}')\n",
    "print(f'alpha: {round(a_hat,2)} with sd: {np.std(a)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated parameters on disclosed data.\n",
      "beta: 1146.96\n",
      "alpha: -1776.36\n"
     ]
    }
   ],
   "source": [
    "b_hat = sum((x-x.mean())*(y-y.mean()))/sum((x-x.mean())**2)\n",
    "a_hat = y.mean() - b_hat*x.mean()\n",
    "print(f'Estimated parameters on disclosed data.')\n",
    "print(f'beta: {round(b_hat,2)}')\n",
    "print(f'alpha: {round(a_hat,2)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ii. Design a two-step protocol that first estimates the means $\\bar{x}$ and $\\bar{y}$, and then uses these to release the terms $(x_i - \\bar{x}) (y_i - \\bar{y})$, $(x_i - \\bar{x})^2$ using the Laplace mechanism. Fit the regression model using the noisy information. Compare with the above. In order to make the results comparable, the total privacy loss should be $\\epsilon=1$. Because you are releasing the composition of four queries, you will need to use $\\epsilon/4$ as the $\\epsilon$ for each query to reach total privacy loss of $\\epsilon$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_ = x.mean()\n",
    "y_ = y.mean()\n",
    "data = data.assign(cvar=lambda x: (x.Age-x_)*(x.Balance-y_), \n",
    "                   xvar=lambda x: (x.Age-x_)**2)\n",
    "cvar = data.values[:,2]\n",
    "xvar = data.values[:,3]\n",
    "eps = 1\n",
    "dF_cvar = cvar.max() - cvar.min()\n",
    "dF_xvar = xvar.max() - xvar.min()\n",
    "#dF_x = x.max() - x.min()\n",
    "#dF_y = y.max() - y.min()\n",
    "dF_x = 50\n",
    "dF_y = 100000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "b, a = [], []\n",
    "np.random.seed(111)\n",
    "for _ in range(100):\n",
    "    cvar_prv = cvar + npr.laplace(0,1,size=N)*(dF_cvar/(eps/4))\n",
    "    xvar_prv = xvar + npr.laplace(0,1,size=N)*(dF_xvar/(eps/4))\n",
    "    x_prv = x + npr.laplace(0,1,size=N)*(dF_x/(eps/4))\n",
    "    y_prv = y + npr.laplace(0,1,size=N)*(dF_y/(eps/4))\n",
    "    y_prv = np.clip(y_prv, -150000, 150000)\n",
    "    b_hat = sum(cvar_prv)/sum(xvar_prv)\n",
    "    a_hat = y_prv.mean() - b_hat*x_prv.mean()\n",
    "    b.append(b_hat)\n",
    "    a.append(a_hat)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta: 825.016762953705 with sd: 6763.371942541809\n",
      "alpha: -12245.576410486377 with sd: 356387.0963511269\n",
      "Empirical estimates\n",
      "beta: 1146.96\n",
      "alpha: -1776.36\n"
     ]
    }
   ],
   "source": [
    "b_hat = np.mean(b)\n",
    "a_hat = np.mean(a)\n",
    "print(f'beta: {b_hat} with sd: {np.std(b)}')\n",
    "print(f'alpha: {a_hat} with sd: {np.std(a)}')\n",
    "print(f'Empirical estimates')\n",
    "print(f'beta: 1146.96')\n",
    "print(f'alpha: -1776.36')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both estimations are far off compared to the original data, yet when adding noise to the Age and Balance directly and doing the estimation on those like in case 1, the estimations on the beta are far too concentrated and most is captured by the alpha.\n",
    "With the 2-step approach, adding the noise to the (co-)variance, both estimates are fairly unconcentrated which at least does not give a false sense of accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "iii. Repeat the above cases with total privacy loss $\\epsilon = 10$ and compare the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "private beta: 911.99 with sd: 59.88864771668498\n",
      "private alpha: 8804.88 with sd: 2777.5948511889064\n",
      "`true` beta: 1146.96\n",
      "`true` alpha: -1776.36\n"
     ]
    }
   ],
   "source": [
    "data = pd.read_csv('../data/balances.csv', header=0, index_col=0)\n",
    "\n",
    "x = data.values[:,0]\n",
    "y = data.values[:,1]\n",
    "N = len(data)\n",
    "eps = 10\n",
    "#df_x = x.max() - x.min()\n",
    "#df_y = y.max() - y.min()\n",
    "dF_x = 25\n",
    "dF_y = 100000\n",
    "b, a = [], []\n",
    "np.random.seed(42)\n",
    "for _ in range(100):\n",
    "    x_prv = x + npr.laplace(0,1,size=N)*(dF_x/(eps/2))\n",
    "    y_prv = y + npr.laplace(0,1,size=N)*(dF_y/(eps/2))\n",
    "    #y_prv = np.clip(y_prv, -150000, 150000)\n",
    "    b_hat = sum((x_prv-x_prv.mean())*(y_prv-y_prv.mean()))/sum((x_prv-x_prv.mean())**2)\n",
    "    a_hat = y_prv.mean() - b_hat*x_prv.mean()\n",
    "    b.append(b_hat)\n",
    "    a.append(a_hat)\n",
    "b_hat = np.mean(b)\n",
    "a_hat = np.mean(a)\n",
    "b_hat, a_hat \n",
    "print(f'private beta: {round(b_hat,2)} with sd: {np.std(b)}')\n",
    "print(f'private alpha: {round(a_hat,2)} with sd: {np.std(a)}')\n",
    "b_hat = sum((x-x.mean())*(y-y.mean()))/sum((x-x.mean())**2)\n",
    "a_hat = y.mean() - b_hat*x.mean()\n",
    "print(f'`true` beta: {round(b_hat,2)}')\n",
    "print(f'`true` alpha: {round(a_hat,2)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta: 1153.3546206290018 with sd: 260.4937937831631\n",
      "alpha: -2122.4210251790846 with sd: 11923.126886969232\n",
      "`true` beta: 1146.96\n",
      "`true` alpha: -1776.36\n"
     ]
    }
   ],
   "source": [
    "eps = 10\n",
    "df_x = 25\n",
    "df_y = 100000\n",
    "dF_cvar = dF_x * df_y\n",
    "dF_xvar = dF_x**2\n",
    "\n",
    "b, a = [], []\n",
    "np.random.seed(1)\n",
    "for _ in range(100):\n",
    "    x_prv = x + npr.laplace(0,1,size=N)*(dF_x/(eps/4))\n",
    "    y_prv = y + npr.laplace(0,1,size=N)*(dF_y/(eps/4))\n",
    "    y_prv = np.clip(y_prv, 0, 100000)\n",
    "    x_ = x_prv.mean()\n",
    "    y_ = y_prv.mean()\n",
    "    cvar = (x-x_)*(y-y_)\n",
    "    xvar = (x-x_)**2\n",
    "    cvar_prv = cvar + npr.laplace(0,1,size=N)*(dF_cvar/(eps/4))\n",
    "    xvar_prv = xvar + npr.laplace(0,1,size=N)*(dF_xvar/(eps/4))\n",
    "    b_hat = sum(cvar_prv)/sum(xvar_prv)\n",
    "    a_hat = y_prv.mean() - b_hat*x_prv.mean()\n",
    "    b.append(b_hat)\n",
    "    a.append(a_hat)\n",
    "\n",
    "b_hat = np.mean(b)\n",
    "a_hat = np.mean(a)\n",
    "print(f'beta: {b_hat} with sd: {np.std(b)}')\n",
    "print(f'alpha: {a_hat} with sd: {np.std(a)}')\n",
    "print(f'`true` beta: 1146.96')\n",
    "print(f'`true` alpha: -1776.36')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With more privacy budget, we can observe two things:\n",
    "a) When adding direct noise prior to the estimation (case 1), most of the noise still ends up in the intercept $\\alpha$ and artifically surpresses a directional effect through $\\beta$. This only slowly dissipates, ie with an almost not privatized data set such as $\\epsilon = 100$, the estimates approach the true values.\n",
    "\n",
    "b) Adding the noise to the (co-)variances (case 2), adding privacy budget translates directly into proportionally better estimates. It does not seem to prefer either intercept or slope right off the bat and seems to be the prefered method, when using regression on privatized data.\n",
    "The primary benefit seems to be less precision."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}