10.3 箱型圖、Q-Q 圖與 ECDF#

import rstatistics as r

g = r.set_seed(11)
a = [round(v, 4) for v in g.rnorm(40, 100.0, 8.0)]
b = [round(v, 4) for v in g.rnorm(40, 108.0, 8.0)]
c = [round(v, 4) for v in g.rnorm(40, 104.0, 14.0)]

for name, v in [("A", a), ("B", b), ("C", c)]:
    s = r.summary(v)
    print(name, "mean =", round(r.mean(v), 3), "sd =", round(r.sd(v), 3),
          "median =", r.median(v), "IQR =", round(r.IQR(v), 3))

print()
for name, v in [("A", a), ("B", b), ("C", c)]:
    bs = r.boxplot_stats(v)
    print(f"{name}: stats = {[round(q,3) for q in bs.stats]}  conf = {[round(q,3) for q in bs.conf]}  out = {len(bs.out)}")

A mean = 97.443 sd = 6.197 median = 97.21555000000001 IQR = 7.215
B mean = 108.292 sd = 7.641 median = 109.26304999999999 IQR = 10.362
C mean = 103.865 sd = 14.17 median = 103.327 IQR = 16.692

A: stats = [87.38, 92.784, 97.216, 100.156, 110.589]  conf = [95.374, 99.057]  out = 1
B: stats = [90.506, 102.192, 109.263, 113.066, 126.718]  conf = [106.547, 111.979]  out = 0
C: stats = [71.815, 95.538, 103.327, 112.277, 135.662]  conf = [99.145, 107.509]  out = 0

import rstatistics as r

g = r.set_seed(11)
a = [round(v, 4) for v in g.rnorm(40, 100.0, 8.0)]
b = [round(v, 4) for v in g.rnorm(40, 108.0, 8.0)]
c = [round(v, 4) for v in g.rnorm(40, 104.0, 14.0)]

grp = [("A 早班", a), ("B 晚班", b), ("C 混合", c)]
stats = [r.boxplot_stats(v) for _, v in grp]

for (name, v), s in zip(grp, stats):
    print(f"{name}: 五數 = {[round(q,3) for q in s.stats]}  鬚線內 = [{s.stats[0]:.3f}, {s.stats[4]:.3f}]  離群 = {[round(o,3) for o in s.out]}")

lo = min(s.stats[0] for s in stats)
hi = max(s.stats[4] for s in stats)

gr = r.Graphics(width=8.5, height=5.5)
gr.plot_window([0.4, 3.6], [lo, hi])
gr.title("三組產量的並排箱型圖", "", "產品重量 (公克)")
xt = [1.0, 2.0, 3.0]
gr.axis(1, xt, ["A 早班", "B 晚班", "C 混合"])
yt = [v for v in r.pretty([lo, hi], 5) if lo <= v <= hi]
gr.axis(2, yt, [f"{v:g}" for v in yt])
gr.box()

for i, ((name, v), s) in enumerate(zip(grp, stats)):
    xc = float(i + 1)
    bw = 0.22
    # 鬚線
    gr.segments([xc], [s.stats[0]], [xc], [s.stats[1]], "black")
    gr.segments([xc], [s.stats[3]], [xc], [s.stats[4]], "black")
    gr.segments([xc - bw / 2], [s.stats[0]], [xc + bw / 2], [s.stats[0]], "black")
    gr.segments([xc - bw / 2], [s.stats[4]], [xc + bw / 2], [s.stats[4]], "black")
    # 箱子
    gr.rect(xc - bw / 2, s.stats[1], xc + bw / 2, s.stats[3], "black", "lightgray")
    # 中位數
    gr.segments([xc - bw / 2], [s.stats[2]], [xc + bw / 2], [s.stats[2]], "black", 2.0)
    # 離群點
    for o in s.out:
        gr.points([xc], [o], pch=1, col="black")
    # 中位數信賴區間
    gr.segments([xc + bw / 2 + 0.04], [s.conf[0]], [xc + bw / 2 + 0.04], [s.conf[1]], "red", 1.5)
gr.save("img/10.3-box-groups.svg")
print("已存 10.3-box-groups.svg")

A 早班: 五數 = [87.38, 92.784, 97.216, 100.156, 110.589]  鬚線內 = [87.380, 110.589]  離群 = [112.364]
B 晚班: 五數 = [90.506, 102.192, 109.263, 113.066, 126.718]  鬚線內 = [90.506, 126.718]  離群 = []
C 混合: 五數 = [71.815, 95.538, 103.327, 112.277, 135.662]  鬚線內 = [71.815, 135.662]  離群 = []
已存 10.3-box-groups.svg

set.seed(11)
a <- round(rnorm(40,100,8),4)
b <- round(rnorm(40,108,8),4)
c <- round(rnorm(40,104,14),4)
for (nm in c("A","B","C")) {
  v <- switch(nm, A=a, B=b, C=c)
  bs <- boxplot.stats(v)
  cat(sprintf("%s: stats = %s  conf = %s  out = %s\n", nm,
      paste(round(bs$stats,3), collapse=" "), paste(round(bs$conf,3), collapse=" "),
      paste(round(bs$out,3), collapse=" ")))
}
cat("\n=== summary ===\n")
for (nm in c("A","B","C")) {
  v <- switch(nm, A=a, B=b, C=c)
  cat(sprintf("%s mean = %.3f sd = %.3f median = %.3f IQR = %.3f\n", nm, mean(v), sd(v), median(v), IQR(v)))
}

A: stats = 87.38 92.784 97.216 100.156 110.589  conf = 95.374 99.057  out = 112.364
B: stats = 90.506 102.192 109.263 113.066 126.718  conf = 106.547 111.979  out = 
C: stats = 71.815 95.538 103.327 112.277 135.662  conf = 99.145 107.509  out = 

=== summary ===
A mean = 97.443 sd = 6.197 median = 97.216 IQR = 7.215
B mean = 108.292 sd = 7.641 median = 109.263 IQR = 10.362
C mean = 103.865 sd = 14.170 median = 103.327 IQR = 16.692

import rstatistics as r

print("=== ppoints ===")
print("ppoints(10) =", [round(v, 4) for v in r.ppoints(10)])
print("ppoints(5)  =", [round(v, 4) for v in r.ppoints(5)])
print("ppoints(60) 的前後三個 =", [round(v, 4) for v in r.ppoints(60)[:3]], [round(v, 4) for v in r.ppoints(60)[-3:]])

=== ppoints ===
ppoints(10) = [0.061, 0.1585, 0.2561, 0.3537, 0.4512, 0.5488, 0.6463, 0.7439, 0.8415, 0.939]
ppoints(5)  = [0.119, 0.3095, 0.5, 0.6905, 0.881]
ppoints(60) 的前後三個 = [0.0083, 0.025, 0.0417] [0.9583, 0.975, 0.9917]

import rstatistics as r

g = r.set_seed(11)
norm = [round(v, 4) for v in g.rnorm(60, 0.0, 1.0)]
heavy = [round(v, 4) for v in g.rt(60, df=3)]

for name, v in [("常態", norm), ("重尾 t(3)", heavy)]:
    t, s = r.qqnorm_coords(v)
    print(f"--- {name} ---")
    print("  理論分位前三:", [round(q, 4) for q in t[:3]], " 樣本分位前三:", [round(q, 4) for q in s[:3]])
    print("  理論分位後三:", [round(q, 4) for q in t[-3:]], " 樣本分位後三:", [round(q, 4) for q in s[-3:]])
    qx = r.quantile(t, [0.25, 0.75])
    qy = r.quantile(s, [0.25, 0.75])
    slope = (qy[1] - qy[0]) / (qx[1] - qx[0])
    icept = qy[0] - slope * qx[0]
    print(f"  qqline: slope = {slope:.4f}, intercept = {icept:.4f}")
    print(f"  樣本 sd = {r.sd(v):.4f}")
    print(f"  左端 (i=1): 理論 {t[0]:.3f} 樣本 {s[0]:.3f} 差距 {s[0]-t[0]:+.3f}")
    print(f"  右端 (i=n): 理論 {t[-1]:.3f} 樣本 {s[-1]:.3f} 差距 {s[-1]-t[-1]:+.3f}")
    hs = r.shapiro_test(v)
    print(f"  Shapiro: W = {hs.statistic:.4f}, p = {hs.p_value:.4g}")

--- 常態 ---
  理論分位前三: [-2.394, -1.96, -1.7317]  樣本分位前三: [-2.1868, -1.6365, -1.5775]
  理論分位後三: [1.7317, 1.96, 2.394]  樣本分位後三: [1.3236, 1.4879, 1.5456]
  qqline: slope = 0.9607, intercept = -0.1854
  樣本 sd = 0.8330
  左端 (i=1): 理論 -2.394 樣本 -2.187 差距 +0.207
  右端 (i=n): 理論 2.394 樣本 1.546 差距 -0.848
  Shapiro: W = 0.9895, p = 0.8884
--- 重尾 t(3) ---
  理論分位前三: [-2.394, -1.96, -1.7317]  樣本分位前三: [-2.455, -2.3613, -1.6047]
  理論分位後三: [1.7317, 1.96, 2.394]  樣本分位後三: [4.4352, 5.4379, 5.4829]
  qqline: slope = 1.1523, intercept = -0.0649
  樣本 sd = 1.5554
  左端 (i=1): 理論 -2.394 樣本 -2.455 差距 -0.061
  右端 (i=n): 理論 2.394 樣本 5.483 差距 +3.089
  Shapiro: W = 0.8543, p = 4.083e-06

import rstatistics as r

g = r.set_seed(11)
norm = [round(v, 4) for v in g.rnorm(60, 0.0, 1.0)]
heavy = [round(v, 4) for v in g.rt(60, df=3)]

gr = r.Graphics(width=10.0, height=5.0, mfrow=[1, 2])

# 面板 1：常態
gr.qqnorm(norm, main="常態資料的 Q-Q 圖", xlab="理論分位數", ylab="樣本分位數")
gr.qqline(norm, col="red")

# 面板 2：重尾
gr.next_panel()
gr.qqnorm(heavy, main="重尾資料 (t, df=3) 的 Q-Q 圖", xlab="理論分位數", ylab="樣本分位數")
gr.qqline(heavy, col="red")

gr.save("img/10.3-qqnorm.svg")
print("已存 10.3-qqnorm.svg")

已存 10.3-qqnorm.svg

set.seed(11)
norm <- round(rnorm(60,0,1),4)
heavy <- round(rt(60,df=3),4)

cat("=== ppoints ===\n")
cat("ppoints(10) =", round(ppoints(10),4), "\n")
cat("ppoints(5)  =", round(ppoints(5),4), "\n")

cat("\n=== qqnorm 的結構 ===\n")
q <- qqnorm(norm, plot.it=FALSE)
cat("names:", names(q), " length(q$x) =", length(q$x), " length(q$y) =", length(q$y), "\n")
cat("sort(q$x) 前三（= rstatistics 的理論分位）:", round(head(sort(q$x),3),4), "\n")
cat("sort(q$y) 前三:", round(head(sort(q$y),3),4), "\n")
cat("sorted y 後三:", round(tail(sort(q$y),3),4), "\n")
cat("heavy sorted y 後三:", round(tail(sort(heavy),4),4), "\n")

cat("\n=== qqline ===\n")
qq_manual <- function(v) {
  n <- length(v); p <- ppoints(n)
  list(x = qnorm(p), y = sort(v))
}
qqline_manual <- function(v) {
  q <- qq_manual(v)
  qx <- quantile(q$x, c(0.25,0.75)); qy <- quantile(q$y, c(0.25,0.75))
  slope <- (qy[2]-qy[1])/(qx[2]-qx[1])
  list(slope = slope, intercept = qy[1] - slope*qx[1])
}
cat("qqline(norm): slope =", qqline_manual(norm)$slope, " intercept =", qqline_manual(norm)$intercept, "\n")
cat("qqline(heavy): slope =", qqline_manual(heavy)$slope, " intercept =", qqline_manual(heavy)$intercept, "\n")

cat("\n=== sd & shapiro ===\n")
cat("sd(norm) =", sd(norm), "\n")
cat("Shapiro(norm): W =", shapiro.test(norm)$statistic, " p =", shapiro.test(norm)$p.value, "\n")
cat("Shapiro(heavy): W =", shapiro.test(heavy)$statistic, " p =", shapiro.test(heavy)$p.value, "\n")

=== ppoints ===
ppoints(10) = 0.061 0.1585 0.2561 0.3537 0.4512 0.5488 0.6463 0.7439 0.8415 0.939
ppoints(5)  = 0.119 0.3095 0.5 0.6905 0.881

=== qqnorm 的結構 ===
names: x y  length(q$x) = 60  length(q$y) = 60
sort(q$x) 前三（= rstatistics 的理論分位）: -2.394 -1.96 -1.7317
sort(q$y) 前三: -2.1868 -1.6365 -1.5775
sorted y 後三: 1.3236 1.4879 1.5456
heavy sorted y 後三: 3.0001 4.4352 5.4379 5.4829

=== qqline ===
qqline(norm): slope = 0.960729  intercept = -0.18535
qqline(heavy): slope = 1.152289  intercept = -0.0648625

=== sd & shapiro ===
sd(norm) = 0.8330054
Shapiro(norm): W = 0.9895462  p = 0.8884449
Shapiro(heavy): W = 0.8542791  p = 4.083211e-06

import rstatistics as r

g = r.set_seed(11)
norm = [round(v, 4) for v in g.rnorm(60, 0.0, 1.0)]

gr = r.Graphics(width=6.0, height=5.0)
e = gr.ecdf(norm, main="g.ecdf 畫出的經驗分布函數", xlab="x", ylab="F(x)")
print("e =", e)
print("e.n =", e.n)
print("e.x[:5] =", [round(v, 4) for v in e.x[:5]])
print("e.y[:5] =", [round(v, 4) for v in e.y[:5]])
print("e.x[-3:] =", [round(v, 4) for v in e.x[-3:]])
print("e.y[-3:] =", [round(v, 4) for v in e.y[-3:]])

e = ECDF of 60 observations
e.n = 60
e.x[:5] = [-2.1868, -1.6365, -1.5775, -1.5383, -1.5166]
e.y[:5] = [0.0167, 0.0333, 0.05, 0.0667, 0.0833]
e.x[-3:] = [1.3236, 1.4879, 1.5456]
e.y[-3:] = [0.9667, 0.9833, 1.0]

import rstatistics as r

g = r.set_seed(11)
norm = [round(v, 4) for v in g.rnorm(60, 0.0, 1.0)]
heavy = [round(v, 4) for v in g.rt(60, df=3)]

def ecdf_poly(v):
    xs = sorted(v)
    n = len(v)
    px = [xs[0]]
    py = [0.0]
    for i in range(n):
        px.append(xs[i]); py.append((i + 1) / n)
        if i + 1 < n:
            px.append(xs[i + 1]); py.append((i + 1) / n)
    return px, py

gr = r.Graphics(width=10.0, height=5.0, mfrow=[1, 2])

for k, (name, v) in enumerate([("常態資料", norm), ("重尾資料 (t, df=3)", heavy)]):
    if k > 0:
        gr.next_panel()
    px, py = ecdf_poly(v)
    lo, hi = min(v), max(v)
    gr.plot_window([lo, hi], [0.0, 1.0])
    gr.title(f"{name} 的 ECDF 與理論常態 CDF", "x", "F(x)")
    xt = [t for t in r.pretty([lo, hi], 5) if lo <= t <= hi]
    gr.axis(1, xt, [f"{t:g}" for t in xt])
    yt = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
    gr.axis(2, yt, [f"{t:g}" for t in yt])
    gr.box()
    # 理論的標準常態 CDF
    ts = [lo + (hi - lo) * i / 200.0 for i in range(201)]
    gr.lines(ts, [r.pnorm(t, 0.0, 1.0) for t in ts], "blue", 1.5)
    # 經驗分布函數的階梯
    gr.lines(px, py, "black", 1.2)
    # KS 統計量：最大的垂直差距
    best = 0.0
    bx = 0.0
    for i, s in enumerate(sorted(v)):
        F = r.pnorm(s, 0.0, 1.0)
        d1 = abs((i + 1) / len(v) - F)
        d2 = abs(i / len(v) - F)
        if d1 > best:
            best = d1; bx = s
        if d2 > best:
            best = d2; bx = s
    Fy = sum(1 for u in v if u <= bx) / len(v)
    Ft = r.pnorm(bx, 0.0, 1.0)
    gr.segments([bx], [min(Fy, Ft)], [bx], [max(Fy, Ft)], "red", 1.5)
    gr.text(bx - 1.9, max(Fy, Ft) + 0.04, f"最大差距 D = {best:.3f}", col="red")
gr.save("img/10.3-ecdf.svg")
print("已存 10.3-ecdf.svg")

已存 10.3-ecdf.svg

set.seed(11)
norm <- round(rnorm(60,0,1),4)
heavy <- round(rt(60,df=3),4)

cat("=== ecdf ===\n")
e <- ecdf(norm)
ee <- environment(e)
xs <- get("x", ee); ys <- get("y", ee)
cat("class:", class(e), " length =", length(xs), "\n")
print(round(head(xs,5),4)); print(round(head(ys,5),4)); print(round(tail(ys,3),4))

cat("\n=== KS（單樣本，對理論分布）===\n")
print(ks.test(norm,"pnorm"))
print(ks.test(heavy,"pnorm"))

=== ecdf ===
class: ecdf stepfun function  length = 60
[1] -2.1868 -1.6365 -1.5775 -1.5383 -1.5166
[1] 0.0167 0.0333 0.0500 0.0667 0.0833
[1] 0.9667 0.9833 1.0000

=== KS（單樣本，對理論分布）===
    Exact one-sample Kolmogorov-Smirnov test

data:  norm
D = 0.15416, p-value = 0.1037
alternative hypothesis: two-sided

    Exact one-sample Kolmogorov-Smirnov test

data:  heavy
D = 0.13856, p-value = 0.1817
alternative hypothesis: two-sided