The Normal Distribution in R


A distribution is the manner in which a set of values are spread across a possible range of values. A common way of visualizing a distribution is a histogram which shows the number of elements, or frequency, within ranges of values:

> x = c(3, 5, 2, 3, 3, 6, 3, 10, 5, 5, 5, 7, 8, 7, 1, 5, 5, 4, 4, 7)
> hist(x)
Example Histogram

The Normal Distribution

One of the most commonly-encountered distributions in statistics is the normal distribution. The distribution of values for a wide variety of measurements of social and environmental phenomena roughly follow a normal distribution. Examples:

Although the formula for the normal distribution is complex, R contains a number of functions that allow analysis of data. The smoothed histogram associated with the normal distribution is popularly known as the bell curve:

> x = seq(-3, 3, 0.1)
> plot(x = x, y = dnorm(x), type="l", bty="n")
The Normal Distribution Bell Curve

The bell curve is a density curve, and the area under the bell curve between a set of values represents the percent of numbers in the distribution between those values.

The middle value of a normal distribution is the mean, and the width of the bell curve is defined by the standard deviation.

Percentages of Values Within A Normal Distribution

The number of standard deviations a value departs from the mean is called the z-score. The Z-score of the mean is zero. For example, if the mean of a distribution is 7 and the standard deviation is 2, a value of 4 has a z-score of -1.5.

> meanval = 7
> stdev = 2
> x = 3:11
> z = (x - meanval) / stdev
> data.frame(x, z)
   x    z
1  3 -2.0
2  4 -1.5
3  5 -1.0
4  6 -0.5
5  7  0.0
6  8  0.5
7  9  1.0
8 10  1.5
9 11  2.0

Normal Distribution Functions

For all major types of distributions, R provides four different functions:

For the normal distribution, these functions are:

> ?dnorm
dnorm(x, mean = 0, sd = 1, log = FALSE)
pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
rnorm(n, mean = 0, sd = 1)

dnorm is the density function for the normal distribution. This was used to draw the figures above.

pnorm is the probability function (the integral of the density function) and returns the percentage of values below the given value:

> z = seq(-3, 3, 0.1)
> plot(z, pnorm(z), type="l")
The Normal Distribution Probability Curve

pnorm() can be used to calculate the percentage of values between a specific range of values.

For example, the normal distribution of height for men in the United States has a mean of 70 inches (5'10") and a standard deviation of four inches. To find the percent of men between 5' (60 inches) and 6' (72 inches) tall in the USA:

> pnorm(72, mean=70, sd=4) - pnorm(60, mean=70, sd=4)
[1] 0.6852528
# approximately 68.5% of American men are between 5' and 6' tall

qnorm() is the quantile function which returns distribution value for a given probability.

Continuing on the example given above, if you want to find the maximum height for 75% of the men in the USA:

> qnorm(0.75, mean=70, sd=4)
[1] 72.69796
# 75% of the men in the USA are 72.7 inches or shorter

rnorm() will generate a pseudo-random sequence of numbers with a normal distribution. This function is useful for creating simulated data for testing. For example, if we wanted to generate 1,000 simulated heights for men in the USA:

> heights = rnorm(1000, mean=70, sd=4)
> mean(heights)
[1] 70.07384
> sd(heights)
[1] 3.880646
> hist(heights)
Example Random Simulation of Heights Using rnorm()

Evaluating Set of Normally Distributed Numbers

When you have a set of numbers, R provides functions for evaluating the distribution of those numbers.

Two fundamental functions are mean() and sd(), which return the mean and standard deviation of a sequence of numbers, respectively. This example is a final point tally for a geography class:

> points = c(452, 437, 542, 447, 411, 484, 392, 317, 502, 526, 455, 
	454, 628, 500, 457, 478, 511, 365, 407, 440, 357, 298, 
	594, 395, 617, 505, 541, 458, 381)
> mean(points)
[1] 460.3793
> sd(points)
[1] 81.33239
> hist(points)
Example Points Histogram

As is obvious from the histogram, the distribution is not a smooth bell. The normal curve is theoretical and complex reality will rarely line up perfectly with that theoretical construct. In this example, there were a few high-performing students, and a number of students who struggled, stretching the left tail of the bell curve out a bit.

The shifting of the mean of a bell curve off-center is called skew. The skewness() function in the e1071 library will evaluate how skewed a curve is. In this example, the positive value for skew fits our perception that the center of the curve is shifted to the right.

> library(e1071)
> skewness(points)
[1] 0.1320943

If you do not have the e1071 library installed on your system, you can install it with:

> install.packages("e1071")

Another way in which a distribution can be distorted away from normal is kurtosis, which indicates whether the middle is stretched up higher or lower than would be in a clean normal distribution. The negative value returned from the kurtosis() function for the distribution above indicates that the distribution is more spread out than it would be if it were perfectly normal.

> kurtosis(points)
[1] -0.4461787

You can also visualize a normal distribution by creating a rank-order bar chart with the point totals displayed from low to high. Since each bar represents 1/count percent of the distribution, the number of bars can be use to create a vector of percentiles, which can then be used with the qnorm() function to generate the expected normal values for each point total.

This type of chart can be useful for showing students how they stand when a class is graded on a normal curve.

> points = points[order(points)]
> x = barplot(points)
> percentile = (1:length(points)) / (length(points) + 1)
> y = qnorm(percentile, mean=mean(points), sd=sd(points))
> lines(x, y, lwd=3, col="red")
Example Rank-Order Bar Plot With A Normal Curve

Tests For Normality

When dealing with a set of numbers where you are not certain if the distribution is actually normal you can use the ks.test() to perform a Kolmogorov-Smirnov test against a normal distribution. In the results a high p-value indicates normality:

> ks.test(points, "pnorm", mean=mean(points), sd=sd(points))

        One-sample Kolmogorov-Smirnov test

data:  points
D = 0.097876, p-value = 0.9186
alternative hypothesis: two-sided

In contrast, a low p-value indicates the distribution is not normal. In this example, the runif() function is used to generate a random sequence of numbers from a uniform distribution, which is not a normal distribution:

> x = runif(1000)
> ks.test(x, "pnorm", mean=mean(x), sd=sd(x))

        One-sample Kolmogorov-Smirnov test

data:  x
D = 0.071493, p-value = 7.269e-05
alternative hypothesis: two-sided
Histogram of a Random Uniform Distribution, Showing What a Non-Normal Distribution Might Look Like

Another test for normality is the Shapiro-Wilks test, and the shapiro.test() function has a simpler set of function parameters than ks.test() if all you want to know is how normal a distribution of numbers is. As with the Kolmogorov-Smirnov test, higher p-values indicate the closer a distribution is to normal.

> shapiro.test(points)

        Shapiro-Wilk normality test

data:  points
W = 0.98271, p-value = 0.9009