This chapter collects three tools that extend the univariate toolkit but sit outside the main path from data to inference. Kernel densities smooth the histogram into a continuous curve. Data transformations reshape a variable before summarizing or modeling it, and Jensen’s inequality says how the mean of transformed data relates to the transformed mean. Inverse sampling generates random draws from either an empirical or a theoretical distribution. Each section stands on its own, so you can read them in any order.
When we want a smoother picture than the jagged steps of a histogram, we replace fixed bins with a moving window centered on each value of \(x\) we want to evaluate.
A kernel density estimates the density at \(x\) by averaging a smooth weight function (the kernel) centered on each data point, with a bandwidth \(h\) that controls how wide the weight spreads. The bandwidth is a half-width, exactly as in the histogram: each kernel below weights observations within \(h\) of \(x\) and gives zero weight beyond. Unlike the histogram, the bins around each \(x\) can overlap, which is why the result is smooth instead of stepped.
Kernel densities are useful when we want a continuous picture of shape: easier to read than a histogram and better suited to overlaying two or three distributions on the same axes, where stacked or overlapping bars become visually muddled. Different kernels weight nearby points in different ways, producing curves with visually distinct (but qualitatively similar) shapes. The simplest kernel function is the “uniform” (or “rectangular”) kernel, which places a weight of \(1/2\) on every point inside the interval \(\left[ x-h, x + h\right]\) and zero outside: \[\begin{aligned} k_{U}\left( \hat{X}_{i}, x, h \right) &= \frac{\mathbf{1}\left(\frac{|\hat{X}_{i}-x|}{h} \leq 1\right)}{2} = \frac{\mathbf{1}\left( \hat{X}_{i} \in \left[ x-h, x + h\right] \right) }{2} \\ \hat{f}_{U}(x) &= \frac{1}{nh} \sum_{i}^{n} k_{U}(\hat{X}_{i}, x, h) = \frac{ \sum_{i}^{n} \mathbf{1}\left( \hat{X}_{i} \in \left[ x-h, x + h\right] \right) }{n 2 h} \end{aligned}\]
For example, take the dataset \(\{3, 3.1, 0.02\}\) with bandwidth \(h=1/2\) and evaluate the density at \(x=3\). The interval \([x-h, x+h]=[2.5, 3.5]\) contains the two points \(3\) and \(3.1\), so \(\sum_{i}\mathbf{1}\left(\hat{X}_{i} \in [2.5,3.5]\right)=2\) and \(\hat{f}_{U}(3)=\frac{2}{n 2h}=\frac{2}{3\times 1}=2/3\).
Notice that the uniform kernel is essentially the histogram but without the restriction that \(x\) must be a midpoint of exclusive bins. Typically, the points \(x\) are chosen to be either the unique observations or some equidistant set of “design points” (e.g., at \(512\) evenly spaced values of \(x\) the dataset, not just the midpoints of exclusive bins).1
# Practical Example
X <- USArrests[, 'Murder']
hist(X,
breaks=seq(0, 20, by=2), # bin width =2, so h=1
freq=FALSE,
border=NA,
main=NA,
xlab='Murder Arrests')
title('Murder Arrests', font.main=1)
# Density Estimate with half-width h=3
# 'density' reads bw as the kernel's standard deviation, so pass bw=h/sqrt(3)
h <- 3
lines( density(X, bw=h/sqrt(3), kernel='rectangular') )
# Raw Observations
rug(X, col=grey(0, .5))
We can also replace the uniform kernel with other kernel functions to create even smoother lines.
There are many kernels, but these are the most intuitive and commonly used.
# Kernel Density Functions
x <- 0 # Design point
X <- seq(-2, 2, length.out=1001) # Where to compute weight
h <- 1 # Half-width: all three kernels are zero beyond |X-x| > h
plot.new()
plot.window(xlim=c(-1.2, 1.2), ylim=c(0, 1))
## Uniform/rectangular
d_unif <- function(X, x, h){
u <- abs(X-x)/h
fu <- 1/2*(u <= 1)
return(fu)
}
kU <- d_unif(X, x, h)
kern_cols <- hcl.colors(3, alpha=.45)
lines( kU ~ X, col=kern_cols[1], lty=1)
## Epanechnikov
d_epan <- function(X, x, h){
u <- abs(X-x)/h
fu <- 3/4*(1-u^2)*(u <= 1)
return(fu)
}
kE <- d_epan(X, x, h)
lines( kE ~ X, col=kern_cols[2], lty=1)
## Try others using the 'density' function
# 'density' reads bw as the kernel's standard deviation, not the half-width h
#kE <- density(x=0, bw=h/sqrt(5), kernel='epanechnikov')
## Tricubic
d_tricub <- function(X, x, h){
u <- abs(X-x)/h
fu <- 70/81*(1-u^3)^3*(u <= 1)
return(fu)
}
kT <- d_tricub(X, x, h)
lines( kT ~ X, col=kern_cols[3], lty=1)
rug(0, lwd=2)
axis(1)
axis(2)
legend('topright', lty=1, col=kern_cols,
legend=c('Uniform (h=1)', 'Epanechnikov (h=1)', 'Tricubic (h=1)'))
Each of the three functions integrates to one over \(u=(\hat{X}_{i}-x)/h\), so the same \(\frac{1}{nh}\) rescaling above turns any of them into a density.
Once we have picked a kernel (which particular one is not particularly important) we can use it to compute density estimates at each design point.
# Practical Example
X <- USArrests[, 'Murder']
hist(X,
breaks=seq(0, 20, by=2), # bin width 2h=2, so h=1
freq=FALSE,
border=NA,
main=NA,
xlab='Murder Arrests')
title('Murder Arrests', font.main=1)
# Same half-width h=3 as the rectangular kernel above, but a smoother weight
# For the Epanechnikov kernel, bw=h/sqrt(5)
h <- 3
lines( density(X, bw=h/sqrt(5), kernel='epanechnikov') )
# Raw Observations
rug(X, col=grey(0, .5))
Transformations can stabilize variance, reduce skewness, and make model errors closer to Gaussian.
Perhaps the most common examples are power transformations: \(y= x^\lambda\), which includes \(\sqrt{x}\) and \(x^2\).
Other examples include the exponential transformation: \(y=\exp(x)\) for any \(x\in (-\infty, \infty)\) and logarithmic transformation: \(y=\log(x)\) for any \(x>0\).
The Box–Cox Transform nests many cases used by statisticians. For \(x>0\) and parameter \(\lambda\), \[ y=\begin{cases} \dfrac{x^\lambda-1}{\lambda}, & \lambda\neq 0,\\ \log\left(x\right) & \lambda=0. \end{cases} \] This function is continuous over \(\lambda\).
# Box–Cox transform and inverse
bc_transform <- function(x, lambda) {
if (any(x <= 0)) stop('Box-Cox requires x > 0')
if (abs(lambda) < 1e-8) log(x) else (x^lambda - 1)/lambda
}
bc_inverse <- function(t, lambda) {
if (abs(lambda) < 1e-8) exp(t) else (lambda*t + 1)^(1/lambda)
}
X <- USArrests[, 'Murder']
hist(X, main=NA, border=NA, freq=FALSE)
par(mfrow=c(1, 3))
for(lambda in c(-1, 0, 1)){
Y <- bc_transform(X, lambda)
hist(Y,
main=NA,
border=NA, freq=FALSE)
title(bquote(paste(lambda, '=', .(lambda))), font.main=1)
}
Be careful about transforming your data, as the interpretation can be harder.
The mean of transformed data is not equivalent to transforming the mean, for example, which is one reason why you want to first summarize your data before transforming it.
We can actually say more about how the mean of transformed data is not equivalent to transforming the mean. To do that, we define two types of functions:
Let \(\hat{Y}_{i}=g( \hat{X}_{i})\), and denote the mean as \(\hat{M}_{Y}\).
If \(g\) is a concave function, then \(g( \hat{M}_{X} ) \geq \hat{M}_{Y}\).
If \(g\) is a convex function, then the inequality reverses: \(g( \hat{M}_{X}) \leq \hat{M}_{Y}\).
mean( exp(x) )
## [1] 10.06429
exp( mean(x) )
## [1] 7.389056To generate a random variable from known distributions, you can use some type of physical machine. E.g., you can roll a fair die to generate Discrete Uniform data or you can roll weighted die to generate Categorical data.
There are also several ways to computationally generate random variables from a probability distribution. Perhaps the most common one is inverse sampling. To generate a random variable using inverse sampling, first sample \(p\) from a uniform distribution and then find the associated quantile via the quantile function \(\hat{F}^{-1}(p)\).2
You can generate a random variable from a known empirical distribution. Inverse sampling randomly selects observations from the dataset with equal probabilities. To implement this, we
Here is an example of generating random murder rates for US states.
# Empirical Distribution
X <- USArrests[, 'Murder']
FX_hat <- ecdf(X)
plot(FX_hat, lwd=2, xlim=c(0, 20),
pch=16, col=grey(0, .5), main=NA)
# Generating random variables via inverse ECDF
p <- runif(3000) ## Multiple Draws
QX_hat <- quantile(FX_hat, p, type=1)
QX_hat[c(1, 2, 3)]
## 52.86985035% 39.33680111% 6.45579686%
## 7.4 6.0 2.2
## Can also do directly from the data
QX_hat <- quantile(X, p, type=1)
QX_hat[c(1, 2, 3)]
## 52.86985035% 39.33680111% 6.45579686%
## 7.4 6.0 2.2If you know the distribution function that generates the data, then you can derive the quantile function and do inverse sampling. That is how computers generate random data from a distribution.
# 4 random data points from 3 different distributions
runif(4)
## [1] 0.8314547 0.2261331 0.7374920 0.8912296
rexp(4)
## [1] 0.2945136 0.2759661 0.4427445 0.3554522
rnorm(4)
## [1] 0.3569255 -0.2070293 -1.1101673 0.2843943Here is an in-depth example of the Dagum distribution. The distribution function is \(F(x)=(1+(x/b)^{-a})^{-c}\). For a given probability \(p\), we can then solve for the quantile as \[F^{-1}(p)=\frac{ b p^{\frac{1}{ac}} }{(1-p^{1/c})^{1/a}}\]. After which, we sample \(p\) from a uniform distribution and then find the associated quantile.
# Theoretical Quantile Function (from VGAM::qdagum)
qdagum <- function(p, scale.b=1, shape1.a, shape2.c) {
# Quantile function (theoretically derived from the CDF)
ans <- scale.b * (expm1(-log(p) / shape2.c))^(-1 / shape1.a)
# Special known cases
ans[p == 0] <- 0
ans[p == 1] <- Inf
# Safety Checks
ans[p < 0] <- NaN
ans[p > 1] <- NaN
if(scale.b <= 0 | shape1.a <= 0 | shape2.c <= 0){ ans <- ans*NaN }
# Return
return(ans)
}
# Generate Random Variables (VGAM::rdagum)
rdagum <-function(n, scale.b=1, shape1.a, shape2.c){
p <- runif(n) # generate random probabilities
x <- qdagum(p, scale.b=scale.b, shape1.a=shape1.a, shape2.c=shape2.c) #find the inverses
return(x)
}
# Example
set.seed(123)
X <- rdagum(3000, 1, 3, 1)
X[c(1, 2, 3)]
## [1] 0.7390476 1.5499868 0.8845006Comment the script you wrote for this chapter, then restart R and check that it runs from a clean session, then check the script with AI as explained in Working with AI. Write three sentences from memory on the main statistical idea of this chapter, and ask the assistant what is wrong, vague, or missing. Finish with your own questions about whatever you found hardest.
Using the dataset \(\{3, 3.1, 0.02\}\) and half-width \(h=1/2\), compute the uniform kernel density \(\hat{f}_{U}(x)\) by hand at \(x=0\), \(x=3\), and \(x=5\). That is, for each \(x\) count the observations within \(h\) of \(x\) and divide by \(n \cdot 2h\). Then write a short R function that does the same count and verify your three answers.
Load USArrests in R and extract the Assault variable. Compute mean(log(X)) and log(mean(X)). Which is larger, and which case of Jensen’s inequality does this illustrate? Then apply the Box–Cox transform with \(\lambda \in \{-1, 0, 1/2, 1\}\) using bc_transform from this chapter, and plot the four histograms side by side. Which \(\lambda\) gives the most symmetric histogram?
The exponential distribution with rate \(\lambda\) has distribution function \(F(x)=1-e^{-\lambda x}\) for \(x \geq 0\). Solve for the quantile function \(F^{-1}(p)\) by hand. Write an R function that draws \(n\) values of \(p\) from a uniform distribution and returns \(F^{-1}(p)\), then generate \(5000\) draws with \(\lambda=2\). Plot a histogram of your draws with freq=FALSE and overlay dexp(x, rate=2) to check that they match.
This chapter collected three tools that extend the univariate toolkit: the kernel density, which smooths a histogram by averaging a weight function centered on each observation with half-width \(h\); data transformations such as the log and the Box–Cox family, together with Jensen’s inequality on how the mean of transformed data compares to the transformed mean; and inverse sampling, which generates random draws by feeding uniform probabilities through a quantile function. The running three-point dataset \(\{3, 3.1, 0.02\}\) pins the first of these down: the uniform kernel density at \(x=3\) with half-width \(h=1/2\) equals \(2/3\), because two points fall inside the interval \([2.5, 3.5]\) and \(2/(3 \times 1)=2/3\). For the second, mean(sqrt(x)) is never larger than sqrt(mean(x)) because the square root is concave, and the inequality flips for a convex function like the exponential. In the next part we extend distributions, statistics, and probability rules to pairs of variables.
One detail matters when moving from the formula to R. The density function scales its bw argument to be the standard deviation of the kernel rather than the half-width, so bw and \(h\) are different numbers. Setting bw from \(h\), rather than reading \(h\) off bw, keeps the code matching the math. The conversion is a constant that depends on the kernel: the half-width is \(\sqrt{3}\) times bw for the rectangular kernel and \(\sqrt{5}\) times bw for the Epanechnikov kernel.↩︎
Drawing random uniform samples with computers is actually quite complex and beyond the scope of this book.↩︎