让我们来看三个概念:总体、样本、模型。
总体是你希望揭示的事实。比如你对高中生的身高体重数据感兴趣,那么地球上所有的学生都是你的总体。如果你对「西瓜到底甜不甜」,那么全宇宙所有的西瓜都是你的总体。想要测量全世界所有高中生的身高体重是一件非常不现实的事情,把全世界的西瓜全都吃一遍也不可能。况且你也吃不到去年的瓜,更吃不到明年的瓜。所以我们通过抽样的方式,从总体中抽出一小部分,来间接性地代表总体。
抽样
这个时候就引出了两个重点,一个是「抽样」这个动词,一个是「代表」这个动词。
先来看抽样:如果抽样方法有问题,你要研究全世界所有的瓜,但是却只吃了楼下菜市场的瓜,那么这就不是一个好的抽样。或者说得更专业一点,这是一个「有偏抽样」。抽样技术有瑕疵的后果就是,你的样本不能很好地代表总体。但你又生成研究是在研究这个总体,因而整个研究的有效性就要打一个大大的问号。其内部逻辑是概念被偷换,这相当危险。在总体中均匀地抽出样本是一个很复杂也很重要的学问,它甚至统计学的命门,我们不应当因为抽样技术里面不涉及公式,就把它和「统计」这个概念分而治之。
再来看代表:理论上讲所有研究都应当在方法小节当中论述清楚,你的样本为什么能代表总体。这点针对任何统计报告都是有效的,无论是商业研究报告,还是学界的论文。但有些可惜的是,依然有大量的研究忽略了这件事。
尽管抽样技术很重要,但与之有同样重要性的是样本量。如果你的样本量过小,那么无论如何都做不到让样本充分代表总体。让我们来做这个实验,构建一个非常随机的总体分布,连续抽 20 次看看样本之间有多大差异,样本分布对总体分布的「代理」有多失真。
下面是本次实验会用到的工具:
# SHARE, NORUN
source("https://kb.not.ci/rand-dist.R")
library(ggplot2)
library(matrixStats)
draw_samples <- function(population, n_samples, sample_size, compute_ks = FALSE) {
mat <- matrix(
sample(population, n_samples * sample_size, replace = TRUE),
nrow = n_samples
)
stats <- data.frame(
sample_id = seq_len(n_samples),
mean = rowMeans(mat),
sd = matrixStats::rowSds(mat)
)
if (compute_ks) {
stats$ks_d <- vapply(seq_len(n_samples), function(i) {
suppressWarnings(ks.test(mat[i, ], sample(population, 1000))$statistic)
}, numeric(1))
}
list(mat = mat, stats = stats)
}
log_sampling <- function(population, result) {
stats <- result$stats
pop_mean <- mean(population)
pop_sd <- sd(population)
cat("总体:\n")
cat(sprintf(" 均值: %.2f 标准差: %.2f\n\n", pop_mean, pop_sd))
cat("样本与总体的趋近程度(KS statistic D,越小越近):\n")
cat(sprintf(" 平均 D: %.3f 最小 D: %.3f 最大 D: %.3f\n\n",
mean(stats$ks_d), min(stats$ks_d), max(stats$ks_d)))
cat("样本间的不一致性:\n")
cat(sprintf(" 样本均值:%.2f ± %.2f(总体均值 %.2f)\n",
mean(stats$mean), sd(stats$mean), pop_mean))
cat(sprintf(" 样本标准差:%.2f ± %.2f(总体标准差 %.2f)\n",
mean(stats$sd), sd(stats$sd), pop_sd))
}
plot_sampling <- function(population, result) {
mat <- result$mat
stats <- result$stats
n_samples <- nrow(mat)
sample_size <- ncol(mat)
samples_df <- data.frame(
x = as.vector(t(mat)),
sample_id = rep(seq_len(n_samples), each = sample_size)
)
samples_df <- merge(samples_df, stats[, c("sample_id", "ks_d")], by = "sample_id")
samples_df <- samples_df[order(samples_df$ks_d), ]
ggplot() +
geom_density(
data = samples_df,
aes(x = x, group = sample_id, color = ks_d),
alpha = 0.5, linewidth = 0.4
) +
geom_density(
data = data.frame(x = population),
aes(x = x),
color = "black", linewidth = 1
) +
scale_color_gradient(low = "#009FFF", high = "#ec2F4B", name = "KS D")
}
population <- random_dist("what ever you like", xmin = 0, xmax = 100, N = 20000)
下面来实际对一个稀奇古怪的分布做一个抽样,先使用江湖传闻的「n=30」进行 20 次抽样技术完美的均匀抽样,看看都能抽出来啥样子的分布。
# GALLERY
result <- draw_samples(population, n_samples = 20, sample_size = 30, compute_ks = TRUE)
log_sampling(population, result)
plot_sampling(population, result)
不难发现,虽然分布的形状有几分相似,但是抽出来的东西真的是各有各的样子,几乎没有什么稳定性可言。
这里的每一根线都是一次「研究的复现」,即一个独立研究团队从总体做一次抽样,得出来的分布。如果你比较倒霉,哪怕用了完美的抽样技术,依然抽出来很扭曲的样本,那事情就会变得难办。
如果样本的分布本身就非常不稳定,那么你应当可以想象最终的统计结论会有多么不稳定。在统计的话语体系下,不稳定就是统计解读失准的风险。我们与风险相伴但总是希望在成本可接受的情况下将它控制在可接受的水平。比如,我们把样本量翻 10 倍,看看样本量为 300 会为样本的稳定性带来什么影响。
result <- draw_samples(population, n_samples = 20, sample_size = 300, compute_ks = TRUE)
log_sampling(population, result)
plot_sampling(population, result)
可以看到形状虽然依然略有失真,但是已经能很好地表达总体的分布,在这个情况下,进行统计分析,结果就更不容易失准。
建模
做完实验,如果你的老板问你实验怎么样,你把成千上万条裸数据直接丢给老板,那你八成要被炒鱿鱼的。没人能从海量的裸数据里面得到有价值的信息。为了让数据变成智慧,我们需要对它建模。也就是对数据做进一步的抽象,揭示数据中关键要素之间的内生关系。
所有模型都要在模型的简洁性和信息解释能力之间做权衡。最复杂的模型就是原始数据,它不丢失任何数据细节,但是对我们来讲没有任何意义。最简单的模型是算一个均值,但对于复杂问题,这种大刀阔斧的简化很有可能让我们挖不出有价值的信息,进而做出偏误的决策、错失机会。
如何建构一个好模型是这个系列文章主要关注的对象。而今天我想和你探讨的是抽样怎么对建模产生影响。
让我们把上面的实验重复更多次,换言之,会有 1000 个团队,组成千军万马一起去世界各地吃瓜,每个团队都算一个甜度均值出来,看看这一大堆均值能构成一个什么样的分布。请注意,每个团队都是独立复现实验,团队和团队之间不做沟通。
# SHARE, NORUN
plot_stat_dist <- function(population, sample_sizes, stat_fn, stat_name = "统计量", n_samples = 1000) {
records <- lapply(sample_sizes, function(sz) {
mat <- draw_samples(population, n_samples, sz)$mat
vals <- apply(mat, 1, stat_fn)
data.frame(stat = vals, sample_size = factor(sz, levels = sample_sizes))
})
df <- do.call(rbind, records)
pop_stat <- stat_fn(population)
ggplot(df, aes(x = stat, fill = sample_size, color = sample_size)) +
geom_density(alpha = 0.25, linewidth = 0.5) +
geom_vline(xintercept = pop_stat, linetype = "dashed", color = "black", linewidth = 0.6) +
scale_fill_brewer(palette = "Set2", name = "样本量") +
scale_color_brewer(palette = "Set2", name = "样本量") +
labs(x = stat_name, y = "密度") +
annotate("text", x = pop_stat, y = Inf,
label = sprintf("总体%s %.2f", stat_name, pop_stat),
hjust = -0.1, vjust = 1.8, size = 3.2, color = "black")
}
log_stat_dist <- function(population, sample_sizes, stat_fn, stat_name = "统计量", n_samples = 1000) {
pop_stat <- stat_fn(population)
cat(sprintf("总体 %s: %.4f\n\n", stat_name, pop_stat))
for (sz in sample_sizes) {
mat <- draw_samples(population, n_samples, sz)$mat
vals <- apply(mat, 1, stat_fn)
bias <- mean(vals) - pop_stat
se <- sd(vals)
rmse <- sqrt(mean((vals - pop_stat)^2))
cv <- se / abs(pop_stat) * 100
ci <- quantile(vals, c(0.025, 0.975))
cat(sprintf("样本量 n = %d\n", sz))
cat(sprintf(" 均值偏离总体: %+.4f(相对偏差 %+.2f%%)\n", bias, bias / abs(pop_stat) * 100))
cat(sprintf(" 标准误: %.4f RMSE: %.4f 变异系数: %.2f%%\n", se, rmse, cv))
cat(sprintf(" 95%% 区间: [%.4f, %.4f] 区间宽度: %.4f\n\n", ci[1], ci[2], ci[2] - ci[1]))
}
}
我们从很小的样本量开始,比如只抽 5 个人(也就是 1000 个团队,每个团队只啃五个瓜),看看执行 1000 次抽样后,我们会得到什么样的均值。
plot_stat_dist(population, 5, mean, "均值")
log_stat_dist(population, 5, mean, "均值")
你应该能够清晰地看出,我们得到的均值分布非常分散,总体的均值是 70.3,但是我们得到的均值从 40 到 90 都存在,这表示不仅样本的分布不能准确地代表总体,根据这这些样本算出来的均值,和总体的均值也有不小的差异。对于一些对分析精确度较高的领域(比如涉及到人生死的医疗领域),这种误差几乎是不可接受的。那么如果每个团队啃 50 个瓜、500 个瓜呢?
plot_stat_dist(population, c(5, 50, 500), mean, "均值")
log_stat_dist(population, c(5, 50, 500), mean, "均值")
你应该可以看得出来,抽样带来的不确定性迅速收窄,整个抽样分布变得又瘦又高。这是每一个数据分析者都希望看到的结果。
在这里请允许我再强调一次,在这个分布中,你执行的实验只是分布中的一个值。对于每次数据分析,你只会做一次「抽样」的动作,剩下的所有分析都基于「样本对总体有良好代表性」的假设之下。我们都不想当那个倒霉鬼抽到特别离谱的分布,得出非常离谱的均值。但这个风险永远都在,为了控制这个风险,增大样本量是行之有效的方法之一(当但,设计好的实验充分控制噪声也是很重要的方法)。
但这是否意味着我们可以无脑堆高样本量?让我们来继续统计学实验。重复执行数次上面的实验,让从 1 依序增加到 500,做一张热力图,看均值分布的范围收窄的趋势,再做一张点图,来看均值分布的尖峰的上涨速度。
# GALLERY
n_sims <- 1000
n_max <- 2000
x_res <- 300
master_mat <- matrix(
sample(population, n_sims * n_max, replace = TRUE),
nrow = n_sims, ncol = n_max
)
cum_mat <- matrixStats::rowCumsums(master_mat)
means_mat <- sweep(cum_mat, 2, seq_len(n_max), "/")
pop_mean <- mean(population)
key_ns <- c(1, 5, 10, 30, 50, 100, 200, 500)
key_means <- means_mat[, key_ns]
summary_df <- data.frame(
n = key_ns,
sample_mean = round(colMeans(key_means), 4),
se = round(matrixStats::colSds(key_means), 4),
rmse = round(sqrt(colMeans((key_means - pop_mean)^2)), 4),
ci_lo = round(apply(key_means, 2, quantile, 0.025), 4),
ci_hi = round(apply(key_means, 2, quantile, 0.975), 4)
)
summary_df$ci_width <- round(summary_df$ci_hi - summary_df$ci_lo, 4)
x_range <- range(means_mat)
x_grid <- seq(x_range[1], x_range[2], length.out = x_res)
density_mat <- t(apply(means_mat, 2, function(v)
density(v, from = x_range[1], to = x_range[2], n = x_res)$y
))
peak_df <- data.frame(
n = seq_len(n_max),
# 原来:apply(density_mat, 1, max) 和 apply(means_mat, 2, sd)
peak = matrixStats::rowMaxs(density_mat),
se = matrixStats::colSds(means_mat)
)
heat_df <- data.frame(
n = rep(seq_len(n_max), times = x_res),
x = rep(x_grid, each = n_max),
dens = as.vector(density_mat)
)
ggplot(heat_df, aes(x = x, y = n, fill = dens)) +
geom_raster(interpolate = TRUE) +
scale_fill_gradientn(
colors = c("#0d0221", "#009FFF", "#ffffff"),
trans = "sqrt",
name = "密度"
) +
theme_minimal() +
xlim(55, 85) +
theme(panel.grid = element_blank())
ggplot(peak_df, aes(x = n, y = peak)) +
geom_point(size = 0.35, alpha = 0.45, color = "#009FFF") +
geom_smooth(method = "loess", span = 0.15, se = FALSE, color = "#ec2F4B", linewidth = 0.8) +
labs(x = "样本量 n", y = "峰值密度") +
theme_minimal()
se_curve <- peak_df$se
se_15 <- se_curve[15]
se_30 <- se_curve[30]
unit_factor <- se_15 / se_30
cat(sprintf("基准:n=15 → 30,追加了 14 个样本,SE 缩小 %.2f 倍\n\n", unit_factor))
cat("要再获得同等幅度的 SE 缩小,后续各步需要追加多少样本:\n\n")
target <- se_30
prev_n <- 30
for (k in seq_len(6)) {
target <- target / unit_factor
n_hit <- which(se_curve <= target)[1]
if (is.na(n_hit)) break
cat(sprintf(" 第 %d 步:n=%4d → %4d,需追加 %3d 个样本,SE 降至 %.2f\n",
k, prev_n, n_hit, n_hit - prev_n, se_curve[n_hit]))
prev_n <- n_hit
}
从上面的运行结果中,你应该能非常直观的看到,不是一个样本就能带来一份确定性。样本量越大,不确定性收束的效率就会越差,如果你只是想要算个均值的话,在样本量超过 500 之后不确定性收束的效率就会变得很不「值当」。但请注意这里不是鼓励你设置一个样本量为 500 的上限,因为你收回来数据不会只算一个均值,随着模型规模的变大,数据的规模要求也会随之增高。像是卡方检验之类的检验对于数据的细节也有要求,有的时候你必须得堆高数据量才能满足统计方法本身对样本特质的要求。此外,你能接受的风险也是一个比需要被考虑到的因素。如果你真的无法接受过大的风险(比如和人的生死有关),那么哪怕堆高样本量带来的模型改进有限,你可能也想付出更多的边际效应税,换得更稳定的结果。
均值分布的钟形曲线有一些很有趣的用法,比如用来判断一个数据是不是界外值。比如有一个研究团队,在一片神秘的沙漠中吃到了一种瓜,他们不知道这是不是西瓜,但是它吃上去像是「只有一点点甜味但又很辣的海绵」。研究团队一致给它的甜度评了 5 分。在满分为 100 分的评定中这是一个非常低的分数。
分布图中,纵轴越高代表出现这个数字的概率越高。而 5 分是一个出现概率非常小的数字。因此只看甜度指标的话,我们可以相对有信心地认为,「这不是西瓜」。
你可能会问,我又没有千军万马,顾不起 1000 个团队就画不出来这个曲线图,没得做分析了么?统计学家针对这个问题给出了很不错的工具。他们用一组公式来描述这个分布的形状,也就是 t 分布的公式了。
但值得注意的是,如果总体的分布就奇形怪状,你的样本量又很小,那么这个公式和真实的抽样分布就会有差异,这是统计结果不准确性的另外一个来源。
# GALLERY
# 计算经验 t 统计量
draw_t_stat <- function(population, sample_size, n_sims = 5000) {
mu <- mean(population)
# 原来:replicate 逐次抽样
mat <- matrix(sample(population, n_sims * sample_size, replace = TRUE),
nrow = n_sims)
t_vals <- (rowMeans(mat) - mu) / (matrixStats::rowSds(mat) / sqrt(sample_size))
data.frame(t = t_vals)
}
# 叠加理论 t 分布曲线
add_t_curve <- function(sample_size, color = "black", linewidth = 0.8) {
stat_function(
fun = dt,
args = list(df = sample_size - 1),
color = color,
linewidth = linewidth
)
}
# 画经验分布 + 理论 t 曲线
plot_t_compare <- function(population, sample_size, n_sims = 5000, xlim = c(-6, 6)) {
df <- draw_t_stat(population, sample_size, n_sims)
ks_test <- ks.test(df$t, "pt", df = sample_size - 1)
cat(sprintf("样本量 n = %d,模拟 t 分布与理论 t 分布的最大差异 = %.2f%%(KS D = %.4f)\n",
sample_size, ks_test$statistic * 100, ks_test$statistic))
ggplot(df, aes(x = t)) +
geom_density(fill = "#009FFF", alpha = 0.20, linewidth = 0.5, color = "#009FFF") +
add_t_curve(sample_size, color = "#ec2F4B", linewidth = 0.9) +
geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.4) +
coord_cartesian(xlim = xlim) +
labs(
x = "t 统计量",
y = "密度",
title = sprintf("样本量 n = %d", sample_size),
subtitle = sprintf("蓝色为真实抽样分布,红色为理论 t 分布,df = %d", sample_size - 1)
) +
theme_minimal()
}
plot_t_compare(population, sample_size = 5, n_sims = 2000, xlim = c(-8, 8))
plot_t_compare(population, sample_size = 50, n_sims = 2000, xlim = c(-5, 5))
应用
不仅一个变量的均值的参数分布遵循 t 分布,一个变量和某常数之差的运算结果也遵循 t 分布(当然)。利用这个性质,我们可以探究一个变量的均值是否显著大于某个特定数字。这就是单样本 t 检验了。
不仅这样,两个变量之差的参数分布也遵循 t 分布。
population_a <- random_dist("One awesome population", xmin = 0, xmax = 100, N = 20000)
population_b <- random_dist("Another awesome population", xmin = 0, xmax = 100, N = 20000)
draw_two_sample_t <- function(pop_a, pop_b, n_a, n_b, n_sims = 5000) {
mu_d <- mean(pop_a) - mean(pop_b)
mat_a <- matrix(sample(pop_a, n_sims * n_a, replace = TRUE), nrow = n_sims)
mat_b <- matrix(sample(pop_b, n_sims * n_b, replace = TRUE), nrow = n_sims)
means_a <- rowMeans(mat_a)
means_b <- rowMeans(mat_b)
vars_a <- matrixStats::rowVars(mat_a)
vars_b <- matrixStats::rowVars(mat_b)
sp2 <- ((n_a - 1) * vars_a + (n_b - 1) * vars_b) / (n_a + n_b - 2)
se <- sqrt(sp2 * (1 / n_a + 1 / n_b))
t_vals <- ((means_a - means_b) - mu_d) / se
data.frame(t = t_vals)
}
plot_two_sample_t_compare <- function(pop_a, pop_b, n_a, n_b, n_sims = 5000, xlim = c(-6, 6)) {
df <- draw_two_sample_t(pop_a, pop_b, n_a, n_b, n_sims)
df_t <- n_a + n_b - 2
ks_test <- ks.test(df$t, "pt", df = df_t)
cat("总体 A:\n")
cat(sprintf(" 均值: %.2f 标准差: %.2f\n", mean(pop_a), sd(pop_a)))
cat("总体 B:\n")
cat(sprintf(" 均值: %.2f 标准差: %.2f\n\n", mean(pop_b), sd(pop_b)))
cat("两个总体的真实均值差:\n")
cat(sprintf(" mu_A - mu_B = %.4f\n\n", mean(pop_a) - mean(pop_b)))
cat("两独立样本均值差的 t 统计量:\n")
cat(sprintf(" 样本量 n_A = %d, n_B = %d, 理论自由度 df = %d\n", n_a, n_b, df_t))
cat(sprintf(" 经验均值: %.4f 经验标准差: %.4f\n", mean(df$t), sd(df$t)))
cat(sprintf(" 95%% 区间: [%.4f, %.4f]\n", quantile(df$t, 0.025), quantile(df$t, 0.975)))
cat(sprintf(" KS D: %.4f\n\n", ks_test$statistic))
ggplot(df, aes(x = t)) +
geom_density(fill = "#009FFF", alpha = 0.20, linewidth = 0.5, color = "#009FFF") +
stat_function(
fun = dt,
args = list(df = df_t),
color = "#ec2F4B",
linewidth = 0.9
) +
geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.4) +
coord_cartesian(xlim = xlim) +
labs(
x = "t 统计量",
y = "密度",
title = sprintf("两独立样本均值之差:n_A = %d, n_B = %d", n_a, n_b),
subtitle = sprintf("蓝色为经验分布,红色为理论 t 分布,df = %d", df_t)
) +
theme_minimal()
}
plot_two_sample_t_compare(population_a, population_b, n_a = 50, n_b = 50, n_sims = 3000, xlim = c(-6, 6))
这是一个很有用的性质,利用它,我们可以探究两个样本之间的差异是否显著大于 0,进而判断两组数据是否来自同一总体。
不仅均值有抽样分布,方差也有抽样分布。具体而言,当你在一个总体中,按照特定样本量进行无数次抽样之后,得到的方差会服从卡方分布。
# GALLERY
sample_sizes <- c(5, 10, 30, 100)
draw_var_stat <- function(population, sample_size, n_sims = 5000) {
sigma2 <- var(population)
mat <- matrix(sample(population, n_sims * sample_size, replace = TRUE),
nrow = n_sims)
s2 <- matrixStats::rowVars(mat)
chi2 <- (sample_size - 1) * s2 / sigma2
data.frame(chi2 = chi2, s2 = s2)
}
# 图一:四种样本量的模拟抽样分布(归一化到 χ²/df)
sim_records <- lapply(sample_sizes, function(sz) {
result <- draw_var_stat(population, sz, n_sims = 5000)
data.frame(
chi2_norm = result$chi2 / (sz - 1),
sample_size = factor(sz, levels = sample_sizes)
)
})
sim_df <- do.call(rbind, sim_records)
ggplot(sim_df, aes(x = chi2_norm, fill = sample_size, color = sample_size)) +
geom_density(alpha = 0.25, linewidth = 0.5) +
scale_fill_brewer(palette = "Set2", name = "样本量") +
scale_color_brewer(palette = "Set2", name = "样本量") +
coord_cartesian(xlim = c(0, 3)) +
labs(x = "χ² / df", y = "密度") +
theme_minimal()
# 图二:四种样本量的理论卡方分布曲线(归一化到 χ²/df)
x_grid <- seq(0.001, 3, length.out = 500)
theory_records <- lapply(sample_sizes, function(sz) {
df_val <- sz - 1
data.frame(
x = x_grid,
y = dchisq(x_grid * df_val, df = df_val) * df_val,
sample_size = factor(sz, levels = sample_sizes)
)
})
theory_df <- do.call(rbind, theory_records)
ggplot(theory_df, aes(x = x, y = y, color = sample_size)) +
geom_line(linewidth = 0.8) +
scale_color_brewer(palette = "Set2", name = "样本量") +
coord_cartesian(xlim = c(0, 3)) +
labs(x = "χ² / df", y = "密度") +
theme_minimal()
cat("左图是模拟抽样的分布,右图是理论分布")
你可能会好奇,这么个分布能有啥用。的确我们平时不太会单独估测总体的方差,但如果我们想要估计一个模型本身的有效性呢?比如,我们给基本款的模型(简洁模型,Compact Model)模型加了几个参数变成扩增模型(扩增模型,Augmented Model),我们想知道模型对误差的解释能力增加了多少。考虑到方差本身就可以加加减减(参考本系列方差对可加性的介绍),你当然也可以把两个方差除一下算个比例出来:PRE = 1 - (简洁模型还剩多少误差 / 扩增模型还剩多少误差),这个分数可以被进一步做线性转换变成 F 分布[1]。
这样我们就得到方差分析啦!
你可以把你手头的模型和一个什么参数都不带,纯算均值的模型比较,也可以看多几个参数是否值得。统计建模的基本工作之一就是在平衡模型的复杂性和模型的解释能力。暴力加参数总是能换得更好的解释能力,F 检验本质上就是在帮我们把关这件事。
严格来讲,在总体非正态的情况下方差的抽样分布和卡方分布存在较大差异。但这个差异的倾向性非常稳定,在你把两个方差除一下之后,它们会在分母和分子上互相抵消,达到一个相对无偏、稳定的统计指标。属于你知不知道都不太会影响你分析结论的豆知识。 ↩︎

Loading comments...