updates based on feedback

Author: sastoudt

README.rst

@@ -55,16 +55,21 @@ Accepts an open-ended set of paired lists, each expected to provide an estimate
 Approximating means
 ~~~~~~~~~~~~~~~~~~~
 
-Estimate a mean and approximate the margin of error. The Census Bureau guidelines do not provide instructions for
-approximating a mean using data from the ACS. They do provide guidance for approximating a mean with data `from the PUMS <https://www2.census.gov/programs-surveys/acs/tech_docs/pums/accuracy/2013_2017AccuracyPUMS.pdf?#>`_.
-Instead, we implement a simulation based approach. First the number of units in each bin is simulated (assuming a normal distribution around the
-estimate). Then for each unit a value within the bin is simulated (assuming a uniform distribution within each bin).
-Note that for quantities, such as income, the Pareto distribution is `often used <https://www2.census.gov/ces/wp/2014/CES-WP-14-21.pdf>`_ instead of the uniform distribution in the upper-most bin. We provide an option to do this (`usePareto = True`).
-We can then calculate the mean directly from the simulated data. We repeat the simulation many times to calculate an empirical margin of error.
+Estimate a mean and approximate the margin of error. 
+    
+The Census Bureau guidelines do not provide instructions for approximating a mean using data from the ACS.
+They do provide guidance for approximating a mean with data `from the PUMS <https://www2.census.gov/programs-surveys/acs/tech_docs/pums/accuracy/2013_2017AccuracyPUMS.pdf?#>`_.
+Instead, we implement a simulation based approach. 
+    
+First the number of units in each bin is simulated (assuming a normal distribution around the estimate). Then for each unit a value within
+the bin is simulated (assuming a uniform distribution within each bin). Note that for quantities, such as income, the Pareto distribution is
+`often used <https://www2.census.gov/ces/wp/2014/CES-WP-14-21.pdf>`_ instead of the uniform distribution in the upper-most bin. We provide an
+option to do this (`pareto = True`). We can then calculate the mean directly from the simulated data. We repeat the simulation many times to
+calculate an empirical margin of error.
+  
 Note that this function assumes you have a lower bound for the smallest bin and an upper bound for the largest bin. We recommend trying different
 lower and upper bounds to assess the sensitivity of the resulting mean to your assumptions.
 
-
 Expects a list of dictionaries that divide the full range of data values into continuous categories. Each dictionary should have four keys:
 
 .. list-table::
@@ -110,7 +115,7 @@ Increasing the number of runs should increase your precision (with diminishing r
  >>> approximate_mean(income)
  (98045.44530685373, 194.54892406267754)
  
- >>> approximate_mean(income, usePareto=True)
+ >>> approximate_mean(income, pareto=True)
  (60364.96525340687, 58.60735554621351)
 
 

census_data_aggregator/__init__.py

@@ -408,12 +408,18 @@ def approximate_percentchange(pair_old, pair_new):
     return percent_change_estimate, percent_change_moe
 
 
-def approximate_mean(range_list, number_replicates=50, usePareto=False):
+def approximate_mean(range_list, simulations=50, pareto=False):
     """
-    Estimate a mean and approximate the margin of error. The Census Bureau guidelines do not provide instructions for
+    Estimate a mean and approximate the margin of error.
+
+    The Census Bureau guidelines do not provide instructions for
     approximating a mean using data from the ACS. They do provide guidance for approximating a mean with data `from the PUMS`_.
-    Instead, we implement a simulation based approach. Due to the stochastic nature of the simulation approach, you will need to set
-    a seed before running this function to ensure replicability. Note that this function assumes you have a lower bound for the smallest
+    Instead, we implement a simulation based approach.
+
+    Due to the stochastic nature of the simulation approach, you will need to set
+    a seed before running this function to ensure replicability.
+
+    Note that this function assumes you have a lower bound for the smallest
     bin and an upper bound for the largest bin. We recommend trying different lower and upper bounds to assess the sensitivity of the
     resulting mean to your assumptions.
 
@@ -424,8 +430,8 @@ def approximate_mean(range_list, number_replicates=50, usePareto=False):
                 * max (int): The maximum value of the range
                 * n (int): The number of people, households or other unit in the range
                 * moe (float): The margin of error for n
-        number_replicates (int): number of replicates for simulation, used to estimate margin of error
-        usePareto (logical): use the Pareto distribution to simulate from in upper bin, otherwise use uniform, Pareto is appropriate for income
+        simulations (int): number of simulations to run, used to estimate margin of error. Defaults to 50.
+        pareto (logical): use the Pareto distribution to simulate from in upper bin, otherwise use uniform. Pareto is appropriate for income. Defaults to False.
 
     Returns:
         A two-item tuple with the mean followed by the approximated margin of error.
@@ -455,7 +461,7 @@ def approximate_mean(range_list, number_replicates=50, usePareto=False):
         ]
         >>> approximate_mean(income)
         (98045.44530685373, 194.54892406267754)
-        >>> approximate_mean(income, usePareto=True)
+        >>> approximate_mean(income, pareto=True)
         (60364.96525340687, 58.60735554621351)
 
     ... _from the PUMS:
@@ -465,34 +471,35 @@ def approximate_mean(range_list, number_replicates=50, usePareto=False):
     # Sort the list
     range_list.sort(key=lambda x: x['min'])
 
-    idx = len(range_list)
-
-    nb1 = range_list[idx - 2]['n']  # number in second to last bin
-    nb = range_list[idx - 1]['n']  # number in last bin
-    lb1 = range_list[idx - 2]['min']  # lower bound of second to last bin
-    lb = range_list[idx - 1]['min']  # lower bound of last bin
+    if pareto:  # need shape parameter if using Pareto distribution
 
-    alpha_hat = (numpy.log(nb1 + nb) - numpy.log(nb)) / (numpy.log(lb) - numpy.log(lb1))
+        nb1 = range_list[-2]['n']  # number in second to last bin
+        nb = range_list[-1]['n']  # number in last bin
+        lb1 = range_list[-2]['min']  # lower bound of second to last bin
+        lb = range_list[-1]['min']  # lower bound of last bin
+        alpha_hat = (numpy.log(nb1 + nb) - numpy.log(nb)) / (numpy.log(lb) - numpy.log(lb1))  # shape parameter for Pareto
 
     simulation_results = []
-    for i in range(number_replicates):
+    for i in range(simulations):
         simulated_values = []
         simulated_n = []
+        # loop through every bin except the last one
         for range_ in range_list[:-1]:
             se = range_['moe'] / 1.645  # convert moe to se
             nn = round(numpy.random.normal(range_['n'], se))  # use moe to introduce randomness into number in bin
             nn = int(nn)  # clean it up
             simulated_values.append(numpy.random.uniform(range_['min'], range_['max'], size=(1, nn)).sum())  # draw random values within the bin, assume uniform
             simulated_n.append(nn)
-        if usePareto:
-            last = range_list[idx - 1]
+        # a special case to handle the last bin
+        if pareto:
+            last = range_list[-1]
             se = last['moe'] / 1.645  # convert moe to se
             nn = round(numpy.random.normal(last['n'], se))  # use moe to introduce randomness into number in bin
             nn = int(nn)  # clean it up
             simulated_values.append(numpy.random.pareto(a=alpha_hat, size=(1, nn)).sum())  # draw random values within the bin, assume uniform
             simulated_n.append(nn)
         else:
-            last = range_list[idx - 1]
+            last = range_list[-1]
             se = last['moe'] / 1.645  # convert moe to se
             nn = round(numpy.random.normal(last['n'], se))  # use moe to introduce randomness into number in bin
             nn = int(nn)  # clean it up

test.py

@@ -248,7 +248,40 @@ def test_mean(self):
 
         numpy.random.seed(711355)
 
-        mean, moe = census_data_aggregator.approximate_mean(range_list, usePareto=True)
+        mean, moe = census_data_aggregator.approximate_mean(range_list, pareto=True)
+
+        self.assertAlmostEqual(mean, 60364.96525340687, places=3)
+        self.assertAlmostEqual(moe, 58.60735554621351, places=3)
+
+    def test_mean_order(self):
+        range_list = [
+            dict(min=50000, max=59999, n=9181800, moe=20965),
+            dict(min=60000, max=74999, n=11818514, moe=30723),
+            dict(min=75000, max=99999, n=14636046, moe=49159),
+            dict(min=100000, max=124999, n=10273788, moe=47842),
+            dict(min=125000, max=149999, n=6428069, moe=37952),
+            dict(min=150000, max=199999, n=6931136, moe=37236),
+            dict(min=200000, max=1000000, n=7465517, moe=42206),
+            dict(min=0, max=9999, n=7942251, moe=17662),
+            dict(min=10000, max=14999, n=5768114, moe=16409),
+            dict(min=15000, max=19999, n=5727180, moe=16801),
+            dict(min=20000, max=24999, n=5910725, moe=17864),
+            dict(min=25000, max=29999, n=5619002, moe=16113),
+            dict(min=30000, max=34999, n=5711286, moe=15891),
+            dict(min=35000, max=39999, n=5332778, moe=16488),
+            dict(min=40000, max=44999, n=5354520, moe=15415),
+            dict(min=45000, max=49999, n=4725195, moe=16890)
+        ]
+        numpy.random.seed(711355)
+        # Calculate the mean and its MOE
+        mean, moe = census_data_aggregator.approximate_mean(range_list)
+
+        self.assertAlmostEqual(mean, 98045.44530685373, places=3)
+        self.assertAlmostEqual(moe, 194.54892406267754, places=3)
+        
+        numpy.random.seed(711355)
+
+        mean, moe = census_data_aggregator.approximate_mean(range_list, pareto=True)
 
         self.assertAlmostEqual(mean, 60364.96525340687, places=3)
         self.assertAlmostEqual(moe, 58.60735554621351, places=3)