{ "cells": [ { "cell_type": "markdown", "id": "8017a9de", "metadata": {}, "source": [ "Title: Sampling from Probability Distributions\n", "Author: Thomas Breuel\n", "Institution: UniKL" ] }, { "cell_type": "code", "execution_count": 1, "id": "9941f709", "metadata": { "collapsed": false }, "outputs": [], "source": [ "\n", "from pylab import *" ] }, { "cell_type": "markdown", "id": "868253f1", "metadata": {}, "source": [ "# Sampling from Distributions" ] }, { "cell_type": "markdown", "id": "11909aaf", "metadata": {}, "source": [ "We can easily sample for distributions in Python.\n", "The basic distributions are `randn` (normal density) and `rand` (uniform density)." ] }, { "cell_type": "code", "execution_count": 4, "id": "5de74510", "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD9CAYAAABDaefJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGq5JREFUeJzt3V1sU+cdx/GfQ9ILRChBbUxlM7lbgoIhYK8oibYymUFA\nLSNKy5YuKTQqME1Fq0Sp1pZudGESJN3KRWFC6gXtgiqVsgsImiAKq+quY9oCJVxMqZpsCmvivAwW\nXhIoCyRnFymOTfyWxI7tk+9HQjrY5/H5G+KfnzznOc+xGIZhCABgChnJLgAAED+EOgCYCKEOACZC\nqAOAiRDqAGAihDoAmEjEUL99+7aKi4vlcrnkdDq1a9cuSVJNTY3sdrvcbrfcbrdOnz7tb1NbW6v8\n/HwVFBSoqakpsdUDAIJYos1Tv3XrlmbPnq27d+/q8ccf11tvvaWPPvpI2dnZ2rlzZ9C+ra2tqqqq\n0rlz5+Tz+bRmzRq1tbUpI4NfCABgOkRN29mzZ0uShoaGNDw8rJycHElSqO+ChoYGVVZWKisrSw6H\nQ3l5eWpubo5zyQCAcKKG+sjIiFwul6xWq1atWqUlS5ZIkg4ePKjly5dr69atunbtmiSpu7tbdrvd\n39Zut8vn8yWodADA/TKj7ZCRkaGLFy/q+vXrWrdunbxer1544QW98cYbkqTdu3fr5Zdf1uHDh0O2\nt1gsMT0GAIgu2souMQ92P/jgg1q/fr3Onz+v3NxcWSwWWSwWbdu2zT/EYrPZ1NnZ6W/T1dUlm80W\ntrBU//OrX/0q6TWYpc50qJE6qTPV/8QiYqhfuXLFP7Ty1Vdf6cyZM3K73ert7fXvc/z4cRUWFkqS\nysrKdPToUQ0NDamjo0Pt7e0qKiqKqRAAwNRFHH7p6elRdXW1RkZGNDIyos2bN2v16tV67rnndPHi\nRVksFj366KN65513JElOp1MVFRVyOp3KzMzUoUOHGGoBgGkUMdQLCwt14cKFcY8fOXIkbJvXX39d\nr7/++tQrSwEejyfZJcQkHepMhxol6ow36px+UeepJ+SgFkvM40MAgFGxZCdXBQGAiRDqAGAihDoA\nmAihDgAmQqgDgIkQ6gBgIoQ6AJgIoY4Zae7c+f71i+bOnZ/scoC44eIjzEijy1fc+xnk5xHpgYuP\nAGCGIdQBwEQIdQAwEUIdAEyEUAcAEyHUAcBECHUAMBFCHQBMhFAHABMh1AHARAh1ADARQh0ATIRQ\nBwATiRjqt2/fVnFxsVwul5xOp3bt2iVJ6u/vV2lpqRYtWqS1a9fq2rVr/ja1tbXKz89XQUGBmpqa\nEls9EGcsyYt0F3Xp3Vu3bmn27Nm6e/euHn/8cb311ls6efKkHnroIb3yyit68803dfXqVdXV1am1\ntVVVVVU6d+6cfD6f1qxZo7a2NmVkBH93sPQuki3c0rssyYtUFpeld2fPni1JGhoa0vDwsHJycnTy\n5ElVV1dLkqqrq3XixAlJUkNDgyorK5WVlSWHw6G8vDw1NzdP9X0AAGKUGW2HkZERffvb39a//vUv\nvfDCC1qyZIn6+vpktVolSVarVX19fZKk7u5ulZSU+Nva7Xb5fL6Qr1tTU+Pf9ng88ng8U3gbAGA+\nXq9XXq93Qm2ihnpGRoYuXryo69eva926dfr444+Dnr83/hhOuOcCQx1ItLlz52tg4GqyywAm5P4O\n7549e6K2iXn2y4MPPqj169frs88+k9VqVW9vrySpp6dHubm5kiSbzabOzk5/m66uLtlstlgPASTM\naKAbAX9ikclJU6SdiKF+5coV/8yWr776SmfOnJHb7VZZWZnq6+slSfX19SovL5cklZWV6ejRoxoa\nGlJHR4fa29tVVFSU4LcAJMpd3fsSoJePdBFx+KWnp0fV1dUaGRnRyMiINm/erNWrV8vtdquiokKH\nDx+Ww+HQsWPHJElOp1MVFRVyOp3KzMzUoUOHIg7NAADiK+qUxoQclCmNmGbBUxUlKbYpjUxvRCqJ\ny5RGAFLg+LrF8kDANuPtSC1RZ78AkMbG16XgHrw0MMAQI1IHPXUAMBFCHQBMhOEX4OvxcsAM6KkD\nAfPRgXRHqAPTgCV9MV0IdZhWYJAmVvTlBAKXKeDqVCQSoY60F64XHLzeSyIFLicwQI8cScUVpUh7\nsd7wItIVpZPfjrwfN99APMWSncx+gckwkwUzG8MvMBlmsmBmo6cOJAy/NWD60VMHEobfGjD9CHUA\nMBFCHQBMhFAHABMh1AHARAh1ADARQh0ATIRQBwATIdSBaRd9VUdgsriiFJh2Yzex5qbViLeIPfXO\nzk6tWrVKS5Ys0dKlS3XgwAFJUk1Njex2u9xut9xut06fPu1vU1tbq/z8fBUUFKipqSmx1QMAgkRc\nere3t1e9vb1yuVwaHBzUY489phMnTujYsWPKzs7Wzp07g/ZvbW1VVVWVzp07J5/PpzVr1qitrU0Z\nGcHfHSy9i3gKv8Ru8pfejWWbzwJiFUt2RuypL1iwQC6XS5I0Z84cLV68WD6fT5JCvnBDQ4MqKyuV\nlZUlh8OhvLw8NTc3T7Z+AMAExTymfunSJbW0tKikpERnz57VwYMHdeTIEa1YsUL79+/XvHnz1N3d\nrZKSEn8bu93u/xK4X01NjX/b4/HI4/FM+k0AgBl5vV55vd4JtYnpzkeDg4PyeDz65S9/qfLycv3n\nP//Rww8/LEnavXu3enp6dPjwYb344osqKSnRs88+K0natm2bnnzyST399NPBB2X4BXHE8AtmiikP\nv0jSnTt3tHHjRm3atEnl5eWSpNzcXP+UrG3btvmHWGw2mzo7O/1tu7q6ZLPZpvIeMMOEu98ogNhE\nDHXDMLR161Y5nU7t2LHD/3hPT49/+/jx4yosLJQklZWV6ejRoxoaGlJHR4fa29tVVFSUoNJhRoE3\nix7dBjAREcfUz549q/fff1/Lli2T2+2WJO3bt08ffPCBLl68KIvFokcffVTvvPOOJMnpdKqiokJO\np1OZmZk6dOgQd34BIgq+O1J2do5u3OhPYj1IdzGNqcf9oIypI4z7x8dj+TlJ9zH1+1+LzwbCicuY\nOpBogePoAKaGUEfSBY6jA5gaQh1JQe8cSAxCHUkx1d45XwpAaKzSiBQWODMkS9Kd+54PPNkIQCLU\nkdLGlqgNPeMEwP0YfgEAEyHUAcBECHUAMBFCHQBMhFAHABMh1AHARAh1ADARQh1IUdwwBJPBxUdA\nSsm8b+mD0QuuBga42AqxoacOpJR7V9GyYiUmh1AHABMh1AHARAh1ADARQh0ATIRQBwATIdQBwEQI\ndQAwkYih3tnZqVWrVmnJkiVaunSpDhw4IEnq7+9XaWmpFi1apLVr1+ratWv+NrW1tcrPz1dBQYGa\nmpoSWz0AIIjFMIywVzn09vaqt7dXLpdLg4ODeuyxx3TixAm99957euihh/TKK6/ozTff1NWrV1VX\nV6fW1lZVVVXp3Llz8vl8WrNmjdra2pSREfzdYbFYFOGwmAFGr5oMdau6WLYn0yae7ZPzWnxmEEt2\nRuypL1iwQC6XS5I0Z84cLV68WD6fTydPnlR1dbUkqbq6WidOnJAkNTQ0qLKyUllZWXI4HMrLy1Nz\nc3M83gsAIAYxr/1y6dIltbS0qLi4WH19fbJarZIkq9Wqvr4+SVJ3d7dKSkr8bex2u3w+X8jXq6mp\n8W97PB55PJ5JlA8A5uX1euX1eifUJqZQHxwc1MaNG/X2228rOzs76Ll7q8iFE+65wFAHAIx3f4d3\nz549UdtEnf1y584dbdy4UZs3b1Z5ebmk0d55b2+vJKmnp0e5ubmSJJvNps7OTn/brq4u2Wy2Cb0J\nmEfg0rEsHwtMj4ihbhiGtm7dKqfTqR07dvgfLysrU319vSSpvr7eH/ZlZWU6evSohoaG1NHRofb2\ndhUVFSWwfKSygYGrGltx0Pj67wASKeLsl7/85S/63ve+p2XLlvmHUWpra1VUVKSKigp9+eWXcjgc\nOnbsmObNmydJ2rdvn959911lZmbq7bff1rp168YflNkvM0LwDBcpcAYHs1+Y/YKJiyU7I4Z6ohDq\nM8P4UM/S6Hrh9xDqhDomYspTGoH44gYQQKIR6kBayOSEM2LCPUqBtHDvtxzuV4rI6KkDgIkQ6gBg\nIoQ6AJgIoQ6kHU6aIjxOlAJph5OmCI+eOgCYCKEOACZCqCOuAldmxHTIZCVMBGFMHXE1tjKjNLp2\nCRJrbHxdYowd9NQBwFQIdQAwEUIdAEyEUAcAEyHUAcBECHUAMBFCHQBMhFAHABMh1AHARAh1ADAR\nQh0ATCRqqG/ZskVWq1WFhYX+x2pqamS32+V2u+V2u3X69Gn/c7W1tcrPz1dBQYGampoSUzUAIKSo\nof7888+rsbEx6DGLxaKdO3eqpaVFLS0teuKJJyRJra2t+vDDD9Xa2qrGxkZt375dIyMjiakcADBO\n1FBfuXKlcnJyxj1uGMa4xxoaGlRZWamsrCw5HA7l5eWpubk5PpUCAKKa9NK7Bw8e1JEjR7RixQrt\n379f8+bNU3d3t0pKSvz72O12+Xy+kO1ramr82x6PRx6PZ7KlAIApeb1eeb3eCbWZVKi/8MILeuON\nNyRJu3fv1ssvv6zDhw+H3DfczRICQx0AMN79Hd49e/ZEbTOp2S+5ubn+O61s27bNP8Ris9nU2dnp\n36+rq0s2m20yh0Aa4W5HQOqYVKj39PT4t48fP+6fGVNWVqajR49qaGhIHR0dam9vV1FRUXwqRcoa\nu9vR+PMsAKZX1OGXyspKffLJJ7py5YoWLlyoPXv2yOv16uLFi7JYLHr00Uf1zjvvSJKcTqcqKirk\ndDqVmZmpQ4cO0XsDgGlkMUJNY0n0QS2WkLNnkJ5Gv7gD70saajvScxPdTnb7VH0tScrS6H1Lpezs\nHN240S+YRyzZyY2nAVMZuxE1N6GemVgmAJgBAk9mz507P9nlIIEIdUwKM17SS+DJ7NFtmBWhjklh\nxguQmgh1ADARTpQCppXJ8NgMRE8dMK17M2EYIptJCHUAMBFCHQBMhFAHABMh1AHARAh1RBR4kZHF\n8gAXHAEpjimNiGjsIiNp/EJSAFINPXUAMBFCHQBMhFAHABMh1AHARAh1ADARQh0ATIRQBwATIdQB\nwEQIdQAwEUIdAEwkaqhv2bJFVqtVhYWF/sf6+/tVWlqqRYsWae3atbp27Zr/udraWuXn56ugoEBN\nTU2JqRrAFGT61/CZO3d+sotBnEUN9eeff16NjY1Bj9XV1am0tFRtbW1avXq16urqJEmtra368MMP\n1draqsbGRm3fvl0jIyOJqRzAJI3dEWl0bR+YSdRQX7lypXJycoIeO3nypKqrqyVJ1dXVOnHihCSp\noaFBlZWVysrKksPhUF5enpqbmxNQNgAglEmt0tjX1yer1SpJslqt6uvrkyR1d3erpKTEv5/dbpfP\n5wv5GjU1Nf5tj8cjj8czmVIAwLS8Xq+8Xu+E2kx56d1o62uHey4w1AEk39y584OGY7Kzc3TjRn8S\nK8L9Hd49e/ZEbTOp2S9Wq1W9vb2SpJ6eHuXm5kqSbDabOjs7/ft1dXXJZrNN5hAAptnY2vmMt6ez\nSYV6WVmZ6uvrJUn19fUqLy/3P3706FENDQ2po6ND7e3tKioqil+1AICIog6/VFZW6pNPPtGVK1e0\ncOFC/frXv9Zrr72miooKHT58WA6HQ8eOHZMkOZ1OVVRUyOl0KjMzU4cOHeLWZwAwjSyGYRjRd4vz\nQS0WJeGwmITRL+Vwt7ObyPZU28fztVKpluS/r3ufxeD/6+DnkBpiyU6uKEWQ4BtN81sWkG648TSC\nBN9oWuIG00B6oacOACZCqAMIgzVi0hHDLwDCuLdGjDQwwDBcuqCnjqCTowDSGz113HdylGCfWTL5\nMjcZeurAjDa2DC/MgVAHABMh1AHARAj1GYqTo4A5EeozVPAyqwDMglAHABMh1AHARAh1ADFgyYB0\nwcVHAGLAkgHpgp76DMKMF8D8CPUZhBkvgPkR6gBgIoQ6AJgIoQ4AJkKoA5ggpjemMqY0Apggpjem\nsimFusPh0Ny5czVr1ixlZWWpublZ/f39euaZZ/Tvf/9bDodDx44d07x58+JVLwAggikNv1gsFnm9\nXrW0tKi5uVmSVFdXp9LSUrW1tWn16tWqq6uLS6EAgOimPKZuGMFznk+ePKnq6mpJUnV1tU6cODHV\nQwAAYjSl4ReLxaI1a9Zo1qxZ+ulPf6qf/OQn6uvrk9VqlSRZrVb19fWFbFtTU+Pf9ng88ng8UykF\nAEzH6/XK6/VOqI3FuL+rPQE9PT165JFHdPnyZZWWlurgwYMqKyvT1atX/fvMnz9f/f39wQe1WMb1\n8JEYc+fO//pK0nsCbzAdbTvW/dLttVKplnR/X1kaPXEqZWfn6MaN4M864iuW7JzS8MsjjzwiSXr4\n4Yf11FNPqbm5WVarVb29vZJGQz83N3cqh0AMAtd0uX+KGUsDILHGblwd2HkI/Jlk6uP0mnSo37p1\nSwMDA5KkmzdvqqmpSYWFhSorK1N9fb0kqb6+XuXl5fGpFGEFBndwrxxIjuDOBD+X02nSY+p9fX16\n6qmnJEl3797Vs88+q7Vr12rFihWqqKjQ4cOH/VMaMZ0yWYURmMGmNKY+6YMyph5XoyFuljFas449\nz4z3de9zHfwzGfwcJi/hY+oAgNRCqAOAiRDqAOIkM6Y7a0WarYWpY0EvAHEyttDX6Fh7aGMzY1gQ\nLBHoqacp7jcKIBRCPU1xURGAUBh+ATANuH5iutBTBzANxpYTQGIR6gCSiFvjxRvDLwCSiFvjxRs9\n9RTHnF4AE0FPPQWFWwOdngyAaAj1FDA+xKVYLuIAgPsR6ikg8Aq7UeGCnGlhACJjTD2tMC0MZhZ+\nJgznlmJHTx1Aigg/E4b1YmJHT30a0dsAkGj01KcRvQ0AiUZPPU7C9cLDr6YY29rTADARhHqcBK6a\nGDg9Mfxqipz0BCZnYksLBHasZsLQJ6GeEPTCgcQZ6xANDAxEDevgjpUR4poQcyHUvzbRk5j3f/sH\noxcOTE1mhM9XoNABb7E8EKGtuRcRS0ioNzY2qqCgQPn5+XrzzTcTcYi4CzV84vV6Y9qf8I6FN9kF\nmIw32QXEyDvJdoEdo1g/X4Ft7kRoO7GefrqJe6gPDw/rZz/7mRobG9Xa2qoPPvhAn3/+ebwPE7NY\nTmCG6wl4vV6mIcaNN9kFmIw32QXEyJvsAqK4F/C/ChqWSefPfdynNDY3NysvL08Oh0OS9OMf/1gN\nDQ1avHjxhF/r5s2bOn36tEZGRiRJNptN3/3udyf0GuGmEcZ6aX5w+yzGyYEZIJ2nH8c91H0+nxYu\nXOj/u91u19///vdJvdann36qH/3oR0GPGcZYEAcuhJWdnaMbN/rHPR4s1rVTwu0X293SAaSjcJ/7\nwMezNDq0o3F/D5dBgY9HMpk2od9FnMXak51sjzdcu4GBqxFeM9bHY9kv3Has+6Xba6VSLbyvMXu+\n/pMKtUTbDlVnqv4bh9vnnjth/x4ugyJnU2iTaXNP3EPdZrOps7PT//fOzk7Z7fagfQJ72wCA+In7\nidIVK1aovb1dly5d0tDQkD788EOVlZXF+zAAgBDi3lPPzMzU7373O61bt07Dw8PaunXrpE6SAgAm\nLiHz1J944gl98cUX+uc//6ldu3aF3W///v3KyMhQf//kTggk2u7du7V8+XK5XC6tXr06aFgplfz8\n5z/X4sWLtXz5cj399NO6fv16sksK6Q9/+IOWLFmiWbNm6cKFC8kuZ5x0uL5iy5YtslqtKiwsTHYp\nEXV2dmrVqlVasmSJli5dqgMHDiS7pHFu376t4uJiuVwuOZ3OiFmVCoaHh+V2u7Vhw4bIOxpJ8uWX\nXxrr1q0zHA6H8d///jdZZUR048YN//aBAweMrVu3JrGa8Jqamozh4WHDMAzj1VdfNV599dUkVxTa\n559/bnzxxReGx+MxPvvss2SXE+Tu3bvGt771LaOjo8MYGhoyli9fbrS2tia7rHH+/Oc/GxcuXDCW\nLl2a7FIi6unpMVpaWgzDMIyBgQFj0aJFKfnvefPmTcMwDOPOnTtGcXGx8emnnya5ovD2799vVFVV\nGRs2bIi4X9KWCdi5c6d+85vfJOvwMcnOzvZvDw4O6qGHHkpiNeGVlpYqI2P0v7K4uFhdXV1Jrii0\ngoICLVq0KNllhBR4fUVWVpb/+opUs3LlSuXk5CS7jKgWLFggl8slSZozZ44WL16s7u7uJFc13uzZ\nsyVJQ0NDGh4e1vz5qXmhUVdXl06dOqVt27ZFnWiSlFBvaGiQ3W7XsmXLknH4CfnFL36hb3zjG6qv\nr9drr72W7HKievfdd/Xkk08mu4y0E+r6Cp/Pl8SKzOPSpUtqaWlRcXFxsksZZ2RkRC6XS1arVatW\nrZLT6Ux2SSG99NJL+u1vf+vvvEWSsJtklJaWqre3d9zje/fuVW1trZqamvyPRfvmSaRwde7bt08b\nNmzQ3r17tXfvXtXV1emll17Se++9l4Qqo9cpjf7bPvDAA6qqqpru8vxiqTMVcaVwYgwODuqHP/yh\n3n77bc2ZMyfZ5YyTkZGhixcv6vr161q3bp28Xq88Hk+yywryxz/+Ubm5uXK73RHXo7onYaF+5syZ\nkI//4x//UEdHh5YvXy5p9NeKxx57TM3NzcrNzU1UOWGFq/N+VVVVSe0BR6vz97//vU6dOqWPPvpo\nmioKLdZ/z1QTy/UVmJg7d+5o48aN2rRpk8rLy5NdTkQPPvig1q9fr/Pnz6dcqP/1r3/VyZMnderU\nKd2+fVs3btzQc889pyNHjoRuMC0j/BGk8onStrY2//aBAweMTZs2JbGa8E6fPm04nU7j8uXLyS4l\nJh6Pxzh//nyyywhy584d45vf/KbR0dFh/O9//0vZE6WGYRgdHR0pf6J0ZGTE2Lx5s7Fjx45klxLW\n5cuXjatXrxqGYRi3bt0yVq5cafzpT39KclWReb1e4wc/+EHEfZK+nnoq/9q7a9cuFRYWyuVyyev1\nav/+/ckuKaQXX3xRg4ODKi0tldvt1vbt25NdUkjHjx/XwoUL9be//U3r16/XE088keyS/AKvr3A6\nnXrmmWdS8vqKyspKfec731FbW5sWLlyYtOHAaM6ePav3339fH3/8sdxut9xutxobG5NdVpCenh59\n//vfl8vlUnFxsTZs2KDVq1cnu6yoomWmxTC4Zh8AzCLpPXUAQPwQ6gBgIoQ6AJgIoQ4AJkKoA4CJ\nEOoAYCL/B73iZ+4f++mAAAAAAElFTkSuQmCC\n" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data = randn(10000)\n", "_=hist(data,bins=100)" ] }, { "cell_type": "markdown", "id": "dd4c6194", "metadata": {}, "source": [ "# Binomial Distributions" ] }, { "cell_type": "markdown", "id": "9b5a6924", "metadata": {}, "source": [ "If we want a binomial variable with $p=0.1$, we can sample as follows:\n", "\n" ] }, { "cell_type": "code", "execution_count": 69, "id": "220e4092", "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "109" ] }, "execution_count": 69, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n = 1000\n", "p = 0.1\n", "sum(rand(n)maxs)*1.0/len(maxs)" ] }, { "cell_type": "markdown", "id": "a55ccc58", "metadata": {}, "source": [ "Interestingly, you can get much greater extreme value differences without any bias at all, simply by having more variance in one population than in another.\n", "\n" ] }, { "cell_type": "code", "execution_count": 104, "id": "d4d44a1c", "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.82399999999999995" ] }, "execution_count": 104, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD9CAYAAACsq4z3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHWRJREFUeJzt3X9sU+f9L/D3SQjxKIwfFXG4Tq4Src1KaCD2WDNNcGua\nJgzzTZSMqBpjYAGZENM2tfAP/AFLriaWqUIVZdWEqqwKQ4JWm2Bpl2WFC55pEIm6xFRdf6Rx4TaJ\nYt9SCGqc+Pe5f+R7XCexY/vYjn1O3i8JYY7Pj8/RA58czvN5nkcQRVEEEREpXk6mAyAiotRgQici\nUgkmdCIilWBCJyJSCSZ0IiKVYEInIlKJuBJ6IBCAXq9HXV0dAODBgweoqalBWVkZamtrMT4+ntYg\niYgotrgS+pkzZ1BeXg5BEAAAbW1tqKmpweDgIKqrq9HW1pbWIImIKLaYCX1kZARdXV1obm6GNAap\ns7MTZrMZAGA2m3HlypX0RklERDHFTOgvvfQSXn75ZeTkfLOr0+mEVqsFAGi1WjidzvRFSEREcVky\n35fvvPMOCgoKoNfrYbFYIu4jCELoVUyk74iIKHFyZmWZ9wn91q1b6OzsRGlpKXbv3o3r169j7969\n0Gq1cDgcAICxsTEUFBTMG5Raf/3mN7/JeAy8N94f7099v+SaN6GfOnUKw8PDuHv3Li5duoTnnnsO\nf/7zn1FfX4+Ojg4AQEdHBxoaGmQHQEREqZFQHbr0CuXYsWO4evUqysrKcP36dRw7diwtwRERUfzm\nfYce7tlnn8Wzzz4LAFizZg2uXbuWtqCUwmg0ZjqEtFHzvQG8P6VT+/3JJYjJvLCJdXJBSOp9EBHR\nYiQ3d3LoPxGRSjChExGpBBM6EZFKMKETEakEEzoRkUrEXbZIlAxTowmOhw4Uri5E1+WuTIdDpEp8\nQqcF4XjogO6QDo6HjkyHQqRaTOhERCrBhE5EpBJM6EREKsGETkSkEkzoREQqwYRORKQSTOhERCrB\nhE5EpBJM6EREKsGh/wpgMu2GwzEBACgsXI6urosZjoiIshETugI4HBPQ6d4GAIyO1mU4GiLKVjFf\nubjdblRVVaGyshLl5eU4fvw4AKClpQVFRUXQ6/XQ6/Xo7u5Oe7BERBRdzCd0jUaDGzduYNmyZfD7\n/diyZQvee+89CIKAI0eO4MiRIwsRJxERxRBXp+iyZcsAAF6vF4FAAKtXrwYALgBNRJRF4krowWAQ\nlZWV0Gq12LZtGzZs2AAAOHv2LDZt2oSDBw9ifHw8rYESEdH84uoUzcnJgc1mw6NHj7B9+3ZYLBYc\nPnwYJ0+eBACcOHECR48eRXt7+5xjW1paQp+NRiOMRmNKAiciUguLxQKLxZL0eRKqclm5ciV27tyJ\n999/f0Zibm5uRl1d5OqL8IRORERzzX7YbW1tlXWemK9c7t+/H3qdMjU1hatXr0Kv18Ph+GblmcuX\nL6OiokJWAERElBoxn9DHxsZgNpsRDAYRDAaxd+9eVFdXY9++fbDZbBAEAaWlpTh37txCxEtERFHE\nTOgVFRXo7++fs/38+fNpCYiIiOThSFGKytRoguOhA4WrC9F1uSvT4RBRDJyci6JyPHRAd0gHx0NH\n7J2JKOOY0ImIVIIJnYhIJZjQiYhUggmdiEglWOWyQOQsUiEdY7ffhU6X2HVSuRDG0OAQDEYDALDi\nhSiL8Ql9gUiLVOh0b4cSe7zH+HzBhK8T7zXi4YcfukM6VrwQZTkmdCIilWBCJyJSCSZ0IiKVYKco\nZRVpugGAHbBEiWJCp6wiTTcAAKPnRjMcDZGy8JULEZFKMKETEakEEzoRkUowoRMRqQQ7RSmlEqlS\nkVPRwioYouiY0CmlEqlSkVPRwioYoujmfeXidrtRVVWFyspKlJeX4/jx4wCABw8eoKamBmVlZait\nrcX4+PiCBEtERNHNm9A1Gg1u3LgBm82GDz74ADdu3MB7772HtrY21NTUYHBwENXV1Whra1uoeImI\nKIqYnaLLli0DAHi9XgQCAaxevRqdnZ0wm80AALPZjCtXrqQ3SiIiiinmO/RgMAiDwQC73Y7Dhw9j\nw4YNcDqd0Gq1AACtVgun0xn1+JaWltBno9EIo9GYdNCUWumYQ52I4mexWGCxWJI+T8yEnpOTA5vN\nhkePHmH79u24cePGjO8FQYAgCFGPD0/olJ2kOdRHR+syHQrRojT7Ybe1tVXWeeKuQ1+5ciV27tyJ\nf//739BqtXA4pkvHxsbGUFBQIOviRESUOvMm9Pv374cqWKampnD16lXo9XrU19ejo6MDANDR0YGG\nhob0R0pERPOa95XL2NgYzGYzgsEggsEg9u7di+rqauj1erzwwgtob29HSUkJ3nrrrYWKl4iIopg3\noVdUVKC/v3/O9jVr1uDatWtpC4qIiBLHkaJZSKo6sdvvQqfLdDSZMzQ4BIPRAIDD/Iniwcm5spBU\ndeLzBTMdSkb54YfukA66Q7rQ/C1EFB0TOhGRSjChExGpBBM6EZFKMKETEakEEzphaKQXVpsBn4xa\nYWo0ZTocIpKJZYsEv+DBiiYd4B6F42NWkxApFZ/QiYhUggmdiEglmNCJiFSC79AzYGjoUxgMdUkt\nKDF7UYpMTBdgajSFRnDGOzQ/fDi/3W6HDroZ28O3EVFi+ISeAX7/Euh0b8PhmJB9Dml6AOkcmZgu\nwPHQkfDQ/PDh/L6Ab8728G1ElBgmdCIilWBCJyJSCSZ0IiKVYEInIlIJVrlkETmVKlLFTCqrW0Jx\n/L97cypOolWpRIwtgX1TSaq+4aIYtNjwCT2LyKlUkSpmUlnd8k0cgbnXi1KlEjG2BPZNJan6hoti\n0GITM6EPDw9j27Zt2LBhA55++mm8+uqrAICWlhYUFRVBr9dDr9eju7s77cESEVF0MV+55OXl4ZVX\nXkFlZSUmJibwve99DzU1NRAEAUeOHMGRI0cWIk4iIoohZkIvLCxEYWEhAGD58uVYv349RkdHAQCi\nKKY3OiIiiltCnaL37t3DwMAAfvCDH6Cnpwdnz57F+fPnsXnzZpw+fRqrVq2ac0xLS0vos9FohNFo\nTDZm1UjFFABEpHwWiwUWiyXp88Sd0CcmJtDU1IQzZ85g+fLlOHz4ME6ePAkAOHHiBI4ePYr29vY5\nx4UndJpJ6tAcHa3LdChElEGzH3ZbW1tlnSeuKhefz4ddu3bhZz/7GRoaGgAABQUFEAQBgiCgubkZ\nfX19sgIgIqLUiJnQRVHEwYMHUV5ejhdffDG0fWxsLPT58uXLqKioSE+EREQUl5ivXHp6enDhwgVs\n3LgRer0eAHDq1ClcvHgRNpsNgiCgtLQU586dS3uwREQUXcyEvmXLFgSDcwet7NixIy0BERGRPBz6\nrzCJDvWXhvEDmFNNs7b4f8DldWHKM4EV6Qo4jcKnFuAwfyIO/VecRIf6S8P4Iy2o4fK6sOKn/wsQ\nlDmeIHxqAQ7zJ2JCJyJSDSZ0IiKVYEInIlIJdoqqmKnRhE9Ge/Dh52vxrRXF8HntMDWaAM9KOBwT\n8Hi9czpDhwaH4A/YYf/SAI/PFfXcvb39cLs96O3tR1WVIb03gm86QBdyXnUipeETuoo5HjqQ92MN\nArkuaJp0yPuxBo6HjlBHaaTJ1fzwI+/HGmiadBARvbPU4/FDEPLh8fjTeQsz4lroedWJlIYJnYhI\nJZjQiYhUggmdiEglmNCJiFSCVS5pIg25XyyLV7AKhSjz+ISeJlIlyezh9mrFKhSizGNCJyJSCSZ0\nIiKVYEInIlIJJnQiIpVglUuSlFTN4nK5cKf/AyD3Y9i/NECEb873QbcHPlcuNJoMBUlEsvEJPUlK\nqmYRRQHIWwr81xJomnRzFrYQRQGCkI+gMte7IFr0Yib04eFhbNu2DRs2bMDTTz+NV199FQDw4MED\n1NTUoKysDLW1tRgfH097sEREFF3MhJ6Xl4dXXnkF//nPf3D79m289tpr+Pjjj9HW1oaamhoMDg6i\nuroabW1tCxEvERFFETOhFxYWorKyEgCwfPlyrF+/HqOjo+js7ITZbAYAmM1mXLlyJb2REhHRvBLq\nFL137x4GBgZQVVUFp9MJrVYLANBqtXA6nRGPaWlpCX02Go0wGo2yg6XM83i8sFr7MOmazHQoCTE1\nmkILSReuLkTX5a4MR0T0DYvFAovFkvR54k7oExMT2LVrF86cOYMVK2aucyMIAgRBiHhceEIn5RNF\nERrNM5iYsGQ6lIQ4HjqgOzQ9x8zoudEMR0M00+yH3dbWVlnniavKxefzYdeuXdi7dy8aGhoATD+V\nOxzTTzxjY2MoKCiQFQAREaVGzIQuiiIOHjyI8vJyvPjii6Ht9fX16OjoAAB0dHSEEj0REWVGzFcu\nPT09uHDhAjZu3Ai9Xg8A+N3vfodjx47hhRdeQHt7O0pKSvDWW2+lPVgiIoouZkLfsmULgsFgxO+u\nXbuW8oCIiEgeDv1PkaGhT2Ew1MmaAkA61m6/Cx3Xhkg7aTGOWNUurIwhpeHQ/xTx+5fIngJAOtbn\ni/w/IUotaTEOKVlHI1XGxLMvUTZgQiciUgkmdCIilWBCJyJSCXaKkipIHZ0AYLfboQN7l2nx4RM6\nqYLU0ak7pIMv4It9AJEKMaETEakEEzoRkUowoRMRqQQTOhGRSrDKRSaTaTccjomsHa4/NPQp/N9y\nwvc/OfqUaLHgE7pMDsdEVg/X9/uXIG/pdxAUMx0JES0UJnQiIpVgQiciUgkmdCIilWBCJyJSCSb0\nFAtfrCLSn1PFI3wBq80A99LP8c+ba+Fe+jmsNgM8whcpu4YoBmG19sHt9qC3tz9l5yWi9GBCT7HZ\ni1Wka/EKMc8PTZMOQl0eArkuCHV50DTpIOb5U3gVARrNMxCEfHg8qTwvEaVDzIR+4MABaLVaVFRU\nhLa1tLSgqKgIer0eer0e3d3daQ2SiIhii5nQ9+/fPydhC4KAI0eOYGBgAAMDA/jRj36UtgCJiCg+\nMRP61q1bsXr16jnbRZEjVoiIsonsof9nz57F+fPnsXnzZpw+fRqrVq2KuF9LS0vos9FohNFolHtJ\nioPo94Y6S72T+ZkOh4jiYLFYYLFYkj6PrIR++PBhnDx5EgBw4sQJHD16FO3t7RH3DU/otADyAE2T\nDhP3P4PY5QOQm+mIiCiG2Q+7ra2tss4jq8qloKAAgiBAEAQ0Nzejr69P1sWJiCh1ZCX0sbGx0OfL\nly/PqIAhIqLMiPnKZffu3fjXv/6F+/fvo7i4GK2trbBYLLDZbBAEAaWlpTh37txCxEpERPOImdAv\nXrw4Z9uBAwfSEgwREcnHBS7CSItWFBYuR1fXxZjbFxOXywWrtU9R5apDg0MwGA0AALvdDh3iW4nE\n1GiC46FjxjHh5ypcXYiuy13pCZooCRz6H0ZatMLhmIhr+2IiitPTACiJH37oDumgO6SDL+CL+zjH\nQ8ecY8LP5XjoSEe4REljQiciUgkmdCIilWBCJyJSCXaKUtbp7e2Hx+MPzcNeVWVIy3XkdnSyg5Sy\nFRM6ZR2Pxw+N5hlMCP8nrfOwSx2dADB6bjTtxxGlG1+5EBGpBBM6EZFKMKETEakEEzoRkUqwUzRL\n9d4xwb30c1htBniELzIdTkwejxdWax88Hu+c73p7+0MVK4uVNJ0AwMoYSh8+oWcpj+iAUJcHTZMO\nYl76Kj1SRRRFaDTPRJzrxePxQxDy01qxku2k6QQ4dQClExM6EZFKMKETEakEEzoRkUqwUzQBQ0Of\nwmCoAwDY7Xehi2967aSJfi+stumh5kroII3G5XIhmObh/AtNmgaAHZ2UDfiEngC/fwl0ureh070N\nny+4cBfOAzRNOsV0kEYjioLqOkelaQDY0UnZgAmdiEglYib0AwcOQKvVoqKiIrTtwYMHqKmpQVlZ\nGWprazE+Pp7WIImIKLaYCX3//v3o7u6esa2trQ01NTUYHBxEdXU12tra0hYgERHFJ2ZC37p1K1av\nXj1jW2dnJ8xmMwDAbDbjypUr6YmOiIjiJqvKxel0QqvVAgC0Wi2cTmfUfVtaWkKfjUYjjEajnEuS\nQohiEFZrH9xTk/C5PbBa+zDpmsx0WERZzWKxwGKxJH2epMsWBUGAIAhRvw9P6LQYCNOLU0xYIAj5\noc9EFN3sh93W1lZZ55FV5aLVauFwTJdpjY2NoaCgQNbFiYgodWQl9Pr6enR0dAAAOjo60NDQkNKg\niIgocTET+u7du/HDH/4Qn376KYqLi/HGG2/g2LFjuHr1KsrKynD9+nUcO3ZsIWIlIqJ5xHyHfvHi\nxYjbr127lvJgiIhIPs7lQklxuVywWvsizoOeKGkhDJ8rFxrNzPNHWjhjMeDCGJQIDv2npIjidFVL\nKkgLYQTDfjZI50/FDwwl4sIYlAgmdCIilWBCJyJSCSZ0IiKVYEKPQFrIwmTanelQ6L9JUwoorXPU\n1GiCwWiA3W7PdCi0CDChRyAtZOFwTGQ6FApRZueo1KnpC/gyHQotAkzoREQqwYRORKQSTOhERCrB\nhE5EpBIc+p9leu+Y4F76ObyT+RG/F/1eWG0GiFBuJ5s0nH9qyg2rtQ/5+Zn/azg0OASD0QAAsNvt\n0EEn+3gO0adM4RN6lvGIDgh1edETdh6gadIBgrKqPcKFTxeg0TwDj8ef4YgAP/yhIfZyKlLCj+cQ\nfcoUJnQiIpVgQiciUgkmdCIilch8b1QWMJl2w+GYgN1+F7qwvjBpCoDZ2ynzpLnTe3v7UVVlSOhY\nj8eb1mkEku1gDZ8DPdnj2UG7uPAJHYDDMQGd7m34fMEZ26UpAGZvp8yT5k6X06EqimJapxFItoM1\nfA70ZI9nB+3iktQTeklJCb797W8jNzcXeXl56OvrS1VcRESUoKQSuiAIsFgsWLNmTariISIimZJ+\n5aK02e+IiNQqqYQuCAKef/55bN68Ga+//nqqYiIiIhmSeuXS09ODdevW4csvv0RNTQ2eeuopbN26\ndcY+LS0toc9GoxFGozGZSyZNqmgBgOHhIRQXP5HxKhaP8AWsNsO8Q/7VzOVyIej2wOfKjbmvtNDF\npGtyxnapciWRaQSkcy301APhVTDDd4dRXFoc+h2QV9lCymaxWGCxWJI+T1J/k9etWwcAWLt2LRob\nG9HX1zdvQs8GUkULAHzySTl0urfxySflGY1JzPND06TDxP3PIHb5AMRObGoiigIEIR/BuN7eTU8b\nMDFhmXWO6coVtzuRjnlBxjHJk6pgAOCTo59Ad0gX+l3aRovL7Ifd1tZWWeeR/cplcnISX3/9NYDp\nJ6x3330XFRUVck9HRERJkv2E7nQ60djYCADw+/3Ys2cPamtrUxYYERElRnZCLy0thc1mS2UsRESU\nBA79X2DSfOe9d0yZDkWVpA7WdA7tzxSpMzXacH5pyP9CdKpK1+LUAtmFQ/8XmDTfuUfkkOx0kDpY\n0zm0P1OkztRow/mlIf9ypgtIlHQtTi2QXZjQiYhUggmdiEglmNCJiFSCCZ2ISCUWTZVLtEUsMsX1\n9RCCS72w2gzRF4SmmFwuF6zWPtV1gM4n0tQBQPxTBnABDPVaNE/o0RaxyBQx1w+hLg+aJh0gLJ5k\nlGqiOD18fzEJX0DD5XMlvBgGF8BQr0WT0ImI1I4JnYhIJZjQiYhUggmdiEglVFPlIlWxFBYuR1fX\nxTnbs6W6hRaOtICFe2oSPrcH7/7TAs23ls2piJHmf+nt7UdVlWHGd729/fB4/FHnhZEW1nCHHR/a\nFnZdX5Tzp1ukipjwaphoFTPS50xVwbASRx7VPKFLVSzSakSzt2dLdQstpOkKGH8AEIR8+AOIWBEj\nzf/i8fjnfOfx+OedF0ZaWCP8eGlb+HWjnT/dIlXEhFfDRKuYkT5nqgqGlTjyqCahExEtdkzoREQq\nobh36MFgEB9++CF8Ph+WL1+O7373uwtyXVEUEQy6MT7+7//+88L/95mIaD6KS+hjY2Nobv7fyM0t\nR07OB+jpuTLj+6GhT2Ew1M3pHE2W3/8I3vwxfDj1KwTdU/Dk/F9YbYaoi1V4hC9C3//z5lr4l34N\nq80A16Q9ZTFRaknTCLjDOjJ9rlxoNN/s09vbH/o+UgdrPOefb+GN2eeX9g3fHq2DVeqMnZpyw2rt\nQ35+7H/esztwZ583kWkGInVkhm/Lpk7XdMtUp67iXrmIoojc3LV4/PEWuN1zOzr9/iURO0dTcV3k\nA0u3Pg7hyaUQ8wLQNOmiLlYh5vlD3wdyXaFh/py3JXtJ0wiEd2QGZ+Vrj8c/bwdrPOef74fA7PNL\n+4Zvj9bBKnXGAtPHxtMJO7sDd/YxiUwzEKkjM3xbNnW6plumOnWTSujd3d146qmn8OSTT+L3v/99\nqmJSDO/IV5kOIW3UfG8A70/pLBZLpkPISrITeiAQwC9/+Ut0d3fjo48+wsWLF/Hxxx+nMras51Px\nPxo13xvA+1M6JvTIZCf0vr4+PPHEEygpKUFeXh5+8pOf4G9/+1sqYyMiogTITuijo6MoLi4O/bmo\nqAijo6MpCWo+giAgELiPr776LfLzF64LQBAEwAt4e+5DHPIAEBbs2kRE8RBEmSsD/PWvf0V3dzde\nf/11AMCFCxfQ29uLs2fPfnNygUmPiEgOOalZdtmiTqfD8PBw6M/Dw8MoKipKOiAiIpJH9juLzZs3\n47PPPsO9e/fg9Xrx5ptvor6+PpWxERFRAmQ/oS9ZsgR/+MMfsH37dgQCARw8eBDr169PZWxERJSA\npHoVd+zYgTt37uDxxx/Hm2++ifLychw/fjzivr/+9a/x5JNPYtOmTRgYGEjmsgvG7XajqqoKlZWV\nUe/NYrFg5cqV0Ov10Ov1+O1vf5uBSJMTCASg1+tRV1cX8Xsltl24+e5P6e1XUlKCjRs3Qq/X45ln\nIg90UnL7xbo/Jbff+Pg4mpqasH79epSXl+P27dtz9km47cQUcLlcoiiKos/nE6uqqsSbN2/O+P7v\nf/+7uGPHDlEURfH27dtiVVVVKi67IGLd240bN8S6urpMhJYyp0+fFn/6059GvA8lt51kvvtTevuV\nlJSIX331VdTvld5+se5Pye23b98+sb29XRTF6fwyPj4+43s5bZeSur9ly5YBALxeLwKBANasWTPj\n+87OTpjNZgBAVVUVxsfH4XQ6U3HptIt1b4CyO39HRkbQ1dWF5ubmiPeh5LYDYt8foOz2A+aPX+nt\nB8RuHyW236NHj3Dz5k0cOHAAwPQr7JUrV87YR07bpSShB4NBVFZWQqvVYtu2bSgvL5/xfaSa9ZGR\nkVRcOu1i3ZsgCLh16xY2bdoEk8mEjz76KEORyvPSSy/h5ZdfRk5O5L8KSm47IPb9Kb39BEHA888/\nj82bN4dKiMMpvf1i3Z9S2+/u3btYu3Yt9u/fD4PBgJ///OeYnJycsY+ctktJQs/JyYHNZsPIyAis\nVmvEYbmzf4oqpUY91r0ZDAYMDw/jzp07+NWvfoWGhobMBCrDO++8g4KCAuj1+nmfcpTadvHcn5Lb\nDwB6enowMDCAf/zjH3jttddw8+bNOfsotf2A2Pen1Pbz+/3o7+/HL37xC/T39+Oxxx5DW1vbnP0S\nbbuUDrVcuXIldu7ciffff3/G9tk16yMjI9ApbIHPaPe2YsWK0GuZHTt2wOfz4cGDB5kIMWG3bt1C\nZ2cnSktLsXv3bly/fh379u2bsY+S2y6e+1Ny+wHAunXrAABr165FY2Mj+vr6Znyv5PYDYt+fUtuv\nqKgIRUVF+P73vw8AaGpqQn9//4x95LRd0gn9/v37GB8fBwBMTU3h6tWr0Ov1M/apr6/H+fPnAQC3\nb9/GqlWroNVqk7102sVzb06nM/RTtK+vD6IoRnzPno1OnTqF4eFh3L17F5cuXcJzzz0XaieJUtsO\niO/+lNx+k5OT+PrrrwFMz7X+7rvvoqKiYsY+Sm6/eO5Pqe1XWFiI4uJiDA4OAgCuXbuGDRs2zNhH\nTtslvcDF2NgYzGYzgsEggsEg9u7di+rqapw7dw4AcOjQIZhMJnR1deGJJ57AY489hjfeeCPZyy6I\neO7tL3/5C/74xz9iyZIlWLZsGS5dupThqOWT/junhraLJNL9Kbn9nE4nGhsbAUz/F37Pnj2ora1V\nTfvFc39Kbr+zZ89iz5498Hq9+M53voM//elPSbed7LlciIgouyhuxSIiIoqMCZ2ISCWY0ImIVIIJ\nnYhIJZjQiYhUggmdiEgl/j9qBdnj41BqNwAAAABJRU5ErkJggg==\n" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "maxs = array([amax(randn(10000)) for i in range(1000)])\n", "max2s = array([amax(1.1*randn(10000)) for i in range(1000)])\n", "_=hist(maxs,bins=100,alpha=0.7)\n", "_=hist(max2s,bins=100,alpha=0.7)\n", "sum(max2s>maxs)*1.0/len(maxs)" ] }, { "cell_type": "markdown", "id": "7de893f1", "metadata": {}, "source": [ "# Conditional Distributions" ] }, { "cell_type": "markdown", "id": "641f47be", "metadata": {}, "source": [ "Let's look at another important example from statistics. Consider a population of size $N = 1000000$:\n", "\n" ] }, { "cell_type": "code", "execution_count": 105, "id": "072924e0", "metadata": { "collapsed": false }, "outputs": [], "source": [ "N = 1000000" ] }, { "cell_type": "markdown", "id": "994cab7b", "metadata": {}, "source": [ "Now assume that about 1:1000 of the population are infected with some disease.\n", "\n" ] }, { "cell_type": "code", "execution_count": 106, "id": "86e52b79", "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1045" ] }, "execution_count": 106, "metadata": {}, "output_type": "execute_result" } ], "source": [ "infected = (rand(N)<1e-3)\n", "sum(infected)" ] }, { "cell_type": "markdown", "id": "a2d3b004", "metadata": {}, "source": [ "Also, assume we have a test that is 99 percent accurate. We give the test to the entire population and measure the results.\n", "\n" ] }, { "cell_type": "code", "execution_count": 107, "id": "79816e05", "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "11120" ] }, "execution_count": 107, "metadata": {}, "output_type": "execute_result" } ], "source": [ "testresult = where(rand(N)<1e-2,1-infected,infected)\n", "sum(testresult)" ] }, { "cell_type": "markdown", "id": "03291ca7", "metadata": {}, "source": [ "Now, let's look at the population of people whose test result is positive.\n", "What fraction of those people are infected?\n", "\n" ] }, { "cell_type": "code", "execution_count": 108, "id": "371eb20b", "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.092985611510791363" ] }, "execution_count": 108, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum(infected[testresult==True])*1.0/sum(testresult==True)" ] }, { "cell_type": "markdown", "id": "c201f550", "metadata": {}, "source": [ "(Implications)\n", "\n", "That's not good: less than 10 percent of the population that tests positive is actually infected. Ninety percent of those who test positive are _false positives_ even though the test is 99 percent accurate.\n", "\n", "This is a common issue in medicine, anti-terrorism, and many other problems:\n", "\n", "- testing everybody for HIV infection\n", "- testing everybody for prostate cancer / breast cancer\n", "- widespread Internet surveillance for \"terrorism\"\n", "\n", "All of these end up generating enormous numbers of false positives." ] }, { "cell_type": "markdown", "id": "4b372949", "metadata": {}, "source": [ "(Bayesian Calculation)\n", "\n", "Mathematically, we write this in terms of _conditional probabilities_:\n", "\n", "prior probability: \n", "\n", "$$P(i)$$\n", "\n", "probability of test: \n", "\n", "$$P(t=0|i=0) = 0.99$$\n", "$$P(t=1|i=1) = 0.99$$\n", "\n", "Bayes formula:\n", "\n", "$$ P(i=1|t=1) = \\frac{P(t=1|i=1) P(i=1)}{P(t=1)} = \\frac{P(t=1|i=1) P(i=1)}{P(t=1|i=0) P(i=0) + P(t=1|i=1) P(i=1) } $$\n", "\n", "$$ ~~~~~~ = \\frac{0.99 \\cdot 0.001}{ 0.01 \\cdot 0.999 + 0.99 \\cdot 0.001 } \\approx \\frac{0.001}{0.01} = 0.1 $$\n", "\n" ] }, { "cell_type": "markdown", "id": "75b760ee", "metadata": {}, "source": [ "(Derivation of Bayes Formula)\n", "\n", "Joint density ($c$ and $x$ could be any variables):\n", "\n", "$$ p(c,x) $$\n", "\n", "Conditional density:\n", "\n", "$$ p(c,x) = P(c|x) ~ p(x) $$\n", "\n", "Bayes formula:\n", "\n", "$$ P(c|x) ~ p(x) ~ = ~ p(x|c) ~ P(c) $$\n", "\n", "$$ P(c|x) = \\frac{p(x|c) P(c)}{p(x)} $$" ] }, { "cell_type": "code", "execution_count": null, "id": "275dfb26", "metadata": { "collapsed": false }, "outputs": [], "source": [] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 5 }