Probability in R Exercises: 15 Real-World Practice Problems

Explanation: sample() draws with replacement from {"H","T"} so each of the 5,000 trials is independent and fair. The expression tosses == "H" is a logical vector and mean() of logicals returns the proportion of TRUE entries, which is exactly the Monte Carlo estimate of P(heads). With 5,000 trials the standard error is roughly 1/sqrt(20000) = 0.007, so landing within 0.01 of 0.5 is exactly what the Law of Large Numbers predicts. A common mistake is forgetting replace = TRUE; without it sample() shuffles a 2-element vector and errors as soon as the third draw.

Explanation: replicate() runs the curly-braced expression 20,000 times and collects the results into a vector. Each trial draws four dice rolls and asks any(rolls == 6), which returns TRUE if at least one six appeared. The theoretical answer is 1 - (5/6)^4 = 0.5177, which matches the simulation to three decimals. This problem is historically famous because gambler Antoine Gombaud (the Chevalier de Mere) lost money betting on the same proposition with three rolls and made money with four, which sparked the Pascal-Fermat correspondence that founded probability theory.

Explanation: rbinom(n, size, prob) returns n independent draws, each of which is the number of successes in size Bernoulli trials with probability prob. Here the simulated mean (2.00) matches n*p = 2.0 exactly to two decimals, and the simulated SD matches sqrt(n*p*(1-p)) = 1.40. The proportion of zero-defect batches (about 0.13) matches dbinom(0, 50, 0.04), which is the analytical probability. QA teams use this kind of simulation to set acceptance-sampling thresholds: a stricter cutoff (zero tolerance) only catches 87 percent of bad batches.

Explanation: pbinom(q, size, prob) returns P(X <= q) by default. To get P(X >= 10), pass q = 9 with lower.tail = FALSE, which returns P(X > 9) and is identical to P(X >= 10) because the distribution is integer-valued. A common mistake is passing q = 10, which excludes the boundary value of interest. The expected count is n*p = 8, so "at least 10" is one standard deviation above the mean and a 25 percent chance is the right order of magnitude. This calculation is the basis of statistical power for one-sample proportion tests.

Explanation: The Poisson distribution models counts of rare independent events over a fixed window. ppois(q, lambda) returns P(X <= q), so to get P(X >= 3) you pass q = 2 with lower.tail = FALSE, which is the cleanest way to compute the upper tail without dealing with 1 - ppois(2, 1.5) and the small floating-point loss that comes with it. The expected count is 1.5 per hour, so seeing three or more is about a 19 percent event. Ops teams calibrate paging thresholds against numbers like this to keep alert noise tolerable.

Explanation: R parameterises the Geometric as the number of failures before the first success, so rgeom(n, p) returns non-negative integers and pgeom(k, p) = 1 - (1-p)^(k+1). All abs_diff values are well under 0.005, which matches the rough Monte Carlo error of sqrt(p*(1-p)/n) for n = 10,000. A common pitfall is mixing this up with the "trials including the success" parameterisation used in some textbooks; in that version the support starts at 1 and pgeom would not match. Comparing empirical to theoretical CDFs is the standard sanity check before using a simulation in production analytics.

Explanation: pnorm(q) defaults to mean = 0, sd = 1 and returns P(Z <= q). The classic "95 percent confidence interval" uses plus and minus 1.96 because that interval traps 95 percent of the standard normal mass (here, 0.9500042 to four decimals of precision). A related shortcut: qnorm(0.975) returns 1.959964, which is where the 1.96 round number comes from. Knowing pnorm, qnorm, dnorm, and rnorm is the four-function toolkit you reach for in every hypothesis-test and simulation problem.

Explanation: The Exponential is memoryless, so the wait time from any moment onward has the same distribution regardless of how long you have already waited. P(wait > 8) for rate 0.2 equals exp(-0.2 * 8) = 0.2019, which matches pexp to six decimals. Note the median (3.47 minutes) is well below the mean (5 minutes), a hallmark of right-skewed distributions. A common modelling mistake is treating the Exponential parameter as the mean rather than the rate; R uses the rate convention (rate = 1/mean), and Python's scipy.stats.expon uses the opposite.

Explanation: The empirical quantile reads the 95th-percentile value directly from the sorted data (here 14.9 mph). The theoretical quantile assumes a normal model with the same mean and SD and asks "where is the 95th percentile of THAT distribution" via qnorm(). The two values disagree by about 0.7 mph, which is a soft signal that the wind distribution is not exactly normal (in fact it is slightly right-skewed). Comparing empirical to model-based quantiles is the cheapest goodness-of-fit check you can run before committing to a parametric model.

Explanation: The analytical answer is 1 - prod((365:341)/365) = 0.5687, so the simulation lands within one Monte Carlo standard error. The result surprises most people because the relevant count is "pairs" not "people": with 25 people there are 300 pairs, and each pair has a 1-in-365 chance of matching. The function duplicated() returns a logical vector flagging the second-and-later occurrence of any repeated value, and any() collapses it to a single TRUE/FALSE per trial. The same pattern catches hash collisions and password reuse, which is why the birthday paradox shows up in cryptography.

Explanation: Switching wins exactly two-thirds of the time because Monty's door-opening rule injects information: he is not allowed to reveal the car, so the door he opens is conditional on the contestant's initial pick. The simulation reproduces the 1/3-versus-2/3 split to two decimals. Notice the if guard around sample(): when only one door is available, sample(x, 1) with a length-1 vector would treat it as sample(1:x, 1) and pick a wrong value, which is one of R's classic foot-guns. Always check length() before calling sample() on an unknown-length pool.

Explanation: Titanic is a 4-dimensional contingency array (Class x Sex x Age x Survived). Slicing [, , "Adult", ] drops the Age dimension while keeping the other three, then apply(..., c("Sex","Survived"), sum) collapses Class to leave a 2x2 table of counts. The conditional probability P(survived | sex) is row-normalised: 75 percent for adult women versus 20 percent for adult men, the textbook "women and children first" effect quantified. Conditional probabilities from contingency tables are the foundation of base-rate reasoning and a standard data-quality audit.

Explanation: The theoretical mean of a fair die is 3.5. The cumulative mean walks toward 3.5 as the sample grows, which is the Law of Large Numbers in action. The shrinkage of the gap is order 1/sqrt(n): from a gap of 0.2 at n = 10 to about 0.002 at n = 2000 is the right ballpark. cumsum() is a vectorised running sum and seq_along(x) is the standard way to get 1:length(x) without breaking on length-zero input. The pair is your go-to idiom for online statistics on streaming data.

Explanation: The Beta is the conjugate prior of the Binomial likelihood: if the prior is Beta(a, b) and you observe k successes in n trials, the posterior is Beta(a + k, b + n - k). With a Beta(2,2) prior and 8 heads in 12 flips, the posterior is Beta(10, 6) and its mean (0.625) is between the prior mean (0.5) and the MLE (0.667), pulled by the prior toward 0.5. The 95 percent credible interval (0.38 to 0.84) covers 0.5, so on this evidence alone you cannot reject the fair-coin hypothesis at the 5 percent level. Conjugacy makes Bayesian inference here a one-line update, no MCMC required.

Explanation: Despite sensitivity of 99 percent, the positive predictive value is only 17 percent. With 1,000 screened people you expect about 10 true positives and 50 false positives (5 percent of 990 healthy people), so among 60 positives only 10 actually have the disease. This base-rate trap is the single most common probabilistic reasoning error in clinical decision-making and machine-learning fraud detection alike. Rule of thumb: when prevalence is much smaller than the false-positive rate, PPV will be poor no matter how high sensitivity is.