HW #4: Transformation
Dr. Teresa Vasquez, EdD, MBA, MSDS (in progress)
October 01, 2026

This report was compiled from my own Jupyter notebook (HW4.ipynb, included in this folder) with a custom Claude Code script built by Dr. Teresa Vasquez. The math, code, and answers are my worked solutions; the script converts the notebook into this R Markdown file and knits it. Every number below is recomputed by R when the file is knit.

Opening

This post is for a senior data scientist and the people in their lab. Labs transform data all the time: concentrations and dose responses get logged, measurements get squared or square-rooted before a model sees them. The question here is what a transformation does to a distribution. Where do the mean and the median land after you transform, and how much does a percentile estimated from a sample of 100 move around?

I work through one example from start to finish. I start with \(X \sim \text{Gamma}(\text{shape} = 5, \text{scale} = 2)\), then look at \(Y = \log(X)\), \(Y = X^2\), and \(Y = \sqrt{X}\). Each answer is found two ways: analytically (with formulas) and by simulation (with random draws in R). When the two agree, I trust the result.

# Setup: seed, plot size, and black-background plot style
set.seed(6110)
options(repr.plot.width = 12, repr.plot.height = 5)
dark <- function() par(bg = "black", col.axis = "white", col.lab = "white",
                       col.main = "white", fg = "white")
k <- 5; theta <- 2   # always pass scale = 2 by name in the d/p/q/r functions

The distribution

A Gamma with shape \(k = 5\) and scale \(\theta = 2\) can be read as a waiting time: how long until the 5th event when events arrive at a rate of \(1/\theta = 1/2\) per unit of time. It takes only positive values and has a long right tail.

The Gamma has two parameters, and the second one can be given as a scale \(\theta\) or as a rate \(\beta = 1/\theta\). R’s dgamma, pgamma, qgamma, and rgamma assume a rate unless told otherwise, so every call here passes scale = 2 by name. Leave that off and you get a different distribution with no error message.

Part 1 builds the probability density function (PDF), the cumulative distribution function (CDF), the mean, and the median, then plots them.

Part 1: Base distribution

“For the distribution below, generate a figure of the probability density function (PDF) and cumulative distribution function (CDF). Mark the mean and median in the figure.”

Distribution: \(X \sim \text{GAMMA}(\text{shape} = 5, \text{scale} = 2)\)

Step 1. Build the normalizing constant \(\Gamma(k)\,\theta^k\).

\[\Gamma(5) = 4! = 4 \cdot 3 \cdot 2 \cdot 1 = 24, \qquad \theta^k = 2^5 = 32\]

\[\Gamma(k)\,\theta^k = 24 \times 32 = 768\]

gamma(5) * 2^5
## [1] 768

Step 2. Write the probability density function (PDF) with numbers in.

Substitute \(k = 5\), \(\theta = 2\) into the Gamma probability density function (PDF):

\[f_X(x) = \frac{x^{5-1} e^{-x/2}}{\Gamma(5)\,2^5} = \frac{x^{4} e^{-x/2}}{768}, \qquad x > 0.\]

Check at \(x = 10\):

\[f_X(10) = \frac{10^4 e^{-5}}{768} = \frac{10000 \times 0.0067379}{768} = \frac{67.379}{768} = 0.08773.\]

R returns dgamma(10, shape = 5, scale = 2) \(= 0.08773\), so the formula is right.

c(hand = 10^4 * exp(-10 / 2) / 768,
  r    = dgamma(10, shape = 5, scale = 2))
##       hand          r 
## 0.08773368 0.08773368

Step 3. The cumulative distribution function (CDF).

The shape is a whole number, so this Gamma has a closed-form cumulative distribution function (CDF). \(X\) is the waiting time until the 5th event of a Poisson process with rate \(1/\theta = 1/2\). The event \(\{X \le x\}\) means at least 5 events by time \(x\), and the count by time \(x\) is Poisson with mean \(x/2\):

\[F_X(x) = P(X \le x) = 1 - P(\text{fewer than 5 events by } x) = 1 - \sum_{j=0}^{4} e^{-x/2}\frac{(x/2)^j}{j!}, \qquad x > 0.\]

Check at \(x = 10\), where \(x/2 = 5\):

\[F_X(10) = 1 - e^{-5}\left(1 + 5 + \frac{25}{2} + \frac{125}{6} + \frac{625}{24}\right) = 1 - 0.0067379 \times (1 + 5 + 12.5 + 20.8333 + 26.0417)\]

\[= 1 - 0.0067379 \times 65.375 = 1 - 0.44049 = 0.5595.\]

R returns pgamma(10, shape = 5, scale = 2) \(= 0.5595\).

j <- 0:4
c(hand = 1 - sum(exp(-10 / 2) * (10 / 2)^j / factorial(j)),
  r    = pgamma(10, shape = 5, scale = 2))
##      hand         r 
## 0.5595067 0.5595067

Step 4. Mean and variance.

\[E[X] = k\theta = 5 \times 2 = 10\]

\[\operatorname{Var}(X) = k\theta^2 = 5 \times 2^2 = 5 \times 4 = 20\]

The variance is needed again in Part 2 (b), because \(E[X^2] = \operatorname{Var}(X) + (E[X])^2\).

Step 5. Median.

The median \(m\) solves \(F_X(m) = 0.5\):

\[1 - \sum_{j=0}^{4} e^{-m/2}\frac{(m/2)^j}{j!} = 0.5.\]

This has no closed-form solution, so R solves it: qgamma(0.5, shape = 5, scale = 2).

\[\operatorname{median}(X) = 9.3418\]

Check: the median \(9.3418\) is below the mean \(10\), which is what a right-skewed distribution does. Also \(F_X(10) = 0.5595 > 0.5\) from Step 3, so the median must be below 10.

med_x <- qgamma(0.5, shape = 5, scale = 2)
mean_x <- k * theta
c(mean = mean_x, median = med_x, var = k * theta^2)
##      mean    median       var 
## 10.000000  9.341818 20.000000

Step 6. Figure: probability density function (PDF) and cumulative distribution function (CDF) of \(X\) with the mean and median marked.

Left panel: \(f_X(x) = x^4 e^{-x/2}/768\) drawn with dgamma. Right panel: \(F_X(x)\) drawn with pgamma. Red solid line at the mean \(10\); green dotted line at the median \(9.3418\). The long right tail pulls the mean to the right of the median. On the cumulative distribution function (CDF) panel the median line crosses the curve at height \(0.5\) and the mean line crosses at \(0.5595\).

Final answer, Part 1: \(f_X(x) = x^4 e^{-x/2}/768\) for \(x > 0\); mean \(= 10\); variance \(= 20\); median \(= 9.3418\).

dark()
par(mfrow = c(1, 2))
curve(dgamma(x, shape = 5, scale = 2), from = 0, to = 30, lwd = 3, col = "deepskyblue",
      main = "Probability density function (PDF) of X", xlab = "x", ylab = "f(x)")
abline(v = mean_x, col = "red", lwd = 2)
abline(v = med_x, col = "limegreen", lwd = 2, lty = 3)
legend(x = 12, y = 0.095, seg.len = 1.5, c("mean = 10", sprintf("median = %.4f", med_x)),
       col = c("red", "limegreen"), lty = c(1, 3), lwd = 2, text.col = "white", bty = "n")
curve(pgamma(x, shape = 5, scale = 2), from = 0, to = 30, lwd = 3, col = "deepskyblue",
      main = "Cumulative distribution function (CDF) of X", xlab = "x", ylab = "F(x)")
abline(v = mean_x, col = "red", lwd = 2)
abline(v = med_x, col = "limegreen", lwd = 2, lty = 3)
abline(h = 0.5, col = "gray", lty = 2)
legend(x = 12, y = 0.45, seg.len = 1.5, c("mean = 10", sprintf("median = %.4f", med_x)),
       col = c("red", "limegreen"), lty = c(1, 3), lwd = 2, text.col = "white", bty = "n")

Reading the figure. The mean is \(10\) and the median is \(9.3418\). The median sits left of the mean because the long right tail pulls the mean toward it. In the cumulative distribution function (CDF) panel, the median is where the curve crosses \(0.5\). At the mean, the curve is already at \(0.5595\), so about 56 percent of values fall below the mean.

Transformations: two methods, and where the mean and median move

Analytic method. All three transformations (\(\log\), square, square root) are increasing on \(x > 0\). That allows two shortcuts. The CDF carries straight over, \(F_Y(y) = F_X(g^{-1}(y))\). The PDF is the old PDF at \(g^{-1}(y)\) times the Jacobian, \(f_Y(y) = f_X(g^{-1}(y))\,\lvert d g^{-1}(y)/dy \rvert\). The Jacobian corrects for the transformation stretching or squeezing the axis, so the total area stays 1.

Simulation method. Draw \(10^5\) values with rgamma(), apply the transformation to each draw, and estimate the PDF with density() and the CDF with ecdf(). Every figure below puts the simulated curve (blue) under the analytic curve (orange, dashed). When the two lines sit on top of each other, both methods agree.

What happens to the mean and median. The median carries straight through: \(\operatorname{median}(Y) = g(\operatorname{median}(X))\), because an increasing \(g\) keeps every value on the same side of the middle. The mean does not. Jensen’s inequality says that for a concave \(g\) (log, square root), \(E[g(X)] < g(E[X])\), and for a convex \(g\) (squaring) the inequality flips.

Transformation Mean of \(Y\) \(g(E[X]) = g(10)\) Median of \(Y\) = \(g(9.3418)\) Shape after
\(\log(X)\) 2.1993 2.3026 2.2345 mean now below median: the right skew is gone and a slight left skew remains
\(X^2\) 120 100 87.2696 much stronger right skew
\(\sqrt{X}\) 3.0843 3.1623 3.0564 still right-skewed, but less

Part 2: Transformations

“For the distribution below, generate a figure of the probability density function (PDF) and cumulative distribution function (CDF) of each transformation: (a) \(Y = \log(X)\), (b) \(Y = X^2\), and (c) \(Y = \sqrt{X}\) random variable. For each transformation, mark the mean and median in the figure. You are required to use both simulation and analytic methods in order to find the probability density function (PDF) and cumulative distribution function (CDF) of the transformation.”

Simulation side, drawn once for all three transformations: x_sim <- rgamma(1e5, shape = 5, scale = 2).

Simulation draw and plotting helper.

One draw of \(10^5\) values is reused for (a), (b), and (c), so the three transformations are applied to the same simulated \(X\) values. The helper overlays the analytic curve (orange dashed) on the simulated curve (blue) so agreement is visible.

Simulated summaries of x_sim: mean \(= 9.9924\) (analytic \(10\)), median \(= 9.3364\) (analytic \(9.3418\)).

x_sim <- rgamma(1e5, shape = 5, scale = 2)
c(sim_mean = mean(x_sim), sim_median = median(x_sim))
##   sim_mean sim_median 
##   9.992433   9.336437
plot_transform <- function(y, dens, cdf, mu, med, from, to, label) {
  dark(); par(mfrow = c(1, 2))
  d <- density(y, from = from, to = to, n = 1024)
  plot(d, col = "deepskyblue", lwd = 3, xlim = c(from, to), ylim = c(0, 1.15 * max(d$y)),
       main = paste("Probability density function (PDF) of", label), xlab = "y", ylab = "f(y)")
  curve(dens(x), from = from, to = to, n = 1001, add = TRUE, col = "orange", lwd = 2, lty = 2)
  abline(v = mu, col = "red", lwd = 2); abline(v = med, col = "limegreen", lwd = 2, lty = 3)
  legend("topright", c("simulated", "analytic", sprintf("mean = %.4f", mu), sprintf("median = %.4f", med)),
         col = c("deepskyblue", "orange", "red", "limegreen"), lty = c(1, 2, 1, 3), lwd = 2,
         text.col = "white", bty = "n", cex = 0.85)
  plot(ecdf(y), do.points = FALSE, verticals = TRUE, col = "deepskyblue", lwd = 3, xlim = c(from, to),
       main = paste("Cumulative distribution function (CDF) of", label), xlab = "y", ylab = "F(y)")
  curve(cdf(x), from = from, to = to, n = 1001, add = TRUE, col = "orange", lwd = 2, lty = 2)
  abline(v = mu, col = "red", lwd = 2); abline(v = med, col = "limegreen", lwd = 2, lty = 3)
  abline(h = 0.5, col = "gray", lty = 2)
}

Part 2 (a): \(Y = \log(X)\)

Step 1. Monotonic?

\(g(x) = \log x\) has derivative \(g'(x) = 1/x > 0\) for every \(x > 0\), so \(g\) is strictly increasing on the support of \(X\). The cumulative distribution function (CDF) shortcut and the median rule are allowed.

Step 2. Invert, \(x = g^{-1}(y)\).

\[y = \log x \;\Longrightarrow\; x = g^{-1}(y) = e^{y}\]

Step 3. Analytic cumulative distribution function (CDF), \(F_Y(y) = F_X(g^{-1}(y))\).

\[F_Y(y) = P(\log X \le y) = P(X \le e^{y}) = F_X(e^{y})\]

Substitute \(e^y\) into the Step 3 formula from Part 1:

\[F_Y(y) = 1 - \sum_{j=0}^{4} e^{-e^{y}/2}\,\frac{(e^{y}/2)^j}{j!}, \qquad -\infty < y < \infty.\]

In R: pgamma(exp(y), shape = 5, scale = 2).

cdf_a <- function(y) pgamma(exp(y), shape = 5, scale = 2)

Step 4. Jacobian, \(\lvert d g^{-1}(y)/dy \rvert\).

\[\frac{dx}{dy} = \frac{d}{dy} e^{y} = e^{y}, \qquad \left\lvert e^{y} \right\rvert = e^{y} \text{ (already positive)}.\]

Step 5. Analytic probability density function (PDF), \(f_Y(y) = f_X(g^{-1}(y)) \cdot \lvert \text{Jacobian} \rvert\).

\[f_Y(y) = f_X(e^{y}) \cdot e^{y} = \frac{(e^{y})^4\, e^{-e^{y}/2}}{768} \cdot e^{y} = \frac{e^{4y}\, e^{y}\, e^{-e^{y}/2}}{768} = \frac{e^{5y}\, e^{-e^{y}/2}}{768}\]

\[\boxed{f_Y(y) = \frac{\exp\!\left(5y - e^{y}/2\right)}{768}, \qquad -\infty < y < \infty}\]

Check: differentiating the Step 3 cumulative distribution function (CDF) gives \(\frac{d}{dy}F_X(e^y) = f_X(e^y)\,e^y\), the same expression. The code cell also integrates \(f_Y\) numerically over the real line and gets 1.

dens_a <- function(y) exp(5 * y - exp(y) / 2) / 768
integrate(dens_a, lower = -Inf, upper = Inf)
## 1 with absolute error < 1.4e-06

Step 6. Support of \(Y\).

\(x \in (0, \infty)\). As \(x \to 0^{+}\), \(\log x \to -\infty\); as \(x \to \infty\), \(\log x \to \infty\). Values of \(x\) in \((0, 1)\) map to negative \(y\). So

\[\text{support of } Y = (-\infty, \infty), \text{ all real numbers.}\]

Step 7. Median mark, \(\operatorname{median}(Y) = g(\operatorname{median}(X))\).

\[\operatorname{median}(Y) = \log(\operatorname{median}(X)) = \log(9.3418) = 2.2345\]

med_a <- log(med_x)
med_a
## [1] 2.234501

Step 8. Mean mark, via a moment tool.

Use the log moment tool with \(k = 5\), \(\theta = 2\):

\[E[Y] = E[\log X] = \psi(5) + \log 2.\]

For a whole number \(k\), \(\psi(k) = -\gamma + \sum_{j=1}^{k-1} \frac{1}{j}\) with Euler’s constant \(\gamma = 0.5772157\):

\[\psi(5) = -0.5772157 + \left(1 + \tfrac{1}{2} + \tfrac{1}{3} + \tfrac{1}{4}\right) = -0.5772157 + 2.0833333 = 1.5061177\]

\[E[Y] = 1.5061177 + 0.6931472 = 2.1992649 \approx 2.1993\]

Contrast: \(\log(E[X]) = \log 10 = 2.3026\), which is \(0.1033\) too high. Since \(\log\) is concave, Jensen’s inequality gives \(E[\log X] < \log E[X]\).

mu_a <- digamma(5) + log(2)
c(moment_tool = mu_a,
  integral = integrate(function(y) y * dens_a(y), -Inf, Inf)$value)
## moment_tool    integral 
##    2.199265    2.199265

Step 9. Figure: analytic and simulated probability density function (PDF) and cumulative distribution function (CDF), mean and median marked.

The simulated curves (blue) and analytic curves (orange dashed) lie on top of each other. The log pulls the long right tail of \(X\) in, so \(Y\) is close to symmetric with a slightly longer left tail, and the mean \(2.1993\) now sits just left of the median \(2.2345\). Simulated check from y: mean \(= 2.1984\), median \(= 2.2339\).

Final answer, Part 2 (a): \(f_Y(y) = \exp(5y - e^y/2)/768\) on all real \(y\); \(F_Y(y) = F_X(e^y)\); median \(= 2.2345\); mean \(= \psi(5) + \log 2 = 2.1993\).

y <- log(x_sim)
c(sim_mean = mean(y), sim_median = median(y), analytic_mean = mu_a, analytic_median = med_a)
##        sim_mean      sim_median   analytic_mean analytic_median 
##        2.198358        2.233925        2.199265        2.234501
plot_transform(y, dens_a, cdf_a, mu_a, med_a, 0.3, 3.8, "Y = log(X)")

Part 2 (b): \(Y = X^2\)

Step 1. Monotonic?

\(g(x) = x^2\) has derivative \(g'(x) = 2x > 0\) for every \(x > 0\). On the support of \(X\) (positive numbers only) \(g\) is strictly increasing. On the whole real line \(x^2\) is not monotonic, but \(X\) never takes negative values, so that does not matter here.

Step 2. Invert, \(x = g^{-1}(y)\).

\[y = x^2 \;\Longrightarrow\; x = g^{-1}(y) = \sqrt{y} \quad (\text{positive root, because } x > 0)\]

Step 3. Analytic cumulative distribution function (CDF), \(F_Y(y) = F_X(g^{-1}(y))\).

\[F_Y(y) = P(X^2 \le y) = P(X \le \sqrt{y}) = F_X(\sqrt{y})\]

\[F_Y(y) = 1 - \sum_{j=0}^{4} e^{-\sqrt{y}/2}\,\frac{(\sqrt{y}/2)^j}{j!}, \qquad y > 0.\]

In R: pgamma(sqrt(y), shape = 5, scale = 2).

cdf_b <- function(y) pgamma(sqrt(y), shape = 5, scale = 2)

Step 4. Jacobian, \(\lvert d g^{-1}(y)/dy \rvert\).

\[\frac{dx}{dy} = \frac{d}{dy} y^{1/2} = \frac{1}{2}\, y^{-1/2} = \frac{1}{2\sqrt{y}} > 0 \text{ for } y > 0.\]

Step 5. Analytic probability density function (PDF), \(f_Y(y) = f_X(g^{-1}(y)) \cdot \lvert \text{Jacobian} \rvert\).

\[f_Y(y) = f_X(\sqrt{y}) \cdot \frac{1}{2\sqrt{y}} = \frac{(\sqrt{y})^4\, e^{-\sqrt{y}/2}}{768} \cdot \frac{1}{2\, y^{1/2}} = \frac{y^{2}\, e^{-\sqrt{y}/2}}{1536\, y^{1/2}}\]

\[\boxed{f_Y(y) = \frac{y^{3/2}\, e^{-\sqrt{y}/2}}{1536}, \qquad y > 0}\]

Check: the code cell integrates \(f_Y\) over \((0, \infty)\) numerically and gets 1.

dens_b <- function(y) y^(3/2) * exp(-sqrt(y) / 2) / 1536
integrate(dens_b, lower = 0, upper = Inf)
## 1 with absolute error < 8.7e-05

Step 6. Support of \(Y\).

\(x \in (0, \infty)\) and \(x^2\) maps \((0, \infty)\) onto \((0, \infty)\), so the support of \(Y\) is \(y > 0\).

Step 7. Median mark, \(\operatorname{median}(Y) = g(\operatorname{median}(X))\).

\[\operatorname{median}(Y) = (\operatorname{median}(X))^2 = 9.3418^2 = 87.2696\]

(Using the unrounded median \(9.341818\): \(9.341818^2 = 87.26956\).)

med_b <- med_x^2
med_b
## [1] 87.26956

Step 8. Mean mark, via a moment tool.

Use the power moment tool with \(r = 2\):

\[E[Y] = E[X^2] = \theta^2\,\frac{\Gamma(k+2)}{\Gamma(k)} = 2^2 \cdot \frac{\Gamma(7)}{\Gamma(5)} = 4 \cdot \frac{6!}{4!} = 4 \cdot \frac{720}{24} = 4 \times 30 = 120\]

Second check with Part 1: \(E[X^2] = \operatorname{Var}(X) + (E[X])^2 = 20 + 10^2 = 20 + 100 = 120\).

Contrast: \((E[X])^2 = 100\), which is \(20\) too low. The gap is exactly \(\operatorname{Var}(X)\).

mu_b <- 2^2 * gamma(7) / gamma(5)
c(moment_tool = mu_b,
  integral = integrate(function(y) y * dens_b(y), 0, Inf)$value)
## moment_tool    integral 
##         120         120

Step 9. Figure: analytic and simulated probability density function (PDF) and cumulative distribution function (CDF), mean and median marked.

Simulated and analytic curves agree. Squaring stretches the right tail a lot, so \(Y\) is strongly right skewed and the mean \(120\) sits well to the right of the median \(87.2696\). Simulated check from y: mean \(= 119.8962\), median \(= 87.1690\). The simulated mean wobbles more here than in (a) because squaring magnifies the largest draws.

Final answer, Part 2 (b): \(f_Y(y) = y^{3/2} e^{-\sqrt{y}/2}/1536\) for \(y > 0\); \(F_Y(y) = F_X(\sqrt{y})\); median \(= 87.2696\); mean \(= 120\).

y <- x_sim^2
c(sim_mean = mean(y), sim_median = median(y), analytic_mean = mu_b, analytic_median = med_b)
##        sim_mean      sim_median   analytic_mean analytic_median 
##       119.89621        87.16905       120.00000        87.26956
plot_transform(y, dens_b, cdf_b, mu_b, med_b, 0, 700, "Y = X^2")

Part 2 (c): \(Y = \sqrt{X}\)

Step 1. Monotonic?

\(g(x) = \sqrt{x} = x^{1/2}\) has derivative \(g'(x) = \frac{1}{2\sqrt{x}} > 0\) for every \(x > 0\), so \(g\) is strictly increasing on the support of \(X\).

Step 2. Invert, \(x = g^{-1}(y)\).

\[y = \sqrt{x} \;\Longrightarrow\; x = g^{-1}(y) = y^{2}\]

Step 3. Analytic cumulative distribution function (CDF), \(F_Y(y) = F_X(g^{-1}(y))\).

\[F_Y(y) = P(\sqrt{X} \le y) = P(X \le y^2) = F_X(y^2)\]

\[F_Y(y) = 1 - \sum_{j=0}^{4} e^{-y^2/2}\,\frac{(y^2/2)^j}{j!}, \qquad y > 0.\]

In R: pgamma(y^2, shape = 5, scale = 2).

cdf_c <- function(y) pgamma(y^2, shape = 5, scale = 2)

Step 4. Jacobian, \(\lvert d g^{-1}(y)/dy \rvert\).

\[\frac{dx}{dy} = \frac{d}{dy} y^{2} = 2y > 0 \text{ for } y > 0.\]

Step 5. Analytic probability density function (PDF), \(f_Y(y) = f_X(g^{-1}(y)) \cdot \lvert \text{Jacobian} \rvert\).

\[f_Y(y) = f_X(y^2) \cdot 2y = \frac{(y^2)^4\, e^{-y^2/2}}{768} \cdot 2y = \frac{2\, y^{8}\, y\, e^{-y^2/2}}{768} = \frac{y^{9}\, e^{-y^2/2}}{384}\]

\[\boxed{f_Y(y) = \frac{y^{9}\, e^{-y^2/2}}{384}, \qquad y > 0}\]

Check: the code cell integrates \(f_Y\) over \((0, \infty)\) numerically and gets 1.

dens_c <- function(y) y^9 * exp(-y^2 / 2) / 384
integrate(dens_c, lower = 0, upper = Inf)
## 1 with absolute error < 3.9e-05

Step 6. Support of \(Y\).

\(x \in (0, \infty)\) and \(\sqrt{x}\) maps \((0, \infty)\) onto \((0, \infty)\), so the support of \(Y\) is \(y > 0\).

Step 7. Median mark, \(\operatorname{median}(Y) = g(\operatorname{median}(X))\).

\[\operatorname{median}(Y) = \sqrt{\operatorname{median}(X)} = \sqrt{9.3418} = 3.0564\]

med_c <- sqrt(med_x)
med_c
## [1] 3.056439

Step 8. Mean mark, via a moment tool.

Use the power moment tool with \(r = 1/2\):

\[E[Y] = E[X^{1/2}] = \theta^{1/2}\,\frac{\Gamma(k + 1/2)}{\Gamma(k)} = \sqrt{2}\,\frac{\Gamma(5.5)}{\Gamma(5)}\]

Build \(\Gamma(5.5)\) with \(\Gamma(z+1) = z\,\Gamma(z)\) and \(\Gamma(1/2) = \sqrt{\pi}\):

\[\Gamma(5.5) = 4.5 \times 3.5 \times 2.5 \times 1.5 \times 0.5 \times \sqrt{\pi} = 29.53125 \times 1.7724539 = 52.34278\]

\[E[Y] = 1.4142136 \times \frac{52.34278}{24} = 1.4142136 \times 2.1809491 = 3.0843\]

Contrast: \(\sqrt{E[X]} = \sqrt{10} = 3.1623\), which is \(0.0780\) too high. The square root is concave, so Jensen’s inequality gives \(E[\sqrt{X}] < \sqrt{E[X]}\).

mu_c <- sqrt(2) * gamma(5.5) / gamma(5)
c(moment_tool = mu_c,
  integral = integrate(function(y) y * dens_c(y), 0, Inf)$value)
## moment_tool    integral 
##    3.084328    3.084328

Step 9. Figure: analytic and simulated probability density function (PDF) and cumulative distribution function (CDF), mean and median marked.

Simulated and analytic curves agree. The square root compresses the right tail, so \(Y\) is only mildly right skewed and the mean \(3.0843\) sits just right of the median \(3.0564\). Simulated check from y: mean \(= 3.0830\), median \(= 3.0556\).

Final answer, Part 2 (c): \(f_Y(y) = y^9 e^{-y^2/2}/384\) for \(y > 0\); \(F_Y(y) = F_X(y^2)\); median \(= 3.0564\); mean \(= \sqrt{2}\,\Gamma(5.5)/\Gamma(5) = 3.0843\).

y <- sqrt(x_sim)
c(sim_mean = mean(y), sim_median = median(y), analytic_mean = mu_c, analytic_median = med_c)
##        sim_mean      sim_median   analytic_mean analytic_median 
##        3.083002        3.055558        3.084328        3.056439
plot_transform(y, dens_c, cdf_c, mu_c, med_c, 0.5, 6, "Y = sqrt(X)")

The 85th percentile and its sampling distribution

A real lab sees one sample, not the whole distribution. From a sample of 100 values you can compute an 85th percentile, but a different sample of 100 gives a different number. The sampling distribution of the 85th percentile is the distribution of that number across many repeated samples. It shows how far one sample’s estimate can land from the truth.

To build it, I draw 1000 samples of size 100 for each transformation, take quantile(..., probs = 0.85) of each sample, and plot the PDF and CDF of those 1000 values. The true value to compare against comes from the same rule as the median: an increasing \(g\) carries every percentile, so \(q_Y(0.85) = g(q_X(0.85))\).

Part 3: Sampling distribution of the 85th percentile

“For each of the transformation, generate 1000 samples of size 100. Then, generate a figure of the probability density function (PDF) and cumulative distribution function (CDF) of the sampling distribution of the 85th percentile.”

Step 1. What varies, which R call, and one assumption.

  1. The random sample of 100 draws from \(X \sim \text{Gamma}(5, \text{scale } 2)\) changes every run, so the 85th percentile computed from it changes too.
  2. quantile(s, probs = 0.85), where s is the transformed sample (default type = 7, linear interpolation between order statistics).
  3. The 100 values in each sample are independent draws from the same \(\text{Gamma}(\text{shape} = 5, \text{scale} = 2)\) distribution, and the 1000 samples are independent of each other.

Step 2. The true 85th percentiles (targets for the figures).

For \(X\): qgamma(0.85, shape = 5, scale = 2) \(= 14.5339\).

Each transformation is increasing, so the percentile rule works the same way the median rule does, \(q_Y(0.85) = g(q_X(0.85))\):

\[\log(14.5339) = 2.6765, \qquad 14.5339^2 = 211.2353, \qquad \sqrt{14.5339} = 3.8123.\]

q85_x <- qgamma(0.85, shape = 5, scale = 2)
truth <- c(log = log(q85_x), square = q85_x^2, sqrt = sqrt(q85_x))
c(x = q85_x, truth)
##          x        log     square       sqrt 
##  14.533936   2.676486 211.235296   3.812340
one_stat <- function(g) quantile(g(rgamma(100, shape = 5, scale = 2)), probs = 0.85)

plot_sampling <- function(vals, true_q, label) {
  dark(); par(mfrow = c(1, 2))
  plot(density(vals), col = "deepskyblue", lwd = 3, xlab = "85th percentile", ylab = "density",
       main = paste("Probability density function (PDF)\n85th percentile of", label))
  abline(v = true_q, col = "red", lwd = 2); abline(v = mean(vals), col = "limegreen", lwd = 2, lty = 3)
  legend("topright", c(sprintf("true = %.4f", true_q), sprintf("mean of 1000 = %.4f", mean(vals))),
         col = c("red", "limegreen"), lty = c(1, 3), lwd = 2, text.col = "white", bty = "n", cex = 0.85)
  plot(ecdf(vals), do.points = FALSE, verticals = TRUE, col = "deepskyblue", lwd = 3, xlab = "85th percentile",
       main = paste("Cumulative distribution function (CDF)\n85th percentile of", label))
  abline(v = true_q, col = "red", lwd = 2); abline(h = 0.5, col = "gray", lty = 2)
}

Step 3. (a) \(Y = \log(X)\): 1000 samples of size 100, 85th percentile of each, figure.

Across the 1000 samples: mean of the 85th percentiles \(= 2.6691\), standard deviation \(= 0.0599\), middle 95 percent from \(2.5540\) to \(2.7903\). True value \(= 2.6765\). The sampling distribution is close to a symmetric bell. The center sits close to the true value; the small gap is the sample-percentile bias at \(n = 100\) plus simulation noise.

vals_log <- replicate(1000, one_stat(log))
c(mean = mean(vals_log), sd = sd(vals_log), quantile(vals_log, c(0.025, 0.975)), truth = truth[["log"]])
##       mean         sd       2.5%      97.5%      truth 
## 2.66906687 0.05990511 2.55401717 2.79033491 2.67648633
plot_sampling(vals_log, truth[["log"]], "Y = log(X)")

Step 4. (b) \(Y = X^2\): 1000 samples of size 100, 85th percentile of each, figure.

Across the 1000 samples: mean of the 85th percentiles \(= 209.5073\), standard deviation \(= 23.9408\), middle 95 percent from \(166.6572\) to \(261.5192\). True value \(= 211.2353\). The sampling distribution is mildly right skewed with a longer upper tail: squaring stretches large values more than small ones, so it is also much wider in absolute terms. The center sits close to the true value; the small gap is the sample-percentile bias at \(n = 100\) plus simulation noise.

vals_square <- replicate(1000, one_stat(function(v) v^2))
c(mean = mean(vals_square), sd = sd(vals_square), quantile(vals_square, c(0.025, 0.975)), truth = truth[["square"]])
##      mean        sd      2.5%     97.5%     truth 
## 209.50727  23.94083 166.65722 261.51915 211.23530
plot_sampling(vals_square, truth[["square"]], "Y = X^2")

Step 5. (c) \(Y = \sqrt{X}\): 1000 samples of size 100, 85th percentile of each, figure.

Across the 1000 samples: mean of the 85th percentiles \(= 3.7991\), standard deviation \(= 0.1176\), middle 95 percent from \(3.5686\) to \(4.0374\). True value \(= 3.8123\). The sampling distribution is close to a symmetric bell. The center sits close to the true value; the small gap is the sample-percentile bias at \(n = 100\) plus simulation noise.

vals_sqrt <- replicate(1000, one_stat(sqrt))
c(mean = mean(vals_sqrt), sd = sd(vals_sqrt), quantile(vals_sqrt, c(0.025, 0.975)), truth = truth[["sqrt"]])
##      mean        sd      2.5%     97.5%     truth 
## 3.7990520 0.1175712 3.5685570 4.0373733 3.8123400
plot_sampling(vals_sqrt, truth[["sqrt"]], "Y = sqrt(X)")

Step 6. Summary table of the three sampling distributions.

Transformation True 85th percentile Mean of 1000 Standard deviation 2.5th percentile 97.5th percentile
\(\log(X)\) 2.6765 2.6691 0.0599 2.5540 2.7903
\(X^2\) 211.2353 209.5073 23.9408 166.6572 261.5192
\(\sqrt{X}\) 3.8123 3.7991 0.1176 3.5686 4.0374

Assumptions: independent and identically distributed draws within and across samples; quantile() type 7; set.seed(6110) at the top with the order log, square, square root. A different seed or order shifts these simulated numbers slightly.

summ <- function(v) c(mean = mean(v), sd = sd(v), quantile(v, c(0.025, 0.975)))
round(rbind(log = c(truth = truth[["log"]], summ(vals_log)),
            square = c(truth = truth[["square"]], summ(vals_square)),
            sqrt = c(truth = truth[["sqrt"]], summ(vals_sqrt))), 4)
##           truth     mean      sd     2.5%    97.5%
## log      2.6765   2.6691  0.0599   2.5540   2.7903
## square 211.2353 209.5073 23.9408 166.6572 261.5192
## sqrt     3.8123   3.7991  0.1176   3.5686   4.0374

Assumptions

What the lab should watch for when transforming data

  1. Check the support first. \(\log\) needs strictly positive values. Zeros or readings below a detection limit have to be handled before you log, not after.
  2. The mean does not transform; the median does. If you average logged concentrations and back-transform, \(e^{2.1993} = 9.02\), not the original mean of \(10\). That is the geometric mean, and it is always smaller. Report it as a geometric mean, or back-transform the median, which does carry over.
  3. Transformations change skew, not just scale. Logging this Gamma removed the right skew and left a slight left skew. Squaring made the tail much heavier. Check the shape after transforming instead of assuming it became normal.
  4. Squaring makes estimates noisier. The 85th percentile of \(X^2\) had a standard deviation of about 11 percent of its value across samples, compared with about 2 percent for \(\log(X)\) and 3 percent for \(\sqrt{X}\). Any statistic that leans on the tail gets less stable after a convex transformation.
  5. Sample percentiles run a little low at \(n = 100\). All three sampling distributions were centered slightly below the true 85th percentile. The gap is small, but it points the same way every time.

What I learned

Part 1: Base distribution

\(X \sim \text{Gamma}(\text{shape } k = 5, \text{scale } \theta = 2)\) has two parameters, and the second can be written as a scale \(\theta\) or a rate \(\beta = 1/\theta\). R defaults to rate, so I pass scale = theta by name every time.

The probability density function (PDF), \(f(x)\), is the curve. Area under it is probability, and the total area is 1. The height by itself is not a probability. The cumulative distribution function (CDF), \(F(x) = P(X \le x)\), is the area to the left of \(x\). It never decreases and stays between 0 and 1.

For this Gamma, the normalizing constant is \(\Gamma(5)\,2^5 = 768\), \(f(10) \approx 0.0877\), and \(F(10) \approx 0.5595\). The mean is \(10\) and the variance is \(20\). The median solves \(F(m) = 0.5\), so \(m = F^{-1}(0.5) = 9.3418\). \(F^{-1}\) is the CDF run backward: probability in, \(x\) out.

Part 2: Transformations

Log, square, and square root are all increasing on \(x > 0\), so the CDF shortcut \(F_Y(y) = F_X(g^{-1}(y))\) and the median rule both apply. The PDF of \(Y\) is \(f_X(g^{-1}(y))\) times the Jacobian \(\lvert dg^{-1}/dy \rvert\).

Support matters. \(\log(X)\) covers all real numbers, while \(X^2\) and \(\sqrt{X}\) stay positive.

The means come from the moment tools: \(E[\log X] = \psi(5) + \log 2 \approx 2.1993\), \(E[X^2] = 120\), and \(E[\sqrt{X}] \approx 3.0843\). Jensen’s inequality explains the direction. Log and square root are concave, so \(E[g(X)] < g(E[X])\); squaring is convex and goes the other way. The medians are just \(g(9.3418)\): \(2.2345\), \(87.2696\), and \(3.0564\).

The simulation (rgamma(), then density() and ecdf()) matches the analytic curves in every figure.

Part 3: Sampling distribution of the 85th percentile

The 85th percentile changes from sample to sample, and 1000 samples show how much. Increasing transformations carry percentiles the same way they carry medians: \(q_Y(0.85) = g(q_X(0.85))\).

The centers of the three sampling distributions sit slightly below the true values. That comes from sample-percentile bias at \(n = 100\) plus simulation noise. \(X^2\) has by far the widest spread, while log and square root come out close to symmetric.

The assumptions are independent and identically distributed (IID) draws, quantile() type 7, and the seed and draw order.

↑ Back to top