分布的机理：四种运算生成一切
2026-07-04
姊妹文档：
exp_normal_boltzmann.Rmd从最大熵原理出发，推导指数、正态、玻尔兹曼三个分布。本文关注生成机制——每个分布是通过什么运算从更简单的分布产生的。
本文的核心论点：统计学中的几十种分布不是彼此独立的发明，而是从一个起点（Uniform）通过四种运算反复生成的。
\[\text{Uniform}(0,1) \xrightarrow{\text{变换 / 求和 / 极值 / 比值}} \text{所有常见分布}\]
par(mar = c(0.5, 0.5, 2, 0.5))
plot(NA, xlim = c(0, 14), ylim = c(0, 10), axes = FALSE, bty = "n",
xlab = "", ylab = "",
main = "Road map: four operations generate all distributions")
# ── colors by operation ──
col_src <- "gray30" # source
col_tf <- "#2E86AB" # transform
col_sum <- "#A23B72" # sum
col_rat <- "#F18F01" # ratio
col_ext <- "#C73E1D" # extreme
draw_box <- function(x, y, label, col, w = 2.2, h = 0.55) {
rect(x - w/2, y - h/2, x + w/2, y + h/2,
col = col, border = "white", lwd = 2)
text(x, y, label, col = "white", font = 2, cex = 0.72)
}
draw_arr <- function(x0, y0, x1, y1, label = "", adj_x = 0, adj_y = 0.15) {
arrows(x0, y0, x1, y1, length = 0.10, lwd = 1.4, col = "gray50")
if (nchar(label) > 0)
text((x0+x1)/2 + adj_x, (y0+y1)/2 + adj_y,
label, cex = 0.58, col = "gray30", font = 3)
}
# ── Row 1: source ──
draw_box(7, 9.3, "Uniform(0,1)", col_src, w = 2.5)
# ── Row 2: transform ──
draw_box(2, 7.5, "Exp(lambda)", col_tf)
draw_box(7, 7.5, "Normal(mu, sigma)", col_tf)
draw_box(12, 7.5, "Weibull(k, lambda)",col_tf)
draw_arr(5.8, 9.05, 2, 7.8, "F^-1(U) = -ln(1-U)/lam")
draw_arr(7, 9.05, 7, 7.8, "F^-1(U) = qnorm(U)")
draw_arr(8.2, 9.05, 12, 7.8, "(-ln(1-U))^(1/k)")
# ── Row 3: sum / transform derivatives ──
draw_box(2, 5.5, "Gamma(k, lambda)", col_sum)
draw_box(5.5, 5.5, "Poisson(lambda*t)", col_sum)
draw_box(9, 5.5, "Chi-squared(k)", col_sum)
draw_box(12, 5.5, "LogNormal(mu, sigma)", col_tf)
draw_arr(2, 7.2, 2, 5.8, "k copies summed")
draw_arr(2, 7.2, 5.5, 5.8, "count in window", adj_x = 0.3)
draw_arr(7, 7.2, 9, 5.8, "k squares summed")
draw_arr(7, 7.2, 12, 5.8, "exp(X)", adj_x = 0.5)
# ── Row 4: ratio ──
draw_box(7, 3.5, "t(k)", col_rat)
draw_box(11, 3.5, "F(m, n)",col_rat)
draw_arr(7, 5.2, 7, 3.8, "Z / sqrt(chi2/k)", adj_x = -1.5)
draw_arr(9, 5.2, 7.8, 3.8, "")
draw_arr(9, 5.2, 11, 3.8, "chi2_m/m / chi2_n/n")
# ── Row 5: extreme ──
draw_box(2, 1.5, "Beta(a, b)", col_ext)
draw_box(5.5, 1.5, "Gumbel", col_ext)
draw_box(9, 1.5, "Frechet", col_ext)
draw_box(12, 1.5, "Pareto(alpha)", col_ext)
draw_arr(7, 9.05, 2, 1.8, "order stats", adj_x = -1.3, adj_y = 0.2)
draw_arr(7, 7.2, 5.5,1.8, "max of n (light tail)")
draw_arr(12, 1.8, 9, 1.8, "max of n", adj_y = 0.2)
# ── legend ──
legend(0, 3.8,
legend = c("Source", "Transform (Sec 2)",
"Sum (Sec 3)", "Ratio (Sec 4)", "Extreme (Sec 5)"),
fill = c(col_src, col_tf, col_sum, col_rat, col_ext),
bty = "n", cex = 0.75, border = "white")
1 起点：Uniform(0,1)
计算机生成的一切随机数都来自 Uniform(0,1)。它的 CDF 是恒等函数 \(F(u) = u\)——这个性质看似平淡，却是逆变换法的根基。
U <- runif(N)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(U, breaks = 60, col = "steelblue", border = "white", freq = FALSE,
main = "Uniform(0,1): perfectly flat", xlab = "u", ylab = "Density")
abline(h = 1, col = "tomato", lwd = 2, lty = 2)
legend("topright", "Theoretical density = 1",
col = "tomato", lwd = 2, lty = 2, bty = "n")
plot(ecdf(U), main = "CDF of Uniform(0,1): F(u) = u",
xlab = "u", ylab = "F(u)", col = "steelblue", lwd = 2)
abline(0, 1, col = "tomato", lwd = 2, lty = 2)
legend("topleft", c("Empirical CDF", "Theoretical F(u) = u"),
col = c("steelblue", "tomato"), lwd = 2, bty = "n")
2 运算一：变换——一个随机变量进，一个随机变量出
2.1 逆 CDF 变换（万能方法）
2.1.1 先说为什么要有这个变换：计算机只会掷”均匀骰子”
这个变换不是凭空的数学游戏，它来自计算机模拟的一个硬需求。
计算机的随机数发生器（RNG）只会产生一种随机：\(U\sim\text{Uniform}(0,1)\)——\([0,1]\) 上的均匀随机数。它天生不会直接吐指数分布、正态分布、泊松分布……那怎么办？
答案：拿唯一能白拿的均匀随机，喂进逆 CDF，把它”掰”成你要的分布。
\[U\sim\text{Uniform}(0,1)\ \xrightarrow{\ X=F^{-1}(U)\ }\ X\sim F.\]
所以这一节的主角其实是”怎么用一个均匀骰子造出任意分布的随机数”。\(U\) 是原料，\(F^{-1}\)（逆 CDF / 分位函数）是成型机。这也是整份文档”四种运算生成一切分布”的第一种运算——一切都从 uniform（和 normal）出发。
定理：若 \(U \sim \text{Uniform}(0,1)\)，\(F\) 是任意分布的 CDF，则 \(X = F^{-1}(U) \sim F\)。
说人话：拿一个均匀随机数 \(U\)，喂进逆 CDF，得到的新随机数 \(X=F^{-1}(U)\) 正好服从你想要的那个分布。“\(X \sim F\)” 读作”\(X\) 服从 \(F\)“，是一句话的缩写：\(P(X\leq x)=F(x)\) 对所有 \(x\) 成立——即 \(X\) 出数的累积规律恰好由 \(F\) 描述。
记号（三个东西类型完全不同，别混）：
- 大写 \(X = F^{-1}(U)\) 是随机变量——会随机出数的”骰子”，每抽一次吐一个具体数（2.46, 0.34, …）。由随机的 \(U\) 生成，本身也随机。目标是证明它服从分布 \(F\)。
- 小写 \(x\) 是一个任意但固定的实数——横轴上一个具体的坐标点。证明中对它取遍所有值，从而通过 \(P(X\leq x)=F(x)\)（对每个 \(x\) 都成立）确定 \(X\) 的分布。
先看直觉版（假设 \(F\) 连续、严格递增，绝大多数常见分布都满足）。核心那一步 \(F^{-1}(U)\leq x \iff U\leq F(x)\) 就是两边同时套 \(F\)：
\[ \begin{aligned} F^{-1}(U)\leq x &\iff F\!\left(F^{-1}(U)\right)\leq F(x) &&\text{（}F\text{ 单调递增，两边套 }F\text{ 保持}\leq\text{）}\\ &\iff U\leq F(x) &&\text{（}F\text{ 与 }F^{-1}\text{ 抵消：}F(F^{-1}(U))=U\text{）}. \end{aligned} \]
第一步：\(F\) 单调递增，所以 \(a\leq b\) 与 \(F(a)\leq F(b)\) 是一回事，两边套 \(F\) 不改变不等号方向。 第二步：\(F\) 和它的逆 \(F^{-1}\) 互相抵消，\(F(F^{-1}(U))\) 就还原成 \(U\)。
再看严格版（把连续、严格递增这两个假设也去掉）。用广义逆
\[F^{-1}(u)=\inf\{t:F(t)\geq u\}\]
表示“CDF 第一次达到高度 \(u\) 时的横坐标”。因此，固定任意 \(x\)，
\[ \begin{aligned} F^{-1}(U)\leq x &\iff \text{CDF 在横坐标 }x\text{ 之前已经达到高度 }U\\ &\iff F(x)\geq U\\ &\iff U\leq F(x). \end{aligned} \]
这两个不等式描述的是同一个随机事件，只是前者从横轴读，后者从纵轴读。于是
\[P(X\leq x)=P\!\left(F^{-1}(U)\leq x\right)=P\!\left(U\leq F(x)\right)=F(x).\]
最后一步使用了 \(U\sim\text{Uniform}(0,1)\)：对任意 \(a\in[0,1]\)，\(P(U\leq a)=a\)。
直觉：\(U\) 在概率轴 \([0,1]\) 上均匀采样，\(F^{-1}\) 把概率值”翻译”回 \(x\) 值。\(F\) 越陡处密度越高，\(F\) 越平处密度越低。
下图固定红色的 \(x\)。此时 CDF 已上升到红色高度 \(F(x)\)：
- 若抽到的 \(U\leq F(x)\)（紫色区域），从该高度向右碰到 CDF，再向下投影，落点一定在 \(x\) 左侧，即 \(F^{-1}(U)\leq x\)。
- 若 \(U>F(x)\)（灰色示例），交点和落点都会在 \(x\) 右侧。
t_grid <- seq(-4, 4, length.out = 1000)
F_grid <- plogis(t_grid)
x_cut <- 0.8
Fx_cut <- plogis(x_cut)
u_demo <- 0.52
x_demo <- qlogis(u_demo)
u_above <- 0.84
x_above <- qlogis(u_above)
op <- par(mar = c(4.8, 5, 3.8, 1))
plot(t_grid, F_grid, type = "n", xlim = c(-4, 4), ylim = c(0, 1),
xlab = "horizontal position t", ylab = "cumulative probability",
main = expression("The same cutoff, read on two axes"),
axes = FALSE, xaxs = "i", yaxs = "i")
# All possible U below F(x): the event U <= F(x).
rect(-4, 0, 4, Fx_cut,
col = adjustcolor("#6F42C1", alpha.f = 0.09), border = NA)
axis(1, at = c(x_demo, x_cut), labels = FALSE, col = "gray35")
axis(2, at = c(0, 0.25, 0.5, 0.75, 1), las = 1,
col.axis = "gray35")
box(col = "gray55")
lines(t_grid, F_grid, col = "#2E5AAC", lwd = 3)
# The fixed cutoff x and its CDF height F(x).
segments(x_cut, 0, x_cut, Fx_cut, col = "#D62728", lwd = 2, lty = 2)
segments(-4, Fx_cut, x_cut, Fx_cut, col = "#D62728", lwd = 2, lty = 2)
points(x_cut, Fx_cut, pch = 21, bg = "#D62728", col = "white", cex = 1.25)
text(x_cut + 0.13, -0.075, expression(x), col = "#D62728",
font = 2, xpd = NA)
text(-3.92, Fx_cut + 0.045, expression(F(x)), col = "#D62728",
font = 2, adj = 0)
# One sampled U <= F(x): horizontal to the CDF, then down to F^-1(U).
segments(-4, u_demo, x_demo, u_demo,
col = "#6F42C1", lwd = 2.5, lty = 3)
segments(x_demo, 0, x_demo, u_demo,
col = "#6F42C1", lwd = 2.5, lty = 3)
points(x_demo, u_demo, pch = 21, bg = "#6F42C1",
col = "white", cex = 1.25)
text(-3.92, u_demo + 0.035, expression(U), col = "#6F42C1",
font = 2, adj = 0)
text(x_demo, -0.075, expression(F^{-1}(U)), col = "#6F42C1",
font = 2, xpd = NA)
# A counterexample to the event: U > F(x) maps to the right of x.
segments(-4, u_above, x_above, u_above,
col = "gray55", lwd = 1.6, lty = 3)
segments(x_above, 0, x_above, u_above,
col = "gray55", lwd = 1.6, lty = 3)
points(x_above, u_above, pch = 21, bg = "gray55",
col = "white", cex = 1.05)
text(x_above + 0.12, u_above + 0.035,
"U > F(x) lands right of x", col = "gray35", adj = 0, cex = 0.78)
# Labels that explicitly connect the two descriptions of the event.
arrows(-3.65, 0.04, -3.65, Fx_cut - 0.04, code = 3,
length = 0.07, col = "#6F42C1", lwd = 1.8)
text(-3.42, Fx_cut / 2, expression(U <= F(x)),
col = "#6F42C1", srt = 90, font = 2)
arrows(-3.8, 0.09, x_cut - 0.08, 0.09, code = 3,
length = 0.07, col = "#6F42C1", lwd = 1.8)
text(-1.45, 0.135, expression(F^{-1}(U) <= x),
col = "#6F42C1", font = 2)
legend("bottomright",
legend = c(expression(F(t)), expression("fixed cutoff " * (x * "," * F(x))),
expression("sampled " * U <= F(x))),
col = c("#2E5AAC", "#D62728", "#6F42C1"),
lwd = c(3, 2, 2.5), lty = c(1, 2, 3),
bty = "n", cex = 0.82)
图中为了直观使用连续 CDF；若 CDF 是阶梯状的，则把 \(F^{-1}(u)\) 理解为“第一次达到或超过 \(u\) 的位置”，上述事件等价仍然成立。
2.1.2 指数分布：手推解析逆 CDF
\(F(x) = 1 - e^{-x}\)。逆 CDF 做的事情是：先给定纵轴上的累计概率 \(u\)，再反过来寻找横轴上的 \(x\)。
逐步反解：
\[ \begin{aligned} u &= 1-e^{-x}\\ 1-u &= e^{-x}\\ \log(1-u) &= -x\\ x &= -\log(1-u) \end{aligned} \]
例如，\(u=0.8\) 时，\(x=-\log(0.2)\approx1.61\)；也就是说，指数分布约有 80% 的概率落在 \([0,1.61]\) 内。
一张图讲完整条链，两列六格，所有面板同尺寸、坐标轴严格对齐。真实数据沿红色编号箭头走一个「П」形：右下（输入）→ 右中（机器）→ 左中（竖着的输出）→ 左下（放倒的输出）。把 Uniform 的 5 段颜色（各 20% 概率）从头贴上去：
- 右上 — CDF \(F(x)=1-e^{-x}\)：横轴 \(x\)、纵轴 \(u\)。它和正下方机器图之间是竖直红色双箭头
x <-> u axes flip——机器图本身就是 inverse CDF，两轴对调直接上下对照。 - 右中 — inverse CDF = 机器（真实数据穿过）：12 个真实 \(U\) 各画一条轨迹——从 \(u\) 轴竖直上升碰到曲线，再水平射向左边。底边小刻度是 400 个真实 \(U\)（均匀散开）。
- 右下 — 输入 \(U\sim\text{Uniform}(0,1)\)：5 个等宽色块 = 等概率。灰色细框是 10 万个真实 \(U\) 的直方图——真数据不完美水平。
- 左中 — 输出（\(x\) 轴竖着）：和左下同一张图、同一批数据、同一套颜色，只是 \(x\) 轴还竖着——竖轴与右边机器图的竖轴完全相同，机器射出的轨迹箭头直接落进来。灰色横条 = 10 万个真实 \(-\log(1-U)\) 的直方图，贴着理论曲线 \(e^{-x}\)。
- 左下 — 输出（放倒）：把左中整张图顺时针转
90° 放平，就是教科书里常见的指数密度图。红色箭头
3. rotate 90 -> lay it flat标出这个动作。
编号红箭头 = 真实数据流水线：
1. feed REAL uniform draws（右下 ↑ 右中）：真实均匀数一个个喂进 \(x=-\log(1-u)\)——一次一个 log，不求导。- 右中 → 左中的无字红箭头：变换后的 \(X\) 落在同一根竖 \(x\) 轴上，直接数直方图——密度是数出来的。
3. rotate 90 -> lay it flat（左中 ↓ 左下）：把竖着的图放倒,得到常规密度图。旋转只是排版动作，数学内容一丝不变。
红点追踪：\(u=0.8\) → \(x=-\log(0.2)=1.61\) → 密度 \(f(1.61)=0.2\)，四张图里同一个红点。右上、左中右侧、右中右侧、左下 \(x\) 轴上的红色刻度 \(0,\,0.22,\,0.51,\,0.92,\,1.61\)（band edges）是同一批 \(x\) 分界。
prob_breaks<-seq(0,1,by=0.2);quantile_breaks<-qexp(prob_breaks,rate=1)
band_cols<-c("#4E79A7","#59A14F","#EDC948","#F28E2B","#E15759")
u_star<-0.8;x_star<- -log(1-u_star);tcol<-function(c,a)adjustcolor(c,alpha.f=a)
RED<-"#D62728";DRED<-"#7a1f1f"
# 3x2 grid, two equal columns. Every plot panel shares MAR -> columns align exactly.
layout(matrix(c(1,2, 3,4, 5,6),3,2,byrow=TRUE))
par(cex.lab=1.35,cex.main=1.15,cex.axis=1.0,mgp=c(3,0.8,0))
MAR <- c(6.2,5.2,3.4,3.6)
X_real <- -log(1-U)
h <- hist(X_real[X_real<=5],breaks=seq(0,5,by=0.1),plot=FALSE)
hs <- h$density*mean(X_real<=5) # histogram heights, scaled to full density
u12 <- sort(runif(12)); x12 <- -log(1-u12)
u400 <- runif(400)
# (1) top-left: pipeline text
par(mar=c(2,2,2,2))
plot(NA,xlim=c(0,1),ylim=c(0,1),xlab="",ylab="",main="",axes=FALSE)
text(0.5,0.5,paste0(
"the real-data pipeline:\n\n",
"input (bottom right): computer draws U\n",
" 0.915, 0.937, 0.286, ...\n",
"machine (middle right): x = -log(1-u)\n",
" 0.915 -> 2.46 0.286 -> 0.34\n",
"output (middle left): histogram of 100,000\n",
" such x's, on the SAME vertical x-axis\n",
"lay it flat (bottom left): rotate 90 degrees\n",
" -> the usual density plot e^-x"),cex=1.25,col="gray20")
# (2) CDF, top right
par(mar=MAR)
plot(NA,xlim=c(0,5),ylim=c(0,1),xlab="",ylab=expression(u==F(x)),main=expression("CDF: "*F(x)==1-e^{-x}),yaxs="i")
for(i in 1:5)rect(0,prob_breaks[i],5,prob_breaks[i+1],col=tcol(band_cols[i],0.30),border=NA)
curve(1-exp(-x),0,5,add=TRUE,lwd=2.6,col="#2E5AAC")
segments(0,u_star,x_star,u_star,lty=2,col=RED);segments(x_star,0,x_star,u_star,lty=2,col=RED);points(x_star,u_star,pch=21,bg=RED,col="white",cex=1.3)
mtext(expression(x),side=1,line=2.4,cex=0.95)
axis(1,at=round(quantile_breaks[1:5],2),labels=round(quantile_breaks[1:5],2),las=2,line=3.0,cex.axis=0.85,col.axis=DRED,col=DRED)
# (3) OUTPUT STOOD UPRIGHT, middle left: same plot as bottom, x-axis VERTICAL
par(mar=MAR)
plot(NA,xlim=c(1.08,0),ylim=c(0,5),xlab="",ylab=expression(x==-log(1-u)),
main=expression("output, x-axis vertical: "*f(x)==e^{-x}),yaxs="i",xaxs="i")
for(i in 1:5){
l<-quantile_breaks[i];r<-min(quantile_breaks[i+1],5)
xx<-seq(l,r,length.out=150)
polygon(c(0,dexp(xx),0),c(l,xx,r),col=tcol(band_cols[i],0.72),border=NA)
}
rect(hs,h$breaks[-length(h$breaks)],0,h$breaks[-1],border="gray30",col=NA,lwd=0.7)
xg<-seq(0,5,length.out=400); lines(dexp(xg),xg,lwd=2,col="gray20")
abline(h=quantile_breaks[2:5],col="white",lwd=1.2)
f_star<-dexp(x_star)
segments(f_star,0,f_star,x_star,lty=2,col=RED);segments(1.08,x_star,f_star,x_star,lty=2,col=RED)
points(f_star,x_star,pch=21,bg=RED,col="white",cex=1.3)
mtext(expression(f(x)),side=1,line=2.4,cex=0.95)
axis(4,at=round(quantile_breaks[1:5],2),las=1,cex.axis=0.8,col.axis=DRED,col=DRED)
mtext("band edges x",side=4,line=2.3,cex=0.80,col=DRED)
text(0.75,4.35,"same colors, same data as below --\njust not laid flat yet.\ngray bars: real histogram",cex=1.05,col="gray25",font=2)
# (4) inverse CDF = MACHINE with real data, middle right
par(mar=MAR)
plot(NA,xlim=c(0,1),ylim=c(0,5),xlab=expression(u),ylab=expression(x==-log(1-u)),
main=expression("inverse CDF: "*x==-log(1-u)*" -- REAL DATA"),yaxs="i",xaxs="i")
for(i in 1:5)rect(prob_breaks[i],0,prob_breaks[i+1],5,col=tcol(band_cols[i],0.30),border=NA)
curve(-log(1-x),0,0.985,add=TRUE,lwd=2.6,col="#6F42C1")
for(i in 2:5)segments(0,quantile_breaks[i],prob_breaks[i],quantile_breaks[i],lty=3,col="gray40")
for(j in 1:12){
segments(u12[j],0,u12[j],x12[j],col="gray35",lwd=0.9)
arrows(u12[j],x12[j],0.006,x12[j],length=0.05,col="gray35",lwd=0.9)
points(u12[j],x12[j],pch=21,bg=band_cols[findInterval(u12[j],prob_breaks,rightmost.closed=TRUE)],col="white",cex=1.0)
}
rug(u400,side=1,ticksize=0.025,col=tcol("gray20",0.35),lwd=0.6)
axis(4,at=round(quantile_breaks[1:5],2),las=1,cex.axis=0.8,col.axis=DRED,col=DRED)
mtext("band edges x",side=4,line=2.3,cex=0.80,col=DRED)
segments(u_star,0,u_star,x_star,lty=2,col=RED);segments(0,x_star,u_star,x_star,lty=2,col=RED);points(u_star,x_star,pch=21,bg=RED,col="white",cex=1.3)
text(0.42,4.55,"bottom ticks: 400 real U's (even)\narrows exit left: same vertical x-axis\nas the plot on the left",cex=1.05,col="gray25",font=2)
# (5) OUTPUT LAID FLAT, bottom left: the usual density plot
par(mar=MAR)
plot(NA,xlim=c(0,5),ylim=c(0,1.08),xlab="",ylab=expression(f(x)),
main=expression("output laid flat: X ~ Exp(1), "*f(x)==e^{-x}),yaxs="i")
for(i in 1:5){l<-quantile_breaks[i];r<-min(quantile_breaks[i+1],5);xx<-seq(l,r,length.out=150);polygon(c(l,xx,r),c(0,dexp(xx),0),col=tcol(band_cols[i],0.72),border=NA)}
rect(h$breaks[-length(h$breaks)],0,h$breaks[-1],hs,border="gray30",col=NA,lwd=0.7)
curve(dexp(x),0,5,add=TRUE,lwd=2,col="gray20")
abline(v=quantile_breaks[2:5],col="white",lwd=1.2)
segments(x_star,0,x_star,f_star,lty=2,col=RED);segments(0,f_star,x_star,f_star,lty=2,col=RED);points(x_star,f_star,pch=21,bg=RED,col="white",cex=1.3)
text(2.6,0.62,"gray bars: histogram of\n100,000 real -log(1-U)\ncurve: theory e^-x",cex=1.1,col="gray25",font=2)
mtext(expression(x),side=1,line=2.4,cex=0.95)
axis(1,at=round(quantile_breaks[1:5],2),labels=round(quantile_breaks[1:5],2),las=2,line=3.0,cex.axis=0.85,col.axis=DRED,col=DRED)
mtext("band edges = inverse-CDF values -log(1-u)",side=1,line=5.4,cex=0.80,col=DRED)
# (6) INPUT uniform, bottom right
par(mar=MAR)
plot(NA,xlim=c(0,1),ylim=c(0,1.15),xlab=expression(u),ylab="density",main="input: U ~ Uniform(0,1)",yaxs="i")
for(i in 1:5){rect(prob_breaks[i],0,prob_breaks[i+1],1,col=tcol(band_cols[i],0.72),border="white",lwd=2);text(mean(prob_breaks[i:(i+1)]),0.5,"20%",col="white",font=2,cex=0.9)}
h_u <- hist(U,breaks=seq(0,1,by=0.04),plot=FALSE)
rect(h_u$breaks[-length(h_u$breaks)],0,h_u$breaks[-1],h_u$density,border="gray30",col=NA,lwd=0.7)
segments(0,1,1,1,lwd=2,col="gray25")
text(0.5,1.09,"gray outline: histogram of 100,000 real U's (not perfectly flat)",cex=1.0,col="gray25",font=2)
# overlays
par(fig=c(0,1,0,1),new=TRUE,mar=c(0,0,0,0));plot(0:1,0:1,type="n",axes=FALSE,xlim=c(0,1),ylim=c(0,1),xaxs="i",yaxs="i")
# flip: CDF <-> machine (right column, rows 1-2)
arrows(0.82,0.702,0.82,0.661,code=3,length=0.10,lwd=2.8,col=RED)
text(0.828,0.684,"x <-> u axes flip",cex=1.3,col=RED,font=2,pos=4)
# 1. uniform -> machine (right column, rows 3->2)
arrows(0.90,0.318,0.90,0.358,length=0.13,lwd=3.4,col=RED)
text(0.895,0.338,"1. feed REAL uniform draws",cex=1.25,col=RED,font=2,pos=2)
# 2. machine -> upright output (middle row, right -> left)
arrows(0.550,0.60,0.490,0.60,length=0.13,lwd=3.4,col=RED)
# 3. upright -> laid flat (left column, rows 2->3)
arrows(0.10,0.358,0.10,0.318,length=0.13,lwd=3.4,col=RED)
text(0.105,0.338,"3. rotate 90 -> lay it flat",cex=1.25,col=RED,font=2,pos=4)
# The hand-rolled inverse CDF and R's qexp are the same formula.
X_exp_hand <- -log(1 - U)
X_exp_qexp <- qexp(U, rate = 1)
cat("Max | -log(1-U) - qexp(U) | =",
max(abs(X_exp_hand - X_exp_qexp)), "\n")## Max | -log(1-U) - qexp(U) | = 0注意两个端点：\(u=0\) 映射到 \(x=0\)；当 \(u\to1\) 时，\(1-u\to0\)，所以 \(-\log(1-u)\to\infty\)。这正是指数分布右侧长尾的来源。qexp(u, rate = 1)
只是同一个公式的 R 实现，不是另一种生成机制。
2.1.3 Weibull：在指数逆 CDF 上再做幂变换
Weibull\((k, \lambda)\) 的 CDF：\(F(x) = 1 - e^{-(x/\lambda)^k}\)，反解得：
\[F^{-1}(u) = \lambda \cdot \bigl(-\log(1-u)\bigr)^{1/k}\]
也就是说：Weibull 的逆 CDF = 指数逆 CDF 的 \(1/k\) 次方。
下面把 exp 那张流水线图原样复刻一遍，只把机器换成 Weibull 的逆 CDF（\(k=2,\ \lambda=1\)）。看点是同一条流水线，只多了一个开方：
- 机器变成 \(x=\sqrt{-\log(1-u)}\) ——就是 exp 机器的输出再开平方。开方把 exp 的大值压扁（\(2.46\to1.57\)），长尾被收紧。
- CDF 变成 S 形（有拐点），exp 的 CDF 没有。
- 密度有峰了（峰在 \(x=1/\sqrt2\approx0.71\)），不再像 exp 那样在 0 处最高、一路单调跌。
- 红点照旧追踪 \(u=0.8\to x=\sqrt{1.61}=1.27\to f(1.27)=0.51\)；band edges \(0,\,0.47,\,0.71,\,0.96,\,1.27\) 正是 exp 那批 \(0,\,0.22,\,0.51,\,0.92,\,1.61\) 的平方根。
k <- 2; lam <- 1
prob_breaks <- seq(0,1,by=0.2); wq <- qweibull(prob_breaks,k,lam)
band_cols <- c("#4E79A7","#59A14F","#EDC948","#F28E2B","#E15759")
u_star <- 0.8; xw_star <- qweibull(u_star,k,lam) # 1.269
tcol <- function(c,a) adjustcolor(c,alpha.f=a)
RED<-"#D62728"; DRED<-"#7a1f1f"; XM <- 3
layout(matrix(c(1,2, 3,4, 5,6),3,2,byrow=TRUE))
par(cex.lab=1.35,cex.main=1.15,cex.axis=1.0,mgp=c(3,0.8,0))
MAR <- c(6.2,5.2,3.4,3.6)
W_real <- qweibull(U,k,lam) # = sqrt(-log(1-U))
hw <- hist(W_real[W_real<=XM],breaks=seq(0,XM,by=0.06),plot=FALSE)
hws <- hw$density*mean(W_real<=XM)
u12 <- sort(runif(12)); w12 <- qweibull(u12,k,lam)
u400 <- runif(400)
# (1) pipeline text
par(mar=c(2,2,2,2))
plot(NA,xlim=c(0,1),ylim=c(0,1),xlab="",ylab="",main="",axes=FALSE)
text(0.5,0.5,paste0(
"same pipeline as Exp -- ONE new move:\n\n",
"Weibull inverse CDF = (Exp inverse CDF)^(1/k)\n",
" x = ( -log(1-u) )^(1/2) (k = 2)\n\n",
"input: U = 0.915, 0.286, ...\n",
"machine: 0.915 -> sqrt(2.46) = 1.57\n",
" 0.286 -> sqrt(0.34) = 0.58\n",
"output: histogram -> f(x) = 2x e^(-x^2)\n\n",
"sqrt squeezes the long Exp tail:\n",
"density now has a PEAK, not a cliff at 0"),cex=1.25,col="gray20")
# (2) CDF, top right
par(mar=MAR)
plot(NA,xlim=c(0,XM),ylim=c(0,1),xlab="",ylab=expression(u==F(x)),
main=expression("CDF: "*F(x)==1-e^{-x^2}),yaxs="i")
for(i in 1:5)rect(0,prob_breaks[i],XM,prob_breaks[i+1],col=tcol(band_cols[i],0.30),border=NA)
curve(pweibull(x,k,lam),0,XM,add=TRUE,lwd=2.6,col="#2E5AAC")
segments(0,u_star,xw_star,u_star,lty=2,col=RED);segments(xw_star,0,xw_star,u_star,lty=2,col=RED)
points(xw_star,u_star,pch=21,bg=RED,col="white",cex=1.3)
mtext(expression(x),side=1,line=2.4,cex=0.95)
axis(1,at=round(wq[1:5],2),labels=round(wq[1:5],2),las=2,line=3.0,cex.axis=0.85,col.axis=DRED,col=DRED)
text(2.0,0.35,"S-shaped now --\nExp's CDF had no\ninflection point",cex=1.05,col="gray25",font=2)
# (3) OUTPUT upright, middle left
par(mar=MAR)
plot(NA,xlim=c(0.95,0),ylim=c(0,XM),xlab="",ylab=expression(x==(-log(1-u))^{1/2}),
main=expression("output, x-axis vertical: "*f(x)==2*x*e^{-x^2}),yaxs="i",xaxs="i")
for(i in 1:5){
l<-wq[i];r<-min(wq[i+1],XM); xx<-seq(l,r,length.out=150)
polygon(c(0,dweibull(xx,k,lam),0),c(l,xx,r),col=tcol(band_cols[i],0.72),border=NA)
}
rect(hws,hw$breaks[-length(hw$breaks)],0,hw$breaks[-1],border="gray30",col=NA,lwd=0.7)
xg<-seq(0,XM,length.out=400); lines(dweibull(xg,k,lam),xg,lwd=2,col="gray20")
abline(h=wq[2:5],col="white",lwd=1.2)
fw_star<-dweibull(xw_star,k,lam)
segments(fw_star,0,fw_star,xw_star,lty=2,col=RED);segments(0.95,xw_star,fw_star,xw_star,lty=2,col=RED)
points(fw_star,xw_star,pch=21,bg=RED,col="white",cex=1.3)
mtext(expression(f(x)),side=1,line=2.4,cex=0.95)
axis(4,at=round(wq[1:5],2),las=1,cex.axis=0.8,col.axis=DRED,col=DRED)
mtext("band edges x",side=4,line=2.3,cex=0.80,col=DRED)
text(0.55,2.55,"same colors, same data\nas below -- not laid\nflat yet",cex=1.05,col="gray25",font=2)
# (4) MACHINE, middle right
par(mar=MAR)
plot(NA,xlim=c(0,1),ylim=c(0,XM),xlab=expression(u),ylab=expression(x==(-log(1-u))^{1/2}),
main=expression("inverse CDF: "*x==(-log(1-u))^{1/2}*" -- REAL DATA"),yaxs="i",xaxs="i")
for(i in 1:5)rect(prob_breaks[i],0,prob_breaks[i+1],XM,col=tcol(band_cols[i],0.30),border=NA)
curve(qweibull(x,k,lam),0,0.9998,add=TRUE,lwd=2.6,col="#6F42C1")
for(i in 2:5)segments(0,wq[i],prob_breaks[i],wq[i],lty=3,col="gray40")
for(j in 1:12){
segments(u12[j],0,u12[j],w12[j],col="gray35",lwd=0.9)
arrows(u12[j],w12[j],0.006,w12[j],length=0.05,col="gray35",lwd=0.9)
points(u12[j],w12[j],pch=21,bg=band_cols[findInterval(u12[j],prob_breaks,rightmost.closed=TRUE)],col="white",cex=1.0)
}
rug(u400,side=1,ticksize=0.025,col=tcol("gray20",0.35),lwd=0.6)
axis(4,at=round(wq[1:5],2),las=1,cex.axis=0.8,col.axis=DRED,col=DRED)
mtext("band edges x",side=4,line=2.3,cex=0.80,col=DRED)
segments(u_star,0,u_star,xw_star,lty=2,col=RED);segments(0,xw_star,u_star,xw_star,lty=2,col=RED)
points(u_star,xw_star,pch=21,bg=RED,col="white",cex=1.3)
text(0.40,2.65,"= the Exp machine, then sqrt:\ntall Exp values get squeezed\n(2.46 -> 1.57), so no long tail",cex=1.05,col="gray25",font=2)
# (5) OUTPUT laid flat, bottom left
par(mar=MAR)
plot(NA,xlim=c(0,XM),ylim=c(0,0.95),xlab="",ylab=expression(f(x)),
main=expression("output laid flat: X ~ Weibull(2,1), "*f(x)==2*x*e^{-x^2}),yaxs="i")
for(i in 1:5){l<-wq[i];r<-min(wq[i+1],XM);xx<-seq(l,r,length.out=150)
polygon(c(l,xx,r),c(0,dweibull(xx,k,lam),0),col=tcol(band_cols[i],0.72),border=NA)}
rect(hw$breaks[-length(hw$breaks)],0,hw$breaks[-1],hws,border="gray30",col=NA,lwd=0.7)
curve(dweibull(x,k,lam),0,XM,add=TRUE,lwd=2,col="gray20")
abline(v=wq[2:5],col="white",lwd=1.2)
segments(xw_star,0,xw_star,fw_star,lty=2,col=RED);segments(0,fw_star,xw_star,fw_star,lty=2,col=RED)
points(xw_star,fw_star,pch=21,bg=RED,col="white",cex=1.3)
text(2.05,0.62,"gray bars: histogram of\n100,000 real sqrt(-log(1-U))\npeak at x = 0.71 -- unlike Exp",cex=1.1,col="gray25",font=2)
mtext(expression(x),side=1,line=2.4,cex=0.95)
axis(1,at=round(wq[1:5],2),labels=round(wq[1:5],2),las=2,line=3.0,cex.axis=0.85,col.axis=DRED,col=DRED)
mtext("band edges = Weibull inverse-CDF values",side=1,line=5.4,cex=0.80,col=DRED)
# (6) INPUT uniform, bottom right
par(mar=MAR)
plot(NA,xlim=c(0,1),ylim=c(0,1.15),xlab=expression(u),ylab="density",main="input: U ~ Uniform(0,1)",yaxs="i")
for(i in 1:5){rect(prob_breaks[i],0,prob_breaks[i+1],1,col=tcol(band_cols[i],0.72),border="white",lwd=2)
text(mean(prob_breaks[i:(i+1)]),0.5,"20%",col="white",font=2,cex=0.9)}
h_u <- hist(U,breaks=seq(0,1,by=0.04),plot=FALSE)
rect(h_u$breaks[-length(h_u$breaks)],0,h_u$breaks[-1],h_u$density,border="gray30",col=NA,lwd=0.7)
segments(0,1,1,1,lwd=2,col="gray25")
text(0.5,1.09,"same 100,000 real U's as the Exp figure",cex=1.0,col="gray25",font=2)
# overlays
par(fig=c(0,1,0,1),new=TRUE,mar=c(0,0,0,0));plot(0:1,0:1,type="n",axes=FALSE,xlim=c(0,1),ylim=c(0,1),xaxs="i",yaxs="i")
arrows(0.82,0.702,0.82,0.661,code=3,length=0.10,lwd=2.8,col=RED)
text(0.828,0.684,"x <-> u axes flip",cex=1.3,col=RED,font=2,pos=4)
arrows(0.90,0.318,0.90,0.358,length=0.13,lwd=3.4,col=RED)
text(0.895,0.338,"1. feed REAL uniform draws",cex=1.25,col=RED,font=2,pos=2)
arrows(0.550,0.60,0.490,0.60,length=0.13,lwd=3.4,col=RED)
arrows(0.10,0.358,0.10,0.318,length=0.13,lwd=3.4,col=RED)
text(0.105,0.338,"3. rotate 90 -> lay it flat",cex=1.25,col=RED,font=2,pos=4)
k <- 2; lambda <- 1
X_wei_hand <- lambda * (-log(1 - U))^(1/k)
X_wei_qwei <- qweibull(U, shape = k, scale = lambda)
cat("Max absolute difference (should be ~0):",
max(abs(X_wei_hand - X_wei_qwei)), "\n")## Max absolute difference (should be ~0): 02.1.4 正态分布：数值逆 CDF
\(\Phi(x) = \int_{-\infty}^x
\frac{1}{\sqrt{2\pi}} e^{-t^2/2} dt\)
没有闭合形式。qnorm(u) 用有理函数近似计算 \(\Phi^{-1}(u)\)。概念完全一样：给定概率
\(u\)，找满足 \(\Phi(x) = u\) 的 \(x\)。
p_in <- c(0.025, 0.16, 0.50, 0.84, 0.975)
x_out <- qnorm(p_in)
data.frame(
p_input = p_in,
"qnorm(p)" = round(x_out, 4),
"pnorm(qnorm(p))" = round(pnorm(x_out), 6)
)## p_input qnorm.p. pnorm.qnorm.p..
## 1 0.025 -1.9600 0.025
## 2 0.160 -0.9945 0.160
## 3 0.500 0.0000 0.500
## 4 0.840 0.9945 0.840
## 5 0.975 1.9600 0.975还是那条流水线，第三遍。这次的新点只有一个：\(\Phi^{-1}\)
没有解析式，机器是数值解——qnorm(u) 在内部解方程
\(\Phi(x)=u\)。流程一根毛都没变：
- CDF 对称 S 形，在 \(x=0\) 处最陡 → 中间那条 20% 带（\(-0.25\sim0.25\)）在 \(x\) 轴上最窄，密度峰就在这。
- 两头都摊向无穷：最外两条带各自伸向 \(\pm\infty\)——正态是双尾，exp/Weibull 只有右尾。
- 红点照旧：\(u=0.8\to x=\Phi^{-1}(0.8)=0.84\to f(0.84)=0.28\)；band edges \(\pm0.84,\ \pm0.25\) 左右对称。
prob_breaks <- seq(0,1,by=0.2); nq <- qnorm(prob_breaks) # -Inf ... Inf
nq_fin <- nq[2:5] # -0.84 -0.25 0.25 0.84
band_cols <- c("#4E79A7","#59A14F","#EDC948","#F28E2B","#E15759")
u_star <- 0.8; xn_star <- qnorm(u_star); fn_star <- dnorm(xn_star)
tcol <- function(c,a) adjustcolor(c,alpha.f=a)
RED<-"#D62728"; DRED<-"#7a1f1f"; XL <- -3; XR <- 3
layout(matrix(c(1,2, 3,4, 5,6),3,2,byrow=TRUE))
par(cex.lab=1.35,cex.main=1.15,cex.axis=1.0,mgp=c(3,0.8,0))
MAR <- c(6.2,5.2,3.4,3.6)
N_real <- qnorm(U)
hn <- hist(N_real[abs(N_real)<=XR],breaks=seq(XL,XR,by=0.12),plot=FALSE)
hns <- hn$density*mean(abs(N_real)<=XR)
u12 <- sort(runif(12)); n12 <- pmin(pmax(qnorm(u12),XL),XR)
u400 <- runif(400)
# (1) pipeline text
par(mar=c(2,2,2,2))
plot(NA,xlim=c(0,1),ylim=c(0,1),xlab="",ylab="",main="",axes=FALSE)
text(0.5,0.5,paste0(
"same pipeline -- new wrinkle:\n\n",
"qnorm(u) has NO closed formula.\n",
"R solves pnorm(x) = u numerically.\n",
"concept unchanged: probability -> position\n\n",
"input: U = 0.915, 0.286, ...\n",
"machine: 0.915 -> qnorm = 1.37\n",
" 0.286 -> qnorm = -0.57\n",
"output: histogram -> bell curve\n\n",
"middle 20% band is NARROWEST (peak);\n",
"outer bands stretch into BOTH tails"),cex=1.25,col="gray20")
# (2) CDF, top right
par(mar=MAR)
plot(NA,xlim=c(XL,XR),ylim=c(0,1),xlab="",ylab=expression(u==Phi(x)),
main=expression("CDF: "*Phi(x)*" (no closed form)"),yaxs="i")
for(i in 1:5)rect(XL,prob_breaks[i],XR,prob_breaks[i+1],col=tcol(band_cols[i],0.30),border=NA)
curve(pnorm(x),XL,XR,add=TRUE,lwd=2.6,col="#2E5AAC")
segments(XL,u_star,xn_star,u_star,lty=2,col=RED);segments(xn_star,0,xn_star,u_star,lty=2,col=RED)
points(xn_star,u_star,pch=21,bg=RED,col="white",cex=1.3)
mtext(expression(x),side=1,line=2.4,cex=0.95)
axis(1,at=round(nq_fin,2),labels=round(nq_fin,2),las=2,line=3.0,cex.axis=0.85,col.axis=DRED,col=DRED)
text(1.7,0.28,"symmetric S;\nsteepest at x = 0",cex=1.05,col="gray25",font=2)
# (3) OUTPUT upright, middle left
par(mar=MAR)
plot(NA,xlim=c(0.46,0),ylim=c(XL,XR),xlab="",ylab=expression(x==qnorm(u)),
main=expression("output, x-axis vertical: bell curve"),yaxs="i",xaxs="i")
for(i in 1:5){
l<-max(nq[i],XL);r<-min(nq[i+1],XR); xx<-seq(l,r,length.out=150)
polygon(c(0,dnorm(xx),0),c(l,xx,r),col=tcol(band_cols[i],0.72),border=NA)
}
rect(hns,hn$breaks[-length(hn$breaks)],0,hn$breaks[-1],border="gray30",col=NA,lwd=0.7)
xg<-seq(XL,XR,length.out=400); lines(dnorm(xg),xg,lwd=2,col="gray20")
abline(h=nq_fin,col="white",lwd=1.2)
segments(fn_star,XL,fn_star,xn_star,lty=2,col=RED);segments(0.46,xn_star,fn_star,xn_star,lty=2,col=RED)
points(fn_star,xn_star,pch=21,bg=RED,col="white",cex=1.3)
mtext(expression(f(x)),side=1,line=2.4,cex=0.95)
axis(4,at=round(nq_fin,2),las=1,cex.axis=0.8,col.axis=DRED,col=DRED)
mtext("band edges x",side=4,line=2.3,cex=0.80,col=DRED)
text(0.30,2.35,"same colors, same data\nas below -- not laid\nflat yet",cex=1.05,col="gray25",font=2)
# (4) MACHINE, middle right
par(mar=MAR)
plot(NA,xlim=c(0,1),ylim=c(XL,XR),xlab=expression(u),ylab=expression(x==qnorm(u)),
main=expression("inverse CDF: "*x==qnorm(u)*" -- REAL DATA"),yaxs="i",xaxs="i")
for(i in 1:5)rect(prob_breaks[i],XL,prob_breaks[i+1],XR,col=tcol(band_cols[i],0.30),border=NA)
curve(qnorm(x),0.0015,0.9985,add=TRUE,lwd=2.6,col="#6F42C1")
for(i in 2:5)segments(0,nq[i],prob_breaks[i],nq[i],lty=3,col="gray40")
for(j in 1:12){
segments(u12[j],XL,u12[j],n12[j],col="gray35",lwd=0.9)
arrows(u12[j],n12[j],0.006,n12[j],length=0.05,col="gray35",lwd=0.9)
points(u12[j],n12[j],pch=21,bg=band_cols[findInterval(u12[j],prob_breaks,rightmost.closed=TRUE)],col="white",cex=1.0)
}
rug(u400,side=1,ticksize=0.025,col=tcol("gray20",0.35),lwd=0.6)
axis(4,at=round(nq_fin,2),las=1,cex.axis=0.8,col.axis=DRED,col=DRED)
mtext("band edges x",side=4,line=2.3,cex=0.80,col=DRED)
segments(u_star,XL,u_star,xn_star,lty=2,col=RED);segments(0,xn_star,u_star,xn_star,lty=2,col=RED)
points(u_star,xn_star,pch=21,bg=RED,col="white",cex=1.3)
text(0.50,2.45,"flat middle = crowded x near 0;\nsteep ends = sparse tails\n(both directions)",cex=1.05,col="gray25",font=2)
# (5) OUTPUT laid flat, bottom left
par(mar=MAR)
plot(NA,xlim=c(XL,XR),ylim=c(0,0.46),xlab="",ylab=expression(f(x)),
main=expression("output laid flat: X ~ N(0,1), bell curve"),yaxs="i")
for(i in 1:5){l<-max(nq[i],XL);r<-min(nq[i+1],XR);xx<-seq(l,r,length.out=150)
polygon(c(l,xx,r),c(0,dnorm(xx),0),col=tcol(band_cols[i],0.72),border=NA)}
rect(hn$breaks[-length(hn$breaks)],0,hn$breaks[-1],hns,border="gray30",col=NA,lwd=0.7)
curve(dnorm(x),XL,XR,add=TRUE,lwd=2,col="gray20")
abline(v=nq_fin,col="white",lwd=1.2)
segments(xn_star,0,xn_star,fn_star,lty=2,col=RED);segments(XL,fn_star,xn_star,fn_star,lty=2,col=RED)
points(xn_star,fn_star,pch=21,bg=RED,col="white",cex=1.3)
text(-1.85,0.36,"gray bars: histogram of\n100,000 real qnorm(U)\npeak f(0) = 0.40",cex=1.1,col="gray25",font=2)
mtext(expression(x),side=1,line=2.4,cex=0.95)
axis(1,at=round(nq_fin,2),labels=round(nq_fin,2),las=2,line=3.0,cex.axis=0.85,col.axis=DRED,col=DRED)
mtext("band edges = qnorm(0.2, 0.4, 0.6, 0.8)",side=1,line=5.4,cex=0.80,col=DRED)
# (6) INPUT uniform, bottom right
par(mar=MAR)
plot(NA,xlim=c(0,1),ylim=c(0,1.15),xlab=expression(u),ylab="density",main="input: U ~ Uniform(0,1)",yaxs="i")
for(i in 1:5){rect(prob_breaks[i],0,prob_breaks[i+1],1,col=tcol(band_cols[i],0.72),border="white",lwd=2)
text(mean(prob_breaks[i:(i+1)]),0.5,"20%",col="white",font=2,cex=0.9)}
h_u <- hist(U,breaks=seq(0,1,by=0.04),plot=FALSE)
rect(h_u$breaks[-length(h_u$breaks)],0,h_u$breaks[-1],h_u$density,border="gray30",col=NA,lwd=0.7)
segments(0,1,1,1,lwd=2,col="gray25")
text(0.5,1.09,"same 100,000 real U's as before",cex=1.0,col="gray25",font=2)
# overlays
par(fig=c(0,1,0,1),new=TRUE,mar=c(0,0,0,0));plot(0:1,0:1,type="n",axes=FALSE,xlim=c(0,1),ylim=c(0,1),xaxs="i",yaxs="i")
arrows(0.82,0.702,0.82,0.661,code=3,length=0.10,lwd=2.8,col=RED)
text(0.828,0.684,"x <-> u axes flip",cex=1.3,col=RED,font=2,pos=4)
arrows(0.90,0.318,0.90,0.358,length=0.13,lwd=3.4,col=RED)
text(0.895,0.338,"1. feed REAL uniform draws",cex=1.25,col=RED,font=2,pos=2)
arrows(0.550,0.60,0.490,0.60,length=0.13,lwd=3.4,col=RED)
arrows(0.10,0.358,0.10,0.318,length=0.13,lwd=3.4,col=RED)
text(0.105,0.338,"3. rotate 90 -> lay it flat",cex=1.25,col=RED,font=2,pos=4)
2.1.5 验证：用逆变换生成三种分布
X_exp_inv <- qexp(U, rate = 1)
X_norm_inv <- qnorm(U)
X_wei_inv <- qweibull(U, shape = 2, scale = 1)
band_cols <- c("#4E79A7","#59A14F","#EDC948","#F28E2B","#E15759")
tcol <- function(c,a) adjustcolor(c,alpha.f=a)
p5 <- seq(0,1,by=0.2)
# One panel: histogram with bars colored by their 20% probability band,
# gray theory curve on top, red dashed marker at F^-1(0.8).
band_hist <- function(x, qfun, dfun, main, xlim, bw){
h <- hist(x[x>=xlim[1] & x<=xlim[2]], breaks=seq(xlim[1],xlim[2],by=bw), plot=FALSE)
qs <- qfun(p5) # band edges in x
mids <- (h$breaks[-length(h$breaks)]+h$breaks[-1])/2
cols <- tcol(band_cols[pmin(pmax(findInterval(mids,qs),1),5)],0.72)
dens <- h$density*mean(x>=xlim[1] & x<=xlim[2])
plot(NA,xlim=xlim,ylim=c(0,max(dens)*1.12),xlab=expression(x),ylab="density",
main=main,yaxs="i")
rect(h$breaks[-length(h$breaks)],0,h$breaks[-1],dens,col=cols,border="white",lwd=0.4)
curve(dfun(x),xlim[1],xlim[2],add=TRUE,lwd=2.2,col="gray20")
xs <- qfun(0.8)
segments(xs,0,xs,dfun(xs),lty=2,col="#D62728",lwd=1.6)
points(xs,dfun(xs),pch=21,bg="#D62728",col="white",cex=1.2)
}
par(mfrow=c(1,3),mar=c(4.6,4.6,3.5,1),cex.lab=1.35,cex.main=1.15,cex.axis=1.0)
band_hist(X_exp_inv, function(p)qexp(p,1), function(x)dexp(x,1),
expression(F^{-1}*(U) %->% "Exp(1)"), c(0,6), 0.12)
band_hist(X_norm_inv, qnorm, dnorm,
expression(F^{-1}*(U) %->% "Normal(0,1)"), c(-3,3), 0.12)
band_hist(X_wei_inv, function(p)qweibull(p,2,1), function(x)dweibull(x,2,1),
expression(F^{-1}*(U) %->% "Weibull(2,1)"), c(0,3), 0.06)
| Distribution | \(F^{-1}(u)\) | Method |
|---|---|---|
| Exp(\(\lambda\)) | \(-\frac{1}{\lambda}\log(1-u)\) | Analytic |
| Weibull(\(k,\lambda\)) | \(\lambda(-\log(1-u))^{1/k}\) | Analytic |
| Normal | No closed form | qnorm (numerical) |
2.2 幂变换：Exp → Weibull
Weibull 不仅可以通过逆 CDF 得到，也可以从 Exp 直接做幂变换：
若 \(X \sim \text{Exp}(1)\)，则 \(Y = X^{1/k} \sim \text{Weibull}(k, 1)\)。
Weibull 的核心概念是风险率（hazard rate）：
\[h(t) = \frac{k}{\lambda}\left(\frac{t}{\lambda}\right)^{k-1}\]
- \(k < 1\)：风险递减 → 早期失效（婴儿死亡率）
- \(k = 1\)：风险恒定 → 指数分布（无记忆）
- \(k > 1\)：风险递增 → 老化磨损
x_seq <- seq(0.01, 3, length.out = 500)
ks <- c(0.5, 1, 2, 3, 5)
cols <- c("purple", "steelblue", "green3", "orange", "red")
par(mfrow = c(1, 3), mar = c(4, 4, 3.5, 1))
# Density
plot(NA, xlim = c(0, 3), ylim = c(0, 2),
main = "Weibull density", xlab = "t", ylab = "Density")
for (i in seq_along(ks))
lines(x_seq, dweibull(x_seq, ks[i], 1), col = cols[i], lwd = 2)
legend("topright", paste("k =", ks), col = cols, lwd = 2, bty = "n")
# Hazard rate
hazard_w <- function(x, k) k * x^(k - 1)
plot(NA, xlim = c(0, 3), ylim = c(0, 8),
main = "Weibull hazard rate h(t)", xlab = "t", ylab = "h(t)")
for (i in seq_along(ks))
lines(x_seq, hazard_w(x_seq, ks[i]), col = cols[i], lwd = 2)
abline(h = 0, col = "gray", lty = 2)
legend("topright", paste("k =", ks), col = cols, lwd = 2, bty = "n")
# Verify: Exp(1)^{1/k} ~ Weibull(k,1)
X_e <- rexp(10000)
k_demo <- 3
Y_w <- X_e^(1/k_demo)
hist(Y_w, breaks = 80, col = "steelblue", freq = FALSE, border = "white",
main = bquote("Exp(1)"^{1/.(k_demo)} ~ "~ Weibull(" ~ .(k_demo) ~ ", 1)"),
xlab = "y")
curve(dweibull(x, shape = k_demo, scale = 1), add = TRUE, col = "tomato", lwd = 2)
2.3 指数/对数变换：Normal ↔︎ LogNormal
若 \(X \sim N(\mu, \sigma^2)\)，则 \(Y = e^X \sim \text{LogNormal}(\mu, \sigma^2)\)。
现实中，当多个乘法效应叠加时（收入、文件大小、生物量），变量呈对数正态分布——因为取对数后就变成加法，CLT 生效。
mu <- 0; sigma <- 1
X_norm <- rnorm(10000, mu, sigma)
Y_lnorm <- exp(X_norm)
par(mfrow = c(1, 3), mar = c(4, 4, 3.5, 1))
hist(X_norm, breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = "X ~ Normal(0,1)", xlab = "x")
curve(dnorm(x, mu, sigma), add = TRUE, col = "tomato", lwd = 2)
hist(Y_lnorm, breaks = 120, col = "steelblue", freq = FALSE, border = "white",
main = "Y = exp(X) ~ LogNormal(0,1)", xlab = "y",
xlim = c(0, quantile(Y_lnorm, 0.98)))
curve(dlnorm(x, mu, sigma), add = TRUE, col = "tomato", lwd = 2)
hist(log(Y_lnorm), breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = "log(Y) restores Normal", xlab = "log(y)")
curve(dnorm(x, mu, sigma), add = TRUE, col = "tomato", lwd = 2)
2.4 变换之间的层级关系
三种变换不是孤立的，它们有层级：
- 逆 CDF 是万能方法（任何分布都可以）
- 幂变换 是逆 CDF 的特例：Weibull 的 \(F^{-1}(u) = (-\log(1-u))^{1/k}\) 就是先做 Exp 的逆 CDF \((-\log(1-u))\)，再做幂 \(({\cdot})^{1/k}\)
- 指数变换 \(e^X\) 把加法世界（Normal, CLT）映射到乘法世界（LogNormal）
\[\text{Uniform} \xrightarrow{-\log(1-U)} \text{Exp} \xrightarrow{(\cdot)^{1/k}} \text{Weibull}\] \[\text{Uniform} \xrightarrow{\texttt{qnorm}} \text{Normal} \xrightarrow{e^{(\cdot)}} \text{LogNormal}\]
2.5 变换小结
par(mar = c(1, 1, 2, 1))
plot(NA, xlim = c(0, 10), ylim = c(0, 4), axes = FALSE, bty = "n",
xlab = "", ylab = "", main = "Transform operations summary")
draw_box <- function(x, y, label, col, w = 2, h = 0.5) {
rect(x-w/2, y-h/2, x+w/2, y+h/2, col = col, border = "white", lwd = 2)
text(x, y, label, col = "white", font = 2, cex = 0.75)
}
arr <- function(x0, y0, x1, y1, lab = "") {
arrows(x0, y0, x1, y1, length = 0.10, lwd = 1.5, col = "gray40")
text((x0+x1)/2, (y0+y1)/2 + 0.2, lab, cex = 0.6, col = "gray30", font = 3)
}
draw_box(1.5, 3, "Uniform(0,1)", "gray30")
draw_box(5, 3, "Exp(1)", "#2E86AB")
draw_box(8.5, 3, "Weibull(k,1)", "#2E86AB")
draw_box(5, 1, "Normal(0,1)", "#2E86AB")
draw_box(8.5, 1, "LogNormal", "#2E86AB")
arr(2.5, 3, 4, 3, "-log(1-U)")
arr(6, 3, 7.5, 3, "X^(1/k)")
arr(1.5, 2.7, 5, 1.3, "qnorm(U)")
arr(6, 1, 7.5, 1, "exp(X)")
3 运算二：求和——多个随机变量叠加
3.1 中心极限定理（CLT）
不管原始分布是什么形状，只要有有限方差，大量独立样本的均值趋向正态。
ns <- c(1, 2, 5, 30)
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
for (n in ns) {
sums <- replicate(5000, sum(runif(n)))
z <- (sums - n * 0.5) / sqrt(n / 12)
hist(z, breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = paste0("Sum of ", n, " Uniform (standardized)"),
xlab = "z", xlim = c(-4, 4), ylim = c(0, 0.46))
curve(dnorm(x), add = TRUE, col = "tomato", lwd = 2)
}
3.1.1 CLT 与原始分布无关
n <- 30; M <- 5000
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
# Exponential (heavy right skew)
z_exp <- scale(replicate(M, mean(rexp(n, 1))))
hist(z_exp, breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = "Mean of 30 Exp(1)", xlab = "z", ylim = c(0, 0.45))
curve(dnorm(x), add = TRUE, col = "tomato", lwd = 2)
# Poisson (discrete)
z_pois <- scale(replicate(M, mean(rpois(n, 3))))
hist(z_pois, breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = "Mean of 30 Poisson(3)", xlab = "z", ylim = c(0, 0.45))
curve(dnorm(x), add = TRUE, col = "tomato", lwd = 2)
# Bernoulli (extreme skew)
z_bern <- scale(replicate(M, mean(rbinom(n, 1, 0.3))))
hist(z_bern, breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = "Mean of 30 Bernoulli(0.3)", xlab = "z", ylim = c(0, 0.45))
curve(dnorm(x), add = TRUE, col = "tomato", lwd = 2)
3.2 泊松过程：Exp 与 Poisson 的双面
同一个物理过程（事件以恒定速率 \(\lambda\) 随机到达）同时产生两种分布：
- 数事件（固定时间窗）→ Poisson(\(\lambda t\))
- 量间隔（相邻事件）→ Exp(\(\lambda\))
两者共享同一个参数 \(\lambda\)。
lambda <- 3
T_total <- 200
n_arr <- rpois(1, lambda * T_total) * 3
inter <- rexp(n_arr, rate = lambda)
arrivals <- cumsum(inter)
arrivals <- arrivals[arrivals <= T_total]
par(mfrow = c(1, 3), mar = c(4, 4, 3.5, 1))
# Inter-arrival times → Exp
hist(inter[1:1000], breaks = 60, col = "steelblue", freq = FALSE,
border = "white",
main = "Inter-arrival time ~ Exp(lambda)", xlab = "Time",
xlim = c(0, quantile(inter, 0.99)))
curve(dexp(x, rate = lambda), add = TRUE, col = "tomato", lwd = 2)
legend("topright", "Exp(3) theory", col = "tomato", lwd = 2, bty = "n")
# Counts per unit time → Poisson
counts <- tabulate(ceiling(arrivals[arrivals > 0]))
counts <- counts[counts > 0]
tbl <- table(counts) / length(counts)
barplot(tbl, col = "steelblue", border = "white",
main = "Events per unit time ~ Poisson(lambda)",
xlab = "Count", ylab = "Proportion")
k_vals <- 0:12
points(k_vals + 0.5, dpois(k_vals, lambda),
col = "tomato", pch = 16, cex = 1.3)
legend("topright", "Poisson(3) theory", col = "tomato", pch = 16, bty = "n")
# Timeline
n_show <- min(50, length(arrivals))
plot(arrivals[1:n_show], rep(1, n_show),
pch = "|", cex = 1.5, col = "steelblue",
main = "Poisson process event timeline", xlab = "Time",
ylab = "", yaxt = "n", ylim = c(0.5, 1.5))
3.3 k 个 Exp 求和 → Gamma
等待第 \(k\) 个事件 = 等 \(k\) 段间隔之和 → \(\text{Gamma}(k, \lambda)\)。
这是 Gamma 分布最直接的物理意义：泊松过程中，第 \(k\) 次事件发生的时刻。
n_sim <- 10000
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
for (k in c(1, 2, 5, 10, 20, 50)) {
sums <- rowSums(matrix(rexp(n_sim * k, rate = 1), nrow = n_sim))
hist(sums, breaks = 80, col = "steelblue", freq = FALSE, border = "white",
main = paste0("Sum of ", k, " Exp(1) = Gamma(", k, ", 1)"),
xlab = "x", ylab = "")
curve(dgamma(x, shape = k, rate = 1), add = TRUE, col = "tomato", lwd = 2)
}
注意 \(k\) 增大时 Gamma 趋向正态——CLT 再次生效。
精确的包含关系：Exp(\(\lambda\)) = Gamma(1, \(\lambda\))。
par(mfrow = c(1, 2), mar = c(4, 4, 3.5, 1))
x_seq <- seq(0.01, 6, length.out = 500)
plot(x_seq, dgamma(x_seq, shape = 1, rate = 1), type = "l",
col = "steelblue", lwd = 3,
main = "Exp(lambda) = Gamma(1, lambda)\n(exact equivalence)",
xlab = "x", ylab = "Density")
lines(x_seq, dexp(x_seq, rate = 1), col = "tomato", lwd = 2, lty = 2)
legend("topright", c("Gamma(1,1)", "Exp(1)"),
col = c("steelblue", "tomato"), lwd = 2, lty = c(1, 2), bty = "n")
# Gamma family
cols_g <- c("steelblue", "green3", "orange", "red", "purple")
plot(NA, xlim = c(0, 12), ylim = c(0, 0.5),
main = "Gamma(k, 1) family", xlab = "x", ylab = "Density")
for (i in seq_along(c(1, 2, 5, 10, 20))) {
k <- c(1, 2, 5, 10, 20)[i]
lines(x_seq, dgamma(seq(0.01, 12, length.out = 500), k, 1),
col = cols_g[i], lwd = 2)
}
legend("topright", paste("shape =", c(1,2,5,10,20)),
col = cols_g, lwd = 2, bty = "n")
3.4 k 个 Normal² 求和 → Chi-squared
若 \(Z_1, \ldots, Z_k \overset{\text{iid}}{\sim} N(0,1)\)，则：
\[\chi^2_k = \sum_{i=1}^k Z_i^2 \sim \chi^2(k) = \text{Gamma}\!\left(\frac{k}{2},\, \frac{1}{2}\right)\]
精确等价：Chi-squared 是 Gamma 分布的特例。这意味着求和这一节的所有分布最终都汇聚到 Gamma：
- Exp = Gamma(1, \(\lambda\))
- \(k\) 个 Exp 之和 = Gamma(\(k\), \(\lambda\))
- Chi-squared(\(k\)) = Gamma(\(k/2\), \(1/2\))
- \(k \to \infty\) 时 Gamma → Normal（CLT 闭环）
n_sim <- 10000
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
for (k in c(1, 2, 3, 5, 10, 30)) {
Z <- matrix(rnorm(n_sim * k), nrow = n_sim)
chi2 <- rowSums(Z^2)
hist(chi2, breaks = 80, col = "steelblue", freq = FALSE, border = "white",
main = paste0("Chi-sq(", k, ") = sum of ", k, " Z^2"),
xlab = "x", ylab = "")
curve(dchisq(x, df = k), add = TRUE, col = "tomato", lwd = 2)
}
3.5 求和小结
par(mar = c(1, 1, 2, 1))
plot(NA, xlim = c(0, 12), ylim = c(0, 3.5), axes = FALSE, bty = "n",
xlab = "", ylab = "", main = "Sum operations: chain of connections")
draw_box <- function(x, y, label, col, w = 2.2, h = 0.5) {
rect(x-w/2, y-h/2, x+w/2, y+h/2, col = col, border = "white", lwd = 2)
text(x, y, label, col = "white", font = 2, cex = 0.7)
}
arr <- function(x0, y0, x1, y1, lab = "") {
arrows(x0, y0, x1, y1, length = 0.10, lwd = 1.5, col = "gray40")
text((x0+x1)/2, (y0+y1)/2 + 0.2, lab, cex = 0.6, col = "gray30", font = 3)
}
draw_box(2, 3, "Any dist", "gray30")
draw_box(6, 3, "Normal(0,1)", "#A23B72")
draw_box(2, 1, "Exp(lambda)", "#2E86AB")
draw_box(6, 1, "Gamma(k, lam)","#A23B72")
draw_box(10, 1, "Chi-sq(k)", "#A23B72")
draw_box(10, 3, "Poisson(lam*t)","#A23B72")
arr(3.1, 3, 4.9, 3, "CLT: mean of n")
arr(2, 2.7, 2, 1.3, "k=1 special case")
arr(3.1, 1, 4.9, 1, "sum of k copies")
arr(6, 2.7, 6, 1.3, "Z^2 summed")
arr(7.1, 1, 8.9, 1, "= Gamma(k/2, 1/2)")
arr(2, 1.3, 10, 2.7, "count events")
4 运算三：比值——用于统计推断
前两节的分布描述”世界本身”。这一节的分布描述“我们对世界的不确定性”——它们不是数据的分布，而是统计量的分布。
当我们用样本估计总体参数时，样本均值服从正态（CLT），但样本方差引入额外随机性。比值运算正是用来处理这个问题。
4.1 t 分布：Normal / sqrt(Chi-squared / k)
\[t_k = \frac{Z}{\sqrt{\chi^2_k / k}}, \quad Z \perp \chi^2_k\]
现实意义：当你用样本均值估计总体均值，但不知道总体方差时，\(\bar{X}\) 的标准化量服从 \(t(n-1)\) 而非 \(N(0,1)\)。
t 分布比正态尾部更重——直觉：分母是一个随机量，有时偏小，就把比值推向极端。自由度 \(k \to \infty\) 时趋向正态（因为样本方差趋向稳定）。
x_seq <- seq(-6, 6, length.out = 600)
par(mfrow = c(1, 2), mar = c(4, 4, 3.5, 1))
# Density comparison
cols_t <- c("red", "orange", "green3", "steelblue")
dfs <- c(1, 2, 5, 30)
plot(x_seq, dnorm(x_seq), type = "l", col = "black", lwd = 3,
main = "t distribution vs Normal", xlab = "x", ylab = "Density",
ylim = c(0, 0.42))
for (i in seq_along(dfs))
lines(x_seq, dt(x_seq, dfs[i]), col = cols_t[i], lwd = 2)
legend("topright", c("Normal", paste0("t(", dfs, ")")),
col = c("black", cols_t), lwd = 2, bty = "n")
# Right tail (log scale)
x_tail <- x_seq[x_seq > 1.5]
plot(x_tail, dnorm(x_tail), type = "l", col = "black", lwd = 3,
main = "Right tail (log scale)", xlab = "x", ylab = "Density (log)", log = "y")
for (i in seq_along(dfs))
lines(x_tail, dt(x_tail, dfs[i]), col = cols_t[i], lwd = 2)
legend("topright", c("Normal", paste0("t(", dfs, ")")),
col = c("black", cols_t), lwd = 2, bty = "n")
4.1.1 从零模拟 t 分布
k <- 5; n_sim <- 10000
Z <- rnorm(n_sim)
chi2 <- rowSums(matrix(rnorm(n_sim * k)^2, nrow = n_sim))
t_sim <- Z / sqrt(chi2 / k)
par(mfrow = c(1, 1), mar = c(4, 4, 3.5, 1))
hist(t_sim, breaks = 100, col = "steelblue", freq = FALSE, border = "white",
main = paste0("Simulated t(", k, ") = Z / sqrt(Chi-sq(", k, ") / ", k, ")"),
xlab = "t", xlim = c(-6, 6))
curve(dt(x, df = k), add = TRUE, col = "tomato", lwd = 2)
legend("topright", c("Simulation", paste0("t(", k, ") theory")),
fill = c("steelblue", NA), border = c("white", NA),
lty = c(NA, 1), col = c(NA, "tomato"), lwd = c(NA, 2), bty = "n")
4.2 F 分布：两个 Chi-squared 之比
\[F_{m,n} = \frac{\chi^2_m / m}{\chi^2_n / n}\]
F 分布是 ANOVA 的基础：当零假设成立时，两个方差估计量的比值服从 F 分布。
n_sim <- 10000
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
for (params in list(c(1,1), c(2,1), c(5,2), c(10,5), c(20,10), c(50,50))) {
m <- params[1]; n <- params[2]
chi2_m <- rowSums(matrix(rnorm(n_sim * m)^2, nrow = n_sim))
chi2_n <- rowSums(matrix(rnorm(n_sim * n)^2, nrow = n_sim))
f_sim <- (chi2_m / m) / (chi2_n / n)
xlim_val <- quantile(f_sim, 0.97)
hist(f_sim[f_sim <= xlim_val], breaks = 80, col = "steelblue",
freq = FALSE, border = "white",
main = paste0("F(", m, ", ", n, ")"), xlab = "x", ylab = "")
curve(df(x, df1 = m, df2 = n), add = TRUE, col = "tomato", lwd = 2)
}
4.3 正态族关系图
par(mar = c(1, 1, 2, 1))
plot(NA, xlim = c(0, 10), ylim = c(0, 6), axes = FALSE, bty = "n",
xlab = "", ylab = "", main = "Normal family: all from Z ~ N(0,1)")
draw_box <- function(x, y, label, col = "#A23B72", w = 2, h = 0.55) {
rect(x-w/2, y-h/2, x+w/2, y+h/2, col = col, border = "white", lwd = 2)
text(x, y, label, col = "white", font = 2, cex = 0.8)
}
arr <- function(x0, y0, x1, y1, lab = "") {
arrows(x0, y0, x1, y1, length = 0.10, lwd = 1.5, col = "gray40")
text((x0+x1)/2 + 0.15, (y0+y1)/2 + 0.15,
lab, cex = 0.65, col = "gray30", font = 3)
}
draw_box(5, 5.3, "Normal(0,1)")
draw_box(5, 3.3, "Chi-squared(k)")
draw_box(2, 1.3, "t(k)")
draw_box(8, 1.3, "F(m, n)")
arr(5, 5, 5, 3.6, "sum of k squares")
arr(3.8, 3.1, 2.5, 1.6, "Z / sqrt(chi2/k)")
arr(6.2, 3.1, 7.5, 1.6, "chi2_m/m / chi2_n/n")
text(5, 4.3, "k -> inf: t(k) -> Normal", cex = 0.65, col = "gray40", font = 3)
5 运算四：取极值——样本的边界行为
5.1 次序统计量 → Beta
若 \(X_1, \ldots, X_n \overset{\text{iid}}{\sim} \text{Uniform}(0,1)\)，则第 \(k\) 小的样本：
\[X_{(k)} \sim \text{Beta}(k,\, n-k+1)\]
推导（以最大值为例）：\(P(X_{(n)} \leq x) = P(\text{all} \leq x) = x^n\)，求导得 \(f(x) = nx^{n-1}\)，即 Beta(n,1)。
n <- 10; n_sim <- 10000
M <- matrix(runif(n_sim * n), nrow = n_sim, ncol = n)
row_min <- apply(M, 1, min)
row_max <- apply(M, 1, max)
par(mfrow = c(1, 2), mar = c(4, 4, 3.5, 1))
hist(row_min, breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = paste0("min of ", n, " Uniform(0,1)"), xlab = "x")
curve(dbeta(x, 1, n), add = TRUE, col = "tomato", lwd = 2)
legend("topright", "Beta(1, n)", col = "tomato", lwd = 2, bty = "n")
hist(row_max, breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = paste0("max of ", n, " Uniform(0,1)"), xlab = "x")
curve(dbeta(x, n, 1), add = TRUE, col = "tomato", lwd = 2)
legend("topleft", "Beta(n, 1)", col = "tomato", lwd = 2, bty = "n")
5.1.1 一般次序统计量
n <- 20; n_sim <- 10000
M_sorted <- t(apply(matrix(runif(n_sim * n), nrow = n_sim), 1, sort))
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
for (k in c(1, 5, 10, 15, 18, 20)) {
hist(M_sorted[, k], breaks = 60, col = "steelblue", freq = FALSE,
border = "white",
main = paste0("X_(", k, ") ~ Beta(", k, ", ", n-k+1, ")"),
xlab = "x", ylab = "")
curve(dbeta(x, k, n - k + 1), add = TRUE, col = "tomato", lwd = 2)
}
5.2 极值分布：大量样本的最大值
5.2.1 Fisher-Tippett-Gnedenko 定理
对 \(n\) 个 i.i.d. 样本的最大值，经适当线性标准化后，只能收敛到三种分布之一：
| Type | Name | Tail of original | Example |
|---|---|---|---|
| I | Gumbel | Light tail (exp decay) | Normal, Gamma, Exp |
| II | Frechet | Heavy tail (power law) | Pareto, t dist |
| III | Weibull (extreme) | Bounded support | Uniform, Beta |
5.2.2 正态最大值 → Gumbel
n_sim <- 5000
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
for (n in c(10, 50, 200, 1000)) {
max_vals <- replicate(n_sim, max(rnorm(n)))
b_n <- qnorm(1 - 1/n)
a_n <- 1 / (n * dnorm(b_n))
z <- (max_vals - b_n) / a_n
hist(z, breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = paste0("max(Normal) n=", n, " -> Gumbel"),
xlab = "z (standardized)", xlim = c(-3, 8))
curve(dgumbel(x), add = TRUE, col = "tomato", lwd = 2)
}
5.2.3 不同尾部 → 不同极值族
n <- 200
n_sim <- 5000
par(mfrow = c(1, 3), mar = c(4, 4, 3.5, 1))
# I. Normal → Gumbel (light tail)
max_norm <- replicate(n_sim, max(rnorm(n)))
b_n <- qnorm(1 - 1/n); a_n <- 1/(n * dnorm(b_n))
z_norm <- (max_norm - b_n) / a_n
hist(z_norm, breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = "max(Normal) -> Gumbel\n(light tail)",
xlab = "z", xlim = c(-3, 8))
curve(dgumbel(x), add = TRUE, col = "tomato", lwd = 2)
# II. Pareto(2) → Frechet (heavy tail)
max_pareto <- replicate(n_sim, max((1-runif(n))^(-1/2)))
a_n_p <- n^(1/2)
z_pareto <- max_pareto / a_n_p
z_pareto_trim <- z_pareto[z_pareto < quantile(z_pareto, 0.98)]
hist(z_pareto_trim, breaks = 60, col = "steelblue", freq = FALSE,
border = "white",
main = "max(Pareto a=2) -> Frechet\n(heavy tail)",
xlab = "z / n^(1/a)")
frechet_density <- function(x, alpha = 2) {
alpha * x^(-(alpha+1)) * exp(-x^(-alpha))
}
curve(frechet_density(x), add = TRUE, col = "tomato", lwd = 2)
# III. Uniform → Weibull (bounded)
max_unif <- replicate(n_sim, max(runif(n)))
z_unif <- n * (max_unif - 1)
hist(z_unif, breaks = 60, col = "steelblue", freq = FALSE, border = "white",
main = "max(Uniform) -> Weibull (EVT)\n(bounded support)",
xlab = "n * (max - 1)")
curve(exp(x) * (x <= 0), add = TRUE, col = "tomato", lwd = 2)
5.3 重尾：Pareto 与幂律
指数尾：\(P(X > x) \sim e^{-\lambda x}\)，下降极快。 幂律尾：\(P(X > x) \sim x^{-\alpha}\)，下降缓慢，极端值远比直觉频繁。
与玻尔兹曼的联系：Boltzmann 分布 \(P \propto e^{-E/kT}\) 是指数尾的物理实例。Arrhenius 反应速率 \(k \propto e^{-E_a/kT}\) 就是这个尾巴的概率——详见
exp_normal_boltzmann.Rmd第 4 节。
在对数-对数坐标下，幂律尾部呈直线（斜率 \(= -\alpha\)）。
n_sim <- 50000
alpha <- 2
x_exp <- rexp(n_sim, 1)
x_lnorm <- rlnorm(n_sim, 0, 1)
x_pareto <- (1 - runif(n_sim))^(-1/alpha)
ecdf_sf <- function(x) {
x_sorted <- sort(x)
sf <- 1 - (1:length(x)) / length(x)
list(x = x_sorted, sf = sf)
}
sf_exp <- ecdf_sf(x_exp)
sf_lnorm <- ecdf_sf(x_lnorm)
sf_pareto <- ecdf_sf(x_pareto)
par(mfrow = c(1, 2), mar = c(4, 4, 3.5, 1))
# Linear scale
plot(sf_exp$x[sf_exp$x < 8], sf_exp$sf[sf_exp$x < 8],
type = "l", col = "steelblue", lwd = 2,
main = "Survival function P(X > x)", xlab = "x", ylab = "P(X > x)")
lines(sf_lnorm$x[sf_lnorm$x < 8], sf_lnorm$sf[sf_lnorm$x < 8],
col = "green3", lwd = 2)
lines(sf_pareto$x[sf_pareto$x < 8], sf_pareto$sf[sf_pareto$x < 8],
col = "red", lwd = 2)
legend("topright", c("Exp(1)", "LogNormal(0,1)", paste0("Pareto(a=", alpha, ")")),
col = c("steelblue","green3","red"), lwd = 2, bty = "n")
# Log-log: Pareto becomes a straight line
x_range <- c(1, 100)
idx_e <- sf_exp$x > 1
idx_l <- sf_lnorm$x > 1
idx_p <- sf_pareto$x > 1
plot(sf_exp$x[idx_e], sf_exp$sf[idx_e],
type = "l", col = "steelblue", lwd = 2, log = "xy",
main = "Log-log: Pareto = straight line",
xlab = "x (log)", ylab = "P(X > x) (log)",
xlim = x_range, ylim = c(1e-4, 1))
lines(sf_lnorm$x[idx_l], sf_lnorm$sf[idx_l], col = "green3", lwd = 2)
lines(sf_pareto$x[idx_p], sf_pareto$sf[idx_p], col = "red", lwd = 2)
abline(a = 0, b = -alpha, col = "red", lwd = 1, lty = 2)
text(10, 0.01, paste("slope = -", alpha), col = "red", cex = 0.8)
legend("topright", c("Exp(1)", "LogNormal(0,1)", paste0("Pareto(a=", alpha, ")")),
col = c("steelblue","green3","red"), lwd = 2, bty = "n")
5.3.1 Hill 估计量：从数据估计尾部指数
\[\hat{\alpha} = \left(\frac{1}{k}\sum_{i=1}^{k} \log X_{(n-i+1)} - \log X_{(n-k)}\right)^{-1}\]
hill_alpha <- function(x, k) {
x_sorted <- sort(x, decreasing = TRUE)
xi_hat <- mean(log(x_sorted[1:k]) - log(x_sorted[k + 1]))
1 / xi_hat
}
alpha_true <- 2
x_pareto <- (1 - runif(10000))^(-1/alpha_true)
k_seq <- 10:800
hill_est <- sapply(k_seq, function(k) hill_alpha(x_pareto, k))
par(mfrow = c(1, 1), mar = c(4, 4, 3.5, 1))
plot(k_seq, hill_est, type = "l", col = "steelblue", lwd = 2,
main = paste0("Hill estimator (Pareto a=", alpha_true, ", n=10000)"),
xlab = "k (top k order statistics used)", ylab = "alpha-hat",
ylim = c(0, 5))
abline(h = alpha_true, col = "tomato", lwd = 2, lty = 2)
legend("topright", c("Hill estimate", paste0("True alpha = ", alpha_true)),
col = c("steelblue", "tomato"), lwd = 2, lty = c(1, 2), bty = "n")
5.3.2 重尾的后果：矩的存在性
对 Pareto(\(\alpha\))，\(k\) 阶矩存在当且仅当 \(k < \alpha\)。\(\alpha \leq 1\) 时连均值都不存在。
alphas <- c(0.8, 1.5, 2.5, 4.0)
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
for (alpha in alphas) {
x_p <- (1 - runif(5000))^(-1/alpha)
cum_mean <- cumsum(x_p) / seq_along(x_p)
plot(cum_mean, type = "l", col = "steelblue", lwd = 1.5,
main = paste0("Pareto(a=", alpha, ") ",
if (alpha > 1) paste0("mean exists (=", round(alpha/(alpha-1), 2), ")")
else "mean = infinity"),
xlab = "Sample size", ylab = "Running mean")
if (alpha > 1) abline(h = alpha/(alpha-1), col = "tomato", lwd = 2, lty = 2)
}
6 综合实战：fitdistrplus 模型选择
把前面的理论付诸实践：给定数据，如何系统选择最佳分布？
决策逻辑：
- 数据特征 → 缩小候选范围
- 非负且右偏？→ Exp, Gamma, Weibull, LogNormal
- 对称或近似对称？→ Normal, t, Logistic
- 有界区间？→ Beta, Uniform
- 极端值异常频繁？→ Pareto, 重尾分布
- Cullen-Frey 图（
descdist）→ 用偏度²和峰度在分布地图上定位 - 竞争拟合（
fitdist）→ MLE 拟合多个候选 - 诊断图 + 统计量（
denscomp/qqcomp+gofstat）→ 选最佳 - 不确定性（
bootdist）→ 参数有多可靠？
工具链：descdist → fitdist →
gofstat → bootdist
6.1 Cullen-Frey 图：快速筛选候选分布
## summary statistics
## ------
## min: 0.1216946 max: 5.114437
## median: 1.258152
## mean: 1.42749
## estimated sd: 0.8154855
## estimated skewness: 0.9041436
## estimated kurtosis: 3.7055156.2 多分布竞争拟合
f_gamma <- fitdist(true_data, "gamma")
f_lnorm <- fitdist(true_data, "lnorm")
f_weibull <- fitdist(true_data, "weibull")
f_exp <- fitdist(true_data, "exp")
plot_labels <- c("Gamma", "LogNormal", "Weibull", "Exp")
par(mfrow = c(2, 2))
denscomp(list(f_gamma, f_lnorm, f_weibull, f_exp), legendtext = plot_labels)
qqcomp (list(f_gamma, f_lnorm, f_weibull, f_exp), legendtext = plot_labels)
cdfcomp (list(f_gamma, f_lnorm, f_weibull, f_exp), legendtext = plot_labels)
ppcomp (list(f_gamma, f_lnorm, f_weibull, f_exp), legendtext = plot_labels)
6.3 拟合优度统计量
## Goodness-of-fit statistics
## Gamma LogNormal Weibull Exp
## Kolmogorov-Smirnov statistic 0.03736550 0.0535318 0.04041843 0.2200214
## Cramer-von Mises statistic 0.09205037 0.3353942 0.16638674 7.5386079
## Anderson-Darling statistic 0.55176491 2.1390725 1.16773243 41.7520511
##
## Goodness-of-fit criteria
## Gamma LogNormal Weibull Exp
## Akaike's Information Criterion 1112.748 1142.514 1121.152 1357.918
## Bayesian Information Criterion 1121.177 1150.943 1129.582 1362.133
7 全景图：四种运算生成所有分布
par(mar = c(1, 1, 2, 1))
plot(NA, xlim = c(0, 14), ylim = c(0, 11), axes = FALSE, bty = "n",
xlab = "", ylab = "",
main = "Complete map: how distributions are generated")
# ── color scheme ──
col_src <- "gray30"
col_tf <- "#2E86AB"
col_sum <- "#A23B72"
col_rat <- "#F18F01"
col_ext <- "#C73E1D"
draw_box <- function(x, y, label, col, w = 2.3, h = 0.6) {
rect(x-w/2, y-h/2, x+w/2, y+h/2, col = col, border = "white", lwd = 2)
text(x, y, label, col = "white", font = 2, cex = 0.72)
}
draw_arr <- function(x0, y0, x1, y1, label = "", adj_x = 0, adj_y = 0.18,
col_line = "gray45") {
arrows(x0, y0, x1, y1, length = 0.10, lwd = 1.3, col = col_line)
if (nchar(label) > 0)
text((x0+x1)/2 + adj_x, (y0+y1)/2 + adj_y,
label, cex = 0.55, col = "gray25", font = 3)
}
# ══════════════════════════════════════════════════════════
# Layer 1: Source
# ══════════════════════════════════════════════════════════
draw_box(7, 10.3, "Uniform(0,1)", col_src, w = 2.8)
# ══════════════════════════════════════════════════════════
# Layer 2: Transform
# ══════════════════════════════════════════════════════════
text(0.8, 8.6, "TRANSFORM", col = col_tf, font = 2, cex = 0.7, srt = 90)
draw_box(3, 8.6, "Exp(lambda)", col_tf)
draw_box(7, 8.6, "Normal(mu, sigma)", col_tf, w = 2.8)
draw_box(11, 8.6, "Weibull(k, lambda)", col_tf, w = 2.6)
draw_arr(5.7, 10.05, 3, 8.95, "-log(1-U)/lam")
draw_arr(7, 10.0, 7, 8.95, "qnorm(U)")
draw_arr(8.3, 10.05, 11, 8.95, "(-log(1-U))^(1/k)")
# LogNormal from Normal
draw_box(11, 7, "LogNormal(mu, sigma)", col_tf, w = 2.8)
draw_arr(8.4, 8.4, 11, 7.35, "exp(X)", adj_x = 0.3)
# Weibull from Exp
draw_arr(4.2, 8.4, 9.7, 8.75, "X^(1/k)", adj_y = 0.22)
# ══════════════════════════════════════════════════════════
# Layer 3: Sum
# ══════════════════════════════════════════════════════════
text(0.8, 6.3, "SUM", col = col_sum, font = 2, cex = 0.7, srt = 90)
draw_box(3, 6.3, "Gamma(k, lambda)", col_sum, w = 2.6)
draw_box(7, 6.3, "Chi-squared(k)", col_sum, w = 2.4)
draw_box(11, 6.3, "Poisson(lambda*t)", col_sum, w = 2.6)
draw_arr(3, 8.3, 3, 6.65, "sum k copies")
draw_arr(7, 8.3, 7, 6.65, "sum k Z^2")
draw_arr(3, 8.3, 11, 6.65, "count in window", adj_x = 1)
# Gamma = Chi-sq
draw_arr(5.3, 6.3, 4.3, 6.3, "= Gamma(k/2, 1/2)", adj_y = -0.2)
# CLT arrow from Any -> Normal
draw_box(13, 8.6, "Any dist\n(finite var)", "gray50", w = 1.8, h = 0.7)
draw_arr(12.1, 8.6, 8.4, 8.6, "CLT: mean of n", adj_y = 0.2)
# ══════════════════════════════════════════════════════════
# Layer 4: Ratio
# ══════════════════════════════════════════════════════════
text(0.8, 4.3, "RATIO", col = col_rat, font = 2, cex = 0.7, srt = 90)
draw_box(5, 4.3, "t(k)", col_rat, w = 2)
draw_box(9, 4.3, "F(m, n)", col_rat, w = 2)
draw_arr(6, 6.0, 5, 4.65, "Z / sqrt(chi2/k)", adj_x = -1.5)
draw_arr(7, 6.0, 5.5, 4.65, "")
draw_arr(8, 6.0, 9, 4.65, "chi2_m/m / chi2_n/n")
# t -> Normal as k -> inf
text(3.5, 4.9, "k -> inf: t -> Normal", cex = 0.55, col = "gray40", font = 3)
# ══════════════════════════════════════════════════════════
# Layer 5: Extreme
# ══════════════════════════════════════════════════════════
text(0.8, 2.0, "EXTREME", col = col_ext, font = 2, cex = 0.7, srt = 90)
draw_box(2, 2, "Beta(a, b)", col_ext)
draw_box(5.5, 2, "Gumbel", col_ext, w = 1.8)
draw_box(8.5, 2, "Frechet", col_ext, w = 1.8)
draw_box(11.5,2, "Pareto(alpha)", col_ext, w = 2.2)
draw_arr(7, 10.0, 2, 2.35, "order stats X_(k)", adj_x = -1.2)
draw_arr(7, 8.3, 5.5, 2.35, "max n (light tail)", adj_x = -0.5)
draw_arr(11.5,2.35, 8.5, 2.35, "max n (heavy tail)", adj_y = 0.25)
draw_arr(7, 10.0, 11.5,2.35, "1/(1-U)^(1/a)", adj_x = 1.2)
# ── legend ──
legend(0.5, 3.8,
legend = c("Source", "Transform (Sec 2)",
"Sum (Sec 3)", "Ratio (Sec 4)", "Extreme (Sec 5)"),
fill = c(col_src, col_tf, col_sum, col_rat, col_ext),
bty = "n", cex = 0.75, border = "white")
小结：
| Operation | Input | Output |
|---|---|---|
| Inverse CDF \(F^{-1}(U)\) | Uniform | Any distribution |
| Power \(X^{1/k}\) | Exp | Weibull |
| Exponential \(e^X\) | Normal | LogNormal |
| Sum of \(k\) copies | Exp | Gamma |
| Sum of \(k\) copies | Any (CLT) | Normal |
| Sum of \(k\) squares | Normal | Chi-squared (= Gamma) |
| Normal / sqrt(Chi-sq) | Normal, Chi-sq | t |
| Chi-sq / Chi-sq | Chi-sq, Chi-sq | F |
| Order statistic \(X_{(k)}\) | Uniform | Beta |
| Max of \(n\) (light tail) | Normal, Exp | Gumbel |
| Max of \(n\) (heavy tail) | Pareto | Frechet |
| Max of \(n\) (bounded) | Uniform | Weibull (EVT) |
8 附录：Gamma 和 Beta 为什么无处不在
这两个分布不是独立的机理，而是前述运算的自然数学延伸。
par(mfrow = c(1, 2), mar = c(4, 4, 3.5, 1))
x_a <- seq(0, 12, length.out = 400)
cols_a <- c("steelblue","tomato","forestgreen","orange","purple")
# Gamma = sum of k Exp
plot(NA, xlim = c(0, 12), ylim = c(0, 0.55),
main = "Gamma(k, 1) = sum of k Exp(1)\n[k=1: Exp; k->inf: Normal by CLT]",
xlab = "x", ylab = "Density")
for (i in seq_along(c(1,2,5,10,20))) {
k <- c(1,2,5,10,20)[i]
lines(x_a, dgamma(x_a, shape = k, rate = 1), col = cols_a[i], lwd = 2)
}
legend("topright", paste("k =", c(1,2,5,10,20)), col = cols_a, lwd = 2, bty = "n")
# Beta = order statistic
x_b <- seq(0, 1, length.out = 400)
params <- list(c(1,1), c(2,9), c(5,6), c(9,2))
plot(NA, xlim = c(0, 1), ylim = c(0, 4),
main = "Beta(k, n-k+1) = k-th order stat of n Uniform\n[a=b=1: Uniform]",
xlab = "x", ylab = "Density")
for (i in seq_along(params))
lines(x_b, dbeta(x_b, params[[i]][1], params[[i]][2]), col = cols_a[i], lwd = 2)
legend("topright",
sapply(params, function(p) paste0("Beta(", p[1], ",", p[2], ")")),
col = cols_a[seq_along(params)], lwd = 2, bty = "n")
- Gamma：\(k\) 个 Exp 之和。\(k=1\) 是 Exp，\(k \to \infty\) 由 CLT 趋向 Normal。Chi-squared 是 Gamma(\(k/2\), \(1/2\)) 的特例。
- Beta：\(n\) 个 Uniform 的第 \(k\) 小次序统计量。\(\alpha=\beta=1\) 退化为 Uniform。
Generated: 2026-07-04 17:57:26.794819