Statistics, read from an analyst's source
In 1973 the University of California at Berkeley had 12,763 applicants to its graduate programmes. Of the 8,442 men who applied, 44 percent were admitted. Of the 4,321 women, 35 percent were. A gap of nine points is too large to be likely chance. Bickel, Hammel and O'Connell published the analysis in Science in 1975: women had tended to apply to the more competitive departments, which admitted a lower share, and men to the less competitive ones, and once the data were corrected for department they showed a small but statistically significant bias in favour of women. The same applicants, counted two ways, gave opposite verdicts. Nothing in the arithmetic was wrong; the question had been put to the wrong total.
The statistics library for DKE, the Deterministic Knowledge Engine, keeps that lesson and a few others in DKE Python, and this essay teaches the subject by reading it. You record three kinds of fact. A study is written as its strata, each with the successes and the trials of two arms, and the library's rules conclude each arm's rate in each stratum and in the pool, and whether the pooled comparison reverses what every stratum shows. A screening test is written as how common a condition is and how often the test is right about those who have it and those who do not, and the rules conclude what a positive result is worth. A significance test is written as its p-value and the level it is held to, and the rules conclude whether it clears. Each answer names the facts it rests on, and correcting a fact withdraws the answers it supported. The second half of the library calculates: means and variances, a line through points, the normal, binomial and Poisson distributions, the birthday problem, Bayes' rule, Welch's t, Pearson's chi-square and a generator of pseudo-random numbers that a run repeats exactly. Every program below runs as shown once the library's source file has been compiled into a store, and every output is copied from a run.
1. Strata, and the reversal
The kidney-stone data come from a paper of 1986 in the British Medical Journal by Charig and his colleagues, which compared open surgical procedures, treatment A, with closed ones, treatment B. For small stones A succeeded in 81 of 87 cases and B in 234 of 270. For large stones A succeeded in 192 of 263 and B in 55 of 80. A did better in each stratum, and B did better overall, 289 of 350 against 273 of 350, because doctors gave the large stones the better treatment, A, more often, so the totals are skewed. Karl Pearson in 1899 and Udny Yule in 1903 had noticed effects of this kind, Edward Simpson described it in a technical paper in 1951, and Colin Blyth gave it the name Simpson's paradox in 1972.
The library writes a stratum as four counts and lets rules fold them into the study and compare:
@rule
def paradox_b_ahead_in_the_pool():
for t in Study:
if t.pooled_b_ahead == True and none(s in Stratum: s.study == t and s.a_wins != True):
t.paradox = TrueThe rule says that a study is a paradox when B is ahead in the pool and there is no stratum in which A is not ahead. Its twin does the same with the arms exchanged. Here is the stone study:
import statistics
statistics.stratum("stones", "small", 81, 87, 234, 270)
statistics.stratum("stones", "large", 192, 263, 55, 80)
print(statistics.paradoxes(), statistics.pooled("stones"))
a = current(statistics.Stratum.stones_small.a_rate)
b = current(statistics.Stratum.stones_small.b_rate)
print(a.value, b.value)
a = current(statistics.Stratum.stones_large.a_rate)
b = current(statistics.Stratum.stones_large.b_rate)
print(a.value, b.value)['stones'] [0.78, 0.8257142857142857]
0.9310344827586207 0.8666666666666667
0.7300380228136882 0.6875The pooled rates are 78 percent for A and 82.6 for B, and the study is a paradox. A leads 93.1 to 86.7 among small stones and 73.0 to 68.8 among large ones. The library does not say which total is the right one, and could not: for Berkeley's departments and these stones the strata may be the better question, and there are cases where the pool is. It states what the numbers do, and a reader who asks why gets the counts that make the reversal:
import statistics
w = why(statistics.Study.stones.paradox)
for c in w:
print(str(c.path) + " = " + str(c.value))statistics.Study.stones.paradox = True
statistics.Study.stones.pooled_b_ahead = True
statistics.Study.stones.b_rate = 0.8257142857142857
statistics.Study.stones.a_rate = 0.78
statistics.Stratum.stones_small.a_rate = 0.9310344827586207
statistics.Stratum.stones_small.b_wins = False
statistics.Stratum.stones_large.a_rate = 0.7300380228136882
statistics.Stratum.stones_large.b_wins = FalseThe answer runs to the study's folded totals, 273 and 289 successes in 350 trials each, and down to the four counts. Now suppose a clerk finds that 234 was a slip for 255, so that B succeeded in 255 of its 270 small-stone cases. The paradox cannot survive that, because B now leads in a stratum, and the rule's answer to why_not says which premise failed:
import statistics
statistics.stratum("recheck", "small", 81, 87, 234, 270)
statistics.stratum("recheck", "large", 192, 263, 55, 80)
print(statistics.paradoxes())
statistics.stratum("recheck", "small", 81, 87, 255, 270)
print(statistics.paradoxes())
reason = why_not(statistics.Study.recheck.paradox)
for line in str(reason).split("\n"):
print(line)['recheck', 'stones']
['stones']
rule `statistics.paradox_a_ahead_in_the_pool` needs statistics.Study[recheck].pooled_a_ahead == True, and statistics.Study[recheck].pooled_a_ahead is False (bool)
rule `statistics.paradox_b_ahead_in_the_pool` needs statistics.Study[recheck].paradox_b_ahead_in_the_pool__none1 to be absent, and it is presentThe study left the list when the count changed, and the two rules that could have concluded the paradox each name the thing that now stands in their way. The 255 is invented for the example; the stone data of the first program is the published one.
2. Describing data, and Anscombe's four
The mean, the median and the sample variance are the first things anyone computes, and the library computes the variance by Welford's single pass, which updates a running mean and a running sum of squared deviations and so does not subtract two large numbers. The line through points is the least-squares line, which Adrien-Marie Legendre published in 1805 and Carl Friedrich Gauss in 1809, claiming to have used it since 1795, which led to a priority dispute. The astronomical problem that gave the method its use was Ceres: Giuseppe Piazzi found the asteroid on 1 January 1801, followed it for forty days until the Sun's glare hid it, and Gauss's calculations were the ones that let Franz Xaver von Zach find it again.
Francis Anscombe made four small data sets in 1973 for a paper in The American Statistician called "Graphs in Statistical Analysis", to counter the impression, as he put it, that numerical calculations are exact but graphs are rough. All four have 11 points, the same mean of x, 9, the same variance of x, 11, the same mean of y, 7.50, the same line y = 3.00 + 0.500x and the same correlation, 0.816. Here are the first three statistics and the line of each set:
import statistics
x = [10.0, 8.0, 13.0, 9.0, 11.0, 14.0, 6.0, 4.0, 12.0, 7.0, 5.0]
x4 = [8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 19.0, 8.0, 8.0, 8.0]
y1 = [8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68]
y2 = [9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74]
y3 = [7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73]
y4 = [6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89]
print(statistics.mean(x), statistics.variance(x), statistics.mean(y1), statistics.variance(y1))
print(statistics.regression(x, y1))
print(statistics.regression(x, y2))
print(statistics.regression(x, y3))
print(statistics.regression(x4, y4))
print(statistics.median(y1), statistics.median(y3), statistics.std_dev(y1))9.0 10.999999999999996 7.500909090909093 4.127269090909092
[3.0000909090909103, 0.5000909090909091, 0.81642051634484]
[3.000909090909091, 0.5, 0.8162365060002428]
[3.0024545454545466, 0.49972727272727274, 0.8162867394895981]
[3.0017272727272726, 0.49990909090909086, 0.8165214368885028]
7.58 7.11 2.031568135925815Each regression line is an intercept, a slope and a correlation, and each is Anscombe's 3.00, 0.500 and 0.816 to the printed digits. The variance of x is 11 up to the last digit the machine keeps, and the variance of y for the first set, 4.127, sits inside the 4.125 plus or minus 0.003 of Anscombe's tables. Set I is a noisy line, set II lies on a curve, set III is a line with one point far off it, and set IV is ten points at the same x with a single point at 19 that fixes the whole slope. The summary numbers do not tell them apart, and the only figures above that differ visibly are the medians of the y values of the first and third sets, 7.58 and 7.11. The origin of Anscombe's sets is unknown; in 2017 Matejka and Fitzmaurice showed a way, by simulated annealing, to make others with identical statistics and arbitrary pictures, one of them, the Datasaurus dozen, with points that trace a dinosaur.
3. What a positive result is worth
A screening test has a sensitivity, the chance of a positive result for someone who has the condition, and a specificity, the chance of a negative for someone who does not. Neither says what a positive result means, which depends on how common the condition is. The posthumous essay of Thomas Bayes, edited by Richard Price and read to the Royal Society on 23 December 1763, is where the rule for turning one into the other starts; Laplace reached the same rule independently and extended it in 1774. The library writes the rule as two computed fields:
ppv: float = prevalence * sensitivity / (prevalence * sensitivity + (1.0 - prevalence) * (1.0 - specificity))
mostly_false: bool = ppv < 0.5The first is the positive predictive value, the chance that a positive result is right. Take three tests. The first is right nine times in ten about both kinds of person, and the condition is in one person in a hundred. The second is right 99 times in 100 about both, and the condition is in 3 in 1,000, which is the figure of a standard example of a drug test. The third is the first test again where the condition is in three people in ten:
import statistics
statistics.screen("rare", 0.01, 0.9, 0.9)
statistics.screen("drug", 0.003, 0.99, 0.99)
statistics.screen("common", 0.3, 0.9, 0.9)
print(statistics.predictive_values("rare"))
print(statistics.predictive_values("drug"))
print(statistics.predictive_values("common"))
print(statistics.bayes(0.01, 0.9, 0.1))[0.08333333333333336, 0.9988789237668162]
[0.22952086553323015, 0.9999696067108382]
[0.7941176470588237, 0.9545454545454545]
0.08333333333333333The first number of each pair is the positive predictive value and the second the negative. A positive result of the first test is right one time in twelve, and a negative result is right 99.9 percent of the time. The test of the drug, called 99 percent accurate, gives a positive that is right 23 percent of the time, and most of its positives are false. The same test where the condition is common is right 79 percent of the time. The last line is the first calculation again by bayes(), which takes a prior and the two likelihoods. This is the lesson of Section 1 with the roles changed. The test is the same and the verdict about a person is not, because the group the person belongs to is not.
4. A p-value against its level
A significance test ends with a p-value, the chance of a result at least as extreme as the one observed if nothing is going on. Karl Pearson introduced the idea in his chi-square test of 1900. The level at which a p-value is held to be small, 0.05, was proposed by Ronald Fisher in 1925 in Statistical Methods for Research Workers, as a two-tailed limit corresponding to about two standard deviations on a normal curve. The rule the library keeps is only the comparison:
significant: bool = p < alphaimport statistics
statistics.significance("t1", 0.049, 0.05)
statistics.significance("t2", 0.051, 0.05)
print(statistics.significant_tests())['t1']A result at 0.049 clears and one at 0.051 does not, though no reasonable person would say that the data of one differ from the other. The American Statistical Association said as much in 2016, in a statement by Wasserstein and Lazar, which says that a p-value does not measure the probability that the hypothesis studied is true, nor the size of an effect or the importance of a result. The library holds the level beside the p-value as a fact in the store, so that a decision names the level it was made at and a changed level changes the decision, and it does nothing to make a threshold wise.
5. The curves
The normal curve gives the chance of a value below a given number of standard deviations from the mean. The library writes the distribution function through the complementary error function of the mathematical functions, 0.5 * math.erfc(-z / math.sqrt(2.0)), which keeps the digits in the tails. About 68 percent of a normal population lies within one standard deviation of the mean, about 95 percent within two and 99.7 within three, and 1.96 is the value with 2.5 percent beyond it:
import statistics
print(statistics.normal_cdf(1.96), statistics.normal_quantile(0.975))
print(statistics.normal_cdf(1.0) - statistics.normal_cdf(-1.0), statistics.normal_cdf(2.0) - statistics.normal_cdf(-2.0), statistics.normal_cdf(3.0) - statistics.normal_cdf(-3.0))
print(statistics.normal_cdf(-5.0), statistics.normal_pdf(0.0, 0.0, 1.0))
print(statistics.z_score(130.0, 100.0, 15.0))0.9750021048517795 1.9599639845400532
0.6826894921370859 0.9544997361036416 0.9973002039367398
2.8665157187919455e-07 0.3989422804014327
2.0The quantile is found by bisection, 200 halvings of an interval from minus 40 to 40, and agrees with the distribution function to the last digit. Five standard deviations below the mean leaves 2.87 in ten million.
The binomial and Poisson distributions count. The binomial gives the chance of exactly k successes in n independent trials, and the library finds it as the coefficient times the two powers while those are ordinary numbers, and past them by Catherine Loader's saddle-point form of 2000, so it has the chance for any n: 550 heads in 1,100 tosses come up about one time in 42. Siméon Denis Poisson published the distribution that bears his name in 1837, in a book on the probability of judgments in criminal and civil matters. Ladislaus Bortkiewicz showed in 1898, in Das Gesetz der kleinen Zahlen, that deaths of Prussian soldiers from horse kicks were well modelled by it, and R. D. Clarke used it in 1946 on the flying bombs that fell on London. The birthday problem is a binomial in disguise: the problem is attributed to Harold Davenport about 1927 and was first published by Richard von Mises in 1939. The chance that two of 23 people share a birthday, with 365 equally likely days and no leap years, is just over half:
import statistics
print(statistics.binom_pmf(10, 3, 0.5), statistics.binom_cdf(10, 3, 0.5), statistics.binom_pmf(100, 50, 0.5))
print(statistics.poisson_pmf(0, 0.61), statistics.poisson_pmf(1, 0.61), statistics.poisson_pmf(2, 0.61), statistics.poisson_pmf(3, 0.61))
print(statistics.birthday(23), statistics.birthday(57), statistics.birthday(70))0.1171875 0.171875 0.07958923738717875
0.5433508690744998 0.3314440301354449 0.10109042919131069 0.020555053935566506
0.5072972343239855 0.9901224593411699 0.9991595759651571Three heads in ten fair tosses is 120 in 1,024, which is 0.1171875, and three or fewer is 0.171875. Fifty heads in a hundred is 8 percent. The Poisson numbers are the chances of 0, 1, 2 and 3 events at a mean of 0.61 a period, and the birthday chance is 50.7 percent at 23 people, 99.0 at 57 and 99.9 at 70. The mean of 0.61 is chosen for the example; the figures show how fast the tail dies away.
6. Two samples, and a table
Student's t-test is named for a pen name. William Sealy Gosset worked at the Guinness brewery in Dublin and wanted an economical way to monitor the quality of stout from small samples. He published the result in Biometrika in 1908 as "The Probable Error of a Mean", under the name Student because, the usual account says, his employer preferred its staff to publish scientific papers under pen names. It was Ronald Fisher's work that made the distribution well known. The statistic for two samples of unequal variance, which the library calls Welch's t, is the difference of the means over the square root of the sum of each variance over its size:
import statistics
a = [5.0, 7.0, 6.0, 9.0, 8.0, 7.0, 6.0, 8.0]
b = [8.0, 9.0, 7.0, 10.0, 9.0, 8.0, 9.0, 10.0]
print(statistics.mean(a), statistics.mean(b), statistics.welch_t(a, b))7.0 8.75 -2.965614910077132The eight values of each list are invented. Turning the statistic into a p-value needs the t distribution and its degrees of freedom, and the library does not carry them: the distribution the library does carry for a test is Pearson's chi-square, which compares counts with the counts a theory predicts, as the sum over the cells of the squared difference over the expected count. The pea crosses of Gregor Mendel gave, in the second generation from round and wrinkled parents, 5,474 round seeds and 1,850 wrinkled, a ratio of 2.96 to 1, where his theory said 3 to 1. With 7,324 seeds the expected counts are 5,493 and 1,831:
import statistics
q = statistics.chi_square([5474.0, 1850.0], [5493.0, 1831.0])
print(q)
print(statistics.chi_square_p_df1(q))0.2628800291279811
0.6081484045220742The statistic is 0.263 on one degree of freedom, and the chance of a departure from 3 to 1 at least that large is 61 percent. The data are no worse than the theory expects, and Fisher, who reanalysed Mendel's second-generation results in 1936, found the opposite complaint: that the ratios were implausibly close to expectation. Later authors disputed his conclusion, and a book of 2008, Ending the Mendel-Fisher controversy, found no basis to say Mendel had fabricated his results. The library does not decide that. It gives the statistic and, for one degree of freedom and for two, the tail.
7. A generator that repeats
The last script in the library is a generator of pseudo-random numbers, because the language has none: a program that gives a different answer on each run cannot be checked. Derrick Henry Lehmer described the multiplicative congruential generator, in which each number is the previous one times a constant, modulo a prime, in a paper of 1949 and a journal version of 1951. Park and Miller in 1988 recommended a multiplier of 16,807 with the modulus 2^31 − 1, and in 1993 Park, Miller and Stockmeyer suggested 48,271, which the library uses and which C++11 calls minstd_rand. The product stays under 2^63, so it never overflows. From a seed of 1:
import statistics
s = 1
for i in range(4):
s = statistics.lehmer_next(s)
print(s, statistics.lehmer_uniform(s))
print(statistics.monte_carlo_pi(1000, 1), statistics.lehmer_next(0))48271 2.2477936010098986e-05
182605794 0.08503244914348818
1291394886 0.6013526053174179
1914720637 0.8916112770753034
3.148 0The first state is the multiplier itself, and a seed of 0 stays 0 for ever, which the last line shows beside an estimate of pi from a thousand points. That estimate is the Monte Carlo method in its smallest form. Stanisław Ulam conceived the method in 1946, recovering from an illness and playing solitaire, when he wondered whether playing many games and counting the wins would be quicker than working out the odds. Nicholas Metropolis proposed the name, after the casino at Monte Carlo, and in the spring of 1948 von Neumann, Metropolis and others programmed the ENIAC to run the first fully automated calculations of the kind. An early variant addressed Buffon's needle problem. A quadrant of a circle inscribed in a square has pi over 4 of its area, so the share of random points of the square that fall inside it, times four, estimates pi. The library draws pairs of points from the generator. With 10,000 of them from seed 1:
import statistics
print(statistics.monte_carlo_pi(10000, 1))3.1344The answer is 3.1344, where pi is 3.14159; the error of such an estimate falls with the square root of the number of points, so a hundred times as many would be needed for another digit. The run repeats: the same seed and the same count give the same number on any machine. The service bounds the memory one run may hold, and 20,000 points exceed it, so the examples stop at 10,000.
8. What the library leaves out
Nothing in the library chooses how a study should be divided into strata, as the stone data show. pooled() and predictive_values() refuse where the figure they answer with was never written, the sample variance needs two values, the mean one, a regression two different x values, and a list that is too small divides by zero. The distribution functions are the normal, the binomial and the Poisson, the tail of the chi-square is for one and two degrees of freedom, Welch's t has no p-value, and fitting beyond a straight line is not here. The generator is a plain congruential one with a short, known structure, and good for lessons and not for secrets. What the library adds is a written reason: the strata under a reversal, the prior under a positive result, the level under a decision.
Closing
Statistics is the discipline of saying what a number means. A stout brewer wanted to know how far to trust a small sample, an astronomer wanted a lost asteroid, a statistician wanted four sets of eleven points that no formula could tell apart, and each left a procedure. The library keeps some of them in a form a reader can ask questions of. It tells a reversal from a trend and names the counts that make it, says what a positive result is worth for the condition it is asked about, and keeps the level beside the p-value that was held to it. And because each answer carries its premises, one wrong count is corrected where it was made, and what stood on it falls with it.