Think Stats: Probability and Statistics for Programmers by Allen Downey - HTML preview

PLEASE NOTE: This is an HTML preview only and some elements such as links or page numbers may be incorrect.
Download the book in PDF, ePub for a complete version.

Chapter 4

Continuous distributions

The distributions we have used so far are called empirical distributions

because they are based on empirical observations, which are necessarily

finite samples.

The alternative is a continuous distribution, which is characterized by a

CDF that is a continuous function (as opposed to a step function). Many

real world phenomena can be approximated by continuous distributions.

4.1

The exponential distribution

I’ll start with the exponential distribution because it is easy to work with. In

the real world, exponential distributions come up when we look at a series

of events and measure the times between events, which are called interar-

rival times. If the events are equally likely to occur at any time, the distri-

bution of interarrival times tends to look like an exponential distribution.

The CDF of the exponential distribution is:

CDF(x) = 1 − e− λ x

The parameter, λ, determines the shape of the distribution. Figure 4.1 shows

what this CDF looks like with λ = 2.

In general, the mean of an exponential distribution is 1/ λ, so the mean of

this distribution is 0.5. The median is log(2)/ λ, which is roughly 0.35.

38

Chapter 4. Continuous distributions

Exponential CDF

1.0

0.8

0.6

CDF 0.4

0.2

0.0

0.0

0.5

1.0

1.5

2.0

2.5

x

Figure 4.1: CDF of exponential distribution.

To see an example of a distribution that is approximately exponential, we

will look at the interarrival time of babies. On December 18, 1997, 44 babies

were born in a hospital in Brisbane, Australia1. The times of birth for all 44

babies were reported in the local paper; you can download the data from

❤tt♣✿✴✴t❤✐♥❦st❛ts✳❝♦♠✴❜❛❜②❜♦♦♠✳❞❛t.

Figure 4.2 shows the CDF of the interarrival times in minutes. It seems to

have the general shape of an exponential distribution, but how can we tell?

One way is to plot the complementary CDF, 1 − CDF(x), on a log-y scale.

For data from an exponential distribution, the result is a straight line. Let’s

see why that works.

If you plot the complementary CDF (CCDF) of a dataset that you think is

exponential, you expect to see a function like:

y ≈ e− λ x

Taking the log of both sides yields:

log y ≈ - λ x

So on a log-y scale the CCDF is a straight line with slope − λ.

Figure 4.3 shows the CCDF of the interarrivals on a log-y scale. It is not

exactly straight, which suggests that the exponential distribution is only

1This example is based on information and data from Dunn, “A Simple Dataset for

Demonstrating Common Distributions,” Journal of Statistics Education v.7, n.3 (1999).

4.1. The exponential distribution

39

1.0

Time between births

0.8

0.6

CDF 0.4

0.2

0

0.0

20

40

60

80

100

120

140

160

minutes

Figure 4.2: CDF of interarrival times

Time between births

100

10-1

Complementary CDF

10-2 0

20

40

60

80

100

120

140

160

minutes

Figure 4.3: CCDF of interarrival times.

40

Chapter 4. Continuous distributions

an approximation. Most likely the underlying assumption—that a birth is

equally likely at any time of day—is not exactly true.

Exercise 4.1 For small values of n, we don’t expect an empirical distribution

to fit a continuous distribution exactly. One way to evaluate the quality of

fit is to generate a sample from a continuous distribution and see how well

it matches the data.

The function ❡①♣♦✈❛r✐❛t❡ in the r❛♥❞♦♠ module generates random values

from an exponential distribution with a given value of λ. Use it to generate

44 values from an exponential distribution with mean 32.6. Plot the CCDF

on a log-y scale and compare it to Figure 4.3.

Hint: You can use the function ♣②♣❧♦t✳②s❝❛❧❡ to plot the y axis on a log

scale.

Or, if you use ♠②♣❧♦t, the ❈❞❢ function takes a boolean option, ❝♦♠♣❧❡♠❡♥t,

that determines whether to plot the CDF or CCDF, and string options,

①s❝❛❧❡ and ②s❝❛❧❡, that transform the axes; to plot a CCDF on a log-y scale:

♠②♣❧♦t✳❈❞❢✭❝❞❢✱ ❝♦♠♣❧❡♠❡♥t❂❚r✉❡✱ ①s❝❛❧❡❂✬❧✐♥❡❛r✬✱ ②s❝❛❧❡❂✬❧♦❣✬✮

Exercise 4.2 Collect the birthdays of the students in your class, sort them,

and compute the interarrival times in days. Plot the CDF of the interar-

rival times and the CCDF on a log-y scale. Does it look like an exponential

distribution?

4.2

The Pareto distribution

The Pareto distribution is named after the economist Vilfredo Pareto, who

used it to describe the distribution of wealth (see ❤tt♣✿✴✴✇✐❦✐♣❡❞✐❛✳♦r❣✴

✇✐❦✐✴P❛r❡t♦❴❞✐str✐❜✉t✐♦♥). Since then, it has been used to describe phe-

nomena in the natural and social sciences including sizes of cities and

towns, sand particles and meteorites, forest fires and earthquakes.

The CDF of the Pareto distribution is:

x

α

CDF(x) = 1 −

xm

The parameters xm and α determine the location and shape of the distri-

bution. xm is the minimum possible value. Figure 4.4 shows the CDF of a

Pareto distribution with parameters xm = 0.5 and α = 1.

4.2. The Pareto distribution

41

1.0

Pareto CDF

0.8

0.6

CDF 0.4

0.2

0

0.0

2

4

6

8

10

x

Figure 4.4: CDF of a Pareto distribution.

The median of this distribution is xm21/ α, which is 1, but the 95th percentile

is 10. By contrast, the exponential distribution with median 1 has 95th per-

centile of only 1.5.

There is a simple visual test that indicates whether an empirical distribution

fits a Pareto distribution: on a log-log scale, the CCDF looks like a straight

line. If you plot the CCDF of a sample from a Pareto distribution on a linear

scale, you expect to see a function like:

x

α

y ≈

xm

Taking the log of both sides yields:

log y ≈− α (log x − log xm)

So if you plot log y versus log x, it should look like a straight line with slope

α and intercept α log xm.

Exercise 4.3 The r❛♥❞♦♠ module provides ♣❛r❡t♦✈❛r✐❛t❡, which generates

random values from a Pareto distribution. It takes a parameter for α, but

not xm. The default value for xm is 1; you can generate a distribution with a

different parameter by multiplying by xm.

Write a wrapper function named ♣❛r❡t♦✈❛r✐❛t❡ that takes α and xm as pa-

rameters and uses r❛♥❞♦♠✳♣❛r❡t♦✈❛r✐❛t❡ to generate values from a two-

parameter Pareto distribution.

Use your function to generate a sample from a Pareto distribution. Com-

pute the CCDF and plot it on a log-log scale. Is it a straight line? What is

42

Chapter 4. Continuous distributions

the slope?

Exercise 4.4 To get a feel for the Pareto distribution, imagine what the world

would be like if the distribution of human height were Pareto. Choosing the

parameters xm = 100 cm and α = 1.7, we get a distribution with a reasonable

minimum, 100 cm, and median, 150 cm.

Generate 6 billion random values from this distribution. What is the mean

of this sample? What fraction of the population is shorter than the mean?

How tall is the tallest person in Pareto World?

Exercise 4.5 Zipf’s law is an observation about how often different words

are used. The most common words have very high frequencies, but there

are many unusual words, like “hapaxlegomenon,” that appear only a few

times. Zipf’s law predicts that in a body of text, called a “corpus,” the dis-

tribution of word frequencies is roughly Pareto.

Find a large corpus, in any language, in electronic format. Count how many

times each word appears. Find the CCDF of the word counts and plot it on a

log-log scale. Does Zipf’s law hold? What is the value of α, approximately?

Exercise 4.6 The Weibull distribution is a generalization of the exponential

distribution that comes up in failure analysis (see ❤tt♣✿✴✴✇✐❦✐♣❡❞✐❛✳♦r❣✴

✇✐❦✐✴❲❡✐❜✉❧❧❴❞✐str✐❜✉t✐♦♥). Its CDF is

CDF(x) = 1 − e−(x/ λ)k

Can you find a transformation that makes a Weibull distribution look like a

straight line? What do the slope and intercept of the line indicate?

Use r❛♥❞♦♠✳✇❡✐❜✉❧❧✈❛r✐❛t❡ to generate a sample from a Weibull distribu-

tion and use it to test your transformation.

4.3

The normal distribution

The normal distribution, also called Gaussian, is the most commonly used

because it describes so many phenomena, at least approximately. It turns

out that there is a good reason for its ubiquity, which we will get to in Sec-

tion 6.6.

The normal distribution has many properties that make it amenable for

analysis, but the CDF is not one of them. Unlike the other distributions

4.3. The normal distribution

43

1.0

Normal CDF

0.8

0.6

CDF 0.4

0.2

0.0

0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

x

Figure 4.5: CDF of a normal distribution.

we have looked at, there is no closed-form expression for the normal CDF;

the most common alternative is to write it in terms of the error function,

which is a special function written erf(x):

1

x − µ

CDF(x) =

1 + erf

2

σ

2

2

x

erf(x) = √

e−t2dt

π

0

The parameters µ and σ determine the mean and standard deviation of the distribution.

If these formulas make your eyes hurt, don’t worry; they are easy to imple-

ment in Python2. There are many fast and accurate ways to approximate

erf(x). You can download one of them from ❤tt♣✿✴✴t❤✐♥❦st❛ts✳❝♦♠✴❡r❢✳

♣②, which provides functions named ❡r❢ and ◆♦r♠❛❧❈❞❢.

Figure 4.5 shows the CDF of the normal distribution with parameters µ = 2.0

and σ = 0.5. The sigmoid shape of this curve is a recognizable characteristic

of a normal distribution.

In the previous chapter we looked at the distribution of birth weights in the

NSFG. Figure 4.6 shows the empirical CDF of weights for all live births and

the CDF of a normal distribution with the same mean and variance.

The normal distribution is a good model for this dataset. A model is a useful

simplification. In this case it is useful because we can summarize the entire

2As of Python 3.2, it is even easier; ❡r❢ is in the ♠❛t❤ module.

44

Chapter 4. Continuous distributions

Birth weights

1.0

model

data

0.8

0.6

CDF 0.4

0.2

0

0.0

50

100

150

200

250

birth weight (oz)

Figure 4.6: CDF of birth weights with a normal model.

distribution with just two numbers, µ = 116.5 and σ = 19.9, and the resulting error (difference between the model and the data) is small.

Below the 10th percentile there is a discrepancy between the data and the

model; there are more light babies than we would expect in a normal distri-

bution. If we are interested in studying preterm babies, it would be impor-

tant to get this part of the distribution right, so it might not be appropriate

to use the normal model.

Exercise 4.7 The Wechsler Adult Intelligence Scale is a test that is intended

to measure intelligence3. Results are transformed so that the distribution of

scores in the general population is normal with µ = 100 and σ = 15.

Use ❡r❢✳◆♦r♠❛❧❈❞❢ to investigate the frequency of rare events in a normal

distribution. What fraction of the population has an IQ greater than the

mean? What fraction is over 115? 130? 145?

A “six-sigma” event is a value that exceeds the mean by 6 standard devia-

tions, so a six-sigma IQ is 190. In a world of 6 billion people, how many do

we expect to have an IQ of 190 or more4?

Exercise 4.8 Plot the CDF of pregnancy lengths for all live births. Does it

look like a normal distribution?

3Whether it does or not is a fascinating controversy that I invite you to investigate at

your leisure.

4On this topic, you might be interested to read ❤tt♣✿✴✴✇✐❦✐♣❡❞✐❛✳♦r❣✴✇✐❦✐✴

❈❤r✐st♦♣❤❡r❴▲❛♥❣❛♥.

4.4. Normal probability plot

45

Compute the mean and variance of the sample and plot the normal distri-

bution with the same parameters. Is the normal distribution a good model

for this data? If you had to summarize this distribution with two statistics,

what statistics would you choose?

4.4

Normal probability plot

For the exponential, Pareto and Weibull distributions, there are simple

transformations we can use to test whether a continuous distribution is a

good model of a dataset.

For the normal distribution there is no such transformation, but there is an

alternative called a normal probability plot. It is based on rankits: if you generate n values from a normal distribution and sort them, the kth rankit

is the mean of the distribution for the kth value.

Exercise 4.9 Write a function called ❙❛♠♣❧❡ that generates 6 samples from a

normal distribution with µ = 0 and σ = 1. Sort and return the values.

Write a function called ❙❛♠♣❧❡s that calls ❙❛♠♣❧❡ 1000 times and returns a

list of 1000 lists.

If you apply ③✐♣ to this list of lists, the result is 6 lists with 1000 values in

each. Compute the mean of each of these lists and print the results. I predict

that you will get something like this:

{−1.2672, −0.6418, −0.2016, 0.2016, 0.6418, 1.2672}

If you increase the number of times you call ❙❛♠♣❧❡, the results should con-

verge on these values.

Computing rankits exactly is moderately difficult, but there are numerical

methods for approximating them. And there is a quick-and-dirty method

that is even easier to implement:

1. From a normal distribution with µ = 0 and σ = 1, generate a sample

with the same size as your dataset and sort it.

2. Sort the values in the dataset.

3. Plot the sorted values from your dataset versus the random values.

46

Chapter 4. Continuous distributions

250

200

150

100

Birth weights (oz)

50

0 4

3

2

1

0

1

2

3

4

Standard normal values

Figure 4.7: Normal probability plot of birth weights.

For large datasets, this method works well. For smaller datasets, you can

improve it by generating m(n+1) − 1 values from a normal distribution,

where n is the size of the dataset and m is a multiplier. Then select every

mth element, starting with the mth.

This method works with other distributions as well, as long as you know

how to generate a random sample.

Figure 4.7 is a quick-and-dirty normal probability plot for the birth weight

data.

The curvature in this plot suggests that there are deviations from a normal

distribution; nevertheless, it is a good (enough) model for many purposes.

Exercise 4.10 Write a function called ◆♦r♠❛❧P❧♦t that takes a sequence of

values and generates a normal probability plot. You can download a solu-

tion from ❤tt♣✿✴✴t❤✐♥❦st❛ts✳❝♦♠✴r❛♥❦✐t✳♣②.

Use the running speeds from r❡❧❛②✳♣② to generate a normal probability

plot. Is the normal distribution a good model for this data? You can down-

load a solution from ❤tt♣✿✴✴t❤✐♥❦st❛ts✳❝♦♠✴r❡❧❛②❴♥♦r♠❛❧✳♣②.

4.5

The lognormal distribution

If the logarithms of a set of values have a normal distribution, the values

have a lognormal distribution. The CDF of the lognormal distribution is

4.5. The lognormal distribution

47

Adult weight

1.0

model

data

0.8

0.6

CDF 0.4

0.2

0.0 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2

adult weight (log kg)

Figure 4.8: CDF of adult weights (log transform).

the same as the CDF of the normal distribution, with log x substituted for

x.

CDFlognormal(x) = CDFnormal(log x)

The parameters of the lognormal distribution are usually denoted µ and σ.

But remember that these parameters are not the mean and standard devia-

tion; the mean of a lognormal distribution is exp(

2

µ + σ /2) and the standard

deviation is ugly5.

It turns out that the distribution of weights for adults is approximately log-

normal6.

The National Center for Chronic Disease Prevention and Health Promotion

conducts an annual survey as part of the Behavioral Risk Factor Surveil-

lance System (BRFSS)7. In 2008, they interviewed 414,509 respondents and

asked about their demographics, health and health risks.

Among the data they collected are the weights in kilograms of 398,484 re-

spondents. Figure 4.8 shows the distribution of log w, where w is weight in

5See ❤tt♣✿✴✴✇✐❦✐♣❡❞✐❛✳♦r❣✴✇✐❦✐✴▲♦❣✲♥♦r♠❛❧❴❞✐str✐❜✉t✐♦♥.

6I was tipped off to this possibility by a comment (without citation) at ❤tt♣✿✴✴

♠❛t❤✇♦r❧❞✳✇♦❧❢r❛♠✳❝♦♠✴▲♦❣◆♦r♠❛❧❉✐str✐❜✉t✐♦♥✳❤t♠❧. Subsequently I found a paper

that proposes the log transform and suggests a cause: Penman and Johnson, “The Chang-

ing Shape of the Body Mass Index Distribution Curve in the Population,” Preventing

Chronic Disease, 2006 July; 3(3): A74. Online at ❤tt♣✿✴✴✇✇✇✳♥❝❜✐✳♥❧♠✳♥✐❤✳❣♦✈✴♣♠❝✴

❛rt✐❝❧❡s✴P▼❈✶✻✸✻✼✵✼.

7Centers for Disease Control and Prevention (CDC). Behavioral Risk Factor Surveillance

System Survey Data. Atlanta, Georgia: U.S. Department of Health and Human Services,

Centers for Disease Control and Prevention, 2008.

48

Chapter 4. Continuous distributions

kilograms, along with a normal model.

The normal model is a good fit for the data, although the highest weights

exceed what we expect from the normal model even after the log transform.

Since the distribution of log w fits a normal distribution, we conclude that

w fits a lognormal distribution.

Exercise 4.11 Download the BRFSS data from ❤tt♣✿✴✴t❤✐♥❦st❛ts✳❝♦♠✴

❈❉❇❘❋❙✵✽✳❆❙❈✳❣③, and my code for reading it from ❤tt♣✿✴✴t❤✐♥❦st❛ts✳

❝♦♠✴❜r❢ss✳♣②. Run ❜r❢ss✳♣② and confirm that it prints summary statistics

for a few of the variables.

Write a program that reads adult weights from the BRFSS and generates

normal probability plots for w and log w. You can download a solution

from ❤tt♣✿✴✴t❤✐♥❦st❛ts✳❝♦♠✴❜r❢ss❴❢✐❣s✳♣②.

Exercise 4.12 The distribution of populations for cities and towns has been

proposed as an example of a real-world phenomenon that can be described

with a Pareto distribution.

The U.S. Census Bureau publishes data on the population of every incorpo-

rated city and town in the United States. I have written a small program

that downloads this data and stores it in a file. You can download it from

❤tt♣✿✴✴t❤✐♥❦st❛ts✳❝♦♠✴♣♦♣✉❧❛t✐♦♥s✳♣②.

1. Read over the program to make sure you know what it does; then run

it to download and process the data.

2. Write a program that computes and plots the distribution of popula-

tions for the 14,593 cities and towns in the dataset.

3. Plot the CDF on linear and log-x scales so you can get a sense of the

shape of the distribution. Then plot the CCDF on a log-log scale to see

if it has the characteristic shape of a Pareto distribution.

4. Try out the other transformations and plots in this chapter to see if

there is a better model for this data.

What conclusion do you draw about the distribution of sizes for cities

and towns? You can download a solution from ❤tt♣✿✴✴t❤✐♥❦st❛ts✳❝♦♠✴

♣♦♣✉❧❛t✐♦♥s❴❝❞❢✳♣②.

Exercise 4.13 The Internal Revenue Service of the United States (IRS) pro-

vides data about income taxes at ❤tt♣✿✴✴✐rs✳❣♦✈✴t❛①st❛ts.

4.6. Why model?

49

One of their files, containing informatio