commit-2018/11/14@15:49CST

analysis.py

@@ -1,19 +1,119 @@
+#Titan Robotics Team 2022: Data Analysis Module
+#Written by Arthur Lu & Jacob Levine
+#Notes:
+#   this should be imported as a python module using 'import analysis'
+#   this should be included in the local directory or environment variable
+#   this module has not been optimized for multhreaded computing
+#Number of easter eggs: 2
 
-#this should be imported as a python module using 'import analysis'
+#setup:
 
+__all__ = [
+    '_init_device',
+    'c_entities',
+    'nc_entities',
+    'obstacles',
+    'objectives',
+    'load_csv',
+    'basic_stats',
+    'z_score',
+    'stdev_z_split',
+    'histo_analysis', #histo_analysis_old is intentionally left out
+    'poly_regression',
+    'r_squared',
+    'rms',
+    'basic_analysis',
+    ]
+
+#now back to your regularly scheduled programming:
+
+
+import warnings
 import statistics
 import math
 import csv
 import functools
+import numpy as np
+import time
+import torch
+import scipy
+import matplotlib
+from sklearn import *
+
+def _init_device (setting, arg): #initiates computation device for ANNs
+    if setting == "cuda":
+        temp = setting + ":" + arg
+        the_device_woman = torch.device(temp if torch.cuda.is_available() else "cpu")
+        return the_device_woman #name that reference
+    elif setting == "cpu":
+        the_device_woman = torch.device("cpu")
+        return the_device_woman #name that reference
+    else:
+        return "error:specified device does not exist"
 
 class c_entities:
 
     c_names = []
     c_ids = []
     c_pos = []
-    c_porperties = []
+    c_properties = []
     c_logic = []
 
+    def debug(self):
+        print("c_entities has attributes names, ids, positions, properties, and logic. __init__ takes self, 1d array of names, 1d array of ids, 2d array of positions, nd array of properties, and nd array of logic")
+        return[self.c_names, self.c_ids, self.c_pos, self.c_properties, self.c_logic]
+    
+    def __init__(self, names, ids, pos, properties, logic):
+        self.c_names = names
+        self.c_ids = ids
+        self.c_pos = pos
+        self.c_properties = properties
+        self.c_logic = logic
+        return None
+        
+
+    def append(self, n_name, n_id, n_pos, n_property, n_logic):
+        self.c_names.append(n_name)
+        self.c_ids.append(n_id)
+        self.c_pos.append(n_pos)
+        self.c_properties.append(n_property)
+        self.c_logic.append(n_logic)
+        return None
+    
+    def edit(self, search, n_name, n_id, n_pos, n_property, n_logic):
+        position = 0
+        for i in range(0, len(self.c_ids), 1):
+            if self.c_ids[i] == search:
+                position = i
+        if n_name != "null":
+            self.c_names[position] = n_name
+    
+        if n_id != "null":
+            self.c_ids[position] = n_id
+    
+        if n_pos != "null":
+            self.c_pos[position] = n_pos
+    
+        if n_property != "null":
+            self.c_properties[position] = n_property
+    
+        if n_logic != "null":
+            self.c_logic[position] = n_logic
+    
+        return None
+    
+    def search(self, search):
+        position = 0
+        for i in range(0, len(self.c_ids), 1):
+            if self.c_ids[i] == search:
+                position = i
+
+        return [self.c_names[position], self.c_ids[position], self.c_pos[position], self.c_properties[position], self.c_logic[position]]    
+
+    def regurgitate(self):
+
+        return[self.c_names, self.c_ids, self.c_pos, self.c_properties, self.c_logic]
+    
 class nc_entities:
 
     c_names = []
@@ -23,7 +123,7 @@ class nc_entities:
     c_effects = []
 
     def debug(self):
-        print ("nc_entities (non-controlable entities) has attributes names, ids, positions, properties, and effects. __init__ takes self, 1d array of names, 1d array of ids, 2d array of psoitions, 2d array of properties, and 2d array of effects.")
+        print ("nc_entities (non-controlable entities) has attributes names, ids, positions, properties, and effects. __init__ takes self, 1d array of names, 1d array of ids, 2d array of positions, 2d array of properties, and 2d array of effects.")
         return[self.c_names, self.c_ids, self.c_pos, self.c_properties, self.c_effects]
 
     def __init__(self, names, ids, pos, properties, effects):
@@ -40,6 +140,8 @@ class nc_entities:
         self.c_pos.append(n_pos)
         self.c_properties.append(n_property)
         self.c_effects.append(n_effect)
+        
+        return None
 
     def edit(self, search, n_name, n_id, n_pos, n_property, n_effect):
         position = 0
@@ -71,6 +173,10 @@ class nc_entities:
 
         return [self.c_names[position], self.c_ids[position], self.c_pos[position], self.c_properties[position], self.c_effects[position]]        
 
+    def regurgitate(self):
+
+        return[self.c_names, self.c_ids, self.c_pos, self.c_properties, self.c_effects]
+
 class obstacles:
 
     c_names = []
@@ -124,6 +230,10 @@ class obstacles:
 
         return [self.c_names[position], self.c_ids[position], self.c_perim[position], self.c_effects[position]]
 
+    def regurgitate(self):
+
+        return[self.c_names, self.c_ids, self.c_perim, self.c_effects]
+
 class objectives:
     
     c_names = []
@@ -178,6 +288,10 @@ class objectives:
 
         return [self.c_names[position], self.c_ids[position], self.c_pos[position], self.c_effects[position]]
 
+    def regurgitate(self):
+
+        return[self.c_names, self.c_ids, self.c_pos, self.c_effects]
+    
 def load_csv(filepath):
     with open(filepath, newline = '') as csvfile:
         file_array = list(csv.reader(csvfile))
@@ -203,8 +317,17 @@ def basic_stats(data, mode, arg): # data=array, mode = ['1d':1d_basic_stats, 'co
             mode = statistics.mode(data_t)
         except:
             mode = None
-        stdev = statistics.stdev(data_t)
+        try:
-        variance = statistics.variance(data_t)
+            stdev = statistics.stdev(data)
+            
+        except:
+            
+            stdev = None
+        
+        try:
+            variance = statistics.variance(data_t)
+        except:
+            variance = None
         
         out = [mean, median, mode, stdev, variance]
 
@@ -216,7 +339,10 @@ def basic_stats(data, mode, arg): # data=array, mode = ['1d':1d_basic_stats, 'co
         c_data_sorted = []
         
         for i in data:
-            c_data.append(float(i[arg]))
+            try:
+                c_data.append(float(i[arg]))
+            except:
+                pass
             
         mean = statistics.mean(c_data)
         median = statistics.median(c_data)
@@ -224,8 +350,14 @@ def basic_stats(data, mode, arg): # data=array, mode = ['1d':1d_basic_stats, 'co
             mode = statistics.mode(c_data)
         except:
             mode = None
-        stdev = statistics.stdev(c_data)
+        try:
-        variance = statistics.variance(c_data)
+            stdev = statistics.stdev(c_data)
+        except:
+            stdev = None
+        try:
+            variance = statistics.variance(c_data)
+        except:
+            variance = None
         
         out = [mean, median, mode, stdev, variance]
 
@@ -244,8 +376,14 @@ def basic_stats(data, mode, arg): # data=array, mode = ['1d':1d_basic_stats, 'co
             mode = statistics.mode(r_data)
         except:
             mode = None
-        stdev = statistics.stdev(r_data)
+        try:
-        variance = statistics.variance(r_data)
+            stdev = statistics.stdev(r_data)
+        except:
+            stdev = None
+        try:
+            variance = statistics.variance(r_data)
+        except:
+            variance = None
         
         out = [mean, median, mode, stdev, variance]
 
@@ -253,11 +391,11 @@ def basic_stats(data, mode, arg): # data=array, mode = ['1d':1d_basic_stats, 'co
     else:
         return ["mode_error", "mode_error"]
     
-def z_score(point, mean, stdev):
+def z_score(point, mean, stdev): #returns z score with inputs of point, mean and standard deviation of spread
     score = (point - mean)/stdev
     return score
 
-def stdev_z_split(mean, stdev, delta, low_bound, high_bound):
+def stdev_z_split(mean, stdev, delta, low_bound, high_bound): #returns n-th percentile of spread given mean, standard deviation, lower z-score, and upper z-score
 
     z_split = []
 
@@ -275,7 +413,7 @@ def stdev_z_split(mean, stdev, delta, low_bound, high_bound):
 
     return z_split
 
-def histo_analysis(hist_data): #note: depreciated
+def histo_analysis_old(hist_data): #note: depreciated
 
     if hist_data == 'debug':
         return['lower estimate (5%)', 'lower middle estimate (25%)', 'middle estimate (50%)', 'higher middle estimate (75%)', 'high estimate (95%)', 'standard deviation', 'note: this has been depreciated']
@@ -286,6 +424,8 @@ def histo_analysis(hist_data): #note: depreciated
         
     derivative_sorted = sorted(derivative, key=int)
     mean_derivative = basic_stats(derivative_sorted, "1d", 0)[0]
+    
+    print(mean_derivative)
     stdev_derivative = basic_stats(derivative_sorted, "1d", 0)[3]
 
     low_bound = mean_derivative + -1.645 * stdev_derivative
@@ -302,10 +442,10 @@ def histo_analysis(hist_data): #note: depreciated
 
     return [low_est, lm_est, mid_est, hm_est, high_est, stdev_derivative]
 
-def histo_analysis_2(hist_data, delta, low_bound, high_bound):
+def histo_analysis(hist_data, delta, low_bound, high_bound):
 
     if hist_data == 'debug':
-        return ('returns list of predicted values based on historical data; input delta for delta step in z-score and lower and igher bounds in number for standard deviations')
+        return ('returns list of predicted values based on historical data; input delta for delta step in z-score and lower and higher bounds in number for standard deviations')
 
     derivative = []
 
@@ -322,15 +462,108 @@ def histo_analysis_2(hist_data, delta, low_bound, high_bound):
     i = low_bound
 
     while True:
+        
+        if i > high_bound:
+            break
 
-        pred_change = mean_derivative + i * stdev_derivative
+        try:
+            pred_change = mean_derivative + i * stdev_derivative
+            
+        except:
+            
+            pred_change = mean_derivative
 
         predictions.append(float(hist_data[-1:][0]) + pred_change)
 
         i = i + delta
 
-        if i > high_bound:
-
-            break
-
     return predictions
+
+def poly_regression(x, y, power):
+
+    if x == "null":
+
+        x = []
+
+        for i in range(len(y)):
+
+            x.append(i)
+
+    reg_eq = scipy.polyfit(x, y, deg = power)
+
+    print(reg_eq)
+
+    eq_str = ""
+
+    for i in range(0, len(reg_eq), 1):
+
+        if i < len(reg_eq)- 1:
+            eq_str = eq_str + str(reg_eq[i]) + "*(z**" + str(len(reg_eq) - i - 1) + ")+"
+        else:
+            eq_str = eq_str + str(reg_eq[i]) + "*(z**" + str(len(reg_eq) - i - 1) + ")"
+
+    vals = []
+
+    for i in range(0, len(x), 1):
+        print(x[i])
+        z = x[i]
+
+        exec("vals.append(" + eq_str + ")")
+
+    print(vals)
+
+    _rms = rms(vals, y)
+
+    r2_d2 = r_squared(vals, y) 
+
+    return [eq_str, _rms, r2_d2]
+
+def r_squared(predictions, targets): # assumes equal size inputs
+
+    out = metrics.r2_score(targets, predictions)
+
+    return out
+
+def rms(predictions, targets): # assumes equal size inputs
+
+    out = 0
+
+    _sum = 0
+
+    avg = 0
+
+    for i in range(0, len(targets), 1):
+
+        _sum = (targets[i] - predictions[i]) ** 2
+
+    avg = _sum/len(targets)
+
+    out = math.sqrt(avg)
+
+    return float(out)
+
+def basic_analysis(filepath): #assumes that rows are the independent variable and columns are the dependant. also assumes that time flows from lowest column to highest column.
+    
+    data = load_csv(filepath)
+    row = len(data)
+    
+    column = []
+    
+    for i in range(0, row, 1):
+        
+        column.append(len(data[i]))
+        
+    column_max = max(column)
+    row_b_stats = []
+    row_histo = []
+    
+    for i in range(0, row, 1):
+        row_b_stats.append(basic_stats(data, "row", i))
+        row_histo.append(histo_analysis(data[i], 0.67449, -0.67449, 0.67449))
+    
+    column_b_stats = []
+    
+    for i in range(0, column_max, 1):
+        column_b_stats.append(basic_stats(data, "column", i))
+    
+    return[row_b_stats, column_b_stats, row_histo]

analysis_benchmark.py

@@ -0,0 +1,8 @@
+import analysis
+import time
+
+start = time.time()
+analysis.basic_analysis("data.txt")
+end = time.time()
+
+print(end - start)

analysis_test.py

@@ -6,24 +6,33 @@ print(analysis.basic_stats(0, 'debug', 0))
 print(analysis.basic_stats(data, "column", 0))
 print(analysis.basic_stats(data, "row", 0))
 print(analysis.z_score(10, analysis.basic_stats(data, "column", 0)[0],analysis.basic_stats(data, "column", 0)[3]))
-print(analysis.histo_analysis(data[0]))
+print(analysis.histo_analysis(data[0], 0.01, -1, 1))
-print(analysis.histo_analysis_2(data[0], 0.01, -1, 1))
 print(analysis.stdev_z_split(3.3, 0.2, 0.1, -5, 5))
 
+game_c_entities = analysis.c_entities(["bot", "bot", "bot"], [0, 1, 2], [[10, 10], [-10, -10], [10, -10]], ["shit", "bad", "worse"], ["triangle", "square", "circle"])
+game_c_entities.append("bot", 3, [-10, 10], "useless", "pentagram")
+game_c_entities.edit(0, "null", "null", "null", "null", "triagon")
+print(game_c_entities.search(0))
+print(game_c_entities.debug())
+print(game_c_entities.regurgitate())
+
 game_nc_entities = analysis.nc_entities(["cube", "cube", "ball"], [0, 1, 2], [[0, 0.5], [1, 1.5], [2, 2]], ["1;1;1;10', '2;1;1;20", "r=0.5, 5"], ["1", "1", "0"])
 game_nc_entities.append("cone", 3, [1, -1], "property", "effect")
 game_nc_entities.edit(2, "sphere", 10, [5, -5], "new prop", "new effect")
 print(game_nc_entities.search(10))
 print(game_nc_entities.debug())
+print(game_nc_entities.regurgitate())
 
 game_obstacles = analysis.obstacles(["wall", "fortress", "castle"], [0, 1, 2],[[[10, 10], [10, 9], [9, 10], [9, 9]], [[-10, 9], [-10, -9], [-9, -10]], [[5, 0], [4, -1], [-4, -1]]] , [0, 0.01, 10])
 game_obstacles.append("bastion", 3, [[50, 50], [49, 50], [50, 49], [49, 49]], 75)
 game_obstacles.edit(0, "motte and bailey", "null", [[10, 10], [9, 10], [10, 9], [9, 9]], 0.01)
 print(game_obstacles.search(0))
 print(game_obstacles.debug())
+print(game_obstacles.regurgitate())
 
 game_objectives = analysis.objectives(["switch", "scale", "climb"], [0,1,2], [[0,0],[1,1],[2,0]], ["0,1", "1,1", "0,5"])
 game_objectives.append("auto", 3, [0, 10], "1, 10")
 game_objectives.edit(3, "null", 4, "null", "null")
 print(game_objectives.search(4))
 print(game_objectives.debug())
+print(game_objectives.regurgitate())

generate_data.py

@@ -4,7 +4,7 @@ def generate(filename, x, y, low, high):
 
     file = open(filename, "w")
 
-    for i in range (0, y - 1, 1):
+    for i in range (0, y, 1):
 
         temp = ""
         

statistics.py

@@ -0,0 +1,669 @@
+"""
+Basic statistics module.
+
+This module provides functions for calculating statistics of data, including
+averages, variance, and standard deviation.
+
+Calculating averages
+--------------------
+
+==================  =============================================
+Function            Description
+==================  =============================================
+mean                Arithmetic mean (average) of data.
+harmonic_mean       Harmonic mean of data.
+median              Median (middle value) of data.
+median_low          Low median of data.
+median_high         High median of data.
+median_grouped      Median, or 50th percentile, of grouped data.
+mode                Mode (most common value) of data.
+==================  =============================================
+
+Calculate the arithmetic mean ("the average") of data:
+
+>>> mean([-1.0, 2.5, 3.25, 5.75])
+2.625
+
+
+Calculate the standard median of discrete data:
+
+>>> median([2, 3, 4, 5])
+3.5
+
+
+Calculate the median, or 50th percentile, of data grouped into class intervals
+centred on the data values provided. E.g. if your data points are rounded to
+the nearest whole number:
+
+>>> median_grouped([2, 2, 3, 3, 3, 4])  #doctest: +ELLIPSIS
+2.8333333333...
+
+This should be interpreted in this way: you have two data points in the class
+interval 1.5-2.5, three data points in the class interval 2.5-3.5, and one in
+the class interval 3.5-4.5. The median of these data points is 2.8333...
+
+
+Calculating variability or spread
+---------------------------------
+
+==================  =============================================
+Function            Description
+==================  =============================================
+pvariance           Population variance of data.
+variance            Sample variance of data.
+pstdev              Population standard deviation of data.
+stdev               Sample standard deviation of data.
+==================  =============================================
+
+Calculate the standard deviation of sample data:
+
+>>> stdev([2.5, 3.25, 5.5, 11.25, 11.75])  #doctest: +ELLIPSIS
+4.38961843444...
+
+If you have previously calculated the mean, you can pass it as the optional
+second argument to the four "spread" functions to avoid recalculating it:
+
+>>> data = [1, 2, 2, 4, 4, 4, 5, 6]
+>>> mu = mean(data)
+>>> pvariance(data, mu)
+2.5
+
+
+Exceptions
+----------
+
+A single exception is defined: StatisticsError is a subclass of ValueError.
+
+"""
+
+__all__ = [ 'StatisticsError',
+            'pstdev', 'pvariance', 'stdev', 'variance',
+            'median',  'median_low', 'median_high', 'median_grouped',
+            'mean', 'mode', 'harmonic_mean',
+          ]
+
+import collections
+import math
+import numbers
+
+from fractions import Fraction
+from decimal import Decimal
+from itertools import groupby
+from bisect import bisect_left, bisect_right
+
+
+
+# === Exceptions ===
+
+class StatisticsError(ValueError):
+    pass
+
+
+# === Private utilities ===
+
+def _sum(data, start=0):
+    """_sum(data [, start]) -> (type, sum, count)
+
+    Return a high-precision sum of the given numeric data as a fraction,
+    together with the type to be converted to and the count of items.
+
+    If optional argument ``start`` is given, it is added to the total.
+    If ``data`` is empty, ``start`` (defaulting to 0) is returned.
+
+
+    Examples
+    --------
+
+    >>> _sum([3, 2.25, 4.5, -0.5, 1.0], 0.75)
+    (<class 'float'>, Fraction(11, 1), 5)
+
+    Some sources of round-off error will be avoided:
+
+    # Built-in sum returns zero.
+    >>> _sum([1e50, 1, -1e50] * 1000)
+    (<class 'float'>, Fraction(1000, 1), 3000)
+
+    Fractions and Decimals are also supported:
+
+    >>> from fractions import Fraction as F
+    >>> _sum([F(2, 3), F(7, 5), F(1, 4), F(5, 6)])
+    (<class 'fractions.Fraction'>, Fraction(63, 20), 4)
+
+    >>> from decimal import Decimal as D
+    >>> data = [D("0.1375"), D("0.2108"), D("0.3061"), D("0.0419")]
+    >>> _sum(data)
+    (<class 'decimal.Decimal'>, Fraction(6963, 10000), 4)
+
+    Mixed types are currently treated as an error, except that int is
+    allowed.
+    """
+    count = 0
+    n, d = _exact_ratio(start)
+    partials = {d: n}
+    partials_get = partials.get
+    T = _coerce(int, type(start))
+    for typ, values in groupby(data, type):
+        T = _coerce(T, typ)  # or raise TypeError
+        for n,d in map(_exact_ratio, values):
+            count += 1
+            partials[d] = partials_get(d, 0) + n
+    if None in partials:
+        # The sum will be a NAN or INF. We can ignore all the finite
+        # partials, and just look at this special one.
+        total = partials[None]
+        assert not _isfinite(total)
+    else:
+        # Sum all the partial sums using builtin sum.
+        # FIXME is this faster if we sum them in order of the denominator?
+        total = sum(Fraction(n, d) for d, n in sorted(partials.items()))
+    return (T, total, count)
+
+
+def _isfinite(x):
+    try:
+        return x.is_finite()  # Likely a Decimal.
+    except AttributeError:
+        return math.isfinite(x)  # Coerces to float first.
+
+
+def _coerce(T, S):
+    """Coerce types T and S to a common type, or raise TypeError.
+
+    Coercion rules are currently an implementation detail. See the CoerceTest
+    test class in test_statistics for details.
+    """
+    # See http://bugs.python.org/issue24068.
+    assert T is not bool, "initial type T is bool"
+    # If the types are the same, no need to coerce anything. Put this
+    # first, so that the usual case (no coercion needed) happens as soon
+    # as possible.
+    if T is S:  return T
+    # Mixed int & other coerce to the other type.
+    if S is int or S is bool:  return T
+    if T is int:  return S
+    # If one is a (strict) subclass of the other, coerce to the subclass.
+    if issubclass(S, T):  return S
+    if issubclass(T, S):  return T
+    # Ints coerce to the other type.
+    if issubclass(T, int):  return S
+    if issubclass(S, int):  return T
+    # Mixed fraction & float coerces to float (or float subclass).
+    if issubclass(T, Fraction) and issubclass(S, float):
+        return S
+    if issubclass(T, float) and issubclass(S, Fraction):
+        return T
+    # Any other combination is disallowed.
+    msg = "don't know how to coerce %s and %s"
+    raise TypeError(msg % (T.__name__, S.__name__))
+
+
+def _exact_ratio(x):
+    """Return Real number x to exact (numerator, denominator) pair.
+
+    >>> _exact_ratio(0.25)
+    (1, 4)
+
+    x is expected to be an int, Fraction, Decimal or float.
+    """
+    try:
+        # Optimise the common case of floats. We expect that the most often
+        # used numeric type will be builtin floats, so try to make this as
+        # fast as possible.
+        if type(x) is float or type(x) is Decimal:
+            return x.as_integer_ratio()
+        try:
+            # x may be an int, Fraction, or Integral ABC.
+            return (x.numerator, x.denominator)
+        except AttributeError:
+            try:
+                # x may be a float or Decimal subclass.
+                return x.as_integer_ratio()
+            except AttributeError:
+                # Just give up?
+                pass
+    except (OverflowError, ValueError):
+        # float NAN or INF.
+        assert not _isfinite(x)
+        return (x, None)
+    msg = "can't convert type '{}' to numerator/denominator"
+    raise TypeError(msg.format(type(x).__name__))
+
+
+def _convert(value, T):
+    """Convert value to given numeric type T."""
+    if type(value) is T:
+        # This covers the cases where T is Fraction, or where value is
+        # a NAN or INF (Decimal or float).
+        return value
+    if issubclass(T, int) and value.denominator != 1:
+        T = float
+    try:
+        # FIXME: what do we do if this overflows?
+        return T(value)
+    except TypeError:
+        if issubclass(T, Decimal):
+            return T(value.numerator)/T(value.denominator)
+        else:
+            raise
+
+
+def _counts(data):
+    # Generate a table of sorted (value, frequency) pairs.
+    table = collections.Counter(iter(data)).most_common()
+    if not table:
+        return table
+    # Extract the values with the highest frequency.
+    maxfreq = table[0][1]
+    for i in range(1, len(table)):
+        if table[i][1] != maxfreq:
+            table = table[:i]
+            break
+    return table
+
+
+def _find_lteq(a, x):
+    'Locate the leftmost value exactly equal to x'
+    i = bisect_left(a, x)
+    if i != len(a) and a[i] == x:
+        return i
+    raise ValueError
+
+
+def _find_rteq(a, l, x):
+    'Locate the rightmost value exactly equal to x'
+    i = bisect_right(a, x, lo=l)
+    if i != (len(a)+1) and a[i-1] == x:
+        return i-1
+    raise ValueError
+
+
+def _fail_neg(values, errmsg='negative value'):
+    """Iterate over values, failing if any are less than zero."""
+    for x in values:
+        if x < 0:
+            raise StatisticsError(errmsg)
+        yield x
+
+
+# === Measures of central tendency (averages) ===
+
+def mean(data):
+    """Return the sample arithmetic mean of data.
+
+    >>> mean([1, 2, 3, 4, 4])
+    2.8
+
+    >>> from fractions import Fraction as F
+    >>> mean([F(3, 7), F(1, 21), F(5, 3), F(1, 3)])
+    Fraction(13, 21)
+
+    >>> from decimal import Decimal as D
+    >>> mean([D("0.5"), D("0.75"), D("0.625"), D("0.375")])
+    Decimal('0.5625')
+
+    If ``data`` is empty, StatisticsError will be raised.
+    """
+    if iter(data) is data:
+        data = list(data)
+    n = len(data)
+    if n < 1:
+        raise StatisticsError('mean requires at least one data point')
+    T, total, count = _sum(data)
+    assert count == n
+    return _convert(total/n, T)
+
+
+def harmonic_mean(data):
+    """Return the harmonic mean of data.
+
+    The harmonic mean, sometimes called the subcontrary mean, is the
+    reciprocal of the arithmetic mean of the reciprocals of the data,
+    and is often appropriate when averaging quantities which are rates
+    or ratios, for example speeds. Example:
+
+    Suppose an investor purchases an equal value of shares in each of
+    three companies, with P/E (price/earning) ratios of 2.5, 3 and 10.
+    What is the average P/E ratio for the investor's portfolio?
+
+    >>> harmonic_mean([2.5, 3, 10])  # For an equal investment portfolio.
+    3.6
+
+    Using the arithmetic mean would give an average of about 5.167, which
+    is too high.
+
+    If ``data`` is empty, or any element is less than zero,
+    ``harmonic_mean`` will raise ``StatisticsError``.
+    """
+    # For a justification for using harmonic mean for P/E ratios, see
+    # http://fixthepitch.pellucid.com/comps-analysis-the-missing-harmony-of-summary-statistics/
+    # http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2621087
+    if iter(data) is data:
+        data = list(data)
+    errmsg = 'harmonic mean does not support negative values'
+    n = len(data)
+    if n < 1:
+        raise StatisticsError('harmonic_mean requires at least one data point')
+    elif n == 1:
+        x = data[0]
+        if isinstance(x, (numbers.Real, Decimal)):
+            if x < 0:
+                raise StatisticsError(errmsg)
+            return x
+        else:
+            raise TypeError('unsupported type')
+    try:
+        T, total, count = _sum(1/x for x in _fail_neg(data, errmsg))
+    except ZeroDivisionError:
+        return 0
+    assert count == n
+    return _convert(n/total, T)
+
+
+# FIXME: investigate ways to calculate medians without sorting? Quickselect?
+def median(data):
+    """Return the median (middle value) of numeric data.
+
+    When the number of data points is odd, return the middle data point.
+    When the number of data points is even, the median is interpolated by
+    taking the average of the two middle values:
+
+    >>> median([1, 3, 5])
+    3
+    >>> median([1, 3, 5, 7])
+    4.0
+
+    """
+    data = sorted(data)
+    n = len(data)
+    if n == 0:
+        raise StatisticsError("no median for empty data")
+    if n%2 == 1:
+        return data[n//2]
+    else:
+        i = n//2
+        return (data[i - 1] + data[i])/2
+
+
+def median_low(data):
+    """Return the low median of numeric data.
+
+    When the number of data points is odd, the middle value is returned.
+    When it is even, the smaller of the two middle values is returned.
+
+    >>> median_low([1, 3, 5])
+    3
+    >>> median_low([1, 3, 5, 7])
+    3
+
+    """
+    data = sorted(data)
+    n = len(data)
+    if n == 0:
+        raise StatisticsError("no median for empty data")
+    if n%2 == 1:
+        return data[n//2]
+    else:
+        return data[n//2 - 1]
+
+
+def median_high(data):
+    """Return the high median of data.
+
+    When the number of data points is odd, the middle value is returned.
+    When it is even, the larger of the two middle values is returned.
+
+    >>> median_high([1, 3, 5])
+    3
+    >>> median_high([1, 3, 5, 7])
+    5
+
+    """
+    data = sorted(data)
+    n = len(data)
+    if n == 0:
+        raise StatisticsError("no median for empty data")
+    return data[n//2]
+
+
+def median_grouped(data, interval=1):
+    """Return the 50th percentile (median) of grouped continuous data.
+
+    >>> median_grouped([1, 2, 2, 3, 4, 4, 4, 4, 4, 5])
+    3.7
+    >>> median_grouped([52, 52, 53, 54])
+    52.5
+
+    This calculates the median as the 50th percentile, and should be
+    used when your data is continuous and grouped. In the above example,
+    the values 1, 2, 3, etc. actually represent the midpoint of classes
+    0.5-1.5, 1.5-2.5, 2.5-3.5, etc. The middle value falls somewhere in
+    class 3.5-4.5, and interpolation is used to estimate it.
+
+    Optional argument ``interval`` represents the class interval, and
+    defaults to 1. Changing the class interval naturally will change the
+    interpolated 50th percentile value:
+
+    >>> median_grouped([1, 3, 3, 5, 7], interval=1)
+    3.25
+    >>> median_grouped([1, 3, 3, 5, 7], interval=2)
+    3.5
+
+    This function does not check whether the data points are at least
+    ``interval`` apart.
+    """
+    data = sorted(data)
+    n = len(data)
+    if n == 0:
+        raise StatisticsError("no median for empty data")
+    elif n == 1:
+        return data[0]
+    # Find the value at the midpoint. Remember this corresponds to the
+    # centre of the class interval.
+    x = data[n//2]
+    for obj in (x, interval):
+        if isinstance(obj, (str, bytes)):
+            raise TypeError('expected number but got %r' % obj)
+    try:
+        L = x - interval/2  # The lower limit of the median interval.
+    except TypeError:
+        # Mixed type. For now we just coerce to float.
+        L = float(x) - float(interval)/2
+
+    # Uses bisection search to search for x in data with log(n) time complexity
+    # Find the position of leftmost occurrence of x in data
+    l1 = _find_lteq(data, x)
+    # Find the position of rightmost occurrence of x in data[l1...len(data)]
+    # Assuming always l1 <= l2
+    l2 = _find_rteq(data, l1, x)
+    cf = l1
+    f = l2 - l1 + 1
+    return L + interval*(n/2 - cf)/f
+
+
+def mode(data):
+    """Return the most common data point from discrete or nominal data.
+
+    ``mode`` assumes discrete data, and returns a single value. This is the
+    standard treatment of the mode as commonly taught in schools:
+
+    >>> mode([1, 1, 2, 3, 3, 3, 3, 4])
+    3
+
+    This also works with nominal (non-numeric) data:
+
+    >>> mode(["red", "blue", "blue", "red", "green", "red", "red"])
+    'red'
+
+    If there is not exactly one most common value, ``mode`` will raise
+    StatisticsError.
+    """
+    # Generate a table of sorted (value, frequency) pairs.
+    table = _counts(data)
+    if len(table) == 1:
+        return table[0][0]
+    elif table:
+        raise StatisticsError(
+                'no unique mode; found %d equally common values' % len(table)
+                )
+    else:
+        raise StatisticsError('no mode for empty data')
+
+
+# === Measures of spread ===
+
+# See http://mathworld.wolfram.com/Variance.html
+#     http://mathworld.wolfram.com/SampleVariance.html
+#     http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
+#
+# Under no circumstances use the so-called "computational formula for
+# variance", as that is only suitable for hand calculations with a small
+# amount of low-precision data. It has terrible numeric properties.
+#
+# See a comparison of three computational methods here:
+# http://www.johndcook.com/blog/2008/09/26/comparing-three-methods-of-computing-standard-deviation/
+
+def _ss(data, c=None):
+    """Return sum of square deviations of sequence data.
+
+    If ``c`` is None, the mean is calculated in one pass, and the deviations
+    from the mean are calculated in a second pass. Otherwise, deviations are
+    calculated from ``c`` as given. Use the second case with care, as it can
+    lead to garbage results.
+    """
+    if c is None:
+        c = mean(data)
+    T, total, count = _sum((x-c)**2 for x in data)
+    # The following sum should mathematically equal zero, but due to rounding
+    # error may not.
+    U, total2, count2 = _sum((x-c) for x in data)
+    assert T == U and count == count2
+    total -=  total2**2/len(data)
+    assert not total < 0, 'negative sum of square deviations: %f' % total
+    return (T, total)
+
+
+def variance(data, xbar=None):
+    """Return the sample variance of data.
+
+    data should be an iterable of Real-valued numbers, with at least two
+    values. The optional argument xbar, if given, should be the mean of
+    the data. If it is missing or None, the mean is automatically calculated.
+
+    Use this function when your data is a sample from a population. To
+    calculate the variance from the entire population, see ``pvariance``.
+
+    Examples:
+
+    >>> data = [2.75, 1.75, 1.25, 0.25, 0.5, 1.25, 3.5]
+    >>> variance(data)
+    1.3720238095238095
+
+    If you have already calculated the mean of your data, you can pass it as
+    the optional second argument ``xbar`` to avoid recalculating it:
+
+    >>> m = mean(data)
+    >>> variance(data, m)
+    1.3720238095238095
+
+    This function does not check that ``xbar`` is actually the mean of
+    ``data``. Giving arbitrary values for ``xbar`` may lead to invalid or
+    impossible results.
+
+    Decimals and Fractions are supported:
+
+    >>> from decimal import Decimal as D
+    >>> variance([D("27.5"), D("30.25"), D("30.25"), D("34.5"), D("41.75")])
+    Decimal('31.01875')
+
+    >>> from fractions import Fraction as F
+    >>> variance([F(1, 6), F(1, 2), F(5, 3)])
+    Fraction(67, 108)
+
+    """
+    if iter(data) is data:
+        data = list(data)
+    n = len(data)
+    if n < 2:
+        raise StatisticsError('variance requires at least two data points')
+    T, ss = _ss(data, xbar)
+    return _convert(ss/(n-1), T)
+
+
+def pvariance(data, mu=None):
+    """Return the population variance of ``data``.
+
+    data should be an iterable of Real-valued numbers, with at least one
+    value. The optional argument mu, if given, should be the mean of
+    the data. If it is missing or None, the mean is automatically calculated.
+
+    Use this function to calculate the variance from the entire population.
+    To estimate the variance from a sample, the ``variance`` function is
+    usually a better choice.
+
+    Examples:
+
+    >>> data = [0.0, 0.25, 0.25, 1.25, 1.5, 1.75, 2.75, 3.25]
+    >>> pvariance(data)
+    1.25
+
+    If you have already calculated the mean of the data, you can pass it as
+    the optional second argument to avoid recalculating it:
+
+    >>> mu = mean(data)
+    >>> pvariance(data, mu)
+    1.25
+
+    This function does not check that ``mu`` is actually the mean of ``data``.
+    Giving arbitrary values for ``mu`` may lead to invalid or impossible
+    results.
+
+    Decimals and Fractions are supported:
+
+    >>> from decimal import Decimal as D
+    >>> pvariance([D("27.5"), D("30.25"), D("30.25"), D("34.5"), D("41.75")])
+    Decimal('24.815')
+
+    >>> from fractions import Fraction as F
+    >>> pvariance([F(1, 4), F(5, 4), F(1, 2)])
+    Fraction(13, 72)
+
+    """
+    if iter(data) is data:
+        data = list(data)
+    n = len(data)
+    if n < 1:
+        raise StatisticsError('pvariance requires at least one data point')
+    T, ss = _ss(data, mu)
+    return _convert(ss/n, T)
+
+
+def stdev(data, xbar=None):
+    """Return the square root of the sample variance.
+
+    See ``variance`` for arguments and other details.
+
+    >>> stdev([1.5, 2.5, 2.5, 2.75, 3.25, 4.75])
+    1.0810874155219827
+
+    """
+    var = variance(data, xbar)
+    try:
+        return var.sqrt()
+    except AttributeError:
+        return math.sqrt(var)
+
+
+def pstdev(data, mu=None):
+    """Return the square root of the population variance.
+
+    See ``pvariance`` for arguments and other details.
+
+    >>> pstdev([1.5, 2.5, 2.5, 2.75, 3.25, 4.75])
+    0.986893273527251
+
+    """
+    var = pvariance(data, mu)
+    try:
+        return var.sqrt()
+    except AttributeError:
+        return math.sqrt(var)