| ##  Module statistics.py | 
 | ## | 
 | ##  Copyright (c) 2013 Steven D'Aprano <steve+python@pearwood.info>. | 
 | ## | 
 | ##  Licensed under the Apache License, Version 2.0 (the "License"); | 
 | ##  you may not use this file except in compliance with the License. | 
 | ##  You may obtain a copy of the License at | 
 | ## | 
 | ##  http://www.apache.org/licenses/LICENSE-2.0 | 
 | ## | 
 | ##  Unless required by applicable law or agreed to in writing, software | 
 | ##  distributed under the License is distributed on an "AS IS" BASIS, | 
 | ##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | 
 | ##  See the License for the specific language governing permissions and | 
 | ##  limitations under the License. | 
 |  | 
 |  | 
 | """ | 
 | 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. | 
 | 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', | 
 |           ] | 
 |  | 
 |  | 
 | import collections | 
 | import math | 
 |  | 
 | from fractions import Fraction | 
 | from decimal import Decimal | 
 | from itertools import groupby | 
 |  | 
 |  | 
 |  | 
 | # === 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: | 
 |  | 
 |     >>> _sum([1e50, 1, -1e50] * 1000)  # Built-in sum returns zero. | 
 |     (<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: | 
 |             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 subclass. | 
 |                 return x.as_integer_ratio() | 
 |             except AttributeError: | 
 |                 try: | 
 |                     # x may be a Decimal. | 
 |                     return _decimal_to_ratio(x) | 
 |                 except AttributeError: | 
 |                     # Just give up? | 
 |                     pass | 
 |     except (OverflowError, ValueError): | 
 |         # float NAN or INF. | 
 |         assert not math.isfinite(x) | 
 |         return (x, None) | 
 |     msg = "can't convert type '{}' to numerator/denominator" | 
 |     raise TypeError(msg.format(type(x).__name__)) | 
 |  | 
 |  | 
 | # FIXME This is faster than Fraction.from_decimal, but still too slow. | 
 | def _decimal_to_ratio(d): | 
 |     """Convert Decimal d to exact integer ratio (numerator, denominator). | 
 |  | 
 |     >>> from decimal import Decimal | 
 |     >>> _decimal_to_ratio(Decimal("2.6")) | 
 |     (26, 10) | 
 |  | 
 |     """ | 
 |     sign, digits, exp = d.as_tuple() | 
 |     if exp in ('F', 'n', 'N'):  # INF, NAN, sNAN | 
 |         assert not d.is_finite() | 
 |         return (d, None) | 
 |     num = 0 | 
 |     for digit in digits: | 
 |         num = num*10 + digit | 
 |     if exp < 0: | 
 |         den = 10**-exp | 
 |     else: | 
 |         num *= 10**exp | 
 |         den = 1 | 
 |     if sign: | 
 |         num = -num | 
 |     return (num, den) | 
 |  | 
 |  | 
 | 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 | 
 |  | 
 |  | 
 | # === 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) | 
 |  | 
 |  | 
 | # 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 | 
 |     cf = data.index(x)  # Number of values below the median interval. | 
 |     # FIXME The following line could be more efficient for big lists. | 
 |     f = data.count(x)  # Number of data points in the median interval. | 
 |     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') | 
 |     ss = _ss(data, mu) | 
 |     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) |