之前我们介绍过，自由度是一种计量独立信息个数的代币。但这种代币不能告诉我们手里的数据包含的的信息量具体有多少。「信息量」是一个很有意思的词，我们该怎么量化它？

为了解答这个问题，我们得理解「我们如何知晓一个信息」。比如，你和你的搭档被分别安置在两个屋子里。你面前有五个格子，你要通过一个老式电话询问你的搭档，这个词的拼写是什么。你被要求不能直接问「这个词」是怎么拼的，只能通过见解揣测的方式，对方只能回答 Yes or No。

如果你是一个算法大师，你可能会用二分法的方法来询问：第一个字母在 N 之前还是 N 之后？如果对方告诉你在 N 之前，你可以继续问，它在 G 之前还是 G 之后？接着，你会问它在 D 之前还是 D 之后？最后你会问它在 B 之前还是 B 之后？问过这几个问题之后，你知道了，第一个字母是 A。

平均而言，为了确认一个字母，你需要问 4.7 个问题，而对于五个字母，你需要问 23.5 个问题。

不过，如果你是一名语言学大师，你可能会发现这么问问题很蠢。因为人们并不是这样理解语言信息的。让我们来看一个例子：「X CXN STXLL RXXD THXS ​XVXN THXXH XLL THX VXWXLS XRX MXSSXNG​」。我猜你大概还是知道这句话的意思为何：「I CAN STILL READ THIS ​EVEN THOUGH ​ALL THE VOWELS ARE MISSING​」。换句话说，你也可以这样问：「它是吃的吗？」、「它是水果吗？」、「它是红色的吗？」，类似这样有实际语义的问题可以极大的提升你的提问效率。

对于英文这个语言，平均对于每个字母你需要问 2.62 个问题。而对于全人类的所有语言，我们通过口语每秒钟传达的信息大概需要问 39 个 Yes 和 No 才可以描绘清楚。

这是一种描绘信息量的方式，问 Yes or No，而这种信息量被称作香农信息量。

对于不同的场景、不同的问题，我们可能需要不同的信息量衡量方法。比如作为一名分析师，你可能想要知道你根据手里的数据分布，到底能够以多大的信心推算出总体信息的分布。如果我们手里的数据包含的信息量足够多，那么我们自然可以有效率地找到总体的分布，如果我们没有足够多的信息，那这件事情就会变得很麻烦。

什么是足够多的信息？样本量大肯定是一个优势。但如果你的数据本身就包含很多噪音，那么样本量大的优势就会被削弱一些，因为噪音意味着不确定性，意味着寻找总体的效率会被折损。再比如，如果你对数据做了一些奇妙处理，那么信息也有可能会有一些折损。像是你测量了全校学生的身高数据，接着把它们分成了「高个子」和「矮个子」。这样我们损失了数据的细节，引入了更大的统计误差，那么找到总体的难度就会陡升^[1]。

如何衡量「容易找到」？让我们来做一个统计学实验：生成一组随机的正态分布样本，然后按照不同的均值、方差生成假定的总体。我们要看这个样本来自不同总体的可能性有多大高，或者它和不同总体的相容程度有多高。

先来定义一些一会将会用到的分析工具，如果你的目的只是为了来看科普，那么可以不用细读其实现细节。

# SHARE
library(MASS)
library(matrixStats)
library(ggplot2)

set.seed(114514)

# ABC 后验：单组正态，汇总统计量为 (样本均值, 样本标准差)
# 返回 data.frame(mu, sig)，每行为一组被接受的参数
abc_normal <- function(
    data,
    mu_range      = c(mean(data) - 3, mean(data) + 3),
    sig_range     = c(0.1, 8),
    n_sim         = 200000,
    n_pilot       = 10000,
    target_accept = 4000
) {
  n <- length(data)
  obs_mean <- mean(data); obs_sd <- sd(data)

  mu_pr  <- runif(n_sim,   mu_range[1],  mu_range[2])
  sig_pr <- runif(n_sim,   sig_range[1], sig_range[2])

  mu_p  <- runif(n_pilot,  mu_range[1],  mu_range[2])
  sig_p <- runif(n_pilot,  sig_range[1], sig_range[2])
  Yp    <- matrix(rnorm(n * n_pilot, rep(mu_p, each=n), rep(sig_p, each=n)), nrow=n)
  s_mean <- sd(colMeans(Yp)); s_sd <- sd(apply(Yp, 2, sd))

  Y    <- matrix(rnorm(n * n_sim, rep(mu_pr, each=n), rep(sig_pr, each=n)), nrow=n)
  dist <- sqrt(((colMeans(Y) - obs_mean) / s_mean)^2 +
               ((colSds(Y)   - obs_sd  ) / s_sd  )^2)

  idx <- order(dist)[1:target_accept]
  data.frame(mu = mu_pr[idx], sig = sig_pr[idx])
}

# ABC 后验：单组正态，汇总统计量为「超过阈值的计数」
# 返回 data.frame(mu, sig)
abc_normal_binary <- function(
    data,
    cut_val,
    mu_range      = c(mean(data) - 3, mean(data) + 3),
    sig_range     = c(0.1, 8),
    n_sim         = 200000,
    n_pilot       = 10000,
    target_accept = 4000
) {
  n <- length(data); k_obs <- sum(data > cut_val)

  mu_pr  <- runif(n_sim,   mu_range[1],  mu_range[2])
  sig_pr <- runif(n_sim,   sig_range[1], sig_range[2])

  mu_p  <- runif(n_pilot,  mu_range[1],  mu_range[2])
  sig_p <- runif(n_pilot,  sig_range[1], sig_range[2])
  Yp    <- matrix(rnorm(n * n_pilot, rep(mu_p, each=n), rep(sig_p, each=n)), nrow=n)
  s_k   <- sd(colSums(Yp > cut_val))

  Y    <- matrix(rnorm(n * n_sim, rep(mu_pr, each=n), rep(sig_pr, each=n)), nrow=n)
  dist <- abs(colSums(Y > cut_val) - k_obs) / s_k

  idx <- order(dist)[1:target_accept]
  data.frame(mu = mu_pr[idx], sig = sig_pr[idx])
}

# ABC 后验：两组正态（共享标准差），汇总统计量为 (均值1, 均值2, 合并标准差)
# 返回 data.frame(mu1, mu2)
abc_two_group <- function(
    x1, x2,
    mu_range      = c(-2, 2.5),
    sig_range     = c(0.5, 2.5),
    n_sim         = 200000,
    n_pilot       = 10000,
    target_accept = 4000
) {
  n1 <- length(x1); n2 <- length(x2)
  obs_m1 <- mean(x1); obs_m2 <- mean(x2)
  obs_sp <- sqrt(((n1-1)*var(x1) + (n2-1)*var(x2)) / (n1+n2-2))

  mu1_pr <- runif(n_sim, mu_range[1], mu_range[2])
  mu2_pr <- runif(n_sim, mu_range[1], mu_range[2])
  sig_pr <- runif(n_sim, sig_range[1], sig_range[2])

  mu1_p <- runif(n_pilot, mu_range[1], mu_range[2])
  mu2_p <- runif(n_pilot, mu_range[1], mu_range[2])
  sig_p <- runif(n_pilot, sig_range[1], sig_range[2])
  Y1p <- matrix(rnorm(n1*n_pilot, rep(mu1_p, each=n1), rep(sig_p, each=n1)), nrow=n1)
  Y2p <- matrix(rnorm(n2*n_pilot, rep(mu2_p, each=n2), rep(sig_p, each=n2)), nrow=n2)
  s_mean <- sd(colMeans(Y1p))
  s_sp   <- sd(sqrt(((n1-1)*colVars(Y1p) + (n2-1)*colVars(Y2p)) / (n1+n2-2)))

  Y1 <- matrix(rnorm(n1*n_sim, rep(mu1_pr, each=n1), rep(sig_pr, each=n1)), nrow=n1)
  Y2 <- matrix(rnorm(n2*n_sim, rep(mu2_pr, each=n2), rep(sig_pr, each=n2)), nrow=n2)
  sim_sp <- sqrt(((n1-1)*colVars(Y1) + (n2-1)*colVars(Y2)) / (n1+n2-2))
  dist <- sqrt(((colMeans(Y1) - obs_m1) / s_mean)^2 +
               ((colMeans(Y2) - obs_m2) / s_mean)^2 +
               ((sim_sp        - obs_sp) / s_sp  )^2)

  idx <- order(dist)[1:target_accept]
  data.frame(mu1 = mu1_pr[idx], mu2 = mu2_pr[idx])
}

# ABC 后验：两组正态（共享标准差），汇总统计量为「各组超过阈值的计数」
# 返回 data.frame(mu1, mu2)
abc_two_group_binary <- function(
    x1, x2,
    cut_val,
    mu_range      = c(-2, 2.5),
    sig_range     = c(0.5, 2.5),
    n_sim         = 200000,
    n_pilot       = 10000,
    target_accept = 4000
) {
  n1 <- length(x1); n2 <- length(x2)
  k1_obs <- sum(x1 > cut_val); k2_obs <- sum(x2 > cut_val)

  mu1_pr <- runif(n_sim, mu_range[1], mu_range[2])
  mu2_pr <- runif(n_sim, mu_range[1], mu_range[2])
  sig_pr <- runif(n_sim, sig_range[1], sig_range[2])

  mu1_p <- runif(n_pilot, mu_range[1], mu_range[2])
  mu2_p <- runif(n_pilot, mu_range[1], mu_range[2])
  sig_p <- runif(n_pilot, sig_range[1], sig_range[2])
  Y1p <- matrix(rnorm(n1*n_pilot, rep(mu1_p, each=n1), rep(sig_p, each=n1)), nrow=n1)
  Y2p <- matrix(rnorm(n2*n_pilot, rep(mu2_p, each=n2), rep(sig_p, each=n2)), nrow=n2)
  s_k1 <- sd(colSums(Y1p > cut_val)); s_k2 <- sd(colSums(Y2p > cut_val))

  Y1 <- matrix(rnorm(n1*n_sim, rep(mu1_pr, each=n1), rep(sig_pr, each=n1)), nrow=n1)
  Y2 <- matrix(rnorm(n2*n_sim, rep(mu2_pr, each=n2), rep(sig_pr, each=n2)), nrow=n2)
  dist <- sqrt(((colSums(Y1 > cut_val) - k1_obs) / s_k1)^2 +
               ((colSums(Y2 > cut_val) - k2_obs) / s_k2)^2)

  idx <- order(dist)[1:target_accept]
  data.frame(mu1 = mu1_pr[idx], mu2 = mu2_pr[idx])
}

# 二维热图：越亮代表后验密度越高，绿色十字为真实参数位置
plot_heatmap <- function(post, true_x, true_y, xlim, ylim,
                         xlab = names(post)[1], ylab = names(post)[2]) {
  x_col <- names(post)[1]; y_col <- names(post)[2]
  ggplot(post, aes(.data[[x_col]], .data[[y_col]])) +
    stat_density_2d(
      geom    = "raster",
      aes(fill = pmax(log(after_stat(ndensity)), -10)),
      contour = FALSE, n = 150,
      h = c(bandwidth.nrd(post[[x_col]]), bandwidth.nrd(post[[y_col]]))
    ) +
    annotate("point", x = true_x, y = true_y,
             shape = 3, color = "#22c55e", size = 4, stroke = 2) +
    scale_fill_gradientn(colors = hcl.colors(256, "inferno"),
                         limits = c(-10, 0), guide = "none") +
    scale_x_continuous(expand = c(0, 0), limits = xlim) +
    scale_y_continuous(expand = c(0, 0), limits = ylim) +
    labs(x = xlab, y = ylab) +
    theme_minimal() + theme(panel.grid = element_blank())
}

# 叠加等高线图：将多个后验分布画在同一坐标系下，用于比较不确定性大小
# posts 为具名 list，每个元素为结构相同的 data.frame
plot_overlay <- function(posts, true_x, true_y, xlim, ylim,
                         xlab = NULL, ylab = NULL, colors = NULL) {
  x_col <- names(posts[[1]])[1]; y_col <- names(posts[[1]])[2]
  if (is.null(xlab)) xlab <- x_col
  if (is.null(ylab)) ylab <- y_col
  if (is.null(colors))
    colors <- setNames(hcl.colors(length(posts), "Dark 3"), names(posts))

  combined <- do.call(rbind, Map(function(df, nm) { df$label <- nm; df },
                                 posts, names(posts)))

  ggplot(combined, aes(.data[[x_col]], .data[[y_col]], color = label)) +
    stat_density_2d(linewidth = 0.75) +
    annotate("point", x = true_x, y = true_y,
             shape = 3, color = "#22c55e", size = 4, stroke = 2) +
    scale_color_manual(values = colors) +
    scale_x_continuous(expand = c(0, 0), limits = xlim) +
    scale_y_continuous(expand = c(0, 0), limits = ylim) +
    labs(x = xlab, y = ylab, color = NULL) +
    theme_minimal() +
    theme(panel.grid = element_blank(),
          legend.position = "inside",
          legend.position.inside = c(0.82, 0.1))
}

post_area <- function(post) {
  sqrt(det(cov(post[, 1:2])))
}

to_prob <- function(rl) {
  d <- exp(rl); d / sum(d)
}

整个实验设计是这样的，我们用蒙特卡洛投点的方式，让机器随便瞎猜可能的总体均值、方差为何。

在接下来的例子中，我们电脑随机生成 20 万个总体分布，涵盖了不同的均值和标准差组合。对于每一个组合，电脑都模拟出一批假数据，然后把假数据和我们手里的真数据比一比，看看它们到底有多相似。最像的那四千个组合会被我们当成合理猜测。

理解了原理后，我们看连续变量的状况如何：

data <- rnorm(100, mean = 0, sd = 1)
post <- abc_normal(data)

plot_heatmap(post, true_x = 0, true_y = 1,
             xlim = c(-4, 4), ylim = c(0.1, 8),
             xlab = expression(mu), ylab = expression(sigma))

cat("我们作为「上帝」，从均值为 0，标准差为 1 的总体中抽出了 100 个数字作为样本。\n")
cat(sprintf("这 100 个数字的样本均值是 %.3f，样本标准差是 %.3f。\n", mean(data), sd(data)))
cat("可以看到，样本统计量和真实值很接近，却不是一模一样（当然）。\n\n")
cat("让我们切换成「凡人」的视角，我们不知道这个样本对应的总体是什么。但我们有伟大的电脑，我们可以作出猜测。为了达到目的，疯狂的科学家试了 20 万组可能的总体均值和标准差组合，每组都生成一批同样大小的模拟数据，然后比一比模拟数据和真实样本的均值、标准差有多接近。越亮代表越接近，越暗代表越不接近。\n")
cat("我们从中留下了最接近的 4000 组参数。\n")
cat("这相当于问：哪些总体最有可能产生我们看到的这批样本？\n\n")
cat("留下来的参数分布（图中的亮色区域）告诉我们：\n")
cat(sprintf("均值的中心估计约 %.3f，\n", median(post$mu)))
cat(sprintf("标准差约 %.3f。\n", median(post$sig)))

我们每生成一个假定总体，就是一次模拟。对于每次模拟，我们都比对总体和样本的参数（均值和标准差）的匹配程度。的只有当一组参数生成的模拟样本，其均值和标准差同时与真实数据接近时，这组参数才会被保留。

由于必须同时满足两个层面的比对要求，符合条件的参数组合数量很少，因此图表呈现出一个小面积、明亮的圆形。

那么，如果这个样本的方差变大，会发生什么？让我们来试试看。

data <- rnorm(100, mean = 0, sd = 5)
post <- abc_normal(data, sig_range = c(0.1, 20))

plot_heatmap(post, true_x = 0, true_y = 5,
             xlim = c(-4, 4), ylim = c(0.1, 8),
             xlab = expression(mu), ylab = expression(sigma))

cat("这「上帝视角」的我们让总体的标准差变成 5，是之前的 5 倍。\n")
cat(sprintf("100 个样本的均值是 %.3f，标准差是 %.3f。\n", mean(data), sd(data)))
cat("从图中可以看出，样本均值理所当然地没有正好落在 0，但圆变大了好多！\n")
cat("方差变大，推测总体时，我们面临的不确定性就会变高。\n\n")
cat("这次，留下来的参数分布（图中的亮色区域）告诉我们：\n")
cat(sprintf("均值的中心估计约 %.3f，\n", median(post$mu)))
cat(sprintf("标准差约 %.3f。\n", median(post$sig)))

很不意外的，圆圈大了不少。

总结而言：如果你能看到又小、又聚集、又亮的点，那么就说明你可以非常有信心地估出总体的参数，信息量很高。而如果图形连成一大片，显得很分散，那么我们的心里可能就没那么有谱。你能看到这张图有两个维度，一个维度是均值，一个维度是方差。从图中你可以直观地看出，无论是推测均值，还是推测方差，不确定性都多了非常多。

样本容量的多寡虽然是一个评估信息总量的有效指标。但根据样本的品质，信息量也可能会出现多少的差异。想要看得更加全面，我们就得用更加明确的信息量大小衡量工具：费舍信息量。

到目前为止，我们都在用连续的均值和标准差来比对。这是最自然的方式，因为我们保留了数据的全部细节。接下来，我们换一个问题：如果我们只告诉机器「有多少个数字大于 0」，少了「均值和标准差是多少」这样的细致信息，机器能猜得同样准吗？

二值化在学界是一种广泛被使用的恶习。如果只模糊地说一种「感觉」，你大概能猜到，把人分成「快乐」、「不快乐」，「血压高」、「血压低」，「高个子」、「矮个子」几类会让数据的细节折损，但今天我们来看一点新鲜东西：直观地感受到信息的折损有多大。下面让我们沿着前例，画出一副「猜总体准确性」的图。

data <- rnorm(100, mean = 0, sd = 1)
post <- abc_normal_binary(data, cut_val = 0)

plot_heatmap(post, true_x = 0, true_y = 1,
             xlim = c(-3, 3), ylim = c(0.1, 8),
             xlab = expression(mu), ylab = expression(sigma))

cat("真实总体依然是均值 0，标准差 1。但这次我们换了一种「观察」方式：\n")
cat(sprintf("不记录具体数值，只看 100 个样本中有 %d 个大于 0。\n", sum(data > 0)))
cat("这种操作相当于把人分成「高个子」和「矮个子」，只留了符号，丢了大小。\n\n")
cat("可以看到，这个操作彻底改变了输出图形的形状，画面变成了喷射状！\n")
cat("仅凭「有几个大于 0」这一条信息，我们很难同时确定均值和标准差，\n")
cat("因为无数种 (μ, σ) 的组合都可能产生差不多的符号比例。\n\n")
cat("即便如此，留下来的参数分布还是给出了一些线索：\n")
cat(sprintf("均值的中心估计约 %.3f，\n", median(post$mu)))
cat(sprintf("标准差约 %.3f。\n", median(post$sig)))

在使用连续数据寻找总体时，样本的均值直接告诉我们总体的均值大概在哪里，样本的标准差直接告诉我们总体的标准差大概在哪里。这两个数字分别对应总体的两个不同方面，合在一起，我们对总体的位置和散布程度都有了具体的了解，图形因此聚成一个小团。

而二值化之后，我们只能知道「有多少个样本超过了 0」。这个数字同时取决于均值和标准差，但混合在一起之后数据细节全然失去，我们失去了估计均值和标准差的能力。均值高一些的总体，超过 0 的样本会更多；标准差大一些的总体，同样会有更多样本跑到 0 以上。因此，均值高一些同时标准差小一些的总体，和均值低一些同时标准差大一些的总体，完全可以产生几乎相同的超阈值个数。这个数字无法区分这两种总体，于是沿着「均值高则标准差小、均值低则标准差大」的方向，一整条带状区域里的总体都同样说得通，图形就沿着这条带扩散成了喷射状^[2]。

很多研究者会有这样的想法：「信息量有折损又没差，把连续变量变二分之后做卡方，的确是能够变得更显著，能算显著不就行了，管那么多干嘛」。为了看看这个问题造成的影响具体有多大，让我们来做下一个实验。

让我们来处理一个更为实用的问题：我们想要探究两个国家人口的身高有没有差异。为此，研究者分别从香蕉国和仙贝国抽取了一批人，测量了他们的身高。

我们随机生成 20 万种对两个国家均值的猜测。对于每一种猜测，电脑按照猜测的均值，分别为两个国家各生成一批符合正态分布的模拟数据，算出假数据里两个国家各自的样本均值，看它们和实际测量数据的两个样本均值差了多远。最接近的四千种猜测被留下来，画在以两个国家均值为坐标轴的图上。圆圈越小，说明我们对两个国家的均值越有把握。

二值化版本里，我们只记录了两个国家各自有多少人的身高超过了 170cm，其余信息全部丢弃。电脑同样生成 20 万种猜测，对每种猜测生成两批假数据，但这次只数假数据里两个国家各自超过 170cm 的人数，看它们和真实数据里的人数差了多远。最接近的四千种猜测同样画在同一张图上，圆圈大小可以直接和连续版本比^[3]。

x1 <- rnorm(60, mean = 0.0, sd = 1)
x2 <- rnorm(60, mean = 0.6, sd = 1)
post <- abc_two_group(x1, x2)

plot_heatmap(post, true_x = 0.0, true_y = 0.6,
             xlim = c(-2, 2.5), ylim = c(-2, 2.5),
             xlab = expression(mu[1]), ylab = expression(mu[2]))

cat(sprintf("真实总体：组1均值 0.0，组2均值 0.6，共同标准差 1.0。\n"))
cat(sprintf("观测均值：组1 %.3f，组2 %.3f，合并标准差 %.3f。\n\n",
            mean(x1), mean(x2),
            sqrt((59*var(x1) + 59*var(x2)) / 118)))
cat("图中横轴是组1的均值，纵轴是组2的均值。\n")
cat("亮色区域代表与观测数据最相容的参数组合。\n\n")
cat(sprintf("组1均值中心估计：%.3f\n", median(post$mu1)))
cat(sprintf("组2均值中心估计：%.3f\n", median(post$mu2)))

看起来不错，圆圆的一团。但如果是卡方风格的检验会给出什么样的结果呢？让我们再来跑一个模拟：

x1 <- rnorm(60, mean = 0.0, sd = 1)
x2 <- rnorm(60, mean = 0.6, sd = 1)
post <- abc_two_group_binary(x1, x2, cut_val = 0.3)

plot_heatmap(post, true_x = 0.0, true_y = 0.6,
             xlim = c(-2, 2.5), ylim = c(-2, 2.5),
             xlab = expression(mu[1]), ylab = expression(mu[2]))

cat(sprintf("真实总体：组1均值 0.0，组2均值 0.6，共同标准差 1.0。\n"))
cat(sprintf("以 0.3 为切割点，组1有 %d 人超过切割点，组2有 %d 人超过切割点。\n\n",
            sum(x1 > 0.3), sum(x2 > 0.3)))
cat("仅凭两组各自超过切割点的人数，后验分布弥散了很多。\n\n")
cat(sprintf("组1均值中心估计：%.3f\n", median(post$mu1)))
cat(sprintf("组2均值中心估计：%.3f\n", median(post$mu2)))

圆圈很明显大了一整圈，让我们把两张图叠加在一起，看看到底差了多少：

x1 <- rnorm(60, mean = 0.0, sd = 1)
x2 <- rnorm(60, mean = 0.6, sd = 1)

posts <- list(
  "t 风格"   = abc_two_group(x1, x2),
  "卡方风格" = abc_two_group_binary(x1, x2, cut_val = 0.3)
)

plot_overlay(
  posts,
  true_x = 0.0, true_y = 0.6,
  xlim   = c(-2, 2.5), ylim = c(-2, 2.5),
  xlab   = expression(mu[1]), ylab = expression(mu[2]),
  colors = c("t 风格" = "#f5a623", "卡方风格" = "#e040fb")
)

post_t   <- posts[["t 风格"]]
post_chi <- posts[["卡方风格"]]

# 后验联合扩散面积（协方差行列式）
vol_t   <- det(cov(post_t))
vol_chi <- det(cov(post_chi))

# 单维不确定性
sd_t_mu1   <- sd(post_t[, "mu1"])
sd_t_mu2   <- sd(post_t[, "mu2"])

sd_chi_mu1 <- sd(post_chi[, "mu1"])
sd_chi_mu2 <- sd(post_chi[, "mu2"])

# 以 t 风格为基准，计算二值化后的信息折损
vol_t   <- det(cov(post_t[, c("mu1", "mu2")]))
vol_chi <- det(cov(post_chi[, c("mu1", "mu2")]))

info_retained <- vol_t / vol_chi
info_lost <- 1 - info_retained
extra_sample_factor <- vol_chi / vol_t
extra_sample_pct <- (extra_sample_factor - 1) * 100

sd_t_mu1   <- sd(post_t$mu1)
sd_t_mu2   <- sd(post_t$mu2)
sd_chi_mu1 <- sd(post_chi$mu1)
sd_chi_mu2 <- sd(post_chi$mu2)

cor_t   <- cor(post_t$mu1, post_t$mu2)
cor_chi <- cor(post_chi$mu1, post_chi$mu2)

cat("轮廓越扩散，表示总体参数越难辨认。\n\n")

cat(sprintf("t 风格后验扩散面积：%.4f\n", vol_t))
cat(sprintf("卡方风格后验扩散面积：%.4f\n", vol_chi))
cat(sprintf("以 t 风格的后验扩散面积为基准（1.00×），卡方风格相当于把不确定性放大到 %.2f 倍。\n", extra_sample_factor))
cat(sprintf("  μ1 不确定性增加 %.0f%%\n",
            (sd_chi_mu1 / sd_t_mu1 - 1) * 100))
cat(sprintf("  μ2 不确定性增加 %.0f%%\n\n",
            (sd_chi_mu2 / sd_t_mu2 - 1) * 100))
cat(sprintf("对应的信息折损约为 %.1f%%。\n", info_lost * 100))
cat(sprintf(
  "若把连续分析当作基准，二值化后大约需要多 %.1f%% 的样本，才能接近同样的辨识能力。\n\n",
  extra_sample_pct
))

你可以从分析日志中看到巨大的信息折损。理由和前文讲得也一样，超过阈值的比例同时受均值和标准差影响，均值高一些、标准差大一些，可以产生和均值低一些、标准差小一些差不多的比例。仅凭计数这一个数字，机器没有办法把这两种情形分开，于是，大量不同的参数组合都同样合理^[4]。

事实上，这还不是最糟糕的，如果我们把切断点设置在偏离两组均值中心的位置，损失的信息量就会变得更多，下面是一个非常触目惊心的例子：

x1 <- rnorm(60, mean = 0.0, sd = 1)
x2 <- rnorm(60, mean = 0.6, sd = 1)

cuts_cmp <- c(-0.6, 0.3, 0.9)

posts_cmp <- lapply(cuts_cmp, function(cv) {
  abc_two_group_binary(
    x1, x2,
    cut_val = cv
  )
})

names(posts_cmp) <- sprintf("cut=%.1f", cuts_cmp)

# t 风格作为基准
post_t <- abc_two_group(x1, x2)

vol_t <- det(cov(post_t))

vol_cmp <- sapply(posts_cmp, function(p) {
  det(cov(p))
})

loss_ratio <- vol_cmp / vol_t

plot_overlay(
  c(list("t 风格" = post_t), posts_cmp),
  true_x = 0.0,
  true_y = 0.6,
  xlim = c(-2, 2.5),
  ylim = c(-2, 2.5),
  xlab = expression(mu[1]),
  ylab = expression(mu[2]),
  colors = c(
    "t 风格"  = "#f5a623",
    "cut=-0.6" = "#e040fb",
    "cut=0.3"  = "#22c55e",
    "cut=0.9"  = "#f87171"
  )
)

cat("这张图里，我们比较了不同 cut 值下的信息折损。\n")
cat("cut 越接近两组均值的中点，判别能力越强。\n\n")

best_cut <- names(loss_ratio)[which.min(loss_ratio)]
worst_cut <- names(loss_ratio)[which.max(loss_ratio)]

cat(sprintf("\n%s 的信息保留最多。\n", best_cut))
cat(sprintf("%s 的带来的不确定性最大。\n", worst_cut))

spread_ratio <- vol_cmp / vol_t
retain_ratio <- 1 / spread_ratio
loss_ratio_pct <- (1 - retain_ratio) * 100

summary_loss <- data.frame(
  cut = names(posts_cmp),
  spread_ratio = as.numeric(spread_ratio),
  retain_pct = 100 * as.numeric(retain_ratio),
  loss_pct = as.numeric(loss_ratio_pct)
)


cat("以 t 风格的后验扩散面积为基准（1.00×）：\n\n")

# 指标计算
spread_ratio <- vol_cmp / vol_t
retain_ratio <- 1 / spread_ratio
loss_ratio   <- 1 - retain_ratio
extra_sample <- spread_ratio - 1

summary_loss <- data.frame(
  cut = names(posts_cmp),

  `后验扩散倍率` =
    sprintf("%.2f×", spread_ratio),

  `信息折损` =
    sprintf("%.1f%%", loss_ratio * 100),

  `额外样本需求` =
    sprintf("+%.1f%%", extra_sample * 100),

  check.names = FALSE
)

print(summary_loss, row.names = FALSE)

你会发现，偏离均值的二分图，其喷射面积不仅更大，方向也偏斜了。这是因为极端的切断点让两组的超标比例都趋向同一个极端值：切在 0.9，两组几乎无人能过，计数双双趋近于零；切在 −0.6，两组几乎人人能过，计数趋近于极大值。在这两种情况下，「组 1 的均值高」和「两组均值一起高」对计数的影响几乎一样，数据再也无法区分这两种情形。于是图形沿着一条对角线弥散，代表所有能产生差不多计数比例的参数组合都一样说得通。

这个时候，有的研究者可能会大舒一口气，原来只要切中间点就好了嘛！但是，在实际的分析场景中，你甚至都不知道自己手里的那个数据是不是真的混杂了两个总体，因此这种切分行为非常鲁莽。哪怕它混杂了两个总体，总体参数也是未知的，我们根本无法在分析前获知所谓的「最优阈值」，所以这件事情从原理上就没有办法「合法地做」。就算你开了天眼找到了最优切断点，得到的也只是「最优损失」，离无损失差得很远，甚至有点对自己的数据暴殄天物。

一些鸡贼玩家可能会觉得，我岂不是也可以先慢慢调，看看哪个断点下面能拿到显著结果。但事实是，你只能知道自己拿到的是阴性结果还是阳性结果，但你不知道自己的阴性和阳性是真的还是假的。That’s p-hacking, bro.

概念阐释

让我们回顾整个讨论。机器一直在做相似的工作：我们先设定一个样本，然后猜测它的总体参数可能为何。对每一组猜测的总体参数，它都会问同样的问题：「手里这批数据，能有多好地支持它？」这个支持程度被称作似然。

机器一开始对参数毫无头绪，把所有参数值视为同等可信，这个看数据之前的参数分布叫先验。似然告诉我们哪些参数值更能解释手里的数据，哪些几乎不能，用似然对先验加权之后得到的参数分布叫后验。先验、似然、后验，这三个词描述的是我们这台模拟机器的工作语言，它们属于贝叶斯推断的框架。

全文每一张热力图画的都是后验：数据越精细，后验越集中，我们对总体参数的把握越强，圆圈又小又亮。数据被粗暴处理之后，后验变得分散，许多截然不同的总体都能同样合理地解释手里的数据。

我们用贝叶斯的方法非常直觉地看到了一个重要的事情：对于同样的样本量，数据的品质不同，推断总体的能力就不同。这个「能力」有一个正式的度量工具，叫做费舍信息量。它衡量的是，样本对总体参数的估计究竟能精确到什么程度。后验越集中，费舍信息量越高；后验扩散成喷射状，费舍信息量就低。

虽然我们没有直接计算它（因为我知道让你推公式是在折磨你的脑子），但全文每一张图里圆圈的大小画的都是它们的倒影。唯一值得注意的是，全文以贝叶斯为手段构造间接统计量让你直观地看到似然这个概念，但在大多数统计实践中，我们都是直接套公式来算的。贝叶斯在似然框架下并不是一个必然的东西，只是一个方便你理解的工具，甚至提出似然概念的费舍老先生还挺讨厌贝叶斯的。

那么最后，恭喜你！你已经一只脚踏入了费舍似然体系的大门！

---

1. 你对数据的任何整形处理都不会使其包含的信息量变多，只会使其包含的信息量变少。 ↩︎

2. 这是一些你可能不感兴趣的技术细节。我们如何衡量一个假想的总体与我们手中真实样本的「相似度」。这个实验通过统计参数的欧式距离来完成衡量。距离越小，代表假想总体与真实数据的匹配度越高。对于使用完整连续数据的分析，这个距离是在一个二维空间中计算的。这两个维度分别由样本的两个核心特征定义：样本均值和样本标准差。在计算距离前，我们必须先解决一个校准问题：均值上的差异和标准差上的差异无法直接比较。因此，我们通过一个「预热」模拟，为这两个维度各自确定了一把「标准尺」，以衡量它们各自的自然随机波动。计算时，我们将每个维度的差异用其对应的「标准尺」进行归一化，然后将这两个归一化后的值代入欧氏距离公式，从而得到最终的匹配度分数。相对地，对于二值化数据的分析，整个空间被压缩成了一个一维空间。这个维度唯一的度量就是「样本中超过阈值的个数」。因此，其距离计算也简化为仅仅是这个计数值在经过「标准尺」归一化后的差值。图形的巨大差异，正源于我们从一个信息丰富的二维空间，退化到了一个信息贫乏的一维空间去寻找答案。 ↩︎

3. 二值化版本类比的是卡方检验，因为卡方检验的操作就是把连续变量变成计数再推断。两者的模型假设完全一样，不同的只是提供给电脑的信息多寡。 ↩︎

4. 严格来讲这块做的不是 t 检验和卡方检验，为了让系统可比我们建构了「t 风格统计量」和「卡方风格统计量」。其核心原理与 t 检验、卡方检验没有区别，但是为了让两个系统可比，额外做了一些处理。整个模拟的基本原理与前一个实验完全相同，同样是计算一个「风味距离」来衡量匹配程度，只是计算距离时所依赖的关键指标变得更多。对于模仿传统 t 检验的「t 风格」分析，我们是在一个三维空间中计算距离。读者可能会疑问为何不是四维（两个均值，两个标准差），这是因为我们的模拟严格遵循了传统 t 检验的一个核心前提：两组数据的总体方差相等。因此，我们只需要一个「合并标准差」即可。这样，我们的三个维度就是：第一组的均值、第二组的均值，以及两组合并的标准差。与前述实验一样，这三个维度在计算距离前都经过了各自「标准尺」的归一化处理。相比之下，模仿卡方检验的「卡方风格」分析，则是在一个更简化的二维空间中进行。这两个维度是：第一组超过切割点的人数和第二组超过切割点的人数。这个过程不仅将每个样本的连续信息压缩成了单一的计数值，还完全忽略了数据的离散程度信息。「t 风格」利用三个信息丰富的维度来精确定位参数，而「卡方风格」只依赖两个粗略的维度。维度数量和质量上的削减，正是造成图中信息量巨大损失的直接原因。整个分析过程中，两种风格的分析除了所依赖的关键指标不同外，其余所有设置，包括先验范围、数据生成模型和筛选规则都完全一致。因此，将它们的后验分布图叠加在一起时，图形轮廓大小的差异，就非常直观地反映了不同分析方法在信息提取能力上的高下之分。 ↩︎