{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Canonical Noise Mechanism"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The canonical noise mechanism ([see AV23](https://projecteuclid.org/journals/annals-of-statistics/volume-51/issue-2/Canonical-noise-distributions-and-private-hypothesis-tests/10.1214/23-AOS2259.short)) can privatize any float-valued statistic with finite sensitivity.\n",
    "\n",
    "Under $(\\epsilon, \\delta)$-DP,\n",
    "the canonical noise distribution follows the Tulap distribution,\n",
    "which is a combination of discrete laplace noise and continuous uniform noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKxlJREFUeJzt3X90VPWd//FXEshAIJOYQBKi4Yc/Skz5IQ0SZmVdKjEhRFaW2K2VYnQ5cOQMLJAWIZaC4I8gtWJlEWzrAlZSLF1/LCBgiCVoCYKhlF8ShcIBDZOwUjL8OCSQ3O8fPczX0YCZZMJ8ZvJ8nHPPmbn3M/e+PyTMvPK5n3snzLIsSwAAAAYJD3QBAAAAX0dAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYp0OgC2iJxsZGVVVVKTo6WmFhYYEuBwAANINlWTp79qySk5MVHn7tMZKgDChVVVVKSUkJdBkAAKAFTpw4oZtuuumabYIyoERHR0v6RwftdnuAqwEAAM3hdruVkpLi+Ry/lqAMKFdO69jtdgIKAABBpjnTM5gkCwAAjENAAQAAxmlVQFm4cKHCwsI0ffp0z7qLFy/K6XQqPj5eXbt2VV5enqqrq71ed/z4ceXm5ioqKkoJCQmaOXOmLl++3JpSAABACGlxQNm1a5deeeUVDRgwwGv9jBkztG7dOq1du1ZlZWWqqqrS2LFjPdsbGhqUm5ur+vp6bd++XatWrdLKlSs1d+7clvcCAACElBYFlHPnzmncuHH6zW9+oxtuuMGzvra2Vq+++qpeeOEF3XPPPUpPT9eKFSu0fft27dixQ5L03nvv6eDBg3r99dd1xx13KCcnR0899ZSWLl2q+vp6//QKAAAEtRYFFKfTqdzcXGVmZnqtr6io0KVLl7zWp6amqmfPniovL5cklZeXq3///kpMTPS0yc7Oltvt1oEDB5o8Xl1dndxut9cCAABCl8+XGa9Zs0a7d+/Wrl27vrHN5XIpMjJSsbGxXusTExPlcrk8bb4aTq5sv7KtKUVFRZo/f76vpQIAgCDl0wjKiRMnNG3aNK1evVqdOnVqq5q+obCwULW1tZ7lxIkT1+3YAADg+vMpoFRUVKimpkbf+9731KFDB3Xo0EFlZWV66aWX1KFDByUmJqq+vl5nzpzxel11dbWSkpIkSUlJSd+4qufK8yttvs5ms3luysbN2QAACH0+BZQRI0Zo37592rNnj2cZPHiwxo0b53ncsWNHlZaWel5TWVmp48ePy+FwSJIcDof27dunmpoaT5uSkhLZ7XalpaX5qVsAACCY+TQHJTo6Wv369fNa16VLF8XHx3vWT5gwQQUFBYqLi5PdbtfUqVPlcDg0dOhQSVJWVpbS0tI0fvx4LVq0SC6XS3PmzJHT6ZTNZvNTtwAAQDDz+3fxLF68WOHh4crLy1NdXZ2ys7P18ssve7ZHRERo/fr1mjx5shwOh7p06aL8/HwtWLDA36UAAIAgFWZZlhXoInzldrsVExOj2tpa5qMAABAkfPn85rt4AACAcfx+igcAes/e4Hl8bGFuACsBEKwYQQEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjMON2gAYixu+Ae0XIygAAMA4BBQAAGAcTvEA8JuvnpIBgNZgBAUAABiHgAIAAIxDQAEAAMYhoAAwRu/ZG5jHAkASAQUAABiIgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA5fFggg4L5+czZu1gaAERQAAGAcAgoAADAOp3gAtBqnZAD4GyMoAADAOAQUAABgHAIKgOuq9+wNnBIC8K0IKAAAwDg+BZRly5ZpwIABstvtstvtcjgc2rhxo2f78OHDFRYW5rU89thjXvs4fvy4cnNzFRUVpYSEBM2cOVOXL1/2T28AAEBI8OkqnptuukkLFy7UbbfdJsuytGrVKt1///36y1/+ou9+97uSpIkTJ2rBggWe10RFRXkeNzQ0KDc3V0lJSdq+fbtOnjyphx9+WB07dtSzzz7rpy4BAIBg51NAGT16tNfzZ555RsuWLdOOHTs8ASUqKkpJSUlNvv69997TwYMHtWXLFiUmJuqOO+7QU089pVmzZunJJ59UZGRkC7sBAABCSYvnoDQ0NGjNmjU6f/68HA6HZ/3q1avVrVs39evXT4WFhbpw4YJnW3l5ufr376/ExETPuuzsbLndbh04cOCqx6qrq5Pb7fZaAABA6PL5Rm379u2Tw+HQxYsX1bVrV7311ltKS0uTJD300EPq1auXkpOTtXfvXs2aNUuVlZV68803JUkul8srnEjyPHe5XFc9ZlFRkebPn+9rqQAAIEj5HFD69u2rPXv2qLa2Vn/84x+Vn5+vsrIypaWladKkSZ52/fv3V48ePTRixAgdOXJEt9xyS4uLLCwsVEFBgee52+1WSkpKi/cHAADM5vMpnsjISN16661KT09XUVGRBg4cqF/96ldNts3IyJAkHT58WJKUlJSk6upqrzZXnl9t3ook2Ww2z5VDVxYAABC6Wn0flMbGRtXV1TW5bc+ePZKkHj16SJIcDof27dunmpoaT5uSkhLZ7XbPaSIAAACfTvEUFhYqJydHPXv21NmzZ1VcXKytW7dq8+bNOnLkiIqLizVq1CjFx8dr7969mjFjhu6++24NGDBAkpSVlaW0tDSNHz9eixYtksvl0pw5c+R0OmWz2dqkgwAAIPj4FFBqamr08MMP6+TJk4qJidGAAQO0efNm3XvvvTpx4oS2bNmiF198UefPn1dKSory8vI0Z84cz+sjIiK0fv16TZ48WQ6HQ126dFF+fr7XfVMAtA/c7h7AtfgUUF599dWrbktJSVFZWdm37qNXr1569913fTksAABoZ/guHgAAYBwCCgAAMA4BBQAAGIeAAqBN9Z69gQmxAHxGQAEAAMYhoAAAAOP4/F08ANASnOYB4AtGUAAEBeayAO0LAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOB0CXQCA4NV79oZAlwAgRDGCAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADj+BRQli1bpgEDBshut8tut8vhcGjjxo2e7RcvXpTT6VR8fLy6du2qvLw8VVdXe+3j+PHjys3NVVRUlBISEjRz5kxdvnzZP70BAAAhwaeActNNN2nhwoWqqKjQxx9/rHvuuUf333+/Dhw4IEmaMWOG1q1bp7Vr16qsrExVVVUaO3as5/UNDQ3Kzc1VfX29tm/frlWrVmnlypWaO3euf3sFAACCWphlWVZrdhAXF6df/OIXeuCBB9S9e3cVFxfrgQcekCQdOnRIt99+u8rLyzV06FBt3LhR9913n6qqqpSYmChJWr58uWbNmqVTp04pMjKyWcd0u92KiYlRbW2t7HZ7a8oH0AqB+DbjYwtzr/sxAfiHL5/fLZ6D0tDQoDVr1uj8+fNyOByqqKjQpUuXlJmZ6WmTmpqqnj17qry8XJJUXl6u/v37e8KJJGVnZ8vtdntGYZpSV1cnt9vttQAAgNDlc0DZt2+funbtKpvNpscee0xvvfWW0tLS5HK5FBkZqdjYWK/2iYmJcrlckiSXy+UVTq5sv7LtaoqKihQTE+NZUlJSfC0bAAAEEZ8DSt++fbVnzx599NFHmjx5svLz83Xw4MG2qM2jsLBQtbW1nuXEiRNtejwAABBYHXx9QWRkpG699VZJUnp6unbt2qVf/epX+uEPf6j6+nqdOXPGaxSlurpaSUlJkqSkpCTt3LnTa39XrvK50qYpNptNNpvN11IBAECQavV9UBobG1VXV6f09HR17NhRpaWlnm2VlZU6fvy4HA6HJMnhcGjfvn2qqanxtCkpKZHdbldaWlprSwEAACHCpxGUwsJC5eTkqGfPnjp79qyKi4u1detWbd68WTExMZowYYIKCgoUFxcnu92uqVOnyuFwaOjQoZKkrKwspaWlafz48Vq0aJFcLpfmzJkjp9PJCAkAAPDwKaDU1NTo4Ycf1smTJxUTE6MBAwZo8+bNuvfeeyVJixcvVnh4uPLy8lRXV6fs7Gy9/PLLntdHRERo/fr1mjx5shwOh7p06aL8/HwtWLDAv70CAABBrdX3QQkE7oMCmIH7oADwxXW5DwoAAEBbIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOP4FFCKiop05513Kjo6WgkJCRozZowqKyu92gwfPlxhYWFey2OPPebV5vjx48rNzVVUVJQSEhI0c+ZMXb58ufW9AQAAIaGDL43LysrkdDp155136vLly3riiSeUlZWlgwcPqkuXLp52EydO1IIFCzzPo6KiPI8bGhqUm5urpKQkbd++XSdPntTDDz+sjh076tlnn/VDlwAAQLDzKaBs2rTJ6/nKlSuVkJCgiooK3X333Z71UVFRSkpKanIf7733ng4ePKgtW7YoMTFRd9xxh5566inNmjVLTz75pCIjI1vQDQAAEEpaNQeltrZWkhQXF+e1fvXq1erWrZv69eunwsJCXbhwwbOtvLxc/fv3V2Jiomdddna23G63Dhw40ORx6urq5Ha7vRYAABC6fBpB+arGxkZNnz5dd911l/r16+dZ/9BDD6lXr15KTk7W3r17NWvWLFVWVurNN9+UJLlcLq9wIsnz3OVyNXmsoqIizZ8/v6WlAgCAINPigOJ0OrV//359+OGHXusnTZrkedy/f3/16NFDI0aM0JEjR3TLLbe06FiFhYUqKCjwPHe73UpJSWlZ4QAAwHgtOsUzZcoUrV+/Xn/605900003XbNtRkaGJOnw4cOSpKSkJFVXV3u1ufL8avNWbDab7Ha71wIAAEKXTwHFsixNmTJFb731lt5//3316dPnW1+zZ88eSVKPHj0kSQ6HQ/v27VNNTY2nTUlJiex2u9LS0nwpBwAAhCifTvE4nU4VFxfrnXfeUXR0tGfOSExMjDp37qwjR46ouLhYo0aNUnx8vPbu3asZM2bo7rvv1oABAyRJWVlZSktL0/jx47Vo0SK5XC7NmTNHTqdTNpvN/z0EAABBx6cRlGXLlqm2tlbDhw9Xjx49PMsbb7whSYqMjNSWLVuUlZWl1NRU/eQnP1FeXp7WrVvn2UdERITWr1+viIgIORwO/fjHP9bDDz/sdd8UAADQvvk0gmJZ1jW3p6SkqKys7Fv306tXL7377ru+HBoAALQjfBcPAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABjHp4BSVFSkO++8U9HR0UpISNCYMWNUWVnp1ebixYtyOp2Kj49X165dlZeXp+rqaq82x48fV25urqKiopSQkKCZM2fq8uXLre8NAAAICT4FlLKyMjmdTu3YsUMlJSW6dOmSsrKydP78eU+bGTNmaN26dVq7dq3KyspUVVWlsWPHerY3NDQoNzdX9fX12r59u1atWqWVK1dq7ty5/usVAAAIamGWZVktffGpU6eUkJCgsrIy3X333aqtrVX37t1VXFysBx54QJJ06NAh3X777SovL9fQoUO1ceNG3XfffaqqqlJiYqIkafny5Zo1a5ZOnTqlyMjIbz2u2+1WTEyMamtrZbfbW1o+gFbqPXvDdT/msYW51/2YAPzDl8/vVs1Bqa2tlSTFxcVJkioqKnTp0iVlZmZ62qSmpqpnz54qLy+XJJWXl6t///6ecCJJ2dnZcrvdOnDgQJPHqaurk9vt9loAAEDoanFAaWxs1PTp03XXXXepX79+kiSXy6XIyEjFxsZ6tU1MTJTL5fK0+Wo4ubL9yramFBUVKSYmxrOkpKS0tGwAABAEWhxQnE6n9u/frzVr1vizniYVFhaqtrbWs5w4caLNjwkAAAKnQ0teNGXKFK1fv17btm3TTTfd5FmflJSk+vp6nTlzxmsUpbq6WklJSZ42O3fu9Nrflat8rrT5OpvNJpvN1pJSAQBAEPJpBMWyLE2ZMkVvvfWW3n//ffXp08dre3p6ujp27KjS0lLPusrKSh0/flwOh0OS5HA4tG/fPtXU1HjalJSUyG63Ky0trTV9AQAAIcKnERSn06ni4mK98847io6O9swZiYmJUefOnRUTE6MJEyaooKBAcXFxstvtmjp1qhwOh4YOHSpJysrKUlpamsaPH69FixbJ5XJpzpw5cjqdjJIAAABJPgaUZcuWSZKGDx/utX7FihV65JFHJEmLFy9WeHi48vLyVFdXp+zsbL388suethEREVq/fr0mT54sh8OhLl26KD8/XwsWLGhdTwAAQMho1X1QAoX7oABm4D4oAHxx3e6DAgAA0BYIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcToEugAAaI3eszd4Hh9bmBvASgD4k88jKNu2bdPo0aOVnJyssLAwvf32217bH3nkEYWFhXktI0eO9Gpz+vRpjRs3Tna7XbGxsZowYYLOnTvXqo4AAIDQ4XNAOX/+vAYOHKilS5detc3IkSN18uRJz/L73//ea/u4ceN04MABlZSUaP369dq2bZsmTZrke/UAACAk+XyKJycnRzk5OddsY7PZlJSU1OS2Tz75RJs2bdKuXbs0ePBgSdKSJUs0atQoPf/880pOTva1JAAAEGLaZJLs1q1blZCQoL59+2ry5Mn68ssvPdvKy8sVGxvrCSeSlJmZqfDwcH300UdN7q+urk5ut9trAYCW6D17g9e8FQBm8ntAGTlypF577TWVlpbqueeeU1lZmXJyctTQ0CBJcrlcSkhI8HpNhw4dFBcXJ5fL1eQ+i4qKFBMT41lSUlL8XTYAADCI36/iefDBBz2P+/fvrwEDBuiWW27R1q1bNWLEiBbts7CwUAUFBZ7nbrebkAK0U1dGP7hiBwhtbX4flJtvvlndunXT4cOHJUlJSUmqqanxanP58mWdPn36qvNWbDab7Ha71wIAAEJXmweUzz//XF9++aV69OghSXI4HDpz5owqKio8bd5//301NjYqIyOjrcsBAABBwOdTPOfOnfOMhkjS0aNHtWfPHsXFxSkuLk7z589XXl6ekpKSdOTIET3++OO69dZblZ2dLUm6/fbbNXLkSE2cOFHLly/XpUuXNGXKFD344INcwQMAACS1YATl448/1qBBgzRo0CBJUkFBgQYNGqS5c+cqIiJCe/fu1b/+67/qO9/5jiZMmKD09HR98MEHstlsnn2sXr1aqampGjFihEaNGqVhw4bp17/+tf96BQAAgprPIyjDhw+XZVlX3b558+Zv3UdcXJyKi4t9PTQAA3CJLoDrgS8LBICr4J4pQODwZYEAPPjiPQCmYAQFAAAYh4ACAACMQ0ABAADGIaAAAADjEFCAdoyrVPg3AExFQAEAAMYhoAAAAONwHxQATfr6aQ/uiwLgemIEBQAAGIcRFCAEXBntuNYoB3eJBRBMGEEBAADGYQQFCCGMkgAIFQQUAEHJl3uXcJ8TIPgQUABAjD4BpmEOCgAAMA4BBQhR3MIdQDDjFA/QDrUkuARD2Pn65dYtrTkY+gqEOgIKgHaJEAKYjVM8AADAOIygACGOkQIAwYiAAiDkEMqA4EdAAcAHOgDjMAcFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIzjc0DZtm2bRo8ereTkZIWFhentt9/22m5ZlubOnasePXqoc+fOyszM1GeffebV5vTp0xo3bpzsdrtiY2M1YcIEnTt3rlUdAQAAocPngHL+/HkNHDhQS5cubXL7okWL9NJLL2n58uX66KOP1KVLF2VnZ+vixYueNuPGjdOBAwdUUlKi9evXa9u2bZo0aVLLewEAAEKKz18WmJOTo5ycnCa3WZalF198UXPmzNH9998vSXrttdeUmJiot99+Ww8++KA++eQTbdq0Sbt27dLgwYMlSUuWLNGoUaP0/PPPKzk5uRXdAQAAocCvc1COHj0ql8ulzMxMz7qYmBhlZGSovLxcklReXq7Y2FhPOJGkzMxMhYeH66OPPmpyv3V1dXK73V4LAAAIXX4NKC6XS5KUmJjotT4xMdGzzeVyKSEhwWt7hw4dFBcX52nzdUVFRYqJifEsKSkp/iwbAAAYJiiu4iksLFRtba1nOXHiRKBLAgAAbcivASUpKUmSVF1d7bW+urrasy0pKUk1NTVe2y9fvqzTp0972nydzWaT3W73WgAAQOjya0Dp06ePkpKSVFpa6lnndrv10UcfyeFwSJIcDofOnDmjiooKT5v3339fjY2NysjI8Gc5AAAgSPl8Fc+5c+d0+PBhz/OjR49qz549iouLU8+ePTV9+nQ9/fTTuu2229SnTx/9/Oc/V3JyssaMGSNJuv322zVy5EhNnDhRy5cv16VLlzRlyhQ9+OCDXMEDAAAktSCgfPzxx/r+97/veV5QUCBJys/P18qVK/X444/r/PnzmjRpks6cOaNhw4Zp06ZN6tSpk+c1q1ev1pQpUzRixAiFh4crLy9PL730kh+6AwAAQkGYZVlWoIvwldvtVkxMjGpra5mPAkjqPXtDoEsIaccW5ga6BCAk+PL5HRRX8QAAgPaFgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMbpEOgCADRP79kbPI+PLcwNYCUA0PYYQQGAAOo9e4NX+ATwDwQUAABgHAIKALQAIx9A2yKgAAAA4xBQAACAcQgoQAA1dZqAUwfm4WcCXH9+DyhPPvmkwsLCvJbU1FTP9osXL8rpdCo+Pl5du3ZVXl6eqqur/V0GAAAIYm0ygvLd735XJ0+e9CwffvihZ9uMGTO0bt06rV27VmVlZaqqqtLYsWPbogwAABCk2uRGbR06dFBSUtI31tfW1urVV19VcXGx7rnnHknSihUrdPvtt2vHjh0aOnRoW5QDGI/TB+3DlZ8zN9oDvl2bjKB89tlnSk5O1s0336xx48bp+PHjkqSKigpdunRJmZmZnrapqanq2bOnysvLr7q/uro6ud1urwUATHRlvgqhE2gdv4+gZGRkaOXKlerbt69Onjyp+fPn65//+Z+1f/9+uVwuRUZGKjY21us1iYmJcrlcV91nUVGR5s+f7+9SgYDhwwsArs3vASUnJ8fzeMCAAcrIyFCvXr30hz/8QZ07d27RPgsLC1VQUOB57na7lZKS0upaAaAtcUoHaLk2/7LA2NhYfec739Hhw4d17733qr6+XmfOnPEaRamurm5yzsoVNptNNputrUsFjMIoS3Dh5wX4V5vfB+XcuXM6cuSIevToofT0dHXs2FGlpaWe7ZWVlTp+/LgcDkdblwJcV8xDaB/4OQNtw+8jKD/96U81evRo9erVS1VVVZo3b54iIiL0ox/9SDExMZowYYIKCgoUFxcnu92uqVOnyuFwcAUP4ANOHQAIdX4PKJ9//rl+9KMf6csvv1T37t01bNgw7dixQ927d5ckLV68WOHh4crLy1NdXZ2ys7P18ssv+7sMwBhf/euaQAEAzeP3gLJmzZprbu/UqZOWLl2qpUuX+vvQABAUmjolxKgY4K3NJ8kCaDvMfQgOvvycCCrAPxBQAD8jNABA6xFQgOuI8AIAzdPmlxkDAAD4ioACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4XGYMAM3EZeLA9cMICgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAF/Te/YGJkMCQIARUAAAgHEIKAAAwDgEFAAAYBxu1Ab44MrclGMLcwNcCeB947grv5P8jiJUEFAAwHDNCR1M7EaoIaAAV9HUX6cAgOuDgAI0A3+dAsD1RUABWoDz/DAdI4AIdgQUtHutCRt8CKCtMGqH9o6Agnbh62/2hAkEs5aGF0b+EEy4DwoAADAOIygIevxVCLQMpyhhMgIKQtr1PI/PnAEEM4I+TENAQbtEmEAw8tfvLb//CAYEFASV5gxJt3YCIRBq+N1GMCKgAAD85mqnipjvAl9xFQ+uu96zNzT7Lzpf2gIAQgcjKAgKzQkpBBmg9ZgsC1MQUGAMXwMGgQRoO74Eldb+3yUMoSkBDShLly7VL37xC7lcLg0cOFBLlizRkCFDAlkSrqE155YJE0Bw4v9uyzEa1ToBm4PyxhtvqKCgQPPmzdPu3bs1cOBAZWdnq6amJlAlQf6bH9Kc/TC/BAhd/P9GawVsBOWFF17QxIkT9eijj0qSli9frg0bNui///u/NXv27ECV1SRTU3BLZ8U39aZxrVGR1mDuCNC++ft7g671vufLsVr6fs7pqesnIAGlvr5eFRUVKiws9KwLDw9XZmamysvLv9G+rq5OdXV1nue1tbWSJLfb3Sb19Zu3ucn11zreldfsn599zXXX2v9XNbWfr29rrLvQrNq+7quvu6LnjLXNfv2VY7V2PwAg/f/3lKbeG6/1nnKt1zX3mF99/ZX3Vl/2d639XHmPbOr9+WqfD9eq51ptv+pqnznX+nxqzn784cq/hWVZ397YCoAvvvjCkmRt377da/3MmTOtIUOGfKP9vHnzLEksLCwsLCwsIbCcOHHiW7NCUFzFU1hYqIKCAs/zxsZGnT59WvHx8QoLC/Nq63a7lZKSohMnTshut1/vUq+b9tJPqf30tb30U2o/faWfoae99LWt+mlZls6ePavk5ORvbRuQgNKtWzdFRESourraa311dbWSkpK+0d5ms8lms3mti42NveYx7HZ7SP/yXNFe+im1n762l35K7aev9DP0tJe+tkU/Y2JimtUuIFfxREZGKj09XaWlpZ51jY2NKi0tlcPhCERJAADAIAE7xVNQUKD8/HwNHjxYQ4YM0Ysvvqjz5897ruoBAADtV8ACyg9/+EOdOnVKc+fOlcvl0h133KFNmzYpMTGxVfu12WyaN2/eN04JhZr20k+p/fS1vfRTaj99pZ+hp7301YR+hllWc671AQAAuH74NmMAAGAcAgoAADAOAQUAABiHgAIAAIwT0gHl008/1f33369u3brJbrdr2LBh+tOf/hTostrMhg0blJGRoc6dO+uGG27QmDFjAl1Sm6mrq9Mdd9yhsLAw7dmzJ9Dl+NWxY8c0YcIE9enTR507d9Ytt9yiefPmqb6+PtCl+cXSpUvVu3dvderUSRkZGdq5c2egS/KroqIi3XnnnYqOjlZCQoLGjBmjysrKQJd1XSxcuFBhYWGaPn16oEvxuy+++EI//vGPFR8fr86dO6t///76+OOPA12WXzU0NOjnP/+513vPU0891bzvzWkDIR1Q7rvvPl2+fFnvv/++KioqNHDgQN13331yuVyBLs3v/ud//kfjx4/Xo48+qr/+9a/685//rIceeijQZbWZxx9/vFm3Sg5Ghw4dUmNjo1555RUdOHBAixcv1vLly/XEE08EurRWe+ONN1RQUKB58+Zp9+7dGjhwoLKzs1VTUxPo0vymrKxMTqdTO3bsUElJiS5duqSsrCydP38+0KW1qV27dumVV17RgAEDAl2K3/3973/XXXfdpY4dO2rjxo06ePCgfvnLX+qGG24IdGl+9dxzz2nZsmX6r//6L33yySd67rnntGjRIi1ZsiQwBfnl2/8MdOrUKUuStW3bNs86t9ttSbJKSkoCWJn/Xbp0ybrxxhut3/72t4Eu5bp49913rdTUVOvAgQOWJOsvf/lLoEtqc4sWLbL69OkT6DJabciQIZbT6fQ8b2hosJKTk62ioqIAVtW2ampqLElWWVlZoEtpM2fPnrVuu+02q6SkxPqXf/kXa9q0aYEuya9mzZplDRs2LNBltLnc3FzrP/7jP7zWjR071ho3blxA6gnZEZT4+Hj17dtXr732ms6fP6/Lly/rlVdeUUJCgtLT0wNdnl/t3r1bX3zxhcLDwzVo0CD16NFDOTk52r9/f6BL87vq6mpNnDhRv/vd7xQVFRXocq6b2tpaxcXFBbqMVqmvr1dFRYUyMzM968LDw5WZmany8vIAVta2amtrJSnof37X4nQ6lZub6/WzDSX/+7//q8GDB+sHP/iBEhISNGjQIP3mN78JdFl+90//9E8qLS3Vp59+Kkn661//qg8//FA5OTkBqScovs24JcLCwrRlyxaNGTNG0dHRCg8PV0JCgjZt2hRyw3J/+9vfJElPPvmkXnjhBfXu3Vu//OUvNXz4cH366ach88ZoWZYeeeQRPfbYYxo8eLCOHTsW6JKui8OHD2vJkiV6/vnnA11Kq/zf//2fGhoavnG36MTERB06dChAVbWtxsZGTZ8+XXfddZf69esX6HLaxJo1a7R7927t2rUr0KW0mb/97W9atmyZCgoK9MQTT2jXrl36z//8T0VGRio/Pz/Q5fnN7Nmz5Xa7lZqaqoiICDU0NOiZZ57RuHHjAlJP0I2gzJ49W2FhYddcDh06JMuy5HQ6lZCQoA8++EA7d+7UmDFjNHr0aJ08eTLQ3WiW5va1sbFRkvSzn/1MeXl5Sk9P14oVKxQWFqa1a9cGuBffrrn9XLJkic6ePavCwsJAl9wize3nV33xxRcaOXKkfvCDH2jixIkBqhwt5XQ6tX//fq1ZsybQpbSJEydOaNq0aVq9erU6deoU6HLaTGNjo773ve/p2Wef1aBBgzRp0iRNnDhRy5cvD3RpfvWHP/xBq1evVnFxsXbv3q1Vq1bp+eef16pVqwJST9Dd6v7UqVP68ssvr9nm5ptv1gcffKCsrCz9/e9/9/qq6Ntuu00TJkzQ7Nmz27rUVmtuX//85z/rnnvu0QcffKBhw4Z5tmVkZCgzM1PPPPNMW5faKs3t57//+79r3bp1CgsL86xvaGhQRESExo0bF7D/RM3V3H5GRkZKkqqqqjR8+HANHTpUK1euVHh40P094aW+vl5RUVH64x//6HWFWX5+vs6cOaN33nkncMW1gSlTpuidd97Rtm3b1KdPn0CX0ybefvtt/du//ZsiIiI86xoaGhQWFqbw8HDV1dV5bQtWvXr10r333qvf/va3nnXLli3T008/rS+++CKAlflXSkqKZs+eLafT6Vn39NNP6/XXXw/IKGfQneLp3r27unfv/q3tLly4IEnfeFMPDw/3jDiYrrl9TU9Pl81mU2VlpSegXLp0SceOHVOvXr3ausxWa24/X3rpJT399NOe51VVVcrOztYbb7yhjIyMtizRL5rbT+kfIyff//73PaNhwR5OJCkyMlLp6ekqLS31BJTGxkaVlpZqypQpgS3OjyzL0tSpU/XWW29p69atIRtOJGnEiBHat2+f17pHH31UqampmjVrVkiEE0m66667vnGp+KeffhoU76++uHDhwjfeayIiIgL3mRmQqbnXwalTp6z4+Hhr7Nix1p49e6zKykrrpz/9qdWxY0drz549gS7P76ZNm2bdeOON1ubNm61Dhw5ZEyZMsBISEqzTp08HurQ2c/To0ZC8iufzzz+3br31VmvEiBHW559/bp08edKzBLs1a9ZYNpvNWrlypXXw4EFr0qRJVmxsrOVyuQJdmt9MnjzZiomJsbZu3er1s7tw4UKgS7suQvEqnp07d1odOnSwnnnmGeuzzz6zVq9ebUVFRVmvv/56oEvzq/z8fOvGG2+01q9fbx09etR68803rW7dulmPP/54QOoJ2YBiWZa1a9cuKysry4qLi7Oio6OtoUOHWu+++26gy2oT9fX11k9+8hMrISHBio6OtjIzM639+/cHuqw2FaoBZcWKFZakJpdQsGTJEqtnz55WZGSkNWTIEGvHjh2BLsmvrvazW7FiRaBLuy5CMaBYlmWtW7fO6tevn2Wz2azU1FTr17/+daBL8ju3221NmzbN6tmzp9WpUyfr5ptvtn72s59ZdXV1Aakn6OagAACA0Bf8J7YBAEDIIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDj/DxEotPOUh32VAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import opendp.prelude as dp\n",
    "dp.enable_features(\"contrib\")\n",
    "\n",
    "space = dp.atom_domain(T=float, nan=False), dp.absolute_distance(T=float)\n",
    "# calibrate noise to satisfy (1.0, 1e-7)-approxDP when sensitivity is 1.0\n",
    "m_cnd = space >> dp.m.then_canonical_noise(d_in=1.0, d_out=(1.0, 1e-7))\n",
    "\n",
    "# plot 10k samples\n",
    "plt.hist([m_cnd(0.) for _ in range(10_000)], bins=200);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also have a suite of inference tools that can be applied under the assumption that your statistic is binomially distributed.\n",
    "An example of binomially-distributed data can be found when counting the number of married individuals in the California demographics dataset,\n",
    "a microdata-level dataset with attributes for age, sex, education, race, income and marriage status."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import opendp.prelude as dp\n",
    "import polars as pl\n",
    "\n",
    "dp.enable_features(\"contrib\")\n",
    "columns = \"age\", \"sex\", \"educ\", \"race\", \"income\", \"married\"\n",
    "married_data = list(pl.read_csv(\n",
    "    dp.examples.get_california_pums_path(), \n",
    "    has_header=False, \n",
    "    new_columns=columns,\n",
    "    ignore_errors=True,\n",
    ")[\"married\"])\n",
    "\n",
    "context = dp.Context.compositor(\n",
    "    data=married_data,\n",
    "    privacy_unit=dp.unit_of(contributions=1),\n",
    "    privacy_loss=dp.loss_of(epsilon=0.2, delta=1e-8),\n",
    "    split_evenly_over=2\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the canonical noise mechanism to conduct statistical inference on the 'married' column,\n",
    "where values are either 0 or 1, indicating marriage status.\n",
    "\n",
    "The following query counts the number of married individuals in the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = context.query().cast_default(float).impute_constant(0.0).clamp((0.0, 1.0)).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The canonical noise mechanism behaves similarly to the Laplace or Gaussian mechanism in that it is an additive noise mechanism.\n",
    "\n",
    "The implementation of the canonical noise mechanism in OpenDP only supports scalar-valued float inputs and can only be calibrated to satisfy $(\\epsilon, \\delta)$-DP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "543.5862745175556"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "estimate: float = query.canonical_noise().release()\n",
    "estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binomial Assumption"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you expect that the statistic you are privatizing follows a binomial distribution,\n",
    "which is commonly the case in counting queries,\n",
    "there is a suite of tools for conducting hypothesis tests and constructing confidence intervals.\n",
    "\n",
    "These utilities can be used standalone via the [`BinomialCND`](../../python/opendp.extras.numpy.canonical.html#opendp.extras.numpy.canonical.BinomialCND) class,\n",
    "or you can simply provide a `binomial_size` prior before making the release and OpenDP will instantiate the class for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BinomialCND(estimate=502.87558450908193, d_in=1.0000000093132269, d_out=(0.1, 5e-09), size=1000)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from opendp.extras.numpy.canonical import BinomialCND\n",
    "\n",
    "estimate: BinomialCND = query.canonical_noise(binomial_size=1_000).release()\n",
    "estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `binomial_size` argument is a prior about the number of bernoulli draws- typically the number of records in the dataset.\n",
    "\n",
    "Let's get familiar with each of the inference tools."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confidence Intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.4608320830961166, 0.5449190859220473)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ci = estimate.confidence_interval(alpha=.05)\n",
    "ci"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a statistical significance of .05 (with 95% confidence), the true proportion of married individuals is within the interval ``ci``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.46834281999150473"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lower_ci = estimate.confidence_interval(alpha=.05, side='lower')\n",
    "lower_ci"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a statistical significance of .05 (with 95% confidence), the true proportion of married individuals is no less than ``lower_ci``."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### P-Values\n",
    "\n",
    "Consider a statistical test where the null hypothesis states that the probability of a success (incrementing the counter) is `theta`.\n",
    "The p-value tells you the probability of observing `estimate` under the assumption that the null hypothesis is true.\n",
    "A very small p-value (less than some statistical significance level) indicates that the null hypothesis should be rejected.\n",
    "\n",
    "For examples in this document, we'll start with a prior assumption that the probability of a person being married is 0.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, test the null hypothesis that the probability of a person being married is 0.5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8861710889637788"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_value_married = estimate.p_value(theta=theta)\n",
    "p_value_married"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Under this null hypothesis, the probability of observing `estimate` is low:\n",
    "this tells us that if the true marriage rate was actually 0.5, then it would be unlikely we would have observed such a high private estimate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also test the null hypothesis that the probability of a person being married is no less than 0.5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5569144555181105"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_value_married = estimate.p_value(theta=theta, tail='left')\n",
    "p_value_married"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Under this hypothesis, it's highly likely to observe an outcome like `estimate`, because the private estimate corroborates the null hypothesis.\n",
    "\n",
    "When the null hypothesis is that the probability of a person being married is no greater than 0.5, the test more strongly indicates that the outcome is unlikely:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4430855444818895"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_value_married = estimate.p_value(theta=theta, tail='right')\n",
    "p_value_married"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the one-sided test is stronger than the two-sided test, which is reflected in the p-value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uniformly-Most-Powerful Tests\n",
    "`one_sided_uniformly_most_powerful_tests` shows how p-values vary over all potential estimates between 0 and the dataset size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAALMpJREFUeJzt3Xt8VPWd//H3TJKZJEASIDDhEgxeKloQFEoM6HZdU7PK0tpt+6DICkutfWixC2ZbJV5grdWw7UrptlQqLXX3t1pQV21XKC4bvJQaRW4qXlAEDAUSiJBMSEImmfn+/khmIBICA0nOfE9ez8djHpIz58x85uB38uZ7OcdjjDECAABwiNfpAgAAQO9GGAEAAI4ijAAAAEcRRgAAgKMIIwAAwFGEEQAA4CjCCAAAcBRhBAAAOCrZ6QLORCQS0f79+9WvXz95PB6nywEAAGfAGKO6ujoNHTpUXu+p+z+sCCP79+9Xbm6u02UAAICzsHfvXg0fPvyUz1sRRvr16yep9cNkZGQ4XA0AADgTwWBQubm5sd/jp2JFGIkOzWRkZBBGAACwzOmmWDCBFQAAOIowAgAAHEUYAQAAjiKMAAAARxFGAACAowgjAADAUYQRAADgKMIIAABwFGEEAAA4Ku4w8uqrr2rq1KkaOnSoPB6Pnn/++dMe8/LLL+uKK66Q3+/XhRdeqMcff/wsSgUAAG4Udxipr6/X2LFjtXTp0jPaf/fu3ZoyZYquueYabdu2TfPmzdO3v/1tvfjii3EXCwAA3Cfue9Ncf/31uv766894/2XLlmnkyJF65JFHJEmXXHKJNmzYoJ/+9KcqKiqK9+0BAIDLdPuN8srLy1VYWNhuW1FRkebNm3fKY5qamtTU1BT7ORgMdld5AHBG6o41qyrYpNrGkGobmxVsbNGx5rCawxGFwkahloiaw62PiDEyRjJS239N64vEtp3iecBB35o8UrkD0h15724PI5WVlQoEAu22BQIBBYNBNTY2Ki0t7aRjSktL9cADD3R3aQDQoYN1x/TnndXatOeItu8PquLTeh1paHa6LKBbTR071L1h5GyUlJSouLg49nMwGFRubq6DFQFwO2OM1n9wUI+/tkd/3lmtSAedFf1Sk5WVnqLMtNZHWkqSUpK8SknyypfsbfuzR96226V7PJJHnrb/6vg2j6f15888f5q7rAPdKpCR6th7d3sYycnJUVVVVbttVVVVysjI6LBXRJL8fr/8fn93lwYAkqSPDx3VXc+8rc2fHIltGzMsUxNHDtAVI/rr/EF9lDsgXX39CfnvN8B63d6yCgoKtGbNmnbb1q1bp4KCgu5+awA4rT++c0B3PrVNx5ojSktJ0sxJ5+mmiSN03sA+TpcG9Bpxh5GjR49q586dsZ93796tbdu2acCAARoxYoRKSkq0b98+/ed//qck6bbbbtMvfvEL3XXXXfrWt76l9evX66mnntLq1au77lMAwFl4dstf9P2n31LESFddmK2ffOMyDcnsuMcWQPeJO4xs2rRJ11xzTezn6NyOWbNm6fHHH9eBAwdUUVERe37kyJFavXq17rzzTv3sZz/T8OHD9etf/5plvQActWnPYd39328rYqRvfiFXD311jJK8TNoAnOAxxiT8mrJgMKjMzEzV1tYqIyPD6XIAWK62oVlf+ukrOljXpBvG5GjpTVfIw+xRoMud6e9v7k0DoNdZvG6HDtY16fxBffSTr48liAAOI4wA6FXePxDU/3v9E0nSj74yWn1YIQM4jjACoFf5+fqPFDHSDWNyNOnCbKfLASDCCIBeZE91vf64vVKSNPfazzlcDYAowgiAXuN3GytkjPTXFw/SxTn9nC4HQBvCCIBeoSUc0bNb90mSpk8c4XA1AE5EGAHQK/z54091qK5JA/v49DejBjtdDoATEEYA9Ar/+27rXJGi0TlKSeKrD0gktEgArheJGK17r/WGnV+6NOBwNQA+izACwPXeOxDUwbom9fEladIFA50uB8BnEEYAuN5rH1dLkvLPHyh/cpLD1QD4LMIIANcr//hTSVLB+fSKAImIMALA1VrCEb2554gkqYAhGiAhEUYAuNqHVUd1tKlFff3JumQId/0GEhFhBICrvf2XGknSmGGZSvJyd14gERFGALjaW3+plSSNzc1ythAAp0QYAeBq0Z6RscMznS0EwCkRRgC4Vqgloh2VdZKkMYQRIGERRgC41u7qerVEjPr5kzUsK83pcgCcAmEEgGvtqGrtFbko0FceD5NXgURFGAHgWh+1hZGLc/o5XAmAzhBGALhWdL7IRYMJI0AiI4wAcK2dB49Kkj4XIIwAiYwwAsCVWsIRVRxukCRdMLiPw9UA6AxhBIAr7a85ppaIkS/Zq0C/VKfLAdAJwggAV/rkcL0k6bwB6fJyGXggoRFGALjSnk9bh2jOG8gQDZDoCCMAXOmT6raekYHpDlcC4HQIIwBcKdozkkcYARIeYQSAK/3lSGsYGT6AMAIkOsIIAFfaX9MoSdyTBrAAYQSA69Q3tSh4rEWSNCSTZb1AoiOMAHCdA7WtvSL9UpPVLzXF4WoAnA5hBIDr7Ks5JkkamskQDWADwggA1znQNl9kSBZDNIANCCMAXGd/bWvPyBB6RgArEEYAuE60Z2Qok1cBKxBGALhOZbC1ZyRAGAGsQBgB4DrVR0OSpMH9/A5XAuBMEEYAuM6huiZJ0iDCCGAFwggAVwlHjA7XE0YAmxBGALjKp/VNihjJ65EG9iGMADYgjABwlegQzYA+fiV5PQ5XA+BMEEYAuEo0jGT39TlcCYAzRRgB4CpMXgXsQxgB4CqHjhJGANsQRgC4SnVd6zVGBvUljAC2IIwAcJXost6BzBkBrEEYAeAqRxqaJUn90wkjgC0IIwBc5UhD6zANYQSwB2EEgKscrm8LI30II4AtCCMAXKUmNkyT4nAlAM4UYQSAa4RaIjra1CJJGkDPCGANwggA16hpmy/i9UgZqfSMALYgjABwjcNtYSQr3Scv96UBrEEYAeAaR+qZLwLYiDACwDVY1gvY6azCyNKlS5WXl6fU1FTl5+dr48aNne6/ZMkSXXzxxUpLS1Nubq7uvPNOHTt27KwKBoBTiYURJq8CVok7jKxatUrFxcVauHChtmzZorFjx6qoqEgHDx7scP8nn3xS8+fP18KFC/X+++/rN7/5jVatWqV77rnnnIsHgBMdiV5jhGEawCpxh5HFixfr1ltv1ezZs3XppZdq2bJlSk9P14oVKzrc/7XXXtPkyZN10003KS8vT9ddd52mT59+2t4UAIhX8Fjrst7MNMIIYJO4wkgoFNLmzZtVWFh4/AW8XhUWFqq8vLzDYyZNmqTNmzfHwseuXbu0Zs0a3XDDDad8n6amJgWDwXYPADidYGPrBFaW9QJ2SY5n5+rqaoXDYQUCgXbbA4GAPvjggw6Puemmm1RdXa2rrrpKxhi1tLTotttu63SYprS0VA888EA8pQGAgsfawgg9I4BVun01zcsvv6yHH35Yv/zlL7VlyxY9++yzWr16tR588MFTHlNSUqLa2trYY+/evd1dJgAXCDYyTAPYKK6ekezsbCUlJamqqqrd9qqqKuXk5HR4zP3336+bb75Z3/72tyVJY8aMUX19vb7zne/o3nvvldd7ch7y+/3y+/3xlAYAJ/SMxPXVBsBhcfWM+Hw+jR8/XmVlZbFtkUhEZWVlKigo6PCYhoaGkwJHUlKSJMkYE2+9AHBKzBkB7BT3Px+Ki4s1a9YsTZgwQRMnTtSSJUtUX1+v2bNnS5JmzpypYcOGqbS0VJI0depULV68WJdffrny8/O1c+dO3X///Zo6dWoslABAV4iupmHOCGCXuMPItGnTdOjQIS1YsECVlZUaN26c1q5dG5vUWlFR0a4n5L777pPH49F9992nffv2adCgQZo6daoeeuihrvsUAHo9Yww9I4ClPMaCsZJgMKjMzEzV1tYqIyPD6XIAJKCGUIsuXfCiJOm9HxYp3ce8EcBpZ/r7m3vTAHCF6EqaZK9HaSkMAQM2IYwAcIUTrzHi8XgcrgZAPAgjAFzh+HwRhmcA2xBGALgCV18F7EUYAeAK0TkjrKQB7EMYAeAKdW09I339DNMAtiGMAHCFo01hSVIfwghgHcIIAFeob2odpunrZ1kvYBvCCABXONoWRugZAexDGAHgCrGeEZb2AtYhjABwhfpQdJiGMALYhjACwBViE1i5Jw1gHcIIAFc42ra0lzkjgH0IIwBcob6tZ4RhGsA+hBEArnB8NQ1LewHbEEYAuAITWAF7EUYAuEI91xkBrEUYAWC9ppawmsNGEmEEsBFhBID1opNXJamPjzkjgG0IIwCsFx2iSU3xKjmJrzXANrRaANY72sTkVcBmhBEA1uMmeYDdCCMArBdbScOl4AErEUYAWK8x1HZfGi54BliJMALAeg1tYSQ1hTAC2IgwAsB6Dc2tYSSdZb2AlQgjAKzX2HYp+HTmjABWIowAsF5jKCJJSqNnBLASYQSA9RqaW3tG0pgzAliJMALAetHVNMwZAexEGAFgvehqGoZpADsRRgBYrzG6moZhGsBKhBEA1mukZwSwGmEEgPUa2pb2prG0F7ASYQSA9WITWBmmAaxEGAFgvUauwApYjTACwHqxe9MQRgArEUYAWI/rjAB2I4wAsF5DbM4IE1gBGxFGAFjNGBObM8LSXsBOhBEAVjvWHIn9mTAC2IkwAsBq0WuMSNwoD7AVYQSA1aJDNP5kr5K8HoerAXA2CCMArMZKGsB+hBEAVovOGUlliAawFmEEgNWaWo4P0wCwE60XgNWaWlp7RvzJ9IwAtiKMALBarGckha8zwFa0XgBWa2qO9ozwdQbYitYLwGoM0wD2I4wAsFqohZ4RwHa0XgBWY84IYD9aLwCrRYdpfEl8nQG2ovUCsBpzRgD7EUYAWK2pmWEawHa0XgBWa2ICK2C9s2q9S5cuVV5enlJTU5Wfn6+NGzd2un9NTY3mzJmjIUOGyO/363Of+5zWrFlzVgUDwIkYpgHslxzvAatWrVJxcbGWLVum/Px8LVmyREVFRdqxY4cGDx580v6hUEhf+tKXNHjwYD3zzDMaNmyYPvnkE2VlZXVF/QB6Oe5NA9gv7jCyePFi3XrrrZo9e7YkadmyZVq9erVWrFih+fPnn7T/ihUrdPjwYb322mtKSUmRJOXl5Z1b1QDQJtYzwpwRwFpxtd5QKKTNmzersLDw+At4vSosLFR5eXmHx/zhD39QQUGB5syZo0AgoNGjR+vhhx9WOBw+5fs0NTUpGAy2ewBARximAewXVxiprq5WOBxWIBBotz0QCKiysrLDY3bt2qVnnnlG4XBYa9as0f33369HHnlEP/rRj075PqWlpcrMzIw9cnNz4ykTQC8SvTeNj2EawFrd3nojkYgGDx6sxx57TOPHj9e0adN07733atmyZac8pqSkRLW1tbHH3r17u7tMAJZizghgv7jmjGRnZyspKUlVVVXttldVVSknJ6fDY4YMGaKUlBQlJR3vQr3kkktUWVmpUCgkn8930jF+v19+vz+e0gD0UgzTAPaL658SPp9P48ePV1lZWWxbJBJRWVmZCgoKOjxm8uTJ2rlzpyKRSGzbhx9+qCFDhnQYRAAgHlxnBLBf3K23uLhYy5cv13/8x3/o/fff1+233676+vrY6pqZM2eqpKQktv/tt9+uw4cPa+7cufrwww+1evVqPfzww5ozZ07XfQoAvRZXYAXsF/fS3mnTpunQoUNasGCBKisrNW7cOK1duzY2qbWiokJe7/EvhdzcXL344ou68847ddlll2nYsGGaO3eu7r777q77FAB6rVCYYRrAdh5jjHG6iNMJBoPKzMxUbW2tMjIynC4HQAKZvGi99tU06vdzJmtsbpbT5QA4wZn+/qZfE4DVonNGWNoL2IvWC8BqLO0F7EfrBWC145eDZ84IYCvCCABrGWMUYmkvYD1aLwBrRXtFJMIIYDNaLwBrRZf1SiztBWxGGAFgrehN8jweKSXJ43A1AM4WYQSAtaIraXxJXnk8hBHAVoQRANbivjSAO9CCAVgrOkzDsl7AboQRANbigmeAO9CCAViLYRrAHWjBAKx1PIwwTAPYjDACwFqxq6+m8FUG2IwWDMBazBkB3IEWDMBa0dU0PoZpAKsRRgBYiwmsgDvQggFYi2EawB1owQCsxWoawB0IIwCsdfwKrHyVATajBQOwVijMMA3gBrRgANaK9YwwTANYjTACwFrROSM+ekYAq9GCAViL1TSAO9CCAViL64wA7kALBmCt46tpmDMC2IwwAsBaDNMA7kALBmCtUJhhGsANaMEArMXSXsAdCCMArMUEVsAdaMEArMWcEcAdaMEArBXrGeHeNIDVaMEArMWcEcAdCCMArMUwDeAOtGAA1gq10DMCuAFhBIC1mDMCuAMtGICVWsIRtUSMJIZpANvRggFYKXr1VUnyEUYAq9GCAVgpupJGknxJfJUBNqMFA7BSdL5IstejZMIIYDVaMAArsawXcA9aMQArHV9Jw7JewHaEEQBWCnGTPMA1aMUArMQwDeAetGIAVoqupmFZL2A/WjEAKzVxKXjANQgjAKzEMA3gHrRiAFbivjSAe9CKAVgpOmeEYRrAfoQRAFZqCrO0F3ALWjEAKzU1t84ZYTUNYD9aMQArHWsLI6kM0wDWI4wAsFJ0AmsqE1gB69GKAVgp1jPCvWkA6xFGAFjpWDM3ygPcgjACwErHe0b4GgNsd1ateOnSpcrLy1Nqaqry8/O1cePGMzpu5cqV8ng8uvHGG8/mbQEg5hiXgwdcI+4wsmrVKhUXF2vhwoXasmWLxo4dq6KiIh08eLDT4/bs2aPvf//7uvrqq8+6WACIomcEcI+4W/HixYt16623avbs2br00ku1bNkypaena8WKFac8JhwOa8aMGXrggQd0/vnnn1PBACCxtBdwk7jCSCgU0ubNm1VYWHj8BbxeFRYWqry8/JTH/fCHP9TgwYN1yy23nNH7NDU1KRgMtnsAwImOL+0ljAC2iyuMVFdXKxwOKxAItNseCARUWVnZ4TEbNmzQb37zGy1fvvyM36e0tFSZmZmxR25ubjxlAugFmhimAVyjW1txXV2dbr75Zi1fvlzZ2dlnfFxJSYlqa2tjj71793ZjlQBsFF3aS88IYL/keHbOzs5WUlKSqqqq2m2vqqpSTk7OSft//PHH2rNnj6ZOnRrbFom0foEkJydrx44duuCCC046zu/3y+/3x1MagF7mWEtrzwg3ygPsF1cr9vl8Gj9+vMrKymLbIpGIysrKVFBQcNL+o0aN0jvvvKNt27bFHl/+8pd1zTXXaNu2bQy/ADhrXIEVcI+4ekYkqbi4WLNmzdKECRM0ceJELVmyRPX19Zo9e7YkaebMmRo2bJhKS0uVmpqq0aNHtzs+KytLkk7aDgDxOD5MQ88IYLu4w8i0adN06NAhLViwQJWVlRo3bpzWrl0bm9RaUVEhr5cvBwDdqyk2TEPPCGA7jzHGOF3E6QSDQWVmZqq2tlYZGRlOlwPAYcYYjSxZI0l6895CDerHHDMgEZ3p72+6MABYJ3qNEYlhGsANaMUArNPUfGIYYZgGsB1hBIB1ost6vR4p2etxuBoA54owAsA6Jy7r9XgII4DtCCMArMPVVwF3IYwAsE50WW8qV18FXIGWDMA69IwA7kIYAWCd6JwRP2EEcAXCCADrxMIIwzSAK9CSAVjnWAv3pQHchJYMwDqNoRZJUrov7ttrAUhAhBEA1mkMtQ7TpPmYMwK4AWEEgHUa2uaMpDOBFXAFwggA69AzArgLYQSAdRoII4CrEEYAWCcaRtJTmMAKuAFhBIB1otcZSadnBHAFwggA6zS0Le1lmAZwB8IIAOvE5oywmgZwBcIIAOtEV9MwTAO4A2EEgHVYTQO4C2EEgHUaYxNYWU0DuAFhBIB1GKYB3IUwAsA60dU0qUxgBVyBMALAOo1cZwRwFcIIAKs0hyNqDhtJhBHALQgjAKwSXUkjsZoGcAvCCACrRC8Fn+T1yJfEVxjgBrRkAFY58eqrHo/H4WoAdAXCCACrcF8awH0IIwCswjVGAPchjACwSn2Iq68CbkMYAWCV+qbWYZq+fnpGALcgjACwytFjrWGkj5+eEcAtCCMArHI01jNCGAHcgjACwCr1hBHAdQgjAKxyNMQwDeA2hBEAVon2jBBGAPcgjACwSnQCK6tpAPcgjACwytGm1uuM9PWnOFwJgK5CGAFglePDNPSMAG5BGAFglfoQq2kAtyGMALDKUSawAq5DGAFgleMTWAkjgFsQRgBYhaW9gPsQRgBYIxIxsbv20jMCuAdhBIA1GprDsT8TRgD3IIwAsEZ0iMbrkVJT+PoC3ILWDMAadceOzxfxeDwOVwOgqxBGAFgjeKxZkpSZxtVXATchjACwRrCxNYxkpBJGADchjACwRrBtmCYjjcmrgJsQRgBYg54RwJ0IIwCsEZ0zksGcEcBVCCMArBFsbBumoWcEcBXCCABrHO8ZYc4I4CZnFUaWLl2qvLw8paamKj8/Xxs3bjzlvsuXL9fVV1+t/v37q3///iosLOx0fwA4FeaMAO4UdxhZtWqViouLtXDhQm3ZskVjx45VUVGRDh482OH+L7/8sqZPn66XXnpJ5eXlys3N1XXXXad9+/adc/EAepfoahquMwK4S9xhZPHixbr11ls1e/ZsXXrppVq2bJnS09O1YsWKDvd/4okn9N3vflfjxo3TqFGj9Otf/1qRSERlZWXnXDyA3iXWM0IYAVwlrjASCoW0efNmFRYWHn8Br1eFhYUqLy8/o9doaGhQc3OzBgwYcMp9mpqaFAwG2z0AIDZnJJU5I4CbxBVGqqurFQ6HFQgE2m0PBAKqrKw8o9e4++67NXTo0HaB5rNKS0uVmZkZe+Tm5sZTJgCXiq2moWcEcJUeXU2zaNEirVy5Us8995xSU1NPuV9JSYlqa2tjj7179/ZglQASkTGGYRrApeLq68zOzlZSUpKqqqraba+qqlJOTk6nx/7bv/2bFi1apP/7v//TZZdd1um+fr9ffr8/ntIAuFxTS0ShcESS1I9hGsBV4uoZ8fl8Gj9+fLvJp9HJqAUFBac87sc//rEefPBBrV27VhMmTDj7agH0WkcaQpKkZK9H/fyEEcBN4m7RxcXFmjVrliZMmKCJEydqyZIlqq+v1+zZsyVJM2fO1LBhw1RaWipJ+td//VctWLBATz75pPLy8mJzS/r27au+fft24UcB4GZH6luHaLLSffJ4PA5XA6ArxR1Gpk2bpkOHDmnBggWqrKzUuHHjtHbt2tik1oqKCnm9xztcHn30UYVCIX39619v9zoLFy7Uv/zLv5xb9QB6jWjPyIA+zBcB3Oas+jrvuOMO3XHHHR0+9/LLL7f7ec+ePWfzFgDQzuH61jCSle5zuBIAXY170wCwQqxnhDACuA5hBIAVonNG+vchjABuQxgBYIVoz0j/dOaMAG5DGAFgheickQH0jACuQxgBYIXjPSOEEcBtCCMArBALIyztBVyHMALACrEJrPSMAK5DGAGQ8Iwx+rS+SRJzRgA3IowASHj1obCONbfeJC+7LzfRBNyGMAIg4R2qa+0V6eNLUh9ukge4DmEEQMKLhpFB/egVAdyIMAIg4RFGAHcjjABIeIfqjkkijABuRRgBkPAOHW3rGWHyKuBKhBEACY9hGsDdCCMAEh5hBHA3wgiAhBcbpiGMAK5EGAGQ8CprW8PI4H6pDlcCoDsQRgAktKaWsKrbekaGZqU5XA2A7kAYAZDQKmtbl/X6k73qn84dewE3IowASGj7a1rDyLCsNHk8HoerAdAdCCMAEtqB2kZJ0pAs5osAbkUYAZDQ9te0hZFM5osAbkUYAZDQ9rfNGWHyKuBehBEACe1AW8/I0EyGaQC3IowASGgVhxskScP60zMCuBVhBEDCCkeM9h5u7RnJG9jH4WoAdBfCCICEVRk8plA4opQkD3NGABcjjABIWJ9U10uScgekK8nLNUYAtyKMAEhYez5tnS/CEA3gboQRAAnrk09be0ZGDEh3uBIA3YkwAiBh7WobpskbSBgB3IwwAiBhfVhVJ0n6XE4/hysB0J0IIwASUkOoJXaNkYsDhBHAzQgjABLSzoNHZYyU3dengX39TpcDoBsRRgAkpB2VbUM09IoArkcYAZCQCCNA70EYAZCQ3v5LrSTp80MzHK4EQHcjjABIOOGI0fb9rWFkXG6Ws8UA6HaEEQAJZ+fBo2oIhdXHl6TzB/V1uhwA3YwwAiDhvLW3RpI0Zngm96QBegHCCICE88buw5KkK0b0d7gSAD2BMAIgoRhj9PquTyVJBRcMdLgaAD2BMAIgoVQcbtC+mkalJHk04bwBTpcDoAcQRgAklD99VC2pdRVNmi/J4WoA9ATCCICE8r/vVUmS/mZUwOFKAPQUwgiAhFF3rFnlH7f2jFz3ecII0FsQRgAkjBffrVJz2OiCQX10AdcXAXoNwgiAhPH0pr2SpK9ePszhSgD0JMIIgITw8aGjemP3YXk80tfGD3e6HAA9iDACICE89souSdK1owIakpnmcDUAehJhBIDj9h5u0HNb90mSbv/r8x2uBkBPI4wAcNzDa95XKBzRVRdmazwXOgN6HcIIAEetfvuA/ri9Ul6PdN/fXeJ0OQAcQBgB4Jh399fq7v9+W5I055oLNSonw+GKADiBMALAEe8fCGr2b9/U0aYWXXn+AH3vby5yuiQADjmrMLJ06VLl5eUpNTVV+fn52rhxY6f7P/300xo1apRSU1M1ZswYrVmz5qyKBWA/Y4ye2rRXX3v0NR2sa9KonH761c0T5Evm30ZAbxV361+1apWKi4u1cOFCbdmyRWPHjlVRUZEOHjzY4f6vvfaapk+frltuuUVbt27VjTfeqBtvvFHbt28/5+IB2CPUEtEf3tqvG5f+WXc987YaQmFNvnCgVn7nSmWmpThdHgAHeYwxJp4D8vPz9YUvfEG/+MUvJEmRSES5ubn63ve+p/nz55+0/7Rp01RfX68XXnghtu3KK6/UuHHjtGzZsjN6z2AwqMzMTNXW1iojgzFlINE1hsKqDB7TJ5/W670DQW2rqNGfd1arPhSWJKWlJGle4UW65aqRSk6iRwRwqzP9/Z0cz4uGQiFt3rxZJSUlsW1er1eFhYUqLy/v8Jjy8nIVFxe321ZUVKTnn3/+lO/T1NSkpqam2M/BYDCeMs/Ybzbs1t7DDWd1bGcZ7nTprrP4Zzo5+nSxsbOnTx85z+F9z/LznP7Y7nvfzp7u/Dye5vN0emynh3bb+57+VJz9331z2KixuUUNobAaQ2E1hMI60hBS3bGWDvcf3M+vm/JHaEb+eRrUz9/5iwPoNeIKI9XV1QqHwwoE2t9NMxAI6IMPPujwmMrKyg73r6ysPOX7lJaW6oEHHointLOy+u392lJR0+3vA/RG6b4kDc1K0yVDMnTpkAxNvnCgRg/NlNfrcbo0AAkmrjDSU0pKStr1pgSDQeXm5nb5+3xt/HBNuiD7lM97TvOd2enTpzm4s2c7O9TT+bue5tjOdXrs6U7GWb6u1PlnOpe/g+58305ft5ODz+nv4Fzet5v+X072epTuS1JaSpLSfclK8yUpMy1ZgYxU9fUnn9P/NwB6j7jCSHZ2tpKSklRVVdVue1VVlXJycjo8JicnJ679Jcnv98vv7/4u3Bn553X7ewAAgM7FNXPM5/Np/PjxKisri22LRCIqKytTQUFBh8cUFBS021+S1q1bd8r9AQBA7xL3ME1xcbFmzZqlCRMmaOLEiVqyZInq6+s1e/ZsSdLMmTM1bNgwlZaWSpLmzp2rL37xi3rkkUc0ZcoUrVy5Ups2bdJjjz3WtZ8EAABYKe4wMm3aNB06dEgLFixQZWWlxo0bp7Vr18YmqVZUVMjrPd7hMmnSJD355JO67777dM899+iiiy7S888/r9GjR3fdpwAAANaK+zojTuA6IwAA2OdMf39ztSEAAOAowggAAHAUYQQAADiKMAIAABxFGAEAAI4ijAAAAEcRRgAAgKMIIwAAwFGEEQAA4Ki4LwfvhOhFYoPBoMOVAACAMxX9vX26i71bEUbq6uokSbm5uQ5XAgAA4lVXV6fMzMxTPm/FvWkikYj279+vfv36yePxdNnrBoNB5ebmau/evdzzphtxnnsO57pncJ57Bue5Z3TneTbGqK6uTkOHDm13E93PsqJnxOv1avjw4d32+hkZGfyP3gM4zz2Hc90zOM89g/PcM7rrPHfWIxLFBFYAAOAowggAAHBUrw4jfr9fCxculN/vd7oUV+M89xzOdc/gPPcMznPPSITzbMUEVgAA4F69umcEAAA4jzACAAAcRRgBAACOIowAAABH9eowsnTpUuXl5Sk1NVX5+fnauHGj0yVZo7S0VF/4whfUr18/DR48WDfeeKN27NjRbp9jx45pzpw5GjhwoPr27auvfe1rqqqqardPRUWFpkyZovT0dA0ePFg/+MEP1NLS0pMfxSqLFi2Sx+PRvHnzYts4z11n3759+od/+AcNHDhQaWlpGjNmjDZt2hR73hijBQsWaMiQIUpLS1NhYaE++uijdq9x+PBhzZgxQxkZGcrKytItt9yio0eP9vRHSVjhcFj333+/Ro4cqbS0NF1wwQV68MEH2927hPMcv1dffVVTp07V0KFD5fF49Pzzz7d7vqvO6dtvv62rr75aqampys3N1Y9//OOu+QCml1q5cqXx+XxmxYoV5t133zW33nqrycrKMlVVVU6XZoWioiLz29/+1mzfvt1s27bN3HDDDWbEiBHm6NGjsX1uu+02k5uba8rKysymTZvMlVdeaSZNmhR7vqWlxYwePdoUFhaarVu3mjVr1pjs7GxTUlLixEdKeBs3bjR5eXnmsssuM3Pnzo1t5zx3jcOHD5vzzjvP/OM//qN54403zK5du8yLL75odu7cGdtn0aJFJjMz0zz//PPmrbfeMl/+8pfNyJEjTWNjY2yfv/3bvzVjx441r7/+uvnTn/5kLrzwQjN9+nQnPlJCeuihh8zAgQPNCy+8YHbv3m2efvpp07dvX/Ozn/0stg/nOX5r1qwx9957r3n22WeNJPPcc8+1e74rzmltba0JBAJmxowZZvv27eZ3v/udSUtLM7/61a/Ouf5eG0YmTpxo5syZE/s5HA6boUOHmtLSUgerstfBgweNJPPKK68YY4ypqakxKSkp5umnn47t8/777xtJpry83BjT2ni8Xq+prKyM7fPoo4+ajIwM09TU1LMfIMHV1dWZiy66yKxbt8588YtfjIURznPXufvuu81VV111yucjkYjJyckxP/nJT2LbampqjN/vN7/73e+MMca89957RpJ58803Y/v88Y9/NB6Px+zbt6/7irfIlClTzLe+9a122/7+7//ezJgxwxjDee4Knw0jXXVOf/nLX5r+/fu3+964++67zcUXX3zONffKYZpQKKTNmzersLAwts3r9aqwsFDl5eUOVmav2tpaSdKAAQMkSZs3b1Zzc3O7czxq1CiNGDEido7Ly8s1ZswYBQKB2D5FRUUKBoN69913e7D6xDdnzhxNmTKl3fmUOM9d6Q9/+IMmTJigb3zjGxo8eLAuv/xyLV++PPb87t27VVlZ2e5cZ2ZmKj8/v925zsrK0oQJE2L7FBYWyuv16o033ui5D5PAJk2apLKyMn344YeSpLfeeksbNmzQ9ddfL4nz3B266pyWl5frr/7qr+Tz+WL7FBUVaceOHTpy5Mg51WjFjfK6WnV1tcLhcLsvZ0kKBAL64IMPHKrKXpFIRPPmzdPkyZM1evRoSVJlZaV8Pp+ysrLa7RsIBFRZWRnbp6O/g+hzaLVy5Upt2bJFb7755knPcZ67zq5du/Too4+quLhY99xzj95880390z/9k3w+n2bNmhU7Vx2dyxPP9eDBg9s9n5ycrAEDBnCu28yfP1/BYFCjRo1SUlKSwuGwHnroIc2YMUOSOM/doKvOaWVlpUaOHHnSa0Sf69+//1nX2CvDCLrWnDlztH37dm3YsMHpUlxn7969mjt3rtatW6fU1FSny3G1SCSiCRMm6OGHH5YkXX755dq+fbuWLVumWbNmOVydezz11FN64okn9OSTT+rzn/+8tm3bpnnz5mno0KGc516sVw7TZGdnKykp6aQVB1VVVcrJyXGoKjvdcccdeuGFF/TSSy9p+PDhse05OTkKhUKqqalpt/+J5zgnJ6fDv4Poc2gdhjl48KCuuOIKJScnKzk5Wa+88or+/d//XcnJyQoEApznLjJkyBBdeuml7bZdcsklqqiokHT8XHX2vZGTk6ODBw+2e76lpUWHDx/mXLf5wQ9+oPnz5+ub3/ymxowZo5tvvll33nmnSktLJXGeu0NXndPu/C7plWHE5/Np/PjxKisri22LRCIqKytTQUGBg5XZwxijO+64Q88995zWr19/Utfd+PHjlZKS0u4c79ixQxUVFbFzXFBQoHfeeaddA1i3bp0yMjJO+qXQW1177bV65513tG3btthjwoQJmjFjRuzPnOeuMXny5JOWp3/44Yc677zzJEkjR45UTk5Ou3MdDAb1xhtvtDvXNTU12rx5c2yf9evXKxKJKD8/vwc+ReJraGiQ19v+V09SUpIikYgkznN36KpzWlBQoFdffVXNzc2xfdatW6eLL774nIZoJPXupb1+v988/vjj5r333jPf+c53TFZWVrsVBzi122+/3WRmZpqXX37ZHDhwIPZoaGiI7XPbbbeZESNGmPXr15tNmzaZgoICU1BQEHs+uuT0uuuuM9u2bTNr1641gwYNYsnpaZy4msYYznNX2bhxo0lOTjYPPfSQ+eijj8wTTzxh0tPTzX/913/F9lm0aJHJysoyv//9783bb79tvvKVr3S4PPLyyy83b7zxhtmwYYO56KKLevWS08+aNWuWGTZsWGxp77PPPmuys7PNXXfdFduH8xy/uro6s3XrVrN161YjySxevNhs3brVfPLJJ8aYrjmnNTU1JhAImJtvvtls377drFy50qSnp7O091z9/Oc/NyNGjDA+n89MnDjRvP76606XZA1JHT5++9vfxvZpbGw03/3ud03//v1Nenq6+epXv2oOHDjQ7nX27Nljrr/+epOWlmays7PNP//zP5vm5uYe/jR2+WwY4Tx3nf/5n/8xo0ePNn6/34waNco89thj7Z6PRCLm/vvvN4FAwPj9fnPttdeaHTt2tNvn008/NdOnTzd9+/Y1GRkZZvbs2aaurq4nP0ZCCwaDZu7cuWbEiBEmNTXVnH/++ebee+9tt1yU8xy/l156qcPv5FmzZhljuu6cvvXWW+aqq64yfr/fDBs2zCxatKhL6vcYc8Jl7wAAAHpYr5wzAgAAEgdhBAAAOIowAgAAHEUYAQAAjiKMAAAARxFGAACAowgjAADAUYQRAADgKMIIAABwFGEEAAA4ijACAAAcRRgBAACO+v8uYKoQ+mjOowAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from opendp.extras.numpy.canonical import one_sided_uniformly_most_powerful_tests\n",
    "\n",
    "ump_tests = one_sided_uniformly_most_powerful_tests(\n",
    "    d_in=estimate.d_in, d_out=estimate.d_out, size=1_000, theta=theta, alpha=0.05, tail=\"left\"\n",
    ")\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(list(range(0, 1001)), ump_tests);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This visual shows how, when the null hypothesis is that the population rate is no smaller than 0.5, the probability of observing an outcome at least as extreme as theta approaches one."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.13.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}