1. Two different questions

Suppose a test score is represented by a random variable X. If we fix a score x, we can calculate the probability that X is less than or equal to it. If we fix a probability p, we can seek the score that reaches that cumulative level. In one case, the threshold is given; in the other, the probability level is given.

The prefixes d, p, q, and r organize many of the probability-distribution families available in R. For discrete variables, d…() functions return probability masses; for continuous variables, they return densities. A density height is not the probability of observing exactly that value: in the continuous case, probabilities are obtained over intervals. See references [1] and [2].

2. The cumulative distribution function

For every real-valued random variable X, the cumulative distribution function is defined as:

This function is nondecreasing and right-continuous; its limits are 0 as x tends to −∞ and 1 as x tends to +∞. It may be continuous, have jumps, or combine both features. Its existence does not depend on our knowing an elementary expression that can be evaluated by hand.

For a standard normal distribution, the cumulative distribution function is denoted by Φ:

This integral is well defined, although it cannot be expressed using elementary functions. In R, pnorm() evaluates it numerically and qnorm() provides its quantiles. Quantiles are useful because of the question they answer, not because the cumulative distribution function is somehow absent.

The cumulative distribution function starts with a value and measures probability; the quantile starts with a probability and locates a value.

3. The quantile function as a generalized inverse

To include distributions with jumps or constant stretches, we define the lower quantile as:

The symbol inf denotes the infimum: the greatest lower bound of the set of values whose cumulative probability reaches or exceeds p. For a real-valued random variable and a level strictly between 0 and 1, this definition selects the first threshold satisfying the condition.

If F is continuous and strictly increasing, the quantile function coincides with its ordinary inverse and F(Q(p)) = p. With jumps, that equality may fail. The relevant general relationship is:

Here, F(q−) = P(X < q) is the left-hand limit, whereas F(q) = P(X ≤ q) includes any probability mass at the threshold. Thus, in general, a quantile does not leave exactly a proportion p below it: it ensures that the cumulative level is reached or exceeded.

Properties that aid interpretation

• Order: as the level p increases, the quantile does not decrease.
• Units: if X measures scores, incomes, or times, its quantiles retain those units.
• Position: Q(0.5) is a choice of median; Q(0.25) and Q(0.75) are the quantiles associated with the first and third quartiles.
• Convention: when a distribution has several medians, the definition above selects one of them. The convention used should be stated.

The terms percent-point function and inverse cumulative distribution function are also used. A percentile expresses the level p on a scale of one hundred: the 95th percentile corresponds to p = 0.95. For finite samples, its calculation also requires specifying a method.

4. The test-score example: one probability and two tails

Suppose test scores follow, as a working model, a normal distribution with mean μ = 72 and standard deviation σ = 15.2:

The second parameter in this notation is the variance, σ². By contrast, the R argument sd takes the standard deviation directly. For this family:

This identity shows how a standard normal quantile is transformed to the original scale by multiplying it by σ and adding μ. The expression qnorm(p) uses μ = 0 and σ = 1 by default; supplying other parameters specifies a different normal distribution rather than the standard normal distribution.

4.1. Lower tail: the threshold is approximately 84

We seek a score whose lower-tail cumulative probability is 0.7850824. The call and its verification are:

target_p <- 0.7850824
score_mean <- 72
score_sd <- 15.2

lower_threshold <- qnorm(target_p, mean = score_mean,
                    sd = score_sd)
lower_threshold
# Approximately 84.00000012

pnorm(lower_threshold, mean = score_mean, sd = score_sd)
# Approximately 0.7850824

Under the model, approximately 78.50824% of scores are less than or equal to 84. The returned value is very close to 84, but it is not identical in exact arithmetic: the supplied decimal probability has been rounded.

4.2. Upper tail: the threshold is approximately 60

We now seek a value above which that same probability remains. The tail changes while the requested probability level stays fixed:

upper_threshold <- qnorm(target_p, mean = score_mean,
                    sd = score_sd, lower.tail = FALSE)
upper_threshold
# Approximately 59.99999988

pnorm(upper_threshold, mean = score_mean,
      sd = score_sd, lower.tail = FALSE)
# Approximately 0.7850824

pnorm(upper_threshold, mean = score_mean, sd = score_sd)
# Approximately 0.2149176

The interpretation is that approximately 78.50824% of the distribution lies above 60. Its lower-tail cumulative probability is approximately 0.2149176. For this continuous distribution, including or excluding the exact threshold does not change the probability, because every individual point has zero probability mass.

For a normal distribution, symmetry explains why both thresholds are equally far from the mean: 84 − 72 = 72 − 60 = 12, up to numerical rounding. This symmetry is a property of this model and should not be assumed for every distribution.

4.3. Is a normal model appropriate for test scores?

The normal distribution has support on the entire real line. If scores are bounded between 0 and 100, the model assigns some probability to impossible values. With these parameters, the probability above 100 is approximately 3.27%: a discrepancy worth examining, especially when high percentiles are of interest.

p_below_0 <- pnorm(0, mean = score_mean, sd = score_sd)
p_above_100 <- pnorm(100, mean = score_mean,
                     sd = score_sd, lower.tail = FALSE)
c(p_below_0 = p_below_0, p_above_100 = p_above_100)

Its use requires examining the shape of the data and the purpose of the analysis. Depending on the situation, empirical quantiles, a truncated distribution, or a model with bounded support may be more appropriate. The example illustrates the mechanism of inversion; it does not automatically validate normally distributed scores.

5. Syntax, arguments, and numerical precision in R

A valid call with all the main arguments specified explicitly is:

qnorm(p = 0.95, mean = 0, sd = 1,
      lower.tail = TRUE, log.p = FALSE)
# Approximately 1.644854

An explanatory phrase such as “TRUE or FALSE” should not be copied literally as an argument: a specific logical value must be chosen in each call. It is preferable to spell out TRUE and FALSE. The abbreviations T and F can be reassigned and may become ambiguous during a working session.

5.1. Complementary tails

For a normal distribution and an interior probability level p, the upper-tail threshold corresponds to the lower-tail quantile at level 1 − p. However, changing lower.tail returns a different value, not the complementary probability as its output. In extreme tails, it is preferable to supply the probability of the tail of interest directly.

qnorm(log(0.95), log.p = TRUE)
# The same quantile as qnorm(0.95)

qnorm(1 - 1e-20)
# Inf: the subtraction rounds to 1

qnorm(1e-20, lower.tail = FALSE)
# Approximately 9.26234

qnorm(-1000, lower.tail = FALSE, log.p = TRUE)
# Finite quantile for upper-tail probability exp(-1000)

Subtracting a tiny probability from 1 can produce exactly 1 in floating-point representation. Thus, qnorm(1 - 1e-20) may return Inf, whereas the upper-tail call retains a finite threshold. With log.p = TRUE, a probability can be represented by its logarithm even when the ordinary probability is too small to be stored with useful precision.

5.2. Vectorization and limiting cases

qnorm(c(0.025, 0.5, 0.975), mean = 72, sd = 15.2)
# Approximately: 42.20855  72.00000  101.79145

qnorm(c(0, 1))
# -Inf  Inf

For the nondegenerate normal interpretation, we assume sd > 0. R also handles sd = 0 as a point mass and flags negative values of sd as invalid. A probability level outside [0, 1] is also invalid when log.p = FALSE. Options and numerical details are documented in [2].

6. When the distribution has jumps

Consider the number of successes in four independent Bernoulli trials, each with success probability 0.5. For X ~ Binomial(4, 0.5), we obtain:

discrete_quantile <- qbinom(0.5, size = 4, prob = 0.5)
discrete_quantile
# 2

pbinom(discrete_quantile - 1, size = 4, prob = 0.5)
# 0.3125

pbinom(discrete_quantile, size = 4, prob = 0.5)
# 0.6875

qbinom(0.5, size = 4, prob = 0.5, lower.tail = FALSE)
# 2
pbinom(2, size = 4, prob = 0.5, lower.tail = FALSE)
# 0.3125

No value has a cumulative probability exactly equal to 0.5. Up to 1 success, the accumulated probability is 0.3125; at 2 successes, it reaches 0.6875. The quantile at level 0.5 is therefore 2: it is the smallest value that reaches or exceeds the requested level. This is the operational definition of qbinom(), described in [3].

In this example, P(X < 2) = 0.3125, P(X ≤ 2) = 0.6875, P(X > 2) = 0.3125, and P(X ≥ 2) = 0.6875. Including the threshold changes the probability because P(X = 2) = 0.375.

7. Theoretical quantiles and sample quantiles

qnorm() calculates quantiles of a specified normal distribution. quantile() works with a sample and a computational method. Even if the data appear normal, these operations involve different constructions: one depends on the model and its parameters; the other depends on the ordered data and the chosen convention.

Consider six illustrative scores, created solely to demonstrate the computational methods:

scores <- c(58, 65, 70, 72, 80, 91)
levels <- c(0.25, 0.5, 0.75)

quantile(scores, probs = levels, type = 1)
# 25% 50% 75%
#  65  70  80

quantile(scores, probs = levels, type = 7)
#  25%   50%   75%
# 66.25 71.00 78.00

empirical_cdf <- ecdf(scores)
empirical_cdf(70)
# 0.5

7.1. The empirical distribution and its inverse

The empirical cumulative distribution function assigns to each threshold the observed proportion of data less than or equal to it:

The expression 𝟙 denotes an indicator function: it equals 1 when the condition holds and 0 otherwise. For 0 < p ≤ 1, the generalized inverse of this function corresponds to type = 1. Writing the ordered observations as x₍₁₎, …, x_(n), it returns x_(⌈np⌉). For example, with six scores and p = 0.5, the index is 3 and the result is 70.

7.2. Default interpolation

The default R method is type = 7. For a finite sample of size at least two, it places the level p at position h = 1 + (n − 1)p and interpolates linearly between adjacent observations when that position is not an integer. Thus, p = 0.25 gives h = 2.25 and the result is 65 + 0.25(70 − 65) = 66.25.

At level 0.5, interpolation falls between 70 and 72 and returns 71. Both methods follow legitimate conventions; they are not resolving an arithmetic disagreement. The interpolated result need not coincide with an observed score.

R provides nine methods for calculating sample quantiles. It is unnecessary to develop all their algorithms here, but the chosen method should be stated, especially when comparing results across programs. The official documentation presents the nine definitions and their relationship to order statistics; see [4].

7.3. Missing values and reproducibility

scores_with_na <- c(scores, NA_real_)
sum(is.na(scores_with_na))
# 1

quantile(scores_with_na, probs = levels,
         type = 7, na.rm = TRUE)
# 66.25 71.00 78.00

Removing missing values allows the calculation to proceed, but does not explain why they are missing. A study should report the sample size used, the policy for handling missing data, the requested probability level, and the quantile method. If sampling weights or a complex design are involved, the procedure must account for them as well.

8. Quantiles, critical values, and confidence intervals

Quantile functions appear in many inferential procedures because they locate critical values of a reference distribution. However, calling qnorm() or qt() does not, by itself, construct a confidence interval. One must specify the parameter being estimated, the statistic being used, and the assumptions under which its distribution is known or approximated.

8.1. An interval for the mean

If X₁, …, X_n are independent observations from the same normal distribution with mean μ and unknown positive variance, the pivot:

follows a Student t distribution with n − 1 degrees of freedom. Inverting the probability inequality for this pivot gives the two-sided confidence interval:

In this formula, the sample mean is the estimator; S/√n is its estimated standard error; and the quantile t_{1−α/2,n−1} is the critical value calculated by qt(). See [5] and [6].

scores <- c(58, 65, 70, 72, 80, 91)
n <- length(scores)
confidence <- 0.95
alpha <- 1 - confidence
sample_mean <- mean(scores)
sample_sd <- sd(scores)
standard_error <- sample_sd / sqrt(n)
critical_t <- qt(1 - alpha / 2, df = n - 1)

mean_ci <- sample_mean + c(-1, 1) *
  critical_t * standard_error
mean_ci
# Approximately: 60.50649  84.82684

automatic_ci <- t.test(scores, conf.level = confidence)$conf.int
unname(automatic_ci)

stopifnot(isTRUE(all.equal(
  as.numeric(mean_ci), as.numeric(automatic_ci),
  tolerance = 1e-10
)))

For these six illustrative scores, the 95% interval is approximately [60.50649, 84.82684]. The calculation demonstrates the procedure; a sample this small requires careful examination of the plausibility of the normality assumption. For nonnormal populations, approximate use of this interval needs justification and can be problematic with pronounced skewness, heavy tails, dependence, or influential outliers.

8.2. Three objects that should be distinguished

qnorm(c(0.025, 0.975), mean = 72, sd = 15.2)
# Approximately: 42.20855  101.79145
# Central interval of the model, not an interval for its mean

For a known normal distribution with these parameters, the central 95% interval is approximately [42.20855, 101.79145]. It is not a confidence interval for μ. The upper limit also illustrates why caution is necessary when using a normal distribution for scores whose actual maximum is 100.

8.3. Critical values and p-values

In a two-sided test concerning a mean, the critical value is obtained by setting a level α before evaluating the statistic. The p-value is calculated from the observed statistic and its distribution under the null hypothesis. The letter p in the arguments of qnorm() or qt() denotes an input probability; it does not automatically mean “the p-value of a test.”

null_mean <- 72
observed_t <- (sample_mean - null_mean) / standard_error
two_sided_p <- 2 * pt(abs(observed_t), df = n - 1,
                      lower.tail = FALSE)

c(observed_t = observed_t,
  critical_t = critical_t, two_sided_p = two_sided_p)

abs(observed_t) > critical_t
# FALSE
two_sided_p < alpha
# FALSE

In this example, both decisions agree: μ = 72 is not rejected at level 0.05. This does not prove that the null hypothesis is true or that a substantive difference is impossible. The calculation uses a symmetric reference distribution; other tests require the appropriate tail, statistic, and distribution to be established.

9. Simulation through the quantile function

The quantile function also allows realizations from a distribution to be constructed from a uniform variable. If U is uniformly distributed on (0, 1), then Q_X(U) has cumulative distribution function F_X. Defining the quantile through a generalized inverse extends this principle to discrete distributions as well.

set.seed(2021)
uniform_values <- runif(1000)
simulated_scores <- qnorm(uniform_values, mean = 72, sd = 15.2)

mean(simulated_scores)
sd(simulated_scores)
# Sample approximations to 72 and 15.2

# Normal generation using the dedicated function:
set.seed(2021)
direct_scores <- rnorm(1000, mean = 72, sd = 15.2)
# Element-by-element equality between the two methods is not expected

This is the inverse transform principle. A seed makes the sequence reproducible within a compatible random-number-generator configuration; it does not eliminate sampling variation or make the observed mean equal to the exact parameter. From the user’s perspective, rnorm() provides direct normal generation, but this example does not describe every internal algorithm through which R can generate it.

This application makes the article’s central relationship visible: a collection of probability levels is transformed into values of a variable. Beyond simulation, the same relationship appears in descriptive percentiles, decision thresholds, and diagnostics comparing empirical quantiles with theoretical quantiles.

10. Reproducing the plots in R

The plots embedded in this article represent the mathematical functions used in the examples. The following block produces their equivalents using base R tools: two density areas, a cumulative distribution curve, a quantile curve, and a discrete cumulative distribution. No additional packages need to be installed.

plot_quantiles <- function() {
  previous_parameters <- par(no.readonly = TRUE)
  on.exit(par(previous_parameters))

  paper <- "#fffdf7"
  ink <- "#2c2418"
  copper <- "#8b4513"
  area <- "#e4d1bb"
  mu <- 72
  sigma <- 15.2
  p <- 0.7850824
  q_left <- qnorm(p, mean = mu, sd = sigma)
  q_right <- qnorm(p, mean = mu, sd = sigma,
                 lower.tail = FALSE)
  grid <- seq(12, 132, length.out = 1201)
  density_values <- dnorm(grid, mean = mu, sd = sigma)

  par(mfrow = c(2, 1), mar = c(4, 4, 3, 1),
      bg = paper, fg = ink, col.axis = ink,
      col.lab = ink, col.main = ink)

  for (tail_choice in c("lower", "upper")) {
    cutoff <- if (tail_choice == "lower") q_left else q_right
    bounds <- if (tail_choice == "lower") {
      c(12, cutoff)
    } else {
      c(cutoff, 132)
    }
    xs <- seq(bounds[1], bounds[2], length.out = 600)
    ys <- dnorm(xs, mean = mu, sd = sigma)
    plot(grid, density_values, type = "n",
         xlab = "Score", ylab = "Density",
         main = paste("Tail", tail_choice, "- probability", p))
    polygon(c(xs[1], xs, tail(xs, 1)), c(0, ys, 0),
            col = area, border = NA)
    lines(grid, density_values, col = copper, lwd = 2)
    segments(cutoff, 0, cutoff,
             dnorm(cutoff, mean = mu, sd = sigma), col = copper)
    abline(v = mu, col = "#c4956a", lty = 3)
    mtext(sprintf("Threshold: %.2f", cutoff), side = 3,
          line = 0.2, cex = 0.8)
  }

  par(mfrow = c(1, 1))
  plot(grid, pnorm(grid, mean = mu, sd = sigma),
       type = "l", col = copper, lwd = 2, ylim = c(0, 1),
       xlab = "Score", ylab = "Cumulative probability",
       main = "From p to threshold: cumulative distribution")
  segments(12, p, q_left, p, lty = 2, col = copper)
  segments(q_left, 0, q_left, p, lty = 2, col = copper)
  points(q_left, p, pch = 19, col = copper)

  plot_levels <- seq(0.005, 0.995, length.out = 1000)
  plot(plot_levels,
       qnorm(plot_levels, mean = mu, sd = sigma),
       type = "l", col = copper, lwd = 2,
       xlab = "Cumulative probability p", ylab = "Quantile",
       main = "Normal quantile function")
  points(c(0.5, p), c(mu, q_left), pch = 19, col = copper)

  plot(NA_real_, NA_real_, type = "n", xlim = c(-0.5, 4.5),
       ylim = c(0, 1), xlab = "Number of successes",
       ylab = "Cumulative probability",
       main = "Binomial(4, 0.5): cumulative distribution with jumps")
  abline(h = 0.5, lty = 2, col = "#c4956a")
  segments(-0.5, 0, 0, 0, col = copper, lwd = 2)
  for (k in 0:4) {
    cumulative <- pbinom(k, size = 4, prob = 0.5)
    endpoint <- min(k + 1, 4.5)
    segments(k, cumulative, endpoint, cumulative,
             col = copper, lwd = 2)
    points(k, cumulative, pch = 19, col = copper)
    if (k < 4) {
      points(k + 1, cumulative, pch = 21,
             col = copper, bg = paper)
    }
  }
  segments(2, 0, 2, pbinom(2, 4, 0.5), lty = 2, col = copper)
  invisible(NULL)
}

plot_quantiles()
# In RStudio, use the Plots pane arrows to view the plots.
# To save all pages in a PDF from R:
# pdf("quantile_plots.pdf", width = 8, height = 7)
# plot_quantiles()
# dev.off()

Shaded regions under the densities correspond to probabilities; the heights are densities. In the cumulative distribution plot, height represents cumulative probability. In the quantile curve, the roles of the axes change: the probability level becomes the input and the variable’s value becomes the output.

11. Common mistakes

Interpreting any result requires recording four elements: the variable and its units, the distribution or sample used, the requested probability level, and the computational convention. If the result is used for inference, the parameter, statistic, and assumptions supporting the procedure must also be made explicit.

A quantile turns a proportion into a position within a distribution. Its interpretation depends on which distribution is being studied and how that position is defined.

12. Complete script and checks

The expandable section collects the examples and plots in execution order. It can be copied to an .R file and run in a fresh session. The six test scores are illustrative and do not come from an empirical study. The last printed decimal places may vary with the version and platform.

# THE QUANTILE FUNCTION IN R
# JOSÉ MAURICIO GÓMEZ JULIÁN
# Illustrative data; base R and stats functions.

# LOWER TAIL AND VERIFICATION
target_p <- 0.7850824
score_mean <- 72
score_sd <- 15.2

lower_threshold <- qnorm(target_p, mean = score_mean,
                    sd = score_sd)
lower_threshold
# Approximately 84.00000012

pnorm(lower_threshold, mean = score_mean, sd = score_sd)
# Approximately 0.7850824

# UPPER TAIL AND VERIFICATION
upper_threshold <- qnorm(target_p, mean = score_mean,
                    sd = score_sd, lower.tail = FALSE)
upper_threshold
# Approximately 59.99999988

pnorm(upper_threshold, mean = score_mean,
      sd = score_sd, lower.tail = FALSE)
# Approximately 0.7850824

pnorm(upper_threshold, mean = score_mean, sd = score_sd)
# Approximately 0.2149176

# PROBABILITY OUTSIDE THE 0–100 SCALE
p_below_0 <- pnorm(0, mean = score_mean, sd = score_sd)
p_above_100 <- pnorm(100, mean = score_mean,
                     sd = score_sd, lower.tail = FALSE)
c(p_below_0 = p_below_0, p_above_100 = p_above_100)

# COMPLETE CALL
qnorm(p = 0.95, mean = 0, sd = 1,
      lower.tail = TRUE, log.p = FALSE)
# Approximately 1.644854

# LOGARITHMS AND EXTREME TAILS
qnorm(log(0.95), log.p = TRUE)
# The same quantile as qnorm(0.95)

qnorm(1 - 1e-20)
# Inf: the subtraction rounds to 1

qnorm(1e-20, lower.tail = FALSE)
# Approximately 9.26234

qnorm(-1000, lower.tail = FALSE, log.p = TRUE)
# Finite quantile for upper-tail probability exp(-1000)

# MULTIPLE LEVELS AND ENDPOINTS
qnorm(c(0.025, 0.5, 0.975), mean = 72, sd = 15.2)
# Approximately: 42.20855  72.00000  101.79145

qnorm(c(0, 1))
# -Inf  Inf

# BINOMIAL QUANTILE AND ATTAINED PROBABILITIES
discrete_quantile <- qbinom(0.5, size = 4, prob = 0.5)
discrete_quantile
# 2

pbinom(discrete_quantile - 1, size = 4, prob = 0.5)
# 0.3125

pbinom(discrete_quantile, size = 4, prob = 0.5)
# 0.6875

qbinom(0.5, size = 4, prob = 0.5, lower.tail = FALSE)
# 2
pbinom(2, size = 4, prob = 0.5, lower.tail = FALSE)
# 0.3125

# COMPARING TWO SAMPLE-QUANTILE METHODS
scores <- c(58, 65, 70, 72, 80, 91)
levels <- c(0.25, 0.5, 0.75)

quantile(scores, probs = levels, type = 1)
# 25% 50% 75%
#  65  70  80

quantile(scores, probs = levels, type = 7)
#  25%   50%   75%
# 66.25 71.00 78.00

empirical_cdf <- ecdf(scores)
empirical_cdf(70)
# 0.5

# MISSING VALUES
scores_with_na <- c(scores, NA_real_)
sum(is.na(scores_with_na))
# 1

quantile(scores_with_na, probs = levels,
         type = 7, na.rm = TRUE)
# 66.25 71.00 78.00

# STUDENT T INTERVAL FOR THE MEAN
scores <- c(58, 65, 70, 72, 80, 91)
n <- length(scores)
confidence <- 0.95
alpha <- 1 - confidence
sample_mean <- mean(scores)
sample_sd <- sd(scores)
standard_error <- sample_sd / sqrt(n)
critical_t <- qt(1 - alpha / 2, df = n - 1)

mean_ci <- sample_mean + c(-1, 1) *
  critical_t * standard_error
mean_ci
# Approximately: 60.50649  84.82684

automatic_ci <- t.test(scores, conf.level = confidence)$conf.int
unname(automatic_ci)

stopifnot(isTRUE(all.equal(
  as.numeric(mean_ci), as.numeric(automatic_ci),
  tolerance = 1e-10
)))

# CENTRAL INTERVAL OF THE NORMAL MODEL
qnorm(c(0.025, 0.975), mean = 72, sd = 15.2)
# Approximately: 42.20855  101.79145
# Central interval of the model, not an interval for its mean

# CRITICAL QUANTILE AND TWO-SIDED P-VALUE
null_mean <- 72
observed_t <- (sample_mean - null_mean) / standard_error
two_sided_p <- 2 * pt(abs(observed_t), df = n - 1,
                      lower.tail = FALSE)

c(observed_t = observed_t,
  critical_t = critical_t, two_sided_p = two_sided_p)

abs(observed_t) > critical_t
# FALSE
two_sided_p < alpha
# FALSE

# INVERSE TRANSFORM
set.seed(2021)
uniform_values <- runif(1000)
simulated_scores <- qnorm(uniform_values, mean = 72, sd = 15.2)

mean(simulated_scores)
sd(simulated_scores)
# Sample approximations to 72 and 15.2

# Normal generation using the dedicated function:
set.seed(2021)
direct_scores <- rnorm(1000, mean = 72, sd = 15.2)
# Element-by-element equality between the two methods is not expected

# PLOTS WITH BASE R
plot_quantiles <- function() {
  previous_parameters <- par(no.readonly = TRUE)
  on.exit(par(previous_parameters))

  paper <- "#fffdf7"
  ink <- "#2c2418"
  copper <- "#8b4513"
  area <- "#e4d1bb"
  mu <- 72
  sigma <- 15.2
  p <- 0.7850824
  q_left <- qnorm(p, mean = mu, sd = sigma)
  q_right <- qnorm(p, mean = mu, sd = sigma,
                 lower.tail = FALSE)
  grid <- seq(12, 132, length.out = 1201)
  density_values <- dnorm(grid, mean = mu, sd = sigma)

  par(mfrow = c(2, 1), mar = c(4, 4, 3, 1),
      bg = paper, fg = ink, col.axis = ink,
      col.lab = ink, col.main = ink)

  for (tail_choice in c("lower", "upper")) {
    cutoff <- if (tail_choice == "lower") q_left else q_right
    bounds <- if (tail_choice == "lower") {
      c(12, cutoff)
    } else {
      c(cutoff, 132)
    }
    xs <- seq(bounds[1], bounds[2], length.out = 600)
    ys <- dnorm(xs, mean = mu, sd = sigma)
    plot(grid, density_values, type = "n",
         xlab = "Score", ylab = "Density",
         main = paste("Tail", tail_choice, "- probability", p))
    polygon(c(xs[1], xs, tail(xs, 1)), c(0, ys, 0),
            col = area, border = NA)
    lines(grid, density_values, col = copper, lwd = 2)
    segments(cutoff, 0, cutoff,
             dnorm(cutoff, mean = mu, sd = sigma), col = copper)
    abline(v = mu, col = "#c4956a", lty = 3)
    mtext(sprintf("Threshold: %.2f", cutoff), side = 3,
          line = 0.2, cex = 0.8)
  }

  par(mfrow = c(1, 1))
  plot(grid, pnorm(grid, mean = mu, sd = sigma),
       type = "l", col = copper, lwd = 2, ylim = c(0, 1),
       xlab = "Score", ylab = "Cumulative probability",
       main = "From p to threshold: cumulative distribution")
  segments(12, p, q_left, p, lty = 2, col = copper)
  segments(q_left, 0, q_left, p, lty = 2, col = copper)
  points(q_left, p, pch = 19, col = copper)

  plot_levels <- seq(0.005, 0.995, length.out = 1000)
  plot(plot_levels,
       qnorm(plot_levels, mean = mu, sd = sigma),
       type = "l", col = copper, lwd = 2,
       xlab = "Cumulative probability p", ylab = "Quantile",
       main = "Normal quantile function")
  points(c(0.5, p), c(mu, q_left), pch = 19, col = copper)

  plot(NA_real_, NA_real_, type = "n", xlim = c(-0.5, 4.5),
       ylim = c(0, 1), xlab = "Number of successes",
       ylab = "Cumulative probability",
       main = "Binomial(4, 0.5): cumulative distribution with jumps")
  abline(h = 0.5, lty = 2, col = "#c4956a")
  segments(-0.5, 0, 0, 0, col = copper, lwd = 2)
  for (k in 0:4) {
    cumulative <- pbinom(k, size = 4, prob = 0.5)
    endpoint <- min(k + 1, 4.5)
    segments(k, cumulative, endpoint, cumulative,
             col = copper, lwd = 2)
    points(k, cumulative, pch = 19, col = copper)
    if (k < 4) {
      points(k + 1, cumulative, pch = 21,
             col = copper, bg = paper)
    }
  }
  segments(2, 0, 2, pbinom(2, 4, 0.5), lty = 2, col = copper)
  invisible(NULL)
}

plot_quantiles()
# In RStudio, use the Plots pane arrows to view the plots.
# To save all pages in a PDF from R:
# pdf("quantile_plots.pdf", width = 8, height = 7)
# plot_quantiles()
# dev.off()

# Numerical checks with tolerance
stopifnot(
  abs(lower_threshold - 84) < 1e-5,
  abs(upper_threshold - 60) < 1e-5,
  abs(pnorm(lower_threshold, score_mean, score_sd) -
        target_p) < 1e-10,
  abs(pnorm(upper_threshold, score_mean, score_sd,
             lower.tail = FALSE) - target_p) < 1e-10,
  discrete_quantile == 2,
  pbinom(discrete_quantile - 1, 4, 0.5) < 0.5,
  pbinom(discrete_quantile, 4, 0.5) >= 0.5,
  isTRUE(all.equal(
    as.numeric(quantile(scores, levels, type = 1)),
    c(65, 70, 80))),
  isTRUE(all.equal(
    as.numeric(quantile(scores, levels, type = 7)),
    c(66.25, 71, 78)))
)

sessionInfo()

13. References and documentation

The following sources document the functions, conventions, and procedures used here. The mathematical explanations and examples are developed in the article itself; these links provide access to implementation details.

1. R Core Team. Distributions in the “stats” Package. Official R manual. The d, p, q, and r families.
2. R Core Team. The Normal Distribution. Official R manual. dnorm, pnorm, qnorm, and rnorm.
3. R Core Team. The Binomial Distribution. Official R manual. The binomial distribution and its quantile definition.
4. R Core Team. Sample Quantiles. Official R manual. quantile and its nine methods. The documentation cites Hyndman, R. J. and Fan, Y. (1996), Sample Quantiles in Statistical Packages, as a reference for this classification.
5. R Core Team. The Student t Distribution. Official R manual. dt, pt, qt, and rt.
6. R Core Team. Student’s t-Test. Official R manual. t.test and confidence intervals.