We use cookies to enhance your experience on the site

Statistics - the mathematics that tells a real track from a coincidence

Last day of the season, @name. Your notebook holds animal counts from ten waterholes, and the reserve office is asking for one number: how many animals there are per waterhole. Ten observations have to shrink down to a single figure in the report, and whether that report tells the truth depends on which figure you pick.

The second question is worse. You weighed lions in the north of the reserve and in the south, eight of each. The north came out twenty one kilograms heavier. Are northern lions genuinely bigger, or did you simply happen to catch eight heavier animals? Your eye cannot settle that. For three seasons I was certain that the cheetahs of the western valley ran faster, and only working it out on paper showed me that the difference sat comfortably inside the ordinary spread of my own measurements.

Statistics is the toolkit for exactly those two questions: how to summarize a pile of numbers into a few meaningful values, and whether what you are looking at is not simply chance. This lesson walks through both, from the mean in your notebook to the straight line that predicts next season.

Where

stats

comes from

Two libraries do the counting. You already know the first one: NumPy keeps measurements in arrays and computes on them fast, and the convention is to shorten it to

np

. The second is SciPy, a large scientific package built on top of NumPy, holding integrals, optimization, signal processing and, the part that interests us today, a statistics module. That module is called

stats

and it lives inside the

scipy

package, exactly the way

pyplot

lives inside Matplotlib.

Since

stats

sits inside, you reach it with the

from ... import ...

form, and the order of those four words is fixed by Python grammar. First

from

, meaning "where we are taking it from". Straight after it the name of the place, the package

scipy

. Then

import

, meaning "pull it into my program". Last the name of the thing being pulled in, the module

stats

. Together:

from

scipy

import

stats

. The reverse order, tempting for anyone who has written JavaScript before, ends in an instant

SyntaxError

- Python knows no "import something from somewhere" spelling and gives up before the first line ever runs.

There is a second correct route:

import scipy.stats as stats

. It does the same job, but it is longer and needs its own alias through

as

. I recommend

from scipy import stats

, because that is how the SciPy documentation writes it and how practically every example online writes it, so other people's code will catch you off guard less often.

1import numpy as np
2from scipy import stats

After those two lines nothing has been computed and nothing has been printed. Two labels are hanging in memory:

np

leads to arrays and basic mathematics,

stats

leads to tests, distributions and correlations. Everything else in this lesson is nothing more than calls made on those two names.

One number instead of ten

Let us start with the reserve office and its question. The most obvious way to summarize a pile of numbers is the arithmetic mean: the sum of all the measurements divided by how many there are. It is computed by the

mean

function from NumPy, which takes an array and hands back a single value. I pack the measurements themselves into an array with

np.array

, the function you met in the NumPy lesson, and I name that array

populations

- it is the same name you will be using in the practical exercise, so it is worth remembering.

1populations = np.array([120, 450, 85, 200, 330, 150, 280, 95, 400, 175])
2
3print(np.mean(populations))
4# 228.5

That gives 228.5 animals per waterhole, and the number is computed flawlessly. The trouble is that it does not describe the reserve. Check it yourself: six waterholes out of ten lie below 228.5 and only four lie above it. The mean landed higher than most of the observations, because two waterholes with crowds of 400 and 450 animals dragged it upwards. The mean takes every single number equally seriously, so every extreme has a real pull on it.

A middle that cannot be bribed

If the mean can be dragged around, we need a measure that looks not at how big the numbers are but at where they stand in the queue. That measure is the median: line the measurements up from smallest to largest and take the one in the middle. When the count is even, as it is with our ten, the median is the average of the two middle values. It is computed by the

median

function from NumPy, called in exactly the same way as

mean

1print(np.median(populations))
2# 187.5

The sorted reserve looks like this: 85, 95, 120, 150, 175, and then 200, 280, 330, 400, 450. The two middle values are 175 and 200, so the median comes out at 187.5. The gap against the mean is forty one animals - and that gap is information about the shape of the data, not a defect. When the mean sits clearly above the median, it tells you that a few very large values are sitting somewhere in the set.

When a wildebeest herd walks through the waterhole

Now let us do what nature does once a season. At the waterhole that usually held 450 animals a migrating herd of wildebeest passes through and you count 4200. One measurement changes, the other nine stay exactly as they were. That is what an outlier looks like: an observation that is not a mistake at all, only a rare event, and that can still blow your statistics apart. Let us compute both middles again.

1with_herd = np.array([120, 4200, 85, 200, 330, 150, 280, 95, 400, 175])
2
3print(np.mean(with_herd))
4# 603.5
5print(np.median(with_herd))
6# 187.5

Stop at those two lines, because this is the most important moment of the lesson. The mean jumped from 228.5 to 603.5, so it nearly tripled on the strength of a single observation. The median did not budge by so much as a hair: 187.5 before the herd and 187.5 after it. The reason is simple - all the median cares about is that 4200 sits somewhere at the right hand end of the queue. Whether that value is 4200, or 450, or four million, the middle pair stays the same. That is why the median is the measure of central tendency that is resistant to outliers, and it is the one that goes into the report whenever the data can turn strange.

Three confusions come back at this exact point, so let us straighten them out right away. The arithmetic mean is a measure of central tendency, but resistant it certainly is not - you have just watched one measurement shift it by 375 units. Variance and standard deviation are not measures of central tendency at all: they do not say where the middle lies, only how widely the data is spread out. Worse still, both are built on squared distances from the mean, so they react to outliers even more violently than the mean itself does. We will compute that in a moment and you will see the scale of it.

The third middle: the mode

There are three measures of central tendency, so let us complete the set. The mode is the value that occurs most often. It earns its keep wherever measurements repeat: the number of patrols per day, the number of cubs in a litter, the species most frequently seen at the waterhole. It is computed by the

mode

function from the

stats

module, which hands back an object with two fields:

mode

is the value itself,

count

is how many times it occurred.

1patrols = np.array([3, 5, 4, 5, 2, 5, 4])
2
3result = stats.mode(patrols)
4print(result.mode)
5# 5
6print(result.count)
7# 3

Five patrols a day was your most frequent workload, and it happened three times. Notice that we did not compute the mode on the

populations

array - and that is not an oversight. Over there all ten numbers are different, so each one occurs exactly once, and the function would politely return the smallest of them with a count of one. That answer is formally correct and completely useless. Compute the mode where something actually repeats.

How widely the data is spread

Two reserves can share an identical mean and have nothing else in common. In the first one every waterhole gathers around two hundred animals; in the second, half of them stand empty while the other half are bursting at the seams. Measures of dispersion exist to describe that difference. The most important of them is the standard deviation: roughly the typical distance of a single measurement from the mean. It is computed by

np.std

. Its close relative is the variance from

np.var

, which is exactly the standard deviation squared.

1print(np.std(populations))
2# 123.06603918222119
3print(np.var(populations))
4# 15145.25

Against a mean of 228.5 a typical waterhole sits about 123 animals away from it - the spread is enormous, and that is the honest portrait of this reserve. The variance of 15145.25 carries the same information, but in squared units, in "square animals", which nobody can picture. That is why reports quote the standard deviation and leave the variance to the formulas working away inside the calculation. Remember as well that

np.std

treats your data as the whole population by default; when your ten waterholes are a sample drawn from a larger reserve, you add the argument

ddof=1

I promised to show you how those two measures react to the wildebeest herd. The standard deviation of

with_herd

comes out at 1202.87 instead of 123.07 - it grew nearly tenfold on one observation, while the median did not move at all. That is the proof that dispersion is the most sensitive part of descriptive statistics, and no kind of shelter from outliers.

Where half the herd lives

There is a third way of looking at the same data, as resistant as the median but richer. A percentile is the value below which a given share of the observations falls: the 25th percentile cuts off the lowest quarter of the measurements, the 95th tells you what almost everything stays under. Three percentiles have names of their own and are called quartiles: Q1 is the 25th percentile, Q2 is the 50th (which is simply the median), Q3 is the 75th. The gap between Q3 and Q1 is called the interquartile range, IQR for short, and it describes the band in which the middle half of the data lives. There is a ready made

iqr

function in the

stats

module for it. For contrast,

np.ptp

gives the plain range, the difference between the largest and the smallest measurement.

1print(np.percentile(populations, 25))
2# 127.5
3print(np.percentile(populations, 50))
4# 187.5
5print(np.percentile(populations, 75))
6# 317.5
7print(stats.iqr(populations))
8# 190.0
9print(np.ptp(populations))
10# 365

The middle half of the waterholes sits between 127.5 and 317.5 animals, and Q2 came out identical to the median from earlier, because those are two names for the same thing. The interesting part happens when you run those same calls on the

with_herd

array. The range shoots up from 365 to 4115, because it depends on nothing but the extreme values. The IQR stays at exactly 190.0, because Q1 and Q3 still land on the same measurements as before. Hence the rule: when rare events may be hiding in your data, describe it with the median and the IQR rather than the mean and the range.

Shape, or which way the distribution leans

You now know where the middle is and how wide the spread is. One piece is missing: whether the data is laid out symmetrically around that middle or leans to one side. That is measured by skewness, computed by the

skew

function. A value near zero means a symmetric distribution. A positive value means a long tail on the right, that is, a handful of very high measurements. A negative value means a long tail on the left.

1print(stats.skew(populations))
2# 0.5269796953903054
3print(stats.skew(with_herd))
4# 2.6335940542443765

The ordinary reserve has a skewness of 0.53, a mild lean to the right - and that is precisely why the mean of 228.5 sat above the median of 187.5. Once the wildebeest herd is added the skewness leaps to 2.63, which is already a shout of "something extreme is sitting in this data". Skewness is therefore a fast detector of outliers, before you even start hunting them down one by one. There is a sister measure too, kurtosis from

stats.kurtosis

, describing how heavy the tails of a distribution are; I mention it by name only, because reading it sensibly takes a lesson of its own and none of today's conclusions rest on it.

The normal distribution, or where practice data comes from

Measurements in nature very often fall into the same shape: most of the results cluster near the middle, and the further you go from it the rarer things get. That shape is called the normal distribution and it looks like a bell. Two numbers describe it: the center of the bell, called

loc

in SciPy and NumPy, and the width of the bell, which is the standard deviation, called

scale

You will not always have real data at hand, so being able to draw some at random is useful. That is the job of

np.random.normal

, which you give

loc

scale

and

size

, the number of measurements. Random drawing has one drawback: every run produces something different.

np.random.seed

fixes that by setting the generator's starting point - once the seed is set, the same program always spits out the same numbers, which is why the results in this lesson match digit for digit what you will see on your own machine. Let us draw the heights of two hundred giraffes in centimeters, and shorten the results themselves with the built in

round

function, which takes a number and a count of decimal places - without it we would be reading a dozen digits of expansion.

1np.random.seed(7)
2heights = np.random.normal(loc=500, scale=40, size=200)
3
4print(round(np.mean(heights), 2))
5# 498.88
6print(round(np.std(heights), 2))
7# 39.62

We asked for a mean of 500 and a deviation of 40, and we got 498.88 and 39.62. That small discrepancy is not an error - it is ordinary sampling noise, the same noise that keeps two teams counting the same reserve from ever coming back with identical numbers. With two hundred giraffes it is slight; with ten it would be far larger. The same mechanism has siblings for other shapes:

np.random.poisson

draws counts of events in a fixed span of time, for instance how many lions you will see within an hour, and

np.random.exponential

draws the gaps in time between events.

Is the data really normal

Plenty of statistical methods assume that the data comes from a normal distribution, so it is only proper to check that rather than guess it off a chart. The Shapiro-Wilk test does the checking, and you call it as

stats.shapiro

. The spelling has six parts in a fixed order: the module name

stats

, a dot, the function name

shapiro

, an opening bracket, the name of the data and a closing bracket. The test returns two values at once, so we catch them into two variables written on the left of the equals sign and separated by a comma: the first is the test statistic, the second is the p-value, the probability of seeing data like this assuming the distribution really is normal. By convention a p-value above 0.05 means "no grounds to reject normality", and below 0.05 means "this data is probably not normal".

1stat, p_value = stats.shapiro(heights)
2
3print(round(stat, 4))
4# 0.9912
5print(round(p_value, 4))
6# 0.2614

The p-value is 0.2614, comfortably above the 0.05 threshold, so the giraffes pass the test - which is what we expected, since we drew them from a normal distribution ourselves. Notice the careful wording: the test has not proved normality, it has merely failed to find a reason to reject it. That is a subtle but important distinction running through the whole of statistics.

Now let us see what the result looks like when the data genuinely is crooked. Run the same test on the reserve with the wildebeest herd, and this time check the threshold directly with the comparison

p_value < 0.05

, which hands back

True

False

on its own.

1stat, p_value = stats.shapiro(with_herd)
2
3print(round(stat, 4))
4# 0.44
5print(p_value < 0.05)
6# True

The statistic dropped from 0.99 to 0.44 and the p-value fell below one in a million, so the condition

p_value < 0.05

is true and we reject normality. No surprise there - this is the same data that skewness flagged at 2.63. Two different tools, one diagnosis: a single observation is pulling the whole set out of shape.

Is the difference between two herds real

Back to the second question from the start of the lesson. You weighed eight lions in each of two parts of the reserve and you want to know whether the north really does grow heavier animals. Let us begin by writing the measurements down and computing both means with a tool you already know.

1north = np.array([190, 205, 178, 212, 195, 188, 201, 199])
2south = np.array([172, 168, 185, 176, 180, 165, 178, 174])
3
4print(np.mean(north))
5# 196.0
6print(np.mean(south))
7# 174.75

The difference is 21.25 kilograms in the north's favor. On its own it still means nothing, because there are only eight measurements per group and lions can differ from one another by a dozen kilograms or more. The question is this: could a difference that size pop out of pure chance, if the two herds were in reality identical?

The answer comes from the t-test, available as

stats.ttest_ind

for two independent groups. The function takes two arrays and returns the t statistic together with a p-value carrying the same meaning as in the normality test: the probability of seeing a difference this large if there were no real difference at all.

1t_stat, p_value = stats.ttest_ind(north, south)
2
3print(round(t_stat, 4))
4# 4.8177
5print(round(p_value, 6))
6# 0.000273

The p-value is 0.000273, that is 0.03 percent. If the two herds were the same, a difference like this would turn up less often than once in three thousand attempts - so we call the difference real and we can put it in the report. Notice that the means themselves did not change; the test recomputes nothing in the data, it only attaches a verdict on how much the data can be trusted.

Since the t-test is a tool people often mix up, let us state plainly what it is not for. It does not normalize data - rescaling measurements is a separate operation in which you subtract the mean and divide by the standard deviation, and you will meet it in feature engineering. It does not calculate correlation -

stats.pearsonr

and

stats.spearmanr

are there for that, and you will meet them in a moment. It does not generate random numbers - that is the work of the

np.random

family, such as the

np.random.normal

used a paragraph ago. The t-test has exactly one job: to compare the means of two groups and rule on whether the difference between them is significant.

When the data are categories, not numbers

The t-test needs numbers, because it computes means. But what if both variables are categories? You note down the time of the expedition (morning or evening) and whether a lion was seen (yes or no), and you want to know whether the time of day has anything to do with it. Data like that is laid out in a contingency table: rows are one trait, columns the other, and the counts of cases sit inside. It is examined by the chi-square test, called as

stats.chi2_contingency

, which returns four things: the chi-square statistic, the p-value, the number of degrees of freedom

dof

, and the

expected

table, the counts you would expect if the two traits were completely independent of one another.

1observed = np.array([[50, 30],
2                     [20, 40]])
3
4chi2, p_value, dof, expected = stats.chi2_contingency(observed)
5
6print(round(chi2, 4))
7# 10.5292
8print(round(p_value, 4))
9# 0.0012
10print(dof)
11# 1

In the mornings a lion was seen 50 times against 30 blank trips, in the evenings the other way round: 20 sightings against 40 blank ones. The p-value of 0.0012 is far below 0.05, so time of day and meeting a lion are not independent - morning expeditions really do have a better hit rate. The

expected

table would hold the values 40, 40, 30 and 30, a picture of a world in which the time of day makes no difference; the further the observations sit from that table, the larger the chi-square statistic.

Do two numbers travel together

The next question from the field sounds different from the previous ones: we are not comparing two groups, we are checking whether two quantities change together. The longer you sit in the hide, the more animals you see - but how strong is that link? It is measured by the Pearson correlation coefficient, computed by

stats.pearsonr

, which takes two arrays of the same length and returns the coefficient together with a p-value.

There are three things you must know about that coefficient. It always lies between -1 and 1. Its sign tells you the direction: positive means "one goes up, the other goes up", negative means "one goes up, the other goes down". Its absolute value tells you the strength: near 1 is a very strong link, near 0 is no link at all. And the crucial caveat: Pearson measures linear dependence only, that is, how well the points line up along a straight line.

1hours = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
2seen = np.array([4, 3, 9, 12, 10, 18, 15, 21, 19, 26])
3
4corr, p_value = stats.pearsonr(hours, seen)
5
6print(round(corr, 4))
7# 0.9526
8print(round(p_value, 6))
9# 2.1e-05

It came out at 0.9526, which rounds to 0.95. We read that off in a single sentence: a strong positive linear correlation - the longer the watch, the more animals, and the points lie close to a straight line. Had it come out at -0.95, the strength would be identical but the direction reversed: a longer watch would mean fewer animals. No correlation lives around zero, say between -0.1 and 0.1, and certainly not at 0.95. One last trap: 0.95 says nothing about a nonlinear relationship, because Pearson simply does not measure such relationships - a high value is exactly the evidence that the link is close to a straight line.

When Pearson misses an obvious link

Let us watch that limit in action. The population of a certain insect species roughly doubles every week, so the link between the week number and the head count is perfect - except that it is exponential rather than straight. Spearman correlation exists for cases like this, and you call it as

stats.spearmanr

. Instead of the values it takes their ranks, their places in the queue, and asks one question only: when one goes up, does the other go up too?

1week = np.array([1, 2, 3, 4, 5, 6, 7, 8])
2insects = np.array([2, 4, 9, 20, 45, 100, 230, 520])
3
4corr, p_value = stats.pearsonr(week, insects)
5print(round(corr, 4))
6# 0.8171
7
8rho, p_value = stats.spearmanr(week, insects)
9print(round(rho, 4))
10# 1.0

The same data, two different answers. Pearson says 0.82, because an exponential curve makes a poor imitation of a straight line. Spearman says 1.0, meaning the link is perfect - every following week gives a higher count than the one before, without a single exception. That does not mean one of the coefficients is lying: they measure different things. Pearson asks "is this a straight line", Spearman asks "is this going up". When you suspect a nonlinear relationship, reach for Spearman.

What your mean is actually worth

Every mean computed from a sample is only an estimate. You measured two hundred giraffes, not every giraffe in the world, so the true mean of the species lies somewhere near yours but not necessarily exactly on it. A confidence interval puts that "somewhere near" into numbers: a 95 percent interval is the range about which we can say, with 95 percent confidence, that it contains the true mean. Two tools are needed:

stats.sem

computes the standard error of the mean, that is, how much a sample mean can wobble, and

stats.t.interval

turns it into an interval. That second call takes the confidence level, the number of degrees of freedom equal to the count of measurements minus one (the built in

len

gives the length of the array), the center

loc

and the scale

scale

. It hands back a pair of numbers: the lower end of the interval sits at index 0, the upper end at index 1.

1mean = np.mean(heights)
2sem = stats.sem(heights)
3ci_95 = stats.t.interval(0.95, len(heights) - 1, loc=mean, scale=sem)
4
5print(round(sem, 3))
6# 2.808
7print(round(ci_95[0], 2), round(ci_95[1], 2))
8# 493.34 504.42

Instead of the dry "the average giraffe height is 498.88 cm" you can now write a much stronger sentence: the true mean lies, with 95 percent confidence, between 493.34 and 504.42 cm. The mean itself has not changed - it has merely gained an honest statement of its own uncertainty. Had you measured ten giraffes instead of two hundred, the interval would be several times wider, because a smaller sample means greater uncertainty. This is the simplest way to stop a report from pretending to a precision it does not actually have.

Linear regression in four steps

One last and most practical question remains: the herd is growing, so how many animals will there be next season? Correlation only says "they rise together", and what you need is a specific number. Linear regression gives it to you: fitting the points with a straight line of the form

y = a * x + b

, where

is the slope, how much

grows when

grows by one, and

is the intercept, the value of

when

is zero. In SciPy one function does the work,

stats.linregress

, but the road to the result always runs through the same four steps and it is worth keeping them in order.

Step one is to prepare the

and
y
data. Without two parallel arrays there is nothing to fit. We type in the season numbers and the head counts recorded for the herd.

1seasons = np.array([1, 2, 3, 4, 5, 6])
2herd = np.array([34, 41, 45, 52, 58, 63])

Nothing has been computed yet, but the data now has the right shape: six seasons and six matching measurements, lined up in pairs in the same order. This is the only moment at which you decide what is the cause and what is the effect -

is the quantity you know,

is the one you want to predict.

Step two is to call

stats.linregress(x, y)

. The function takes both arrays and returns a single object holding the full set of results.

1result = stats.linregress(seasons, herd)

Silence again - nothing was printed, because the result landed inside the variable

result

. Five values are sitting in there at once:

slope

intercept

rvalue

pvalue

and

stderr

. You can also catch them the older way, unpacking all five into separate variables on one line, but I recommend the

result

object with a dot, because with five values the order is easy to muddle, while

result.slope

reads unambiguously.

Step three is to read the

slope

and the
intercept
, that is, to pull the equation of the line out of the result.

1print(round(result.slope, 2))
2# 5.8
3print(round(result.intercept, 2))
4# 28.53

The line has the equation

y = 5.8 * x + 28.53

. In plain words: the herd grows by an average of 5.8 animals per season, and before the first season the fitted starting point sat at 28.53. Substituting the seventh season gives you a prediction of 5.8 times 7 plus 28.53, so about 69 animals.

Step four is to check the R-squared, because a straight line can be fitted even to a cloud of entirely random points and nobody will warn you about it.

rvalue

is the correlation coefficient you already know, and its square - computed in Python with the exponentiation operator

**

- is called R-squared and tells you what share of the variation in

your line explains. One is a perfect fit, zero is a useless line. While we are here we also glance at

pvalue

, which answers the question of whether the slope of the line differs from zero at all.

1print(round(result.rvalue ** 2, 4))
2# 0.9964
3print(round(result.pvalue, 6))
4# 5e-06

The R-squared is 0.9964, so the line explains over 99 percent of the variation in the herd count - the prediction of 69 animals can be trusted. Had it come out at 0.24, as it would for data jumping up and down with no trend, the

slope

would still have been computed and would still have looked serious, and it would have meant nothing. That is why the order of these four steps is not a matter of taste: without data there is nothing to call, without the call there is nowhere to read

slope

and

intercept

from, and without the R-squared you do not know whether the line you read off describes reality at all.

What you take back to camp

The measures of central tendency are the mean, the median and the mode. The one resistant to outliers is the median - a single measurement can shift the mean by hundreds of units, while variance and standard deviation do not describe the middle at all, only the dispersion, and they react to extremes even more strongly.
You describe dispersion with
```
np.std
```
and
```
np.var
```
, and more robustly with quartiles and
```
stats.iqr
```
. The plain range
```
np.ptp
```
dies at the first outlier.
You import the statistics module in the order
```
from
```
,
```
scipy
```
,
```
import
```
,
```
stats
```
. Written the other way round it is a
```
SyntaxError
```
.
The t-test (
stats.ttest_ind
) is for comparing the means of two groups. It does not normalize data, it does not calculate correlation and it does not generate random numbers.
The chi-square test checks for a link between categorical traits from a contingency table.
A Pearson correlation of 0.95 means a strong positive linear correlation. A strong negative one is -0.95, no correlation is around zero, and a nonlinear relationship is something Pearson does not measure at all - Spearman is there for that.
You call the normality test from six parts:
```
stats
```
,
```
.
```
,
```
shapiro
```
,
```
(
```
,
```
data
```
,
```
)
```
, and you read the answer from the p-value against the 0.05 threshold.
A confidence interval from
```
stats.sem
```
and
```
stats.t.interval
```
adds an honest measure of uncertainty to your mean.
Linear regression has four steps in a fixed order: prepare the
```
x
```
and
```
y
```
data, call
```
stats.linregress(x, y)
```
, read the
```
slope
```
and
```
intercept
```
, check the R-squared.

Statistics does not add a single new observation to your notebook, @name - all it gives you is the right to stand up at the briefing and say that the track is real, rather than that it merely looked that way to you out in the field.

Go to CodeWorlds→