logo

oasis-root

Compiled tree of Oasis Linux based on own branch at <https://hacktivis.me/git/oasis/> git clone https://anongit.hacktivis.me/git/oasis-root.git

statistics.py (43160B)


  1. """
  2. Basic statistics module.
  3. This module provides functions for calculating statistics of data, including
  4. averages, variance, and standard deviation.
  5. Calculating averages
  6. --------------------
  7. ================== ==================================================
  8. Function Description
  9. ================== ==================================================
  10. mean Arithmetic mean (average) of data.
  11. fmean Fast, floating point arithmetic mean.
  12. geometric_mean Geometric mean of data.
  13. harmonic_mean Harmonic mean of data.
  14. median Median (middle value) of data.
  15. median_low Low median of data.
  16. median_high High median of data.
  17. median_grouped Median, or 50th percentile, of grouped data.
  18. mode Mode (most common value) of data.
  19. multimode List of modes (most common values of data).
  20. quantiles Divide data into intervals with equal probability.
  21. ================== ==================================================
  22. Calculate the arithmetic mean ("the average") of data:
  23. >>> mean([-1.0, 2.5, 3.25, 5.75])
  24. 2.625
  25. Calculate the standard median of discrete data:
  26. >>> median([2, 3, 4, 5])
  27. 3.5
  28. Calculate the median, or 50th percentile, of data grouped into class intervals
  29. centred on the data values provided. E.g. if your data points are rounded to
  30. the nearest whole number:
  31. >>> median_grouped([2, 2, 3, 3, 3, 4]) #doctest: +ELLIPSIS
  32. 2.8333333333...
  33. This should be interpreted in this way: you have two data points in the class
  34. interval 1.5-2.5, three data points in the class interval 2.5-3.5, and one in
  35. the class interval 3.5-4.5. The median of these data points is 2.8333...
  36. Calculating variability or spread
  37. ---------------------------------
  38. ================== =============================================
  39. Function Description
  40. ================== =============================================
  41. pvariance Population variance of data.
  42. variance Sample variance of data.
  43. pstdev Population standard deviation of data.
  44. stdev Sample standard deviation of data.
  45. ================== =============================================
  46. Calculate the standard deviation of sample data:
  47. >>> stdev([2.5, 3.25, 5.5, 11.25, 11.75]) #doctest: +ELLIPSIS
  48. 4.38961843444...
  49. If you have previously calculated the mean, you can pass it as the optional
  50. second argument to the four "spread" functions to avoid recalculating it:
  51. >>> data = [1, 2, 2, 4, 4, 4, 5, 6]
  52. >>> mu = mean(data)
  53. >>> pvariance(data, mu)
  54. 2.5
  55. Statistics for relations between two inputs
  56. -------------------------------------------
  57. ================== ====================================================
  58. Function Description
  59. ================== ====================================================
  60. covariance Sample covariance for two variables.
  61. correlation Pearson's correlation coefficient for two variables.
  62. linear_regression Intercept and slope for simple linear regression.
  63. ================== ====================================================
  64. Calculate covariance, Pearson's correlation, and simple linear regression
  65. for two inputs:
  66. >>> x = [1, 2, 3, 4, 5, 6, 7, 8, 9]
  67. >>> y = [1, 2, 3, 1, 2, 3, 1, 2, 3]
  68. >>> covariance(x, y)
  69. 0.75
  70. >>> correlation(x, y) #doctest: +ELLIPSIS
  71. 0.31622776601...
  72. >>> linear_regression(x, y) #doctest:
  73. LinearRegression(slope=0.1, intercept=1.5)
  74. Exceptions
  75. ----------
  76. A single exception is defined: StatisticsError is a subclass of ValueError.
  77. """
  78. __all__ = [
  79. 'NormalDist',
  80. 'StatisticsError',
  81. 'correlation',
  82. 'covariance',
  83. 'fmean',
  84. 'geometric_mean',
  85. 'harmonic_mean',
  86. 'linear_regression',
  87. 'mean',
  88. 'median',
  89. 'median_grouped',
  90. 'median_high',
  91. 'median_low',
  92. 'mode',
  93. 'multimode',
  94. 'pstdev',
  95. 'pvariance',
  96. 'quantiles',
  97. 'stdev',
  98. 'variance',
  99. ]
  100. import math
  101. import numbers
  102. import random
  103. from fractions import Fraction
  104. from decimal import Decimal
  105. from itertools import groupby, repeat
  106. from bisect import bisect_left, bisect_right
  107. from math import hypot, sqrt, fabs, exp, erf, tau, log, fsum
  108. from operator import itemgetter
  109. from collections import Counter, namedtuple
  110. # === Exceptions ===
  111. class StatisticsError(ValueError):
  112. pass
  113. # === Private utilities ===
  114. def _sum(data, start=0):
  115. """_sum(data [, start]) -> (type, sum, count)
  116. Return a high-precision sum of the given numeric data as a fraction,
  117. together with the type to be converted to and the count of items.
  118. If optional argument ``start`` is given, it is added to the total.
  119. If ``data`` is empty, ``start`` (defaulting to 0) is returned.
  120. Examples
  121. --------
  122. >>> _sum([3, 2.25, 4.5, -0.5, 1.0], 0.75)
  123. (<class 'float'>, Fraction(11, 1), 5)
  124. Some sources of round-off error will be avoided:
  125. # Built-in sum returns zero.
  126. >>> _sum([1e50, 1, -1e50] * 1000)
  127. (<class 'float'>, Fraction(1000, 1), 3000)
  128. Fractions and Decimals are also supported:
  129. >>> from fractions import Fraction as F
  130. >>> _sum([F(2, 3), F(7, 5), F(1, 4), F(5, 6)])
  131. (<class 'fractions.Fraction'>, Fraction(63, 20), 4)
  132. >>> from decimal import Decimal as D
  133. >>> data = [D("0.1375"), D("0.2108"), D("0.3061"), D("0.0419")]
  134. >>> _sum(data)
  135. (<class 'decimal.Decimal'>, Fraction(6963, 10000), 4)
  136. Mixed types are currently treated as an error, except that int is
  137. allowed.
  138. """
  139. count = 0
  140. n, d = _exact_ratio(start)
  141. partials = {d: n}
  142. partials_get = partials.get
  143. T = _coerce(int, type(start))
  144. for typ, values in groupby(data, type):
  145. T = _coerce(T, typ) # or raise TypeError
  146. for n, d in map(_exact_ratio, values):
  147. count += 1
  148. partials[d] = partials_get(d, 0) + n
  149. if None in partials:
  150. # The sum will be a NAN or INF. We can ignore all the finite
  151. # partials, and just look at this special one.
  152. total = partials[None]
  153. assert not _isfinite(total)
  154. else:
  155. # Sum all the partial sums using builtin sum.
  156. # FIXME is this faster if we sum them in order of the denominator?
  157. total = sum(Fraction(n, d) for d, n in sorted(partials.items()))
  158. return (T, total, count)
  159. def _isfinite(x):
  160. try:
  161. return x.is_finite() # Likely a Decimal.
  162. except AttributeError:
  163. return math.isfinite(x) # Coerces to float first.
  164. def _coerce(T, S):
  165. """Coerce types T and S to a common type, or raise TypeError.
  166. Coercion rules are currently an implementation detail. See the CoerceTest
  167. test class in test_statistics for details.
  168. """
  169. # See http://bugs.python.org/issue24068.
  170. assert T is not bool, "initial type T is bool"
  171. # If the types are the same, no need to coerce anything. Put this
  172. # first, so that the usual case (no coercion needed) happens as soon
  173. # as possible.
  174. if T is S: return T
  175. # Mixed int & other coerce to the other type.
  176. if S is int or S is bool: return T
  177. if T is int: return S
  178. # If one is a (strict) subclass of the other, coerce to the subclass.
  179. if issubclass(S, T): return S
  180. if issubclass(T, S): return T
  181. # Ints coerce to the other type.
  182. if issubclass(T, int): return S
  183. if issubclass(S, int): return T
  184. # Mixed fraction & float coerces to float (or float subclass).
  185. if issubclass(T, Fraction) and issubclass(S, float):
  186. return S
  187. if issubclass(T, float) and issubclass(S, Fraction):
  188. return T
  189. # Any other combination is disallowed.
  190. msg = "don't know how to coerce %s and %s"
  191. raise TypeError(msg % (T.__name__, S.__name__))
  192. def _exact_ratio(x):
  193. """Return Real number x to exact (numerator, denominator) pair.
  194. >>> _exact_ratio(0.25)
  195. (1, 4)
  196. x is expected to be an int, Fraction, Decimal or float.
  197. """
  198. try:
  199. # Optimise the common case of floats. We expect that the most often
  200. # used numeric type will be builtin floats, so try to make this as
  201. # fast as possible.
  202. if type(x) is float or type(x) is Decimal:
  203. return x.as_integer_ratio()
  204. try:
  205. # x may be an int, Fraction, or Integral ABC.
  206. return (x.numerator, x.denominator)
  207. except AttributeError:
  208. try:
  209. # x may be a float or Decimal subclass.
  210. return x.as_integer_ratio()
  211. except AttributeError:
  212. # Just give up?
  213. pass
  214. except (OverflowError, ValueError):
  215. # float NAN or INF.
  216. assert not _isfinite(x)
  217. return (x, None)
  218. msg = "can't convert type '{}' to numerator/denominator"
  219. raise TypeError(msg.format(type(x).__name__))
  220. def _convert(value, T):
  221. """Convert value to given numeric type T."""
  222. if type(value) is T:
  223. # This covers the cases where T is Fraction, or where value is
  224. # a NAN or INF (Decimal or float).
  225. return value
  226. if issubclass(T, int) and value.denominator != 1:
  227. T = float
  228. try:
  229. # FIXME: what do we do if this overflows?
  230. return T(value)
  231. except TypeError:
  232. if issubclass(T, Decimal):
  233. return T(value.numerator) / T(value.denominator)
  234. else:
  235. raise
  236. def _find_lteq(a, x):
  237. 'Locate the leftmost value exactly equal to x'
  238. i = bisect_left(a, x)
  239. if i != len(a) and a[i] == x:
  240. return i
  241. raise ValueError
  242. def _find_rteq(a, l, x):
  243. 'Locate the rightmost value exactly equal to x'
  244. i = bisect_right(a, x, lo=l)
  245. if i != (len(a) + 1) and a[i - 1] == x:
  246. return i - 1
  247. raise ValueError
  248. def _fail_neg(values, errmsg='negative value'):
  249. """Iterate over values, failing if any are less than zero."""
  250. for x in values:
  251. if x < 0:
  252. raise StatisticsError(errmsg)
  253. yield x
  254. # === Measures of central tendency (averages) ===
  255. def mean(data):
  256. """Return the sample arithmetic mean of data.
  257. >>> mean([1, 2, 3, 4, 4])
  258. 2.8
  259. >>> from fractions import Fraction as F
  260. >>> mean([F(3, 7), F(1, 21), F(5, 3), F(1, 3)])
  261. Fraction(13, 21)
  262. >>> from decimal import Decimal as D
  263. >>> mean([D("0.5"), D("0.75"), D("0.625"), D("0.375")])
  264. Decimal('0.5625')
  265. If ``data`` is empty, StatisticsError will be raised.
  266. """
  267. if iter(data) is data:
  268. data = list(data)
  269. n = len(data)
  270. if n < 1:
  271. raise StatisticsError('mean requires at least one data point')
  272. T, total, count = _sum(data)
  273. assert count == n
  274. return _convert(total / n, T)
  275. def fmean(data):
  276. """Convert data to floats and compute the arithmetic mean.
  277. This runs faster than the mean() function and it always returns a float.
  278. If the input dataset is empty, it raises a StatisticsError.
  279. >>> fmean([3.5, 4.0, 5.25])
  280. 4.25
  281. """
  282. try:
  283. n = len(data)
  284. except TypeError:
  285. # Handle iterators that do not define __len__().
  286. n = 0
  287. def count(iterable):
  288. nonlocal n
  289. for n, x in enumerate(iterable, start=1):
  290. yield x
  291. total = fsum(count(data))
  292. else:
  293. total = fsum(data)
  294. try:
  295. return total / n
  296. except ZeroDivisionError:
  297. raise StatisticsError('fmean requires at least one data point') from None
  298. def geometric_mean(data):
  299. """Convert data to floats and compute the geometric mean.
  300. Raises a StatisticsError if the input dataset is empty,
  301. if it contains a zero, or if it contains a negative value.
  302. No special efforts are made to achieve exact results.
  303. (However, this may change in the future.)
  304. >>> round(geometric_mean([54, 24, 36]), 9)
  305. 36.0
  306. """
  307. try:
  308. return exp(fmean(map(log, data)))
  309. except ValueError:
  310. raise StatisticsError('geometric mean requires a non-empty dataset '
  311. ' containing positive numbers') from None
  312. def harmonic_mean(data, weights=None):
  313. """Return the harmonic mean of data.
  314. The harmonic mean is the reciprocal of the arithmetic mean of the
  315. reciprocals of the data. It can be used for averaging ratios or
  316. rates, for example speeds.
  317. Suppose a car travels 40 km/hr for 5 km and then speeds-up to
  318. 60 km/hr for another 5 km. What is the average speed?
  319. >>> harmonic_mean([40, 60])
  320. 48.0
  321. Suppose a car travels 40 km/hr for 5 km, and when traffic clears,
  322. speeds-up to 60 km/hr for the remaining 30 km of the journey. What
  323. is the average speed?
  324. >>> harmonic_mean([40, 60], weights=[5, 30])
  325. 56.0
  326. If ``data`` is empty, or any element is less than zero,
  327. ``harmonic_mean`` will raise ``StatisticsError``.
  328. """
  329. if iter(data) is data:
  330. data = list(data)
  331. errmsg = 'harmonic mean does not support negative values'
  332. n = len(data)
  333. if n < 1:
  334. raise StatisticsError('harmonic_mean requires at least one data point')
  335. elif n == 1 and weights is None:
  336. x = data[0]
  337. if isinstance(x, (numbers.Real, Decimal)):
  338. if x < 0:
  339. raise StatisticsError(errmsg)
  340. return x
  341. else:
  342. raise TypeError('unsupported type')
  343. if weights is None:
  344. weights = repeat(1, n)
  345. sum_weights = n
  346. else:
  347. if iter(weights) is weights:
  348. weights = list(weights)
  349. if len(weights) != n:
  350. raise StatisticsError('Number of weights does not match data size')
  351. _, sum_weights, _ = _sum(w for w in _fail_neg(weights, errmsg))
  352. try:
  353. data = _fail_neg(data, errmsg)
  354. T, total, count = _sum(w / x if w else 0 for w, x in zip(weights, data))
  355. except ZeroDivisionError:
  356. return 0
  357. if total <= 0:
  358. raise StatisticsError('Weighted sum must be positive')
  359. return _convert(sum_weights / total, T)
  360. # FIXME: investigate ways to calculate medians without sorting? Quickselect?
  361. def median(data):
  362. """Return the median (middle value) of numeric data.
  363. When the number of data points is odd, return the middle data point.
  364. When the number of data points is even, the median is interpolated by
  365. taking the average of the two middle values:
  366. >>> median([1, 3, 5])
  367. 3
  368. >>> median([1, 3, 5, 7])
  369. 4.0
  370. """
  371. data = sorted(data)
  372. n = len(data)
  373. if n == 0:
  374. raise StatisticsError("no median for empty data")
  375. if n % 2 == 1:
  376. return data[n // 2]
  377. else:
  378. i = n // 2
  379. return (data[i - 1] + data[i]) / 2
  380. def median_low(data):
  381. """Return the low median of numeric data.
  382. When the number of data points is odd, the middle value is returned.
  383. When it is even, the smaller of the two middle values is returned.
  384. >>> median_low([1, 3, 5])
  385. 3
  386. >>> median_low([1, 3, 5, 7])
  387. 3
  388. """
  389. data = sorted(data)
  390. n = len(data)
  391. if n == 0:
  392. raise StatisticsError("no median for empty data")
  393. if n % 2 == 1:
  394. return data[n // 2]
  395. else:
  396. return data[n // 2 - 1]
  397. def median_high(data):
  398. """Return the high median of data.
  399. When the number of data points is odd, the middle value is returned.
  400. When it is even, the larger of the two middle values is returned.
  401. >>> median_high([1, 3, 5])
  402. 3
  403. >>> median_high([1, 3, 5, 7])
  404. 5
  405. """
  406. data = sorted(data)
  407. n = len(data)
  408. if n == 0:
  409. raise StatisticsError("no median for empty data")
  410. return data[n // 2]
  411. def median_grouped(data, interval=1):
  412. """Return the 50th percentile (median) of grouped continuous data.
  413. >>> median_grouped([1, 2, 2, 3, 4, 4, 4, 4, 4, 5])
  414. 3.7
  415. >>> median_grouped([52, 52, 53, 54])
  416. 52.5
  417. This calculates the median as the 50th percentile, and should be
  418. used when your data is continuous and grouped. In the above example,
  419. the values 1, 2, 3, etc. actually represent the midpoint of classes
  420. 0.5-1.5, 1.5-2.5, 2.5-3.5, etc. The middle value falls somewhere in
  421. class 3.5-4.5, and interpolation is used to estimate it.
  422. Optional argument ``interval`` represents the class interval, and
  423. defaults to 1. Changing the class interval naturally will change the
  424. interpolated 50th percentile value:
  425. >>> median_grouped([1, 3, 3, 5, 7], interval=1)
  426. 3.25
  427. >>> median_grouped([1, 3, 3, 5, 7], interval=2)
  428. 3.5
  429. This function does not check whether the data points are at least
  430. ``interval`` apart.
  431. """
  432. data = sorted(data)
  433. n = len(data)
  434. if n == 0:
  435. raise StatisticsError("no median for empty data")
  436. elif n == 1:
  437. return data[0]
  438. # Find the value at the midpoint. Remember this corresponds to the
  439. # centre of the class interval.
  440. x = data[n // 2]
  441. for obj in (x, interval):
  442. if isinstance(obj, (str, bytes)):
  443. raise TypeError('expected number but got %r' % obj)
  444. try:
  445. L = x - interval / 2 # The lower limit of the median interval.
  446. except TypeError:
  447. # Mixed type. For now we just coerce to float.
  448. L = float(x) - float(interval) / 2
  449. # Uses bisection search to search for x in data with log(n) time complexity
  450. # Find the position of leftmost occurrence of x in data
  451. l1 = _find_lteq(data, x)
  452. # Find the position of rightmost occurrence of x in data[l1...len(data)]
  453. # Assuming always l1 <= l2
  454. l2 = _find_rteq(data, l1, x)
  455. cf = l1
  456. f = l2 - l1 + 1
  457. return L + interval * (n / 2 - cf) / f
  458. def mode(data):
  459. """Return the most common data point from discrete or nominal data.
  460. ``mode`` assumes discrete data, and returns a single value. This is the
  461. standard treatment of the mode as commonly taught in schools:
  462. >>> mode([1, 1, 2, 3, 3, 3, 3, 4])
  463. 3
  464. This also works with nominal (non-numeric) data:
  465. >>> mode(["red", "blue", "blue", "red", "green", "red", "red"])
  466. 'red'
  467. If there are multiple modes with same frequency, return the first one
  468. encountered:
  469. >>> mode(['red', 'red', 'green', 'blue', 'blue'])
  470. 'red'
  471. If *data* is empty, ``mode``, raises StatisticsError.
  472. """
  473. pairs = Counter(iter(data)).most_common(1)
  474. try:
  475. return pairs[0][0]
  476. except IndexError:
  477. raise StatisticsError('no mode for empty data') from None
  478. def multimode(data):
  479. """Return a list of the most frequently occurring values.
  480. Will return more than one result if there are multiple modes
  481. or an empty list if *data* is empty.
  482. >>> multimode('aabbbbbbbbcc')
  483. ['b']
  484. >>> multimode('aabbbbccddddeeffffgg')
  485. ['b', 'd', 'f']
  486. >>> multimode('')
  487. []
  488. """
  489. counts = Counter(iter(data)).most_common()
  490. maxcount, mode_items = next(groupby(counts, key=itemgetter(1)), (0, []))
  491. return list(map(itemgetter(0), mode_items))
  492. # Notes on methods for computing quantiles
  493. # ----------------------------------------
  494. #
  495. # There is no one perfect way to compute quantiles. Here we offer
  496. # two methods that serve common needs. Most other packages
  497. # surveyed offered at least one or both of these two, making them
  498. # "standard" in the sense of "widely-adopted and reproducible".
  499. # They are also easy to explain, easy to compute manually, and have
  500. # straight-forward interpretations that aren't surprising.
  501. # The default method is known as "R6", "PERCENTILE.EXC", or "expected
  502. # value of rank order statistics". The alternative method is known as
  503. # "R7", "PERCENTILE.INC", or "mode of rank order statistics".
  504. # For sample data where there is a positive probability for values
  505. # beyond the range of the data, the R6 exclusive method is a
  506. # reasonable choice. Consider a random sample of nine values from a
  507. # population with a uniform distribution from 0.0 to 1.0. The
  508. # distribution of the third ranked sample point is described by
  509. # betavariate(alpha=3, beta=7) which has mode=0.250, median=0.286, and
  510. # mean=0.300. Only the latter (which corresponds with R6) gives the
  511. # desired cut point with 30% of the population falling below that
  512. # value, making it comparable to a result from an inv_cdf() function.
  513. # The R6 exclusive method is also idempotent.
  514. # For describing population data where the end points are known to
  515. # be included in the data, the R7 inclusive method is a reasonable
  516. # choice. Instead of the mean, it uses the mode of the beta
  517. # distribution for the interior points. Per Hyndman & Fan, "One nice
  518. # property is that the vertices of Q7(p) divide the range into n - 1
  519. # intervals, and exactly 100p% of the intervals lie to the left of
  520. # Q7(p) and 100(1 - p)% of the intervals lie to the right of Q7(p)."
  521. # If needed, other methods could be added. However, for now, the
  522. # position is that fewer options make for easier choices and that
  523. # external packages can be used for anything more advanced.
  524. def quantiles(data, *, n=4, method='exclusive'):
  525. """Divide *data* into *n* continuous intervals with equal probability.
  526. Returns a list of (n - 1) cut points separating the intervals.
  527. Set *n* to 4 for quartiles (the default). Set *n* to 10 for deciles.
  528. Set *n* to 100 for percentiles which gives the 99 cuts points that
  529. separate *data* in to 100 equal sized groups.
  530. The *data* can be any iterable containing sample.
  531. The cut points are linearly interpolated between data points.
  532. If *method* is set to *inclusive*, *data* is treated as population
  533. data. The minimum value is treated as the 0th percentile and the
  534. maximum value is treated as the 100th percentile.
  535. """
  536. if n < 1:
  537. raise StatisticsError('n must be at least 1')
  538. data = sorted(data)
  539. ld = len(data)
  540. if ld < 2:
  541. raise StatisticsError('must have at least two data points')
  542. if method == 'inclusive':
  543. m = ld - 1
  544. result = []
  545. for i in range(1, n):
  546. j, delta = divmod(i * m, n)
  547. interpolated = (data[j] * (n - delta) + data[j + 1] * delta) / n
  548. result.append(interpolated)
  549. return result
  550. if method == 'exclusive':
  551. m = ld + 1
  552. result = []
  553. for i in range(1, n):
  554. j = i * m // n # rescale i to m/n
  555. j = 1 if j < 1 else ld-1 if j > ld-1 else j # clamp to 1 .. ld-1
  556. delta = i*m - j*n # exact integer math
  557. interpolated = (data[j - 1] * (n - delta) + data[j] * delta) / n
  558. result.append(interpolated)
  559. return result
  560. raise ValueError(f'Unknown method: {method!r}')
  561. # === Measures of spread ===
  562. # See http://mathworld.wolfram.com/Variance.html
  563. # http://mathworld.wolfram.com/SampleVariance.html
  564. # http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
  565. #
  566. # Under no circumstances use the so-called "computational formula for
  567. # variance", as that is only suitable for hand calculations with a small
  568. # amount of low-precision data. It has terrible numeric properties.
  569. #
  570. # See a comparison of three computational methods here:
  571. # http://www.johndcook.com/blog/2008/09/26/comparing-three-methods-of-computing-standard-deviation/
  572. def _ss(data, c=None):
  573. """Return sum of square deviations of sequence data.
  574. If ``c`` is None, the mean is calculated in one pass, and the deviations
  575. from the mean are calculated in a second pass. Otherwise, deviations are
  576. calculated from ``c`` as given. Use the second case with care, as it can
  577. lead to garbage results.
  578. """
  579. if c is not None:
  580. T, total, count = _sum((x-c)**2 for x in data)
  581. return (T, total)
  582. c = mean(data)
  583. T, total, count = _sum((x-c)**2 for x in data)
  584. # The following sum should mathematically equal zero, but due to rounding
  585. # error may not.
  586. U, total2, count2 = _sum((x - c) for x in data)
  587. assert T == U and count == count2
  588. total -= total2 ** 2 / len(data)
  589. assert not total < 0, 'negative sum of square deviations: %f' % total
  590. return (T, total)
  591. def variance(data, xbar=None):
  592. """Return the sample variance of data.
  593. data should be an iterable of Real-valued numbers, with at least two
  594. values. The optional argument xbar, if given, should be the mean of
  595. the data. If it is missing or None, the mean is automatically calculated.
  596. Use this function when your data is a sample from a population. To
  597. calculate the variance from the entire population, see ``pvariance``.
  598. Examples:
  599. >>> data = [2.75, 1.75, 1.25, 0.25, 0.5, 1.25, 3.5]
  600. >>> variance(data)
  601. 1.3720238095238095
  602. If you have already calculated the mean of your data, you can pass it as
  603. the optional second argument ``xbar`` to avoid recalculating it:
  604. >>> m = mean(data)
  605. >>> variance(data, m)
  606. 1.3720238095238095
  607. This function does not check that ``xbar`` is actually the mean of
  608. ``data``. Giving arbitrary values for ``xbar`` may lead to invalid or
  609. impossible results.
  610. Decimals and Fractions are supported:
  611. >>> from decimal import Decimal as D
  612. >>> variance([D("27.5"), D("30.25"), D("30.25"), D("34.5"), D("41.75")])
  613. Decimal('31.01875')
  614. >>> from fractions import Fraction as F
  615. >>> variance([F(1, 6), F(1, 2), F(5, 3)])
  616. Fraction(67, 108)
  617. """
  618. if iter(data) is data:
  619. data = list(data)
  620. n = len(data)
  621. if n < 2:
  622. raise StatisticsError('variance requires at least two data points')
  623. T, ss = _ss(data, xbar)
  624. return _convert(ss / (n - 1), T)
  625. def pvariance(data, mu=None):
  626. """Return the population variance of ``data``.
  627. data should be a sequence or iterable of Real-valued numbers, with at least one
  628. value. The optional argument mu, if given, should be the mean of
  629. the data. If it is missing or None, the mean is automatically calculated.
  630. Use this function to calculate the variance from the entire population.
  631. To estimate the variance from a sample, the ``variance`` function is
  632. usually a better choice.
  633. Examples:
  634. >>> data = [0.0, 0.25, 0.25, 1.25, 1.5, 1.75, 2.75, 3.25]
  635. >>> pvariance(data)
  636. 1.25
  637. If you have already calculated the mean of the data, you can pass it as
  638. the optional second argument to avoid recalculating it:
  639. >>> mu = mean(data)
  640. >>> pvariance(data, mu)
  641. 1.25
  642. Decimals and Fractions are supported:
  643. >>> from decimal import Decimal as D
  644. >>> pvariance([D("27.5"), D("30.25"), D("30.25"), D("34.5"), D("41.75")])
  645. Decimal('24.815')
  646. >>> from fractions import Fraction as F
  647. >>> pvariance([F(1, 4), F(5, 4), F(1, 2)])
  648. Fraction(13, 72)
  649. """
  650. if iter(data) is data:
  651. data = list(data)
  652. n = len(data)
  653. if n < 1:
  654. raise StatisticsError('pvariance requires at least one data point')
  655. T, ss = _ss(data, mu)
  656. return _convert(ss / n, T)
  657. def stdev(data, xbar=None):
  658. """Return the square root of the sample variance.
  659. See ``variance`` for arguments and other details.
  660. >>> stdev([1.5, 2.5, 2.5, 2.75, 3.25, 4.75])
  661. 1.0810874155219827
  662. """
  663. var = variance(data, xbar)
  664. try:
  665. return var.sqrt()
  666. except AttributeError:
  667. return math.sqrt(var)
  668. def pstdev(data, mu=None):
  669. """Return the square root of the population variance.
  670. See ``pvariance`` for arguments and other details.
  671. >>> pstdev([1.5, 2.5, 2.5, 2.75, 3.25, 4.75])
  672. 0.986893273527251
  673. """
  674. var = pvariance(data, mu)
  675. try:
  676. return var.sqrt()
  677. except AttributeError:
  678. return math.sqrt(var)
  679. # === Statistics for relations between two inputs ===
  680. # See https://en.wikipedia.org/wiki/Covariance
  681. # https://en.wikipedia.org/wiki/Pearson_correlation_coefficient
  682. # https://en.wikipedia.org/wiki/Simple_linear_regression
  683. def covariance(x, y, /):
  684. """Covariance
  685. Return the sample covariance of two inputs *x* and *y*. Covariance
  686. is a measure of the joint variability of two inputs.
  687. >>> x = [1, 2, 3, 4, 5, 6, 7, 8, 9]
  688. >>> y = [1, 2, 3, 1, 2, 3, 1, 2, 3]
  689. >>> covariance(x, y)
  690. 0.75
  691. >>> z = [9, 8, 7, 6, 5, 4, 3, 2, 1]
  692. >>> covariance(x, z)
  693. -7.5
  694. >>> covariance(z, x)
  695. -7.5
  696. """
  697. n = len(x)
  698. if len(y) != n:
  699. raise StatisticsError('covariance requires that both inputs have same number of data points')
  700. if n < 2:
  701. raise StatisticsError('covariance requires at least two data points')
  702. xbar = fsum(x) / n
  703. ybar = fsum(y) / n
  704. sxy = fsum((xi - xbar) * (yi - ybar) for xi, yi in zip(x, y))
  705. return sxy / (n - 1)
  706. def correlation(x, y, /):
  707. """Pearson's correlation coefficient
  708. Return the Pearson's correlation coefficient for two inputs. Pearson's
  709. correlation coefficient *r* takes values between -1 and +1. It measures the
  710. strength and direction of the linear relationship, where +1 means very
  711. strong, positive linear relationship, -1 very strong, negative linear
  712. relationship, and 0 no linear relationship.
  713. >>> x = [1, 2, 3, 4, 5, 6, 7, 8, 9]
  714. >>> y = [9, 8, 7, 6, 5, 4, 3, 2, 1]
  715. >>> correlation(x, x)
  716. 1.0
  717. >>> correlation(x, y)
  718. -1.0
  719. """
  720. n = len(x)
  721. if len(y) != n:
  722. raise StatisticsError('correlation requires that both inputs have same number of data points')
  723. if n < 2:
  724. raise StatisticsError('correlation requires at least two data points')
  725. xbar = fsum(x) / n
  726. ybar = fsum(y) / n
  727. sxy = fsum((xi - xbar) * (yi - ybar) for xi, yi in zip(x, y))
  728. sxx = fsum((xi - xbar) ** 2.0 for xi in x)
  729. syy = fsum((yi - ybar) ** 2.0 for yi in y)
  730. try:
  731. return sxy / sqrt(sxx * syy)
  732. except ZeroDivisionError:
  733. raise StatisticsError('at least one of the inputs is constant')
  734. LinearRegression = namedtuple('LinearRegression', ('slope', 'intercept'))
  735. def linear_regression(x, y, /):
  736. """Slope and intercept for simple linear regression.
  737. Return the slope and intercept of simple linear regression
  738. parameters estimated using ordinary least squares. Simple linear
  739. regression describes relationship between an independent variable
  740. *x* and a dependent variable *y* in terms of linear function:
  741. y = slope * x + intercept + noise
  742. where *slope* and *intercept* are the regression parameters that are
  743. estimated, and noise represents the variability of the data that was
  744. not explained by the linear regression (it is equal to the
  745. difference between predicted and actual values of the dependent
  746. variable).
  747. The parameters are returned as a named tuple.
  748. >>> x = [1, 2, 3, 4, 5]
  749. >>> noise = NormalDist().samples(5, seed=42)
  750. >>> y = [3 * x[i] + 2 + noise[i] for i in range(5)]
  751. >>> linear_regression(x, y) #doctest: +ELLIPSIS
  752. LinearRegression(slope=3.09078914170..., intercept=1.75684970486...)
  753. """
  754. n = len(x)
  755. if len(y) != n:
  756. raise StatisticsError('linear regression requires that both inputs have same number of data points')
  757. if n < 2:
  758. raise StatisticsError('linear regression requires at least two data points')
  759. xbar = fsum(x) / n
  760. ybar = fsum(y) / n
  761. sxy = fsum((xi - xbar) * (yi - ybar) for xi, yi in zip(x, y))
  762. sxx = fsum((xi - xbar) ** 2.0 for xi in x)
  763. try:
  764. slope = sxy / sxx # equivalent to: covariance(x, y) / variance(x)
  765. except ZeroDivisionError:
  766. raise StatisticsError('x is constant')
  767. intercept = ybar - slope * xbar
  768. return LinearRegression(slope=slope, intercept=intercept)
  769. ## Normal Distribution #####################################################
  770. def _normal_dist_inv_cdf(p, mu, sigma):
  771. # There is no closed-form solution to the inverse CDF for the normal
  772. # distribution, so we use a rational approximation instead:
  773. # Wichura, M.J. (1988). "Algorithm AS241: The Percentage Points of the
  774. # Normal Distribution". Applied Statistics. Blackwell Publishing. 37
  775. # (3): 477–484. doi:10.2307/2347330. JSTOR 2347330.
  776. q = p - 0.5
  777. if fabs(q) <= 0.425:
  778. r = 0.180625 - q * q
  779. # Hash sum: 55.88319_28806_14901_4439
  780. num = (((((((2.50908_09287_30122_6727e+3 * r +
  781. 3.34305_75583_58812_8105e+4) * r +
  782. 6.72657_70927_00870_0853e+4) * r +
  783. 4.59219_53931_54987_1457e+4) * r +
  784. 1.37316_93765_50946_1125e+4) * r +
  785. 1.97159_09503_06551_4427e+3) * r +
  786. 1.33141_66789_17843_7745e+2) * r +
  787. 3.38713_28727_96366_6080e+0) * q
  788. den = (((((((5.22649_52788_52854_5610e+3 * r +
  789. 2.87290_85735_72194_2674e+4) * r +
  790. 3.93078_95800_09271_0610e+4) * r +
  791. 2.12137_94301_58659_5867e+4) * r +
  792. 5.39419_60214_24751_1077e+3) * r +
  793. 6.87187_00749_20579_0830e+2) * r +
  794. 4.23133_30701_60091_1252e+1) * r +
  795. 1.0)
  796. x = num / den
  797. return mu + (x * sigma)
  798. r = p if q <= 0.0 else 1.0 - p
  799. r = sqrt(-log(r))
  800. if r <= 5.0:
  801. r = r - 1.6
  802. # Hash sum: 49.33206_50330_16102_89036
  803. num = (((((((7.74545_01427_83414_07640e-4 * r +
  804. 2.27238_44989_26918_45833e-2) * r +
  805. 2.41780_72517_74506_11770e-1) * r +
  806. 1.27045_82524_52368_38258e+0) * r +
  807. 3.64784_83247_63204_60504e+0) * r +
  808. 5.76949_72214_60691_40550e+0) * r +
  809. 4.63033_78461_56545_29590e+0) * r +
  810. 1.42343_71107_49683_57734e+0)
  811. den = (((((((1.05075_00716_44416_84324e-9 * r +
  812. 5.47593_80849_95344_94600e-4) * r +
  813. 1.51986_66563_61645_71966e-2) * r +
  814. 1.48103_97642_74800_74590e-1) * r +
  815. 6.89767_33498_51000_04550e-1) * r +
  816. 1.67638_48301_83803_84940e+0) * r +
  817. 2.05319_16266_37758_82187e+0) * r +
  818. 1.0)
  819. else:
  820. r = r - 5.0
  821. # Hash sum: 47.52583_31754_92896_71629
  822. num = (((((((2.01033_43992_92288_13265e-7 * r +
  823. 2.71155_55687_43487_57815e-5) * r +
  824. 1.24266_09473_88078_43860e-3) * r +
  825. 2.65321_89526_57612_30930e-2) * r +
  826. 2.96560_57182_85048_91230e-1) * r +
  827. 1.78482_65399_17291_33580e+0) * r +
  828. 5.46378_49111_64114_36990e+0) * r +
  829. 6.65790_46435_01103_77720e+0)
  830. den = (((((((2.04426_31033_89939_78564e-15 * r +
  831. 1.42151_17583_16445_88870e-7) * r +
  832. 1.84631_83175_10054_68180e-5) * r +
  833. 7.86869_13114_56132_59100e-4) * r +
  834. 1.48753_61290_85061_48525e-2) * r +
  835. 1.36929_88092_27358_05310e-1) * r +
  836. 5.99832_20655_58879_37690e-1) * r +
  837. 1.0)
  838. x = num / den
  839. if q < 0.0:
  840. x = -x
  841. return mu + (x * sigma)
  842. # If available, use C implementation
  843. try:
  844. from _statistics import _normal_dist_inv_cdf
  845. except ImportError:
  846. pass
  847. class NormalDist:
  848. "Normal distribution of a random variable"
  849. # https://en.wikipedia.org/wiki/Normal_distribution
  850. # https://en.wikipedia.org/wiki/Variance#Properties
  851. __slots__ = {
  852. '_mu': 'Arithmetic mean of a normal distribution',
  853. '_sigma': 'Standard deviation of a normal distribution',
  854. }
  855. def __init__(self, mu=0.0, sigma=1.0):
  856. "NormalDist where mu is the mean and sigma is the standard deviation."
  857. if sigma < 0.0:
  858. raise StatisticsError('sigma must be non-negative')
  859. self._mu = float(mu)
  860. self._sigma = float(sigma)
  861. @classmethod
  862. def from_samples(cls, data):
  863. "Make a normal distribution instance from sample data."
  864. if not isinstance(data, (list, tuple)):
  865. data = list(data)
  866. xbar = fmean(data)
  867. return cls(xbar, stdev(data, xbar))
  868. def samples(self, n, *, seed=None):
  869. "Generate *n* samples for a given mean and standard deviation."
  870. gauss = random.gauss if seed is None else random.Random(seed).gauss
  871. mu, sigma = self._mu, self._sigma
  872. return [gauss(mu, sigma) for i in range(n)]
  873. def pdf(self, x):
  874. "Probability density function. P(x <= X < x+dx) / dx"
  875. variance = self._sigma ** 2.0
  876. if not variance:
  877. raise StatisticsError('pdf() not defined when sigma is zero')
  878. return exp((x - self._mu)**2.0 / (-2.0*variance)) / sqrt(tau*variance)
  879. def cdf(self, x):
  880. "Cumulative distribution function. P(X <= x)"
  881. if not self._sigma:
  882. raise StatisticsError('cdf() not defined when sigma is zero')
  883. return 0.5 * (1.0 + erf((x - self._mu) / (self._sigma * sqrt(2.0))))
  884. def inv_cdf(self, p):
  885. """Inverse cumulative distribution function. x : P(X <= x) = p
  886. Finds the value of the random variable such that the probability of
  887. the variable being less than or equal to that value equals the given
  888. probability.
  889. This function is also called the percent point function or quantile
  890. function.
  891. """
  892. if p <= 0.0 or p >= 1.0:
  893. raise StatisticsError('p must be in the range 0.0 < p < 1.0')
  894. if self._sigma <= 0.0:
  895. raise StatisticsError('cdf() not defined when sigma at or below zero')
  896. return _normal_dist_inv_cdf(p, self._mu, self._sigma)
  897. def quantiles(self, n=4):
  898. """Divide into *n* continuous intervals with equal probability.
  899. Returns a list of (n - 1) cut points separating the intervals.
  900. Set *n* to 4 for quartiles (the default). Set *n* to 10 for deciles.
  901. Set *n* to 100 for percentiles which gives the 99 cuts points that
  902. separate the normal distribution in to 100 equal sized groups.
  903. """
  904. return [self.inv_cdf(i / n) for i in range(1, n)]
  905. def overlap(self, other):
  906. """Compute the overlapping coefficient (OVL) between two normal distributions.
  907. Measures the agreement between two normal probability distributions.
  908. Returns a value between 0.0 and 1.0 giving the overlapping area in
  909. the two underlying probability density functions.
  910. >>> N1 = NormalDist(2.4, 1.6)
  911. >>> N2 = NormalDist(3.2, 2.0)
  912. >>> N1.overlap(N2)
  913. 0.8035050657330205
  914. """
  915. # See: "The overlapping coefficient as a measure of agreement between
  916. # probability distributions and point estimation of the overlap of two
  917. # normal densities" -- Henry F. Inman and Edwin L. Bradley Jr
  918. # http://dx.doi.org/10.1080/03610928908830127
  919. if not isinstance(other, NormalDist):
  920. raise TypeError('Expected another NormalDist instance')
  921. X, Y = self, other
  922. if (Y._sigma, Y._mu) < (X._sigma, X._mu): # sort to assure commutativity
  923. X, Y = Y, X
  924. X_var, Y_var = X.variance, Y.variance
  925. if not X_var or not Y_var:
  926. raise StatisticsError('overlap() not defined when sigma is zero')
  927. dv = Y_var - X_var
  928. dm = fabs(Y._mu - X._mu)
  929. if not dv:
  930. return 1.0 - erf(dm / (2.0 * X._sigma * sqrt(2.0)))
  931. a = X._mu * Y_var - Y._mu * X_var
  932. b = X._sigma * Y._sigma * sqrt(dm**2.0 + dv * log(Y_var / X_var))
  933. x1 = (a + b) / dv
  934. x2 = (a - b) / dv
  935. return 1.0 - (fabs(Y.cdf(x1) - X.cdf(x1)) + fabs(Y.cdf(x2) - X.cdf(x2)))
  936. def zscore(self, x):
  937. """Compute the Standard Score. (x - mean) / stdev
  938. Describes *x* in terms of the number of standard deviations
  939. above or below the mean of the normal distribution.
  940. """
  941. # https://www.statisticshowto.com/probability-and-statistics/z-score/
  942. if not self._sigma:
  943. raise StatisticsError('zscore() not defined when sigma is zero')
  944. return (x - self._mu) / self._sigma
  945. @property
  946. def mean(self):
  947. "Arithmetic mean of the normal distribution."
  948. return self._mu
  949. @property
  950. def median(self):
  951. "Return the median of the normal distribution"
  952. return self._mu
  953. @property
  954. def mode(self):
  955. """Return the mode of the normal distribution
  956. The mode is the value x where which the probability density
  957. function (pdf) takes its maximum value.
  958. """
  959. return self._mu
  960. @property
  961. def stdev(self):
  962. "Standard deviation of the normal distribution."
  963. return self._sigma
  964. @property
  965. def variance(self):
  966. "Square of the standard deviation."
  967. return self._sigma ** 2.0
  968. def __add__(x1, x2):
  969. """Add a constant or another NormalDist instance.
  970. If *other* is a constant, translate mu by the constant,
  971. leaving sigma unchanged.
  972. If *other* is a NormalDist, add both the means and the variances.
  973. Mathematically, this works only if the two distributions are
  974. independent or if they are jointly normally distributed.
  975. """
  976. if isinstance(x2, NormalDist):
  977. return NormalDist(x1._mu + x2._mu, hypot(x1._sigma, x2._sigma))
  978. return NormalDist(x1._mu + x2, x1._sigma)
  979. def __sub__(x1, x2):
  980. """Subtract a constant or another NormalDist instance.
  981. If *other* is a constant, translate by the constant mu,
  982. leaving sigma unchanged.
  983. If *other* is a NormalDist, subtract the means and add the variances.
  984. Mathematically, this works only if the two distributions are
  985. independent or if they are jointly normally distributed.
  986. """
  987. if isinstance(x2, NormalDist):
  988. return NormalDist(x1._mu - x2._mu, hypot(x1._sigma, x2._sigma))
  989. return NormalDist(x1._mu - x2, x1._sigma)
  990. def __mul__(x1, x2):
  991. """Multiply both mu and sigma by a constant.
  992. Used for rescaling, perhaps to change measurement units.
  993. Sigma is scaled with the absolute value of the constant.
  994. """
  995. return NormalDist(x1._mu * x2, x1._sigma * fabs(x2))
  996. def __truediv__(x1, x2):
  997. """Divide both mu and sigma by a constant.
  998. Used for rescaling, perhaps to change measurement units.
  999. Sigma is scaled with the absolute value of the constant.
  1000. """
  1001. return NormalDist(x1._mu / x2, x1._sigma / fabs(x2))
  1002. def __pos__(x1):
  1003. "Return a copy of the instance."
  1004. return NormalDist(x1._mu, x1._sigma)
  1005. def __neg__(x1):
  1006. "Negates mu while keeping sigma the same."
  1007. return NormalDist(-x1._mu, x1._sigma)
  1008. __radd__ = __add__
  1009. def __rsub__(x1, x2):
  1010. "Subtract a NormalDist from a constant or another NormalDist."
  1011. return -(x1 - x2)
  1012. __rmul__ = __mul__
  1013. def __eq__(x1, x2):
  1014. "Two NormalDist objects are equal if their mu and sigma are both equal."
  1015. if not isinstance(x2, NormalDist):
  1016. return NotImplemented
  1017. return x1._mu == x2._mu and x1._sigma == x2._sigma
  1018. def __hash__(self):
  1019. "NormalDist objects hash equal if their mu and sigma are both equal."
  1020. return hash((self._mu, self._sigma))
  1021. def __repr__(self):
  1022. return f'{type(self).__name__}(mu={self._mu!r}, sigma={self._sigma!r})'