Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

1Objectives

By the end of this lesson you will be able to:

  1. Find a Normal-model probability (an area under the curve) with xpnorm(), and read the picture it draws.

  2. Find a percentile/cutoff value with xqnorm().

  3. Work with a simple discrete probability distribution (like dbinom() for the binomial model).

  4. Run a small simulation with do(n) * to see a long-run probability, instead of only computing it.

Everything below runs the same way on your own laptop or on the CSUB JupyterHub (https://csub.jupyter.cal-icor.org/) — same functions, same output, same pictures.

2Setup

library(mosaic)
library(BSDA)

31. xpnorm(): area under the Normal curve, with a picture

Your calculator’s normalcdf finds the area under a Normal curve between two values. xpnorm() does the same job — and draws the curve with the area shaded in, so you see the answer, not just compute it.

Suppose IQ scores follow a Normal model with mean 100 and standard deviation 15 (a standard textbook setup). What fraction of people score 115 or below?

xpnorm(115, mean = 100, sd = 15)

If X ~ N(100, 15), then 

	P(X <= 115) = P(Z <= 1) = 0.8413
	P(X >  115) = P(Z >  1) = 0.1587

[1] 0.8413447

What you’ll see: a bell curve for N(100, 15) with the region to the left of 115 shaded, a vertical line at 115 labeled z = 1, and a small unshaded region in the right tail. The printed text gives you both halves for free — the area below 115 (0.8413) and the area above it (0.1587) — and the number returned, 0.8413447, is the same value as P(X <= 115), just with more decimal places.

Try a value below the mean — the score 80:

xpnorm(80, mean = 100, sd = 15)

If X ~ N(100, 15), then 

	P(X <= 80) = P(Z <= -1.333) = 0.09121
	P(X >  80) = P(Z >  -1.333) = 0.9088

[1] 0.09121122

About 9.1% of people score 80 or below — xpnorm() handles negative z-scores exactly the same way, no special steps.

xpnorm() also accepts two values at once, splitting the curve into three shaded regions — handy for a “between” question, like what fraction of people score between 85 and 115?

xpnorm(c(85, 115), mean = 100, sd = 15)

If X ~ N(100, 15), then 

	P(X <=  85) = P(Z <= -1) = 0.1587	P(X <= 115) = P(Z <=  1) = 0.8413
	P(X >   85) = P(Z >  -1) = 0.8413	P(X >  115) = P(Z >   1) = 0.1587

[1] 0.1586553 0.8413447

What you’ll see: the curve now shaded in three colors — a left tail below 85, a middle region between 85 and 115, and a right tail above 115, labeled with their probabilities (0.1587, 0.6827, 0.1587). That middle probability, 0.6827, should look familiar: it’s the “68” in the 68–95–99.7 rule, because 85 and 115 are exactly one SD below and above the mean of 100.

42. xqnorm(): going the other direction

xpnorm() starts with a value and finds a probability. xqnorm() (your calculator’s invNorm) starts with a probability and finds the value — the cutoff score for a given percentile.

What IQ score marks the 90th percentile (the score that beats 90% of people)?

xqnorm(0.90, mean = 100, sd = 15)

If X ~ N(100, 15), then 

	P(X <= 119.2233) = 0.9
	P(X >  119.2233) = 0.1

[1] 119.2233

What you’ll see: the same shaded-curve picture as xpnorm(), but now the vertical line is drawn at the cutoff value (119.22) instead of at a value you supplied — with 90% of the area shaded below it and 10% above. A score of about 119.2 sits right at the 90th percentile.

53. The binomial model: dbinom() and pbinom()

Not every random process is continuous like the Normal model. The binomial model counts successes in a fixed number of yes/no trials — like guessing on a 10-question true/false quiz, each guess independent, each with a 50% chance of being right.

dbinom(x, size, prob) gives the probability of exactly x successes. What’s the probability of guessing exactly 7 out of 10 correct?

dbinom(7, size = 10, prob = 0.5)
[1] 0.1171875

About a 11.7% chance of exactly 7 correct. You can get the whole distribution at once by feeding dbinom() a vector of every possible outcome, 0 through 10:

round(dbinom(0:10, size = 10, prob = 0.5), 4)
 [1] 0.0010 0.0098 0.0439 0.1172 0.2051 0.2461 0.2051 0.1172 0.0439 0.0098
[11] 0.0010

Eleven numbers (for 0 through 10 correct), symmetric around the middle, peaking at 5 correct (24.6%) — exactly what you’d expect when guessing is a coin flip.

pbinom(x, size, prob) gives the cumulative probability — x or fewer successes — matching your calculator’s binomcdf:

pbinom(6, size = 10, prob = 0.5)
[1] 0.828125

About an 82.8% chance of 6 or fewer correct. To flip that into “7 or more correct” (the upper tail), subtract from 1:

1 - pbinom(6, size = 10, prob = 0.5)
[1] 0.171875

64. Simulation: watching probability happen with do(n) *

Every probability above was computed with a formula. mosaic also lets you watch probability play out, by repeating a random process many times and tallying the results. rflip(10) flips a fair coin 10 times; do(1000) * repeats that entire experiment 1000 times and stacks up the results into a data table:

set.seed(1209)
sims <- do(1000) * rflip(10)
head(sims)
   n heads tails prop
1 10     5     5  0.5
2 10     3     7  0.3
3 10     4     6  0.4
4 10     4     6  0.4
5 10     7     3  0.7
6 10     6     4  0.6

set.seed(1209) makes the “random” simulation exactly reproducible — anyone who runs this same code gets these same 1000 rows. Each row is one simulated set of 10 flips; heads counts how many came up heads. Tally up the heads column across all 1000 simulated repetitions:

tally(~ heads, data = sims)
heads
  1   2   3   4   5   6   7   8   9 
  9  54 108 208 232 223 108  48  10 
gf_histogram(~ heads, data = sims, binwidth = 1, fill = "#0072B2", color = "white",
             xlab = "Number of heads in 10 flips",
             ylab = "Number of simulated reps (of 1000)",
             title = "Simulated distribution of heads in 10 coin flips")

What you’ll see: a bell-shaped histogram, peaking at 5 heads, falling off symmetrically toward 0 and 10 — the simulation rediscovers the same shape you computed exactly with dbinom(0:10, size = 10, prob = 0.5) above, just from random repetition instead of a formula.

Now check the simulation against the exact math from Section 3. What fraction of the 1000 simulated repetitions landed on 7 or more heads?

prop(~ heads >= 7, data = sims)
prop_TRUE 
    0.166

The simulation says 16.6%; the exact binomial calculation, 1 - pbinom(6, 10, 0.5), said 17.2%. Close, but not identical — that gap is simulation noise, and it would shrink if you simulated more than 1000 reps. This is the core idea behind simulation-based statistics: random repetition approximates the exact probability, and the more repetitions you run, the closer the approximation gets.

7Summary

8Check your understanding

  1. Adult male height is roughly Normal with mean 70 inches and SD 3 inches. Write the xpnorm() call that finds the probability a random adult male is shorter than 65 inches.

  2. Using the same height model, write the xqnorm() call for the height that marks the top 5% (the 95th percentile).

  3. A multiple-choice quiz has 8 questions, each with a 25% guessing chance (4 choices). Write the dbinom() call for the probability of guessing exactly 3 correct, and the pbinom() call for the probability of 3 or fewer correct.

  4. In your own words, why did the simulated proportion in Section 4 (16.6%) not land on exactly the same number as the exact binomial calculation (17.2%)? What would you expect to happen to that gap if you changed do(1000) to do(100000)?

R Help
L7 · Exploring data
R Help
L9 · CIs & tests