{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Dollar-unit sampling and \"taint\"\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Dollar-unit sampling\n", "\n", "Dollar-unit sampling (DUS) is a form of probability proportional to size (PPS) sampling that originated in statistical financial auditing. In this case, \"size\" is an upper bound on the value of the item; DUS is sometimes called sampling with probability proportional to an upper bound (PPUB).\n", "\n", "The idea is that each item in a population has a reported value and a true value; the true value is known only if the item is audited—which is expensive. There is a possibility that the reported value is too high. We sample units with probability proportional to their reported values.\n", "\n", "Intuitively, it makes sense to focus attention on the items with larger reported values, because \"that's where the money is.\" Errors in the values of those items potentially have a larger effect on the error in the total reported value for the population.\n", "\n", "### The math\n", "\n", "There are $N$ items, $\\{x_j\\}_{j=1}^N$. There are $N$ known numbers $\\{u_j\\}_{j=1}^N$ such that\n", "$$\n", " 0 \\le x_j \\le u_j, \\; j = 1, \\ldots, N.\n", "$$\n", "Define $u \\equiv \\sum_{j=1}^N u_j$, $x \\equiv \\sum_{j=1}^N x_j$, and $\\mu \\equiv x/N = \\frac{1}{N} \\sum_{j=1}^N x_j.$\n", "\n", "We will consider only sampling with replacement (although this can be extended to sampling without replacement).\n", "\n", "In the $i$th draw, we select item $j$ with probability $\\pi_j \\equiv u_j/u$. \n", "\n", "Let $X_i$ be the true value of the item selected in draw $i$, and let $U_i$ be the reported value of the item selected in draw $i$.\n", "\n", "If we sample with replacement, $\\{X_i\\}$ are iid, as are $\\{U_i \\}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Taint\n", "\n", "Define $t_j \\equiv x_j/u_j$. This value, the \"taint,\" is the fraction of the reported value that the true value represents.\n", "\n", "Define $t \\equiv x/u$. (This is _not_ parallel to the definitions of $x$ and $u$.)\n", "\n", "Define $T_i \\equiv X_i/U_i$, the taint of the item selected in the $i$th draw. Then $\\{T_i \\}$ are iid.\n", "\n", "Calculate:\n", "$$\n", " \\mathbb E T_i = \\sum_{j=1}^N t_j \\pi_j = \\sum_{j=1}^N (x_j/u_j) (u_j/u) = \\frac{1}{u}\\sum_{j=1}^N x_j = x/u = t.\n", "$$\n", "Thus the expected value of the taint for any draw, times $u$, is the population total $x$.\n", "\n", "Since $u$ is known, we can translate an upper confidence bound for $\\mathbb E T_i$ into an upper confidence bound for $x$ or for $\\mu$ by multiplication.\n", "\n", "This is the relationship we will exploit to connect auditing problems to the problem of finding a nonparametric confidence bound for the mean of a non-negative or bounded population.\n", "\n", "The main point is that since $\\mathbb P \\{T_i \\in [0, 1]\\} = 1$, we can use any of the methods we've developed for finding confidence bounds for the mean of nonnegative or bounded populations (Binomial with thresholding, Checychev's inequality, Hoeffding's inequality, Markov's inequality, the MDKW inequality, and the Kaplan-Markov method) to make confidence bonds for $\\mathbb E T_i$. In turn, any confidence bound for $t$ can be translated easily into a rigorous confidence bound for the population mean $\\mu$ and population total $x$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Some examples\n", "\n", "### Litigation and elections\n", "\n", "A plaintiff might want to find a confidence bound in one direction; a defendant might want a confidence bound in the opposite direction.\n", "\n", "+ Employment litigation\n", " - wage and hour: overtime is nonnegative, and has an upper bound based on the number of days an employee worked\n", " - missed meal or rest breaks: nonnegative and has an upper bound based on days worked\n", "+ Toxic tort class action\n", " - damages are nonnegative\n", "+ Construction defects, product defects\n", " - damages nonnegative and bounded (e.g., by replacement cost)\n", "+ Illegal charges for loan origination\n", " - nonnegative and bounded by total origination fee\n", "+ Patent and intellectual property infringement\n", " - damages nonnegative; upper bound from royalty rate and number of potentially infringing items\n", "+ Healthcare fraud\n", " - damages nonnegative and bounded by billed amount\n", "+ Tax fraud\n", " - taxes paid is an upper bound\n", "+ Under-refunding of security deposits\n", " - amount not refunded is an upper bound\n", "+ Election integrity\n", " - upper and lower bounds on the error in counting the votes on a ballot depend on the election rules; I'll talk about this in my plenary lecture" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's code to draw a weighted random sample in Python:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# This is the first cell with code: set up the Python environment\n", "%matplotlib inline\n", "from __future__ import print_function, division\n", "import matplotlib.pyplot as plt\n", "import math\n", "import numpy as np\n", "import scipy as sp\n", "import scipy.stats\n", "from scipy.stats import binom\n", "import pandas as pd\n", "from ipywidgets import interact, interactive, fixed\n", "import ipywidgets as widgets\n", "from IPython.display import clear_output, display, HTML" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def weightedRandomSample(n, weights):\n", " '''\n", " Weighted random sample of size n drawn with replacement.\n", " Returns indices of the selected items, and the raw uniform values used to \n", " select them.\n", " '''\n", " if any(weights < 0):\n", " print('negative weight in weightedRandomSample')\n", " return float('NaN')\n", " else:\n", " wc = np.cumsum(weights, dtype=float)/np.sum(weights, dtype=float) # ensure weights sum to 1\n", " theSam = np.random.random_sample((n,))\n", " return np.array(wc).searchsorted(theSam), theSam" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 0 1 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5\n", " 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8\n", " 8 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9 9 9 9 9 9 9 9 9]\n" ] }, { "data": { "text/plain": [ "(array([ 2., 5., 4., 7., 7., 14., 14., 7., 26.]),\n", " array([-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5]),\n", " )" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAC4hJREFUeJzt3W2IZQUdx/HfL9coH6iVnZbNh8ZAjCVwNwaxjLDWwjRS\n34RCsoSwvjBbQ4jNN/lyAx/qRQiru7mQGeEDSi7WtgkihDSri+6Dotiqu63uiJTWG1v99WKOOSMz\nc+/ce2fO3f98PzDMveee2fPnMPPl7Jl7zjiJAADHv4+1PQAAYDAIOgAUQdABoAiCDgBFEHQAKIKg\nA0ARBB0AiiDoAFAEQQeAIpYt5sZWrFiR0dHRxdwkABz3du/e/WaSkU7rLWrQR0dHNT4+vpibBIDj\nnu1XulmPUy4AUARBB4AiCDoAFEHQAaAIgg4ARRB0ACiCoANAEQQdAIog6ABQxKJeKQoAi21006Nt\njyBJOrj5sgXfBkfoAFAEQQeAIgg6ABTRMei2z7T9uO39tvfZ3tgsv8X2Ydt7mo9LF35cAMBsuvml\n6DFJNyV52vapknbb3tm8dkeSWxduPABAtzoGPckRSUeax+/YPiDp9IUeDAAwP/M6h257VNJaSU81\ni26w/aztbbaXD3g2AMA8dB1026dIekDSjUnelnSnpM9LWqPJI/jbZvm6DbbHbY9PTEwMYGQAwEy6\nCrrtEzUZ83uTPChJSd5I8l6S9yXdJen8mb42yZYkY0nGRkY6/kk8AECPunmXiyVtlXQgye1Tlq+a\nstqVkvYOfjwAQLe6eZfLhZKukfSc7T3NspslXW17jaRIOijpugWZEADQlW7e5fKkJM/w0o7BjwMA\n6BVXigJAEQQdAIog6ABQBEEHgCIIOgAUQdABoAiCDgBFEHQAKIKgA0ARBB0AiiDoAFAEQQeAIgg6\nABRB0AGgCIIOAEUQdAAogqADQBEEHQCKIOgAUARBB4AiCDoAFEHQAaAIgg4ARRB0ACiCoANAEQQd\nAIog6ABQBEEHgCIIOgAUQdABoAiCDgBFdAy67TNtP257v+19tjc2y0+zvdP2i83n5Qs/LgBgNt0c\noR+TdFOS1ZIukHS97dWSNknaleQcSbua5wCAlnQMepIjSZ5uHr8j6YCk0yVdLml7s9p2SVcs1JAA\ngM7mdQ7d9qiktZKekrQyyZHmpdclrRzoZACAeek66LZPkfSApBuTvD31tSSRlFm+boPtcdvjExMT\nfQ0LAJhdV0G3faImY35vkgebxW/YXtW8vkrS0Zm+NsmWJGNJxkZGRgYxMwBgBt28y8WStko6kOT2\nKS89Iml983i9pIcHPx4AoFvLuljnQknXSHrO9p5m2c2SNkv6ve1rJb0i6XsLMyIAoBsdg57kSUme\n5eV1gx0HANArrhQFgCIIOgAUQdABoAiCDgBFEHQAKIKgA0ARBB0AiiDoAFAEQQeAIgg6ABRB0AGg\nCIIOAEUQdAAogqADQBEEHQCKIOgAUARBB4AiCDoAFEHQAaAIgg4ARRB0ACiCoANAEQQdAIog6ABQ\nBEEHgCIIOgAUQdABoAiCDgBFEHQAKIKgA0ARBB0AiugYdNvbbB+1vXfKsltsH7a9p/m4dGHHBAB0\n0s0R+j2SLplh+R1J1jQfOwY7FgBgvjoGPckTkt5ahFkAAH3o5xz6DbafbU7JLJ9tJdsbbI/bHp+Y\nmOhjcwCAufQa9DslfV7SGklHJN0224pJtiQZSzI2MjLS4+YAAJ30FPQkbyR5L8n7ku6SdP5gxwIA\nzFdPQbe9asrTKyXtnW1dAMDiWNZpBdv3SbpI0grbhyT9TNJFttdIiqSDkq5bwBkBAF3oGPQkV8+w\neOsCzAIA6ANXigJAEQQdAIroeMoFQHdGNz3a9ghD5eDmy9oeYcnhCB0AiiDoAFAEQQeAIgg6ABRB\n0AGgCIIOAEUQdAAogqADQBEEHQCKIOgAUARBB4AiCDoAFEHQAaAIgg4ARRB0ACiCoANAEQQdAIog\n6ABQBEEHgCIIOgAUQdABoAiCDgBFEHQAKIKgA0ARBB0AiiDoAFAEQQeAIjoG3fY220dt752y7DTb\nO22/2HxevrBjAgA66eYI/R5Jl3xk2SZJu5KcI2lX8xwA0KKOQU/yhKS3PrL4cknbm8fbJV0x4LkA\nAPPU6zn0lUmONI9fl7RyQPMAAHrU9y9Fk0RSZnvd9gbb47bHJyYm+t0cAGAWvQb9DdurJKn5fHS2\nFZNsSTKWZGxkZKTHzQEAOuk16I9IWt88Xi/p4cGMAwDoVTdvW7xP0l8lnWv7kO1rJW2W9E3bL0q6\nuHkOAGjRsk4rJLl6lpfWDXgWAEAfuFIUAIog6ABQBEEHgCIIOgAUQdABoAiCDgBFEHQAKIKgA0AR\nBB0AiiDoAFAEQQeAIgg6ABRB0AGgCIIOAEUQdAAogqADQBEEHQCKIOgAUARBB4AiOv5NUWAuo5se\nbXsEDCm+NxYfR+gAUARBB4AiCDoAFEHQAaAIgg4ARRB0ACiCoANAEQQdAIog6ABQBEEHgCIIOgAU\n0de9XGwflPSOpPckHUsyNoihAADzN4ibc309yZsD+HcAAH3glAsAFNFv0CPpz7Z3294wiIEAAL3p\n95TLV5Mctv0ZSTttP5/kiakrNKHfIElnnXVWn5tr37Dc4/ng5svaHgHAkOnrCD3J4ebzUUkPSTp/\nhnW2JBlLMjYyMtLP5gAAc+g56LZPtn3qB48lfUvS3kENBgCYn35OuayU9JDtD/6d3yZ5bCBTAQDm\nreegJ3lZ0nkDnAUA0AfetggARRB0ACiCoANAEYO49B8tGJb3wwMYHhyhA0ARBB0AiiDoAFAEQQeA\nIgg6ABRB0AGgCIIOAEUQdAAogqADQBEEHQCKIOgAUARBB4AiCDoAFEHQAaAIgg4ARRw390Pn/t8A\nMDeO0AGgCIIOAEUQdAAogqADQBEEHQCKIOgAUARBB4AiCDoAFEHQAaAIgg4ARRB0ACiir6DbvsT2\nC7Zfsr1pUEMBAOav56DbPkHSryR9W9JqSVfbXj2owQAA89PPEfr5kl5K8nKSdyX9TtLlgxkLADBf\n/QT9dEmvTXl+qFkGAGjBgt8P3fYGSRuap/+2/cJCb3MOKyS92eL2hw3740Psi+nYH9P1vT/88762\n/7luVuon6IclnTnl+RnNsmmSbJG0pY/tDIzt8SRjbc8xLNgfH2JfTMf+mO542R/9nHL5m6RzbJ9t\n++OSrpL0yGDGAgDMV89H6EmO2f6hpD9KOkHStiT7BjYZAGBe+jqHnmSHpB0DmmUxDMWpnyHC/vgQ\n+2I69sd0x8X+cJK2ZwAADACX/gNAEUsm6NymYJLtM20/bnu/7X22N7Y90zCwfYLtZ2z/oe1Z2mb7\n07bvt/287QO2v9z2TG2x/ePm52Sv7ftsf6LtmeayJILObQqmOSbppiSrJV0g6folvC+m2ijpQNtD\nDIlfSnosyRcknaclul9sny7pR5LGknxRk2/+uKrdqea2JIIublPwf0mOJHm6efyOJn9Yl/QVvrbP\nkHSZpLvbnqVttj8l6WuStkpSkneT/LPdqVq1TNInbS+TdJKkf7Q8z5yWStC5TcEMbI9KWivpqXYn\nad0vJP1E0vttDzIEzpY0IenXzSmou22f3PZQbUhyWNKtkl6VdETSv5L8qd2p5rZUgo6PsH2KpAck\n3Zjk7bbnaYvt70g6mmR327MMiWWSviTpziRrJf1H0pL8nZPt5Zr8n/zZkj4r6WTb3293qrktlaB3\ndZuCpcL2iZqM+b1JHmx7npZdKOm7tg9q8lTcN2z/pt2RWnVI0qEkH/yv7X5NBn4puljS35NMJPmv\npAclfaXlmea0VILObQoatq3J86MHktze9jxtS/LTJGckGdXk98Vfkgz1UdhCSvK6pNdsn9ssWidp\nf4sjtelVSRfYPqn5uVmnIf8F8YLfbXEYcJuCaS6UdI2k52zvaZbd3Fz1C0jSDZLubQ5+Xpb0g5bn\naUWSp2zfL+lpTb477BkN+RWjXCkKAEUslVMuAFAeQQeAIgg6ABRB0AGgCIIOAEUQdAAogqADQBEE\nHQCK+B8jcHsC8OLlsAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# illustrate the random sampling code\n", "vals = 10\n", "n = 100\n", "w = np.arange(vals)+1 # linearly increasing weights\n", "wrs, raw = weightedRandomSample(n, w)\n", "print(np.sort(wrs))\n", "fig, ax = plt.subplots(nrows=1, ncols=1)\n", "bins = np.arange(np.min(wrs)-0.5, np.max(wrs)+0.5, 1)\n", "ax.hist(wrs, bins=bins)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What's next?\n", "- [Next: Penny Sampling](pennySampling.ipynb)\n", "- [Previous: The Kaplan-Wald Confidence Bound for a Nonnegative Mean](kaplanWald.ipynb)\n", "- [Index](index.ipynb)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "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.1" } }, "nbformat": 4, "nbformat_minor": 1 }