R Probability Distributions Exercises: 20 Real-World Practice Problems

Explanation: dnorm(x, mean, sd) returns the height of the bell curve at x, which is a density (per inch), not a probability. The value 0.1065 does not mean "10.65% of men are 72 inches tall"; it means the density there is 0.1065 per inch. To get an actual probability you integrate the density over an interval, which is what pnorm() does in the next exercise.

Explanation: pnorm() is the cumulative distribution function: pnorm(q) returns the area under the curve to the left of q. The probability of landing inside [a, b] is therefore pnorm(b) - pnorm(a). Both spec limits are 1.33 standard deviations from the mean, so this is the classic pnorm(1.33) - pnorm(-1.33) calculation; about 82% of widgets are in spec.

Explanation: qnorm(p, mean, sd) is the inverse of pnorm(); it returns the value below which a fraction p of the distribution falls. The 95th percentile of N(0, 1) is the famous 1.6449, and the answer here is just 40 + 4 * 1.6449. This same call is how confidence intervals get their z critical values: qnorm(0.975) gives the 1.96 you see in every textbook 95% CI formula.

Explanation: rnorm(n) with no mean or sd defaults to N(0, 1), the standard normal. set.seed() fixes R's pseudo-random generator state so the same seed always returns the same draws; this matters enormously for reproducible simulations and unit tests. Note that calling rnorm() again without re-seeding will return the next 10 numbers in the stream, not the same 10.

Explanation: dbinom(x, size, prob) returns the probability mass at exactly x successes in size trials. The formula is choose(20, 7) * 0.3^7 * 0.7^13, but you should never compute it by hand because dbinom handles the large factorials in log space and stays accurate. For "exactly k" questions use dbinom; for "at most k" or "more than k" questions use pbinom, shown in the next exercise.

Explanation: pbinom(q, size, prob) returns P(X <= q), the cumulative probability up to and including q successes. About 86% of shifts come in at or below the 3-complaint target. A common slip is treating "no more than 3" as pbinom(2, ...); the upper bound is inclusive, so 3 is pbinom(3, ...). For a strict less-than-3 you would use pbinom(2, ...).

Explanation: pbinom(q, ..., lower.tail = FALSE) returns P(X > q), the right tail strictly above q. To get P(X >= 65) you ask for the tail above 64, not above 65, because the bound is exclusive when lower.tail = FALSE. Writing it this way avoids the rounding error you would get from 1 - pbinom(64, 100, 0.6); the lower.tail = FALSE path is numerically stable in the far tails too.

Explanation: dpois(x, lambda) returns the Poisson probability mass at exactly x events when the rate parameter is lambda. The mean and variance of a Poisson are both lambda, so a sample of 8 lives a bit above the mean of 5 and gets about a 6.5% probability. Poisson is the natural model for counts in a fixed window when events are independent and arrive at a constant average rate.

Explanation: ppois(q, lambda, lower.tail = FALSE) returns P(X > q). So "more than 15" is exactly ppois(15, 12, lower.tail = FALSE). About 16% of hours will see more than 15 customers; that headroom drives the staffing buffer. If the question were "15 or more", you would query ppois(14, 12, lower.tail = FALSE) since the bound is strict.

Explanation: A Poisson rate is additive over time: 0.2 per day for 30 days gives lambda = 6 for the month. The window scaling is the most common Poisson gotcha; if you keep lambda = 0.2 and ask about 3 events, you get an absurdly tiny probability because 0.2 is a per-day rate, not a per-month rate. Always confirm the rate and the question are on the same time scale.

Explanation: qt(p, df) returns the quantile of the t distribution. For a two-sided 0.05-level test, the upper rejection cutoff sits at the 0.975 quantile because we split alpha evenly into both tails. Note 2.064 > 1.96; the t critical value is always larger than the normal one because the t has heavier tails. As df grows, qt(0.975, df) converges to 1.96 from above.

Explanation: The two-sided p-value is 2 * P(T <= -|t|), i.e. the lower-tail area at -|t| doubled. Using -abs(t_stat) makes the formula work whether the input statistic is positive or negative. Don't write 2 * (1 - pt(t_stat, df)); that gives the right answer for positive t only and is numerically worse in the tails than the -abs(t) form.

Explanation: Chi-squared tests are one-sided in the upper tail; the entire alpha goes to the right, so the cutoff is qchisq(1 - alpha, df). You will reject the null whenever the observed test statistic exceeds 11.07. Degrees of freedom for a goodness-of-fit test is (number of categories) - 1 - (parameters estimated); with 6 categories and no estimated parameters, df = 5.

Explanation: Like chi-squared, the F test is one-sided in the upper tail because the F statistic is a ratio of variances and large values argue against the null. With p ~ 0.012 you would reject the null of equal group means at alpha = 0.05. The two df arguments are not interchangeable: df1 is the numerator (between-group) df and df2 is the denominator (within-group) df. Swapping them gives a different and wrong p-value.

Explanation: The Beta distribution lives on [0, 1] and is the natural choice for proportions and probabilities. Its mean is shape1 / (shape1 + shape2), so Beta(2, 5) has mean 2/7. With 1000 draws the sample mean is within ~0.003 of the theoretical mean, illustrating the law of large numbers. Increase n to 100000 and the sample mean tightens further around 0.2857.

Explanation: The CLT says the sampling distribution of the mean converges to a normal even when the parent distribution is skewed. Exponential(1) is heavily right-skewed with mean 1 and sd 1, yet the distribution of 5000 sample means (each from n = 30) is nearly normal with mean 1 and sd 1/sqrt(30) = 0.183. replicate() is the idiomatic way to run many independent simulations; for very large simulations swap to vapply(seq_len(M), ..., numeric(1)) for type safety.

Explanation: The empirical proportion mean(x > 2) is a Monte Carlo estimate of P(X > 2). With 100000 draws the Monte Carlo standard error is roughly sqrt(p * (1 - p) / n) = sqrt(0.0227 * 0.9773 / 1e5) = 0.00047, so agreement to about three decimals is exactly what you expect. This kind of cross-check is gold for catching code bugs: if the simulation disagreed with pnorm by more than a few standard errors, the simulation is wrong, not pnorm.

Explanation: For "at least 48" with lower.tail = FALSE, you query the right tail above 47, since P(X >= 48) = P(X > 47). About 22% of batches will hit the 48-out-of-50 bar at a 92% per-wafer yield. To raise that probability you either improve per-wafer yield or accept a more lenient threshold; both can be re-evaluated with one tweak to this call.

Explanation: The Wald formula plugs the sample proportion straight into a normal approximation. The 1.96 multiplier is qnorm(0.975), the 97.5th percentile of N(0, 1); using qnorm instead of a hard-coded 1.96 lets you change the confidence level (say to 90% or 99%) by editing one number. For very small samples or proportions near 0 or 1, the Wald interval underperforms; consider Wilson or Clopper-Pearson (binom.test()), but the Wald formula is the workhorse you will see every day.

Explanation: The defect probability is the sum of the two tail areas outside the spec, not the area inside; resist the urge to write 1 - pnorm(hi, mu, sd) because that ignores the lower tail. Note the asymmetry: with the process mean offset to 25.05 (5 hundredths above target), most of the 0.73% defect probability comes from the upper tail. Centering the process at 25.00 would push the defect rate below 0.5%, which is the kind of result that justifies a recalibration.