1  Variables and Distributions

Time to play in data land! We began our journey with a promise to help you understand and quantify your uncertainty. Now, we shift to the heart and soul of data - the distribution. Every variable has a distribution, and understanding that distribution is key to understanding the data. We will explore the most common distributions and how to use them to make sense of your data. Let’s get started.

1.1 Learning Objectives

  1. Understand the concept of a variable and how they may be represented in data.
  2. Learn about the most common distributions and how to use them.
  3. Understand the concept of a probability density function and how to use it.
  4. Learn how to use distributions to make predictions and quantify uncertainty.

1.2 Variables

Everything we wish to understand about people or things comes in the form of a variable. If I want to know why people like chocolate ice cream, I must assign some form of numerical value to “like” such that I can learn from my observations. Why? Learning via these methods require data. Each variable is a set of observations. We assign numbers to observations and these observations serve as our data. Each variable is unique. Some variables contain data that represent a singular characteristic of a person or thing; other variables contain data that represent multiple characteristics. These characteristics we often call constructs in psychological science. From this point forward, we shall refer to the underlying characteristic that we measure as the construct.

Let’s examine this a bit more closely. Consider two variables height and weight. Height is a variable that contains data representing a singular characteristic of a person - the length of the person’s body when standing. Changes in height measures a singular characteristic of a person that can be altered only by skeletal growth/decay. Weight, on the other hand, contains data representing a single characteristic (mass) but determined by many contributing factors (e.g., water, fat, muscle, bone, etc.). Combined, height and weight make a composite variable called Body Mass Index (BMI) - an often used but derided predictor of health and well-being.

Singular Construct Multiple Constructs
Singular Cause IDEAL (Height) Less IDEAL (Weight)
Multiple Cause Most Common (g) REALLY BAD (BMI)

In the table above, we see that the ideal variable is one that contains data representing a singular characteristic of a person and is determined by a singular cause. Similarly, we see BMI labeled as “REALLY BAD.” These are abstractions right now but they become important as we delve into psychological science and the uncertainty around our constructs.

1.2.1 Types of Variables

We have many types of variables to choose from when we begin our journey into data land. Each type of variable has its own unique characteristics and uses. Let’s explore the most common types of variables. Decades ago, psychological scientists were hamstrung by W.W. Stevens levels of measurement. Stevens argued that there were four levels of measurement: nominal, ordinal, interval, and ratio. These levels of measurement were used to determine the permissible analyses based solely on the logic of arithmetic. Times have changed but arithmetic remains the same. The argument, however, that the level of measurement dictates the analysis is simply untrue. We may do whatever analysis we choose. Which ones make sense depends upon a whole host of factors that we address in this book. So, keep reading. Below is a much easier taxonomy for understanding variables.

1.2.1.1 Continuous Variables

Continuous variables are those that can take on any value within a range. These variables are often measured and can be divided into smaller and smaller units. Examples of continuous variables include height, weight, and temperature. We care not only about the integer values (e.g., 5 feet) but also the decimal values as well (e.g., 5.5 feet). Continuous variables are the most commonly assumed type of variable in psychological science. We emphasize assumed here because we see many instances where continuity, that is, the continuous nature of a variable, makes no sense and yet, in those instances, is treated as such. This is a problem we will return to. For now, a continuous variable is a variable containing any value between a range of possible values including negative and positive infinity.

1.2.1.2 Discrete Variables

Discrete variables are the countable cousins of continuous variables. They take on separated, distinct values - usually whole numbers - with nothing in between. The number of children in a family, the number of times you checked your phone today, the number of correct answers on a quiz: all discrete. You can have 2 children or 3 children, but never 2.4 of them. The gap matters. A Likert scale is the classic troublemaker: a 0-to-5 item offers only the integers 0, 1, 2, 3, 4, 5, so it is technically discrete, yet we routinely treat it as continuous - a convenience we return to later. Much of the mischief in applied statistics comes from treating a discrete variable as if it were continuous (or the reverse), so always ask yourself: can the thing between two values actually exist?

1.2.1.3 Categorical Variables

Categorical variables sort observations into buckets that have no inherent numerical order - think of pet (dog, cat, mouse, fish), eye color, political party, or ice-cream flavor. Because the values are non-numeric, they must be converted into numbers at some point before we can compute with them; we can code them (1 = chocolate, 2 = vanilla) for convenience, but those numbers are labels, not quantities - vanilla is not “twice” chocolate. A categorical variable with exactly two buckets (yes/no, treatment/control, alive/dead) is called dichotomous or binary, and it turns out to be one of the most useful variables in all of statistics, as we will see when point-biserial correlation and ANOVA arrive.

1.2.1.4 Ordinal Variables

Ordinal variables are categories with a rank order but without guaranteed equal spacing between ranks. Think of a race (1st, 2nd, 3rd) or a survey item (strongly disagree → strongly agree). We know 1st beat 2nd, but the ordering does not tell us by how much - the gap between 1st and 2nd may be a mile and the gap between 2nd and 3rd a whisker. Ordinal data live in the uneasy space between categorical and continuous, and reasonable people argue about how to analyze them. We will side-step the dogma: what you can do depends on the question, not on a label Stevens handed down.

1.2.2 Representing Variables Numerically

To compute anything, we must turn observations into numbers. For continuous and discrete variables that is natural - the number is the measurement. For categorical variables we assign codes, but we must remember that the software does not know those codes are arbitrary; it will happily average “flavor” and hand you 1.7 if you let it. In R, we protect ourselves by storing categories as factors, which keeps their labels attached and their order (if any) explicit:

NoteWorking in SPSS, Julia, or Python?

The code tabs below assume this chapter’s data is already loaded. Grab the one-file setup for your language from Getting the Book’s Data, run it once, then load what you need by name - this chapter uses ch01-heights, ch01-skewed, ch01-shapes. For example, book_data("ch01-heights") in R, Julia, or Python, or !bookdata name = "ch01-heights". in SPSS. Every language reads the same shipped files, so your numbers will match the ones printed here exactly.

library(tidyverse)    # dplyr, ggplot2, purrr, tibble, readr, stringr, forcats
source("_common.R")   # book-wide helpers: round2(), fmt_p(), tidy2()

flavor <- fct(c("choc", "van", "choc", "straw", "van"),
              levels = c("choc", "van", "straw"))
as.integer(flavor)                 # the codes R uses under the hood
[1] 1 2 1 3 2
tibble(flavor) |> count(flavor)    # what we actually care about: the counts
flavor n
choc 2
van 2
straw 1
* Categories live as string values with VALUE LABELS attached; the labels
* travel with the data the way an R factor's levels do.
DATA LIST LIST /flavor (A6).
BEGIN DATA
choc
van
choc
straw
van
END DATA.

* AUTORECODE exposes the integer codes sitting under the labels.
AUTORECODE VARIABLES=flavor /INTO flavor_code.
* What we actually care about: the counts.
FREQUENCIES VARIABLES=flavor.
Reading free-form data from INLINE.
+--------+------+
|Variable|Format|
+--------+------+
|flavor  |A6    |
+--------+------+

                             flavor
+-----------+---------+-------+-------------+------------------+
|           |Frequency|Percent|Valid Percent|Cumulative Percent|
+-----------+---------+-------+-------------+------------------+
|Valid choc |        2|  40.0%|        40.0%|             40.0%|
|      straw|        1|  20.0%|        20.0%|             60.0%|
|      van  |        2|  40.0%|        40.0%|            100.0%|
+-----------+---------+-------+-------------+------------------+
|Total      |        5| 100.0%|             |                  |
+-----------+---------+-------+-------------+------------------+
using CategoricalArrays, DataFrames
flavor = categorical(["choc", "van", "choc", "straw", "van"],
                     levels = ["choc", "van", "straw"])
levelcode.(flavor)                  # the codes Julia uses under the hood
combine(groupby(DataFrame(flavor = flavor), :flavor), nrow => :n)
5-element Vector{Int64}:
 1
 2
 1
 3
 2
3×2 DataFrame
 Row │ flavor  n
     │ Cat…    Int64
─────┼───────────────
   1 │ choc        2
   2 │ van         2
   3 │ straw       1
import pandas as pd
flavor = pd.Categorical(["choc", "van", "choc", "straw", "van"],
                        categories=["choc", "van", "straw"])
flavor.codes                        # the codes pandas uses under the hood
pd.Series(flavor).value_counts()    # what we actually care about: the counts
array([0, 1, 0, 2, 1], dtype=int8)
choc     2
van      2
straw    1
Name: count, dtype: int64

1.2.3 Plotting Variables

Before you compute a single statistic, look at your data. A picture reveals shape, spread, gaps, and outliers that a lone number will hide from you. The workhorse for a single variable is the histogram.

1.2.3.1 Histograms

A histogram slices the range of a variable into bins and counts how many observations fall in each. The result is a picture of the variable’s distribution - the pattern of how often each value shows up.

set.seed(20)
heights <- tibble(height = rnorm(500, mean = 68, sd = 4))  # adult heights (inches)

ggplot(heights, aes(x = height)) +
  # y = density, not count, so the bars and the curve can share one axis
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 fill = "steelblue", colour = "white") +
  # the theoretical curve, drawn straight from the parameters we simulated with
  stat_function(fun = dnorm, args = list(mean = 68, sd = 4),
                colour = "firebrick", linewidth = 1) +
  labs(x = "height (inches)", y = "density") +
  theme_book()
Figure 1.1: Five hundred simulated adult heights, with the normal curve they were drawn from. The bars are the sample; the curve is the model.
* The book's own 500 heights, then the histogram with a normal curve on it.
* (NORMAL) asks SPSS to overlay the fitted normal, the same curve the R tab
* draws from the parameters.
INSERT FILE='data/sim/ch01-heights.sps'.

GRAPH /HISTOGRAM(NORMAL) = height.

using CSV, DataFrames, Plots, Distributions
d = CSV.read("data/sim/ch01-heights.csv", DataFrame)   # the book's heights

histogram(d.height, bins = 30, normalize = :pdf, legend = false,
          xlabel = "height (inches)", ylabel = "density")
plot!(x -> pdf(Normal(68, 4), x), 54, 82, linewidth = 2, color = :firebrick)

import numpy as np
import matplotlib.pyplot as plt

import pandas as pd
from scipy.stats import norm

d = pd.read_csv("data/sim/ch01-heights.csv")   # the book's heights

plt.hist(d.height, bins=30, density=True, color="steelblue", edgecolor="white")
grid = np.linspace(54, 82, 400)
plt.plot(grid, norm.pdf(grid, 68, 4), color="firebrick", linewidth=2)
plt.xlabel("height (inches)"); plt.ylabel("density")
plt.show()
(array([0.00257652, 0.00772956, 0.00515304, 0.00772956, 0.00772956,
       0.01545911, 0.02318867, 0.04637733, 0.03864778, 0.06183644,
       0.06698948, 0.11336681, 0.0979077 , 0.06698948, 0.08502511,
       0.10821377, 0.11336681, 0.08502511, 0.04637733, 0.06441296,
       0.06183644, 0.05153037, 0.03864778, 0.02318867, 0.01803563,
       0.01288259, 0.00772956, 0.00257652, 0.00257652, 0.00515304]), array([56.44112954, 57.21737083, 57.99361212, 58.76985341, 59.5460947 ,
       60.32233599, 61.09857728, 61.87481857, 62.65105986, 63.42730114,
       64.20354243, 64.97978372, 65.75602501, 66.5322663 , 67.30850759,
       68.08474888, 68.86099017, 69.63723146, 70.41347275, 71.18971404,
       71.96595532, 72.74219661, 73.5184379 , 74.29467919, 75.07092048,
       75.84716177, 76.62340306, 77.39964435, 78.17588564, 78.95212693,
       79.72836822]), <BarContainer object of 30 artists>)
[<matplotlib.lines.Line2D object at 0x7f25f289e180>]
Text(0.5, 0, 'height (inches)')
Text(0, 0.5, 'density')

That mound-shaped pile is the thing we spend the rest of the book reasoning about. Notice we drew two things: blue bars and a red curve. They are not the same kind of object, and confusing them causes more trouble than almost anything else at this stage, so let us be blunt about it.

ImportantThe bars are your data. The curve is your model.

A histogram is a picture of the data you actually collected. Slice the range into bins, count how many observations land in each. Collect 500 more people and the bars change. Choose different bins and the bars change again. The histogram belongs to your sample.

A density curve is a picture of the idealized distribution you think produced that data. It is a smooth mathematical function fixed entirely by its parameters - here a mean of 68 and a standard deviation of 4 - and it would be the same curve whether you collected 500 people or 5 million. The curve belongs to the model.

They can be drawn on the same axes only if they use the same vertical scale. That is why we plotted y = density rather than the default count: on the density scale, the total area of the bars is 1, and the total area under the curve is 1 as well. Now they are comparable, and the question “does my data look like the model?” becomes a question you can answer by eye.

One more consequence worth holding onto: on a density scale the height of the curve is not a probability. Probability is area, never height. We will lean on that hard in a moment.

Throughout this book, when you see a histogram you will almost always see its curve too. The curve is drawn from theory, not fitted by eye, and it carries no error band - a theoretical distribution is exactly what its parameters say it is. (Later, when we fit models to data rather than draw distributions, error bands come back and they matter enormously.)

1.3 Shapes

Distributions have shapes, and the shape tells a story. Three features do most of the work:

  • Central tendency - where the pile sits (the subject of the next chapter).
  • Spread - how wide the pile is (the chapter after that).
  • Symmetry - whether the pile leans. A distribution with a long right tail is right- (positively) skewed; a long left tail is left- (negatively) skewed. Income is famously right-skewed: most people cluster low, a few soar high, and the mean gets dragged toward the millionaires.

A skewed distribution is not a different species of thing. It is a normal curve that leans. We can build one out of nothing but ordinary normal draws, and doing so keeps us inside a single family for the whole book - one shape, one pair of familiar parameters, plus one knob for the lean:

\[X = \delta\,|Z_0| + \sqrt{1-\delta^2}\;Z_1, \qquad \delta = \frac{\alpha}{\sqrt{1+\alpha^2}}\]

Both \(Z_0\) and \(Z_1\) are standard normals. Folding one of them to be positive (\(|Z_0|\)) is what creates the lean; \(\alpha\) controls how much. Set \(\alpha = 0\) and the second term takes over completely, leaving you with exactly the normal curve you already know. This is the skew-normal, and rsnorm() in the book’s _common.R is those two lines and nothing else.

set.seed(21)
# alpha = 6 leans hard to the right; alpha = 0 would give the plain normal
skewed <- tibble(value = rsnorm(1000, mean = 0, sd = 1, alpha = 6))

ggplot(skewed, aes(x = value)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40,
                 fill = "darkorange", colour = "white") +
  stat_function(fun = dsnorm, args = list(mean = 0, sd = 1, alpha = 6),
                colour = "firebrick", linewidth = 1) +
  labs(x = "value", y = "density") +
  theme_book()
Figure 1.2: A skew-normal with alpha = 6 - still a bell, but leaning right. Set alpha to zero and the red curve becomes the ordinary normal.
* The book's skew-normal draws. The construction - two normal draws, one of
* them folded positive - is spelled out in the R tab; here we read the same
* 1000 values so the picture matches.
INSERT FILE='data/sim/ch01-skewed.sps'.

GRAPH /HISTOGRAM(NORMAL) = value.

using CSV, DataFrames, Plots, Distributions
d = CSV.read("data/sim/ch01-skewed.csv", DataFrame)   # the book's skew-normal
alpha = 6                                             # the lean it was built with

histogram(d.value, bins = 40, normalize = :pdf, legend = false,
          xlabel = "value", ylabel = "density")
# Distributions.jl ships the skew-normal, so the curve is one call
plot!(x -> pdf(SkewNormal(0, 1, alpha), x), -2, 4,
      linewidth = 2, color = :firebrick)

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import skewnorm

import pandas as pd

d = pd.read_csv("data/sim/ch01-skewed.csv")   # the book's skew-normal
value = d.value
alpha = 6                                    # the lean it was built with

plt.hist(value, bins=40, density=True, color="darkorange", edgecolor="white")
grid = np.linspace(-2, 4, 400)
plt.plot(grid, skewnorm.pdf(grid, alpha), color="firebrick", linewidth=2)
plt.title("A right-skewed normal (long tail to the right)")
plt.xlabel("value"); plt.ylabel("density")
plt.show()
(array([0.04422955, 0.07740171, 0.11057388, 0.26537731, 0.44229551,
       0.60815633, 0.63027111, 0.80718931, 0.86247625, 0.77401715,
       0.67450066, 0.65238588, 0.48652507, 0.42018074, 0.46441029,
       0.44229551, 0.44229551, 0.42018074, 0.46441029, 0.23220514,
       0.22114776, 0.24326253, 0.2764347 , 0.13268865, 0.14374604,
       0.18797559, 0.09951649, 0.06634433, 0.05528694, 0.0884591 ,
       0.04422955, 0.03317216, 0.02211478, 0.        , 0.02211478,
       0.05528694, 0.03317216, 0.        , 0.        , 0.01105739]), array([-0.3766044 , -0.28616713, -0.19572986, -0.10529258, -0.01485531,
        0.07558196,  0.16601923,  0.25645651,  0.34689378,  0.43733105,
        0.52776832,  0.6182056 ,  0.70864287,  0.79908014,  0.88951742,
        0.97995469,  1.07039196,  1.16082923,  1.25126651,  1.34170378,
        1.43214105,  1.52257832,  1.6130156 ,  1.70345287,  1.79389014,
        1.88432742,  1.97476469,  2.06520196,  2.15563923,  2.24607651,
        2.33651378,  2.42695105,  2.51738832,  2.6078256 ,  2.69826287,
        2.78870014,  2.87913741,  2.96957469,  3.06001196,  3.15044923,
        3.24088651]), <BarContainer object of 40 artists>)
[<matplotlib.lines.Line2D object at 0x7fb30e6c7470>]
Text(0.5, 1.0, 'A right-skewed normal (long tail to the right)')
Text(0.5, 0, 'value')
Text(0, 0.5, 'density')

1.4 Distributions

A distribution is the full accounting of how often each value (or range of values) occurs. Some distributions are so common - they arise from the same kinds of processes again and again - that we have given them names, formulas, and R functions. Knowing a handful of them is like knowing a handful of chords: most songs are built from them.

1.4.1 Names for Shapes

  • Normal (Gaussian) - the bell curve. Arises whenever many small, independent influences add up (height, measurement error). rnorm().
  • Uniform - every value in a range equally likely. runif().
  • Binomial - the count of “successes” in a fixed number of yes/no trials (heads in 10 coin flips). rbinom().
  • Poisson - the count of rare events in a fixed window (typos per page, emails per hour). rpois().

Each panel below carries its theoretical curve in red, alongside the 1000 values we drew. Look at the red marks before the grey bars: the curve is the shape itself, and the bars are one sample’s noisy attempt at it.

Two of these shapes are continuous (Normal, Uniform) - a value can land anywhere, so the curve is an unbroken line. Two are discrete (Binomial, Poisson) - you can flip 6 heads or 7, never 6.4 - so their “curve” is really a row of points sitting on the whole numbers, with nothing in between. That is why we drew them differently. It matters: for a continuous shape, probability is the area under a stretch of the curve; for a discrete one, probability sits right on each point and you simply add the points up.

set.seed(22)
# draw 1000 values from each of the four named shapes, then stack them
# into one long tibble so ggplot can facet them side by side
draws <- tibble(
  Normal   = rnorm(1000),
  Uniform  = runif(1000),
  Binomial = rbinom(1000, 10, .5),
  Poisson  = rpois(1000, 3)
) |>
  pivot_longer(everything(), names_to = "shape", values_to = "value")

# Each shape's theoretical curve, computed from its parameters - not fitted.
# The two continuous shapes get a smooth line; the two discrete ones get
# points, because their probability sits only on whole numbers.
smooth_curves <- bind_rows(
  tibble(shape = "Normal",  value = seq(-4, 4, length.out = 300)) |>
    mutate(density = dnorm(value)),
  tibble(shape = "Uniform", value = seq(0, 1, length.out = 300)) |>
    mutate(density = dunif(value))
)
point_curves <- bind_rows(
  tibble(shape = "Binomial", value = 0:10) |> mutate(density = dbinom(value, 10, .5)),
  tibble(shape = "Poisson",  value = 0:12) |> mutate(density = dpois(value, 3))
)

ggplot(draws, aes(x = value)) +
  # continuous shapes: narrow bins
  geom_histogram(data = \(d) filter(d, shape %in% c("Normal", "Uniform")),
                 aes(y = after_stat(density)), bins = 30,
                 fill = "grey85", colour = "white") +
  # discrete shapes: one bin per whole number, or the density scale misleads
  geom_histogram(data = \(d) filter(d, shape %in% c("Binomial", "Poisson")),
                 aes(y = after_stat(density)), binwidth = 1,
                 fill = "grey85", colour = "white") +
  geom_line(data = smooth_curves, aes(y = density),
            colour = "firebrick", linewidth = 1) +
  geom_point(data = point_curves, aes(y = density),
             colour = "firebrick", size = 1.8) +
  facet_wrap(~ shape, scales = "free") +
  labs(x = NULL, y = "density") +
  theme_book()
Figure 1.3: Four named distributions, each carrying its theoretical curve. Normal and Uniform are continuous, so the curve is an unbroken line; Binomial and Poisson are discrete, so their probability sits on the whole numbers as points.
* The book's own 1000 draws from each shape. SPSS has no faceting, so each
* shape gets its own GRAPH command.
INSERT FILE='data/sim/ch01-shapes.sps'.

GRAPH /HISTOGRAM(NORMAL) = Normal.
GRAPH /HISTOGRAM = Uniform.
GRAPH /HISTOGRAM = Binomial.
GRAPH /HISTOGRAM = Poisson.

using CSV, DataFrames, Plots
d = CSV.read("data/sim/ch01-shapes.csv", DataFrame)   # the book's draws

shapes = [n => d[!, n] for n in ["Normal", "Uniform", "Binomial", "Poisson"]]

# one panel per shape, laid out 2 x 2
plot([histogram(v, bins = 30, title = k, legend = false) for (k, v) in shapes]...,
     layout = (2, 2))

import numpy as np
import matplotlib.pyplot as plt

import pandas as pd

d = pd.read_csv("data/sim/ch01-shapes.csv")   # the book's draws
shapes = {n: d[n] for n in ["Normal", "Uniform", "Binomial", "Poisson"]}

fig, axes = plt.subplots(2, 2)             # one panel per shape
for ax, (name, v) in zip(axes.flat, shapes.items()):
    ax.hist(v, bins=30, color="lightgrey", edgecolor="white")
    ax.set_title(name)
fig.tight_layout(); plt.show()
(array([ 1.,  0.,  0.,  0.,  1.,  6., 11., 14., 30., 33., 44., 54., 81.,
       75., 94., 87., 86., 82., 83., 53., 46., 44., 23., 25.,  9.,  3.,
        5.,  2.,  4.,  4.]), array([-3.66705338, -3.43451698, -3.20198059, -2.96944419, -2.7369078 ,
       -2.50437141, -2.27183501, -2.03929862, -1.80676222, -1.57422583,
       -1.34168944, -1.10915304, -0.87661665, -0.64408025, -0.41154386,
       -0.17900746,  0.05352893,  0.28606532,  0.51860172,  0.75113811,
        0.98367451,  1.2162109 ,  1.4487473 ,  1.68128369,  1.91382008,
        2.14635648,  2.37889287,  2.61142927,  2.84396566,  3.07650205,
        3.30903845]), <BarContainer object of 30 artists>)
Text(0.5, 1.0, 'Normal')
(array([45., 36., 33., 31., 32., 29., 28., 35., 25., 37., 46., 38., 42.,
       31., 33., 28., 23., 35., 31., 34., 30., 33., 35., 39., 45., 29.,
       24., 35., 32., 26.]), array([0.00293194, 0.03611783, 0.06930371, 0.10248959, 0.13567548,
       0.16886136, 0.20204724, 0.23523313, 0.26841901, 0.30160489,
       0.33479078, 0.36797666, 0.40116254, 0.43434843, 0.46753431,
       0.5007202 , 0.53390608, 0.56709196, 0.60027785, 0.63346373,
       0.66664961, 0.6998355 , 0.73302138, 0.76620726, 0.79939315,
       0.83257903, 0.86576491, 0.8989508 , 0.93213668, 0.96532256,
       0.99850845]), <BarContainer object of 30 artists>)
Text(0.5, 1.0, 'Uniform')
(array([  2.,   0.,   0.,  12.,   0.,   0.,  38.,   0.,   0., 120.,   0.,
         0., 204.,   0.,   0., 239.,   0.,   0., 218.,   0.,   0., 115.,
         0.,   0.,  43.,   0.,   0.,   7.,   0.,   2.]), array([ 0.        ,  0.33333333,  0.66666667,  1.        ,  1.33333333,
        1.66666667,  2.        ,  2.33333333,  2.66666667,  3.        ,
        3.33333333,  3.66666667,  4.        ,  4.33333333,  4.66666667,
        5.        ,  5.33333333,  5.66666667,  6.        ,  6.33333333,
        6.66666667,  7.        ,  7.33333333,  7.66666667,  8.        ,
        8.33333333,  8.66666667,  9.        ,  9.33333333,  9.66666667,
       10.        ]), <BarContainer object of 30 artists>)
Text(0.5, 1.0, 'Binomial')
(array([ 42.,   0.,   0., 156.,   0.,   0., 224.,   0.,   0., 243.,   0.,
         0., 170.,   0.,   0.,  94.,   0.,   0.,  41.,   0.,   0.,  17.,
         0.,   0.,  11.,   0.,   0.,   1.,   0.,   1.]), array([ 0.        ,  0.33333333,  0.66666667,  1.        ,  1.33333333,
        1.66666667,  2.        ,  2.33333333,  2.66666667,  3.        ,
        3.33333333,  3.66666667,  4.        ,  4.33333333,  4.66666667,
        5.        ,  5.33333333,  5.66666667,  6.        ,  6.33333333,
        6.66666667,  7.        ,  7.33333333,  7.66666667,  8.        ,
        8.33333333,  8.66666667,  9.        ,  9.33333333,  9.66666667,
       10.        ]), <BarContainer object of 30 artists>)
Text(0.5, 1.0, 'Poisson')

1.4.2 Probability Density Functions

For a continuous variable, asking “what is the probability of exactly 68.0000 inches?” is a trap - the answer is zero, because there are infinitely many possible values. Instead we describe a continuous distribution with a probability density function (PDF): a curve whose area between two points gives the probability of landing in that range. The whole area under the curve is exactly 1 (something must happen). The normal PDF is the famous bell:

tibble(z = seq(-4, 4, length.out = 400)) |>
  mutate(density = dnorm(z)) |>
  ggplot(aes(x = z, y = density)) +
  geom_line(colour = "steelblue", linewidth = 1) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(x = "z", y = "density") +
  theme_book()
Figure 1.4: The standard normal density. The height of the curve is not a probability - probability is the area under a stretch of it.
* Valid SPSS, but shown rather than run: PSPP (the free engine this book uses
* to execute its SPSS tabs) has not implemented GRAPH /LINE.
*
* The standard normal density across a grid of z values from -4 to 4.
INPUT PROGRAM.
LOOP #i = -400 TO 400.
  COMPUTE z = #i / 100.
  COMPUTE density = PDF.NORMAL(z, 0, 1).
  END CASE.
END LOOP.
END FILE.
END INPUT PROGRAM.
EXECUTE.

GRAPH /LINE(SIMPLE) = VALUE(density) BY z.
using Distributions, Plots
z = range(-4, 4, length = 400)

plot(z, pdf.(Normal(0, 1), z), linewidth = 2, legend = false,
     title = "The standard normal PDF", xlabel = "z", ylabel = "density")
vline!([0], linestyle = :dash, color = :grey)

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

z = np.linspace(-4, 4, 400)

plt.plot(z, norm.pdf(z), color="steelblue", linewidth=2)
plt.axvline(0, linestyle="--", color="grey")
plt.title("The standard normal PDF")
plt.xlabel("z"); plt.ylabel("density")
plt.show()
[<matplotlib.lines.Line2D object at 0x7f46abe33530>]
<matplotlib.lines.Line2D object at 0x7f46aa5c7320>
Text(0.5, 1.0, 'The standard normal PDF')
Text(0.5, 0, 'z')
Text(0, 0.5, 'density')

The height of the curve is density, not probability; probability is the shaded area you sweep out between two z-values. That single idea, probability as area under a density, is the engine of everything inferential to come, starting with the z-distribution. But first we need to pin down the two summaries every distribution demands: its center and its spread.