第八章 参数估计

8.1 什么是参数估计

在统计学中，参数（Parameter）是描述总体分布特征的数值指标。它是一个未知的常数，决定了总体的概率分布。常见的总体参数包括：均值：\(\mu = E(X)\)，描述数据的中心位置；方差：\(\sigma^2 = Var(X) = E[(X-\mu)^2]\)，描述数据的离散程度；比例：\(p\)，描述总体中具有某种特征的比例；相关系数：\(\rho\)，描述两个变量间的相关程度；分布中的其他参数：如泊松分布的\(\lambda\)等。

参数总体的固有属性，通常是未知的，确定的常数。而估计（Estimation）就是指根据样本观测值对总体未知参数进行推断的过程。这一过程可以描述为：从总体中抽取样本 → 计算样本统计量 → 用这个统计量推断总体参数。例如：通过100个人的身高（样本），计算平均身高（样本统计量），从而推断人群的平均身高（参数）。在此过程中，我们将估计总体参数的统计量叫做估计量（Estimator）。

我们给出估计量的规范定义：设总体\(X\)的分布函数为\(F(x;\theta)\)，其中\(\theta\)是未知参数，\(X_1, X_2, ..., X_n\)是来自总体\(X\)的一个样本。构造一个样本函数\(\hat{\theta} = T(X_1, X_2, ..., X_n)\)来估计\(\theta\)，则称\(\hat{\theta}\)为\(\theta\)的估计量。

可以看出，估计量是一个随机变量（因为样本是随机的）。它是样本的函数：\(\hat{\theta} = T(X_1, X_2, ..., X_n)\)。用估计量的一个具体观测值称为估计值。例如，样本均值\(\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i\)是总体均值\(\mu\)的估计量；样本方差\(S^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X})^2\)是总体方差\(\sigma^2\)的估计量。

参数估计主要有两种主要类型：点估计（Point Estimation）和区间估计（Interval Estimation）。点估计是用单个数值来估计总体参数。例如：用样本均值\(\bar{x} = 172.5cm\)估计总体平均身高\(\mu\)。点估计容易理解和计算，但不能给出估计的精度或可靠程度。而区间估计是用一个区间范围来估计总体参数，同时给出该区间包含参数真值的可信程度（置信水平），例如：\(\mu\)的95%置信区间为\([170.2, 174.8]\)。

8.2 点估计

8.2.1 矩估计与最大似然估计

矩估计法（Method of Moments），最大似然估计法（Maximum Likelihood Estimation），贝叶斯估计法（Bayesian Estimation）和最小二乘估计法（Least Squares Estimation）都是常用的点估计方法。因为本书以“频率学派统计学”为主，我们介绍前两种方法。贝叶斯估计法更具备贝叶斯统计学的特征，而最小二乘估计法更像是数学方法而非统计学方法。

a) 矩估计（Method of Moments, MM）

矩（Moment）是描述随机变量分布特征的一系列数字特征。这个词是卡尔皮尔逊从物理学类比过来的。物理学中“质量矩”可以描述质量分布如何决定系统围绕参考点旋转、平衡、扭转的力学趋势。于是，卡尔皮尔逊创造了“概率矩”来描述概率质量分布如何决定分布围绕参考点集中、分散、偏斜的形状趋势。

k阶原点矩(可以对比一下\(mr^k\)的质量矩公式)： \[ \mu_k = E(X^k)=\sum p_i x_i^k, \quad k=1,2,3,... \]

k阶中心矩： \[ \nu_k = E[(X - E(X))^k], \quad k=2,3,4,... \]

由定义可知：一阶原点矩就是期望（均值）：\(\mu_1 = E(X)\)，二阶中心矩是方差：\(\nu_2 = E[(X - \mu)^2]\)，三阶标准化中心矩是偏度（Skewness），四阶标准化中心矩是峰度（Kurtosis）。

矩估计的本质就是用样本矩作为相应总体矩的估计，用样本矩的函数作为总体矩的函数的估计。因为大数定律，当样本容量\(n\)足够大时，样本矩依概率收敛于总体矩。

于是，一般的矩估计的方法步骤：

• 1，计算总体的前\(k\)阶矩（通常用原点矩）
• 2，计算样本的对应矩
• 3，令总体矩等于样本矩，建立方程组
• 4，解方程组得到参数的矩估计

例8-1：

举个正态分布矩估计的例子：假定总体 $ X N(, ^2) $，参数为 \(\mu\) 和 \(\sigma^2\)，做矩估计。

1. 总体矩（因为有两个参数 \(\mu\) 和 \(\sigma^2\)，所以我们只需要前两阶矩）： \[ \mu_1 = E(X) = \mu \] \[ \mu_2 = E(X^2) = \sigma^2 + \mu^2 \]

2. 样本矩： \[ A_1 = \bar{X} \] \[ A_2 = \frac{1}{n}\sum_{i=1}^n X_i^2 \]

3. 建立方程组： \[ \begin{cases} \mu = A_1 \\ \sigma^2 + \mu^2 = A_2 \end{cases} \]

4. 解得： \[ \hat{\mu} = \bar{X} \] \[ \hat{\sigma}^2 = A_2 - \bar{X}^2 = \frac{1}{n}\sum_{i=1}^n (X_i - \bar{X})^2 \]

接下来我们用R语言进行模拟和求解。

# 生成模拟数据 - 正态分布 N(5, 2^2)
set.seed(123)
n <- 1000
true_mu <- 5
true_sigma <- 2
data <- rnorm(n, mean = true_mu, sd = true_sigma)

{
cat("=== 模拟数据 ===\n")
cat(sprintf("真实参数: μ = %.2f, σ = %.2f\n", true_mu, true_sigma))
cat(sprintf("样本量: n = %d\n", n))
cat(sprintf("样本数据 (前10个):\n"))
print(round(data[1:10], 3))

# 矩估计函数
cat("\n=== 矩估计法 ===\n")
  
  # Step 2: 计算样本矩
  n <- length(data)
  A1 <- mean(data)  # 一阶样本矩
  A2 <- sum(data^2) / n  # 二阶样本矩
  
  cat(sprintf("一阶样本矩 (A1): %.4f\n", A1))
  cat(sprintf("二阶样本矩 (A2): %.4f\n", A2))
  
  mu_hat_mm <- A1  # μ的矩估计
  sigma2_hat_mm <- A2 - A1^2  # σ^2的矩估计
  sigma_hat_mm <- sqrt(sigma2_hat_mm)  # σ的矩估计

  # 定义方程组
  equations <- function(vars) {
    mu <- vars[1]
    sigma2 <- vars[2]
    c(
      mu - A1,  
      sigma2 + mu^2 - A2  
    )
  }
  
  # 解析解
  cat("\n解析解方程组:\n")
  cat(sprintf("μ = A1 = %.4f\n",A1))
  cat(sprintf("σ² = A2 - A1² = %.4f - (%.4f)² = %.4f\n", A2, A1, sigma2_hat_mm))
} 

## === 模拟数据 ===
## 真实参数: μ = 5.00, σ = 2.00
## 样本量: n = 1000
## 样本数据 (前10个):
##  [1] 3.879 4.540 8.117 5.141 5.259 8.430 5.922 2.470 3.626 4.109
## 
## === 矩估计法 ===
## 一阶样本矩 (A1): 5.0323
## 二阶样本矩 (A2): 29.2535
## 
## 解析解方程组:
## μ = A1 = 5.0323
## σ² = A2 - A1² = 29.2535 - (5.0323)² = 3.9299

值得注意的是，我们这样得到的 \(\hat{\sigma}^2_M\) 是有偏估计（期望为 \(\frac{n-1}{n} \sigma^2\)），但在矩估计中它就是这样被构造出来的，因为矩估计法不要求无偏性，只要求矩相等。(即，矩估计不保证无偏性)

b) 最大似然估计（Maximum Likelihood Estimation, MLE）

似然（Likelihood），意为“好像是这样”，在概率中是指参数取某个值时，观察到当前样本数据的可能性或概率。（当参数是某个值的时候，样本好像是观察到的这样。）所以，似然函数会以每一个样本出现可能性的连乘形式表达：

设\(X_1, X_2, ..., X_n\)是来自总体\(f(x;\theta)\)的样本，\(\theta\)为未知参数。

• 对于离散分布：似然函数为 \[ L(\theta) = P(X_1=x_1, X_2=x_2, ..., X_n=x_n|\theta) = \prod_{i=1}^n P(X_i=x_i|\theta) \]

• 对于连续分布：似然函数为 \[ L(\theta) = f(x_1,x_2,...,x_n|\theta) = \prod_{i=1}^n f(x_i|\theta) \]

而最大似然估计的本质就是选择使样本观测值出现概率最大的参数值作为参数的估计值。其背后的思想是”已经发生的事情应该是概率最大的事情”。如果某个参数值使得我们观察到的样本出现的可能性最大，那么这个值就是该参数最合理的估计。

于是，最大似然估计就可以表达为：在参数空间\(\Theta\)中寻找\(\hat{\theta}\)，使得 \[ L(\hat{\theta}) = \max_{\theta \in \Theta} L(\theta) \] 或者等价地： \[ \ln L(\hat{\theta}) = \max_{\theta \in \Theta} \ln L(\theta) \]

于是，最大似然估计的方法步骤一般为：

1. 写出似然函数

• 离散总体：\(L(\theta) = \prod_{i=1}^n P(X_i=x_i|\theta)\)
• 连续总体：\(L(\theta) = \prod_{i=1}^n f(x_i|\theta)\) 通常情况下我们会取对数，写成对数似然函数（连加比连乘容易计算）： \[ \ell(\theta) = \ln L(\theta) = \sum_{i=1}^n \ln f(x_i|\theta) \]

2. 对参数求导并令导数为零 \[ \frac{\partial \ell(\theta_1,...,\theta_k)}{\partial \theta_i} = 0 \quad (i=1,...,k) \quad \]

3. 解方程得到参数的极大似然估计\(\hat{\theta}\)（或\(\hat{\theta}_1,...,\hat{\theta}_k\)）

例8-2：

我们还是以正态分布的参数\(\mu\)和\(\sigma^2\)的估计来举例：

1. 对数似然函数： \[ \ell(\mu, \sigma^2) = -\frac{n}{2}\ln(2\pi) - \frac{n}{2}\ln(\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n (x_i-\mu)^2 \]

2. 求导建立方程组： \[ \begin{cases} \frac{\partial \ell}{\partial \mu} = \frac{1}{\sigma^2}\sum_{i=1}^n (x_i-\mu) = 0 \\ \frac{\partial \ell}{\partial \sigma^2} = -\frac{n}{2\sigma^2} + \frac{1}{2(\sigma^2)^2}\sum_{i=1}^n (x_i-\mu)^2 = 0 \end{cases} \]

3. 解方程得： \[ \hat{\mu} = \bar{x} = \frac{1}{n}\sum_{i=1}^n x_i \] \[ \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2 \]

值得注意的是，\(\hat{\sigma}^2\)是\(\sigma^2\)的有偏估计，期望为\(E(\hat{\sigma}^2) = \frac{n-1}{n}\sigma^2\)。(最大似然估计也不保证无偏性。)

接下来我们用R语言进行模拟和求解。

  n <- length(data)
  
  # 写出正态分布的对数似然函数
  log_likelihood <- function(params) {
    mu <- params[1]
    sigma <- exp(params[2])  # 使用log(sigma)确保sigma>0
    # ℓ(μ, σ²) = -n/2 ln(2π) - n/2 ln(σ²) - 1/(2σ²) Σ(x_i-μ)²
    -(n/2)*log(2*pi) - n*log(sigma) - sum((data - mu)^2)/(2*sigma^2)
  }
  
  # 负对数似然函数（用于最小化）
  neg_log_lik <- function(params) {
    -log_likelihood(params)
  }
  
  # 通过数值优化求解
  initial_params <- c(mean(data), log(sd(data)))
  
  # 使用优化算法
  mle_result <- optim(
    par = initial_params,
    fn = neg_log_lik,
    method = "BFGS",
    hessian = TRUE,
    control = list(trace = 0)  # 不显示优化过程
  )
  
  # 参数估计
  mu_hat_mle <- mle_result$par[1]
  sigma_hat_mle <- exp(mle_result$par[2])
  sigma2_hat_mle <- sigma_hat_mle^2
  
  {
  # 展示似然函数表达式
  cat("对数似然函数表达式:\n")
  cat("ℓ(μ, σ²) = -n/2 ln(2π) - n/2 ln(σ²) - 1/(2σ²) Σ(x_i-μ)²\n\n")
  cat("优化结果:\n")
  cat(sprintf("μ的MLE估计: %.4f\n", mu_hat_mle))
  cat(sprintf("σ的MLE估计: %.4f\n", sigma_hat_mle))
  cat(sprintf("σ²的MLE估计: %.4f\n", sigma2_hat_mle))
  # 计算似然值
  log_lik_value <- log_likelihood(c(mu_hat_mle, log(sigma_hat_mle)))
  cat(sprintf("最大对数似然值: %.4f\n", log_lik_value))
  }

## 对数似然函数表达式:
## ℓ(μ, σ²) = -n/2 ln(2π) - n/2 ln(σ²) - 1/(2σ²) Σ(x_i-μ)²
## 
## 优化结果:
## μ的MLE估计: 5.0323
## σ的MLE估计: 1.9824
## σ²的MLE估计: 3.9299
## 最大对数似然值: -2103.2458

8.2.2 估计量的评价标准

根据不同估计方法得到可用于估计总体参数的统计量。但既然是“估计”，就需要对其优劣进行评价。最直观的评价标准是：将样本数据代入估计量计算得到的估计值，是否恰好等于待估参数的真实值。然而，这一标准在实际应用时面临两个根本性的问题：第一，在绝大多数实际问题中，我们并不知道总体参数的真实值，否则就无需进行估计。第二，由于抽样具有随机性，即使采用同一个估计量，基于不同样本计算得到的估计值也会波动，单次估计结果很可能因抽样误差而偏离参数真值。

因为这两个问题，现代统计学并不以“单次估计值完全等于真值”作为判定标准，而是转向对估计量本身统计性质的系统性评价:

对于“真值不可知”的困境，统计推断通常在假设总体分布已知的情形下进行分析.尽管实际应用中该假设未必成立，但这样的理论设定允许我们先在“参数真值可知”的框架下，通过数学推导或数值模拟，研究估计量的各类性质。一旦某种估计量在理论或模拟中被论证具有良好的统计表现，它便可以在实际问题中被推广使用。

针对“抽样随机性导致的估计波动”，统计学建立起一系列评价准则，主要包括无偏性、有效性和一致性等。这些准则并不要求某次估计完全准确，而是从期望、方差以及大样本下的收敛趋势等角度，对估计量的长期表现或渐近行为提出规范，从而为比较和选择估计量提供了系统、可靠的理论依据。

a). 无偏性（Unbiasedness）

如果估计量\(\hat{\theta}\)的期望等于参数真值\(\theta\)，即 \[ E(\hat{\theta}) = \theta \] 则称\(\hat{\theta}\)是\(\theta\)的无偏估计量。

无偏性意味着在重复抽样下，用该估计量进行多次估计，所得估计值的平均值会趋于参数真值 \(\theta\)。这说明了估计量不存在系统性的高估或低估倾向，即其“平均而言”是正确的。可以看出，无偏性是对估计量长期表现的衡量，而不是对某一次抽样结果的评价。某次具体估计的数值可能偏离真值，但大量独立估计的平均值会收敛到真值。无偏性并不保证估计量的方差小，也不保证某一次估计比较精确。无偏估计量不一定唯一，同一个参数可能存在多个无偏估计量。

根据无偏性的定义，我们可以得知：

1. 样本均值是总体均值的无偏估计： \[ E(\bar{X}) = E\left(\frac{1}{n}\sum_{i=1}^n X_i\right) = \frac{1}{n}\sum_{i=1}^n E(X_i) = \mu \]

2. 样本二阶中心矩（除以 \(n\) 的样本方差）是总体方差的有偏估计量（Biased estimator）

\[S^2_{n} = \frac{1}{n} \sum_{i=1}^{n}(X_i - \bar{X})^2, \quad \text{其中} \; \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i\]

通过数理统计推导可得： \[E[S^2_{n}] = \sigma^2 - \frac{\sigma^2}{n} = \sigma^2 \cdot \left( 1 - \frac{1}{n} \right)\]

显然，\(\mathbb{E}[S^2_{biased}] \neq \sigma^2\)（除非 \(n \to \infty\)）。
所以 \(S^2_{biased}\) 是一个有偏估计量，其偏差是 \(-\sigma^2 / n\)，会系统性地低估总体方差。

3. 除以 \(n-1\) 的样本方差是总体方差的有偏估计量无偏估计量（Unbiased estimator）

\[S^2_{n-1} = \frac{1}{n-1} \sum_{i=1}^{n}(X_i - \bar{X})^2\]

可以推导： \[E[S^2_{n-1}] = \sigma^2\]

因此 \(S^2_{n-1}\) 是总体方差 \(\sigma^2\) 的无偏估计量。

例8-3：

设总体 \(X\) 服从正态分布 \(N(5, 2^2)\)，即总体方差真值 \(\sigma^2 = 4\)。

请通过蒙特卡洛模拟验证\(S_n^2\)和\(S_{n-1}^2\)，哪个哪个估计量满足“无偏性”。

# 设置参数
library(tidyverse)
set.seed(123)
mu <- 5
sigma <- 2
n <- 10
simulations <- 1000

# 模拟估计量
biased_variance <- numeric(simulations)
unbiased_variance <- numeric(simulations)

for (i in 1:simulations) {
  # 从正态分布生成样本
  sample_data <- rnorm(n, mean = mu, sd = sigma)
  
  # 有偏估计（除以n）
  biased_variance[i] <- mean((sample_data - mean(sample_data))^2)
  
  # 无偏估计（除以n-1）
  unbiased_variance[i] <- var(sample_data)
}

{
# 计算期望值
cat("理论总体方差:", sigma^2, "\n")
cat("有偏估计的期望值:", mean(biased_variance), "\n")
cat("无偏估计的期望值:", mean(unbiased_variance), "\n")
}

## 理论总体方差: 4 
## 有偏估计的期望值: 3.594174 
## 无偏估计的期望值: 3.993527

# 可视化比较
variance_data <- data.frame(
  Biased = biased_variance,
  Unbiased = unbiased_variance
) %>% 
  pivot_longer(cols = everything(), names_to = "Estimator", values_to = "Variance")

# 计算各估计量的均值
mean_biased <- mean(biased_variance)
mean_unbiased <- mean(unbiased_variance)

# 图像比较
ggplot(variance_data, aes(x = Variance, fill = Estimator)) +
  geom_density(alpha = 0.5) +
  geom_vline(xintercept = sigma^2, linetype = "dashed", color = "red", size = 1.2) +
  geom_vline(xintercept = mean_biased, linetype = "solid", 
             color = "orange", size = 1.2, alpha = 0.8) +
  geom_vline(xintercept = mean_unbiased, linetype = "solid", 
             color = "lightblue", size = 1.2, alpha = 0.8) +
  annotate("text", x = sigma^2, y = 0.08, 
           label = paste0("真值\n", sigma^2),
           color = "red", fontface = "bold", size = 4) +
  labs(subtitle = paste("方差估计量的无偏性比较:","n =", n, "，模拟次数 =", simulations),
       x = "方差估计值",
       y = "密度",
       caption = "红色虚线：总体方差真值 | 彩色实线：各估计量均值") +
  theme_minimal() +
  theme( plot.subtitle = element_text( size = 12),
         legend.position = "top"
  ) +
  scale_fill_manual(
    values = c("Biased" = "orange", "Unbiased" = "lightblue"),
    labels = c("Biased" = paste0("有偏估计 (均值=", round(mean_biased, 3), ")"),
               "Unbiased" = paste0("无偏估计 (均值=", round(mean_unbiased, 3), ")"))
  )

b). 有效性（Efficiency）

前面讨论的无偏性反映了估计量没有系统性偏差。然而，即使同样是“平均正确”的无偏估计，我们也希望它尽可能稳定，这就是有效性准则。有效性基于一个直观的想法：估计量的方差越小，其结果的波动性就越弱，用单次样本得出的估计值也就越有可能接近真值。因此，有效性准则旨在比较不同无偏估计量的精度高低。

设 \(\hat{\theta}_1\) 和 \(\hat{\theta}_2\) 都是参数 \(\theta\) 的无偏估计，如果它们的方差满足 \[ \operatorname{Var}\left(\hat{\theta}_1\right) < \operatorname{Var}\left(\hat{\theta}_2\right), \] 则称 \(\hat{\theta}_1\) 比 \(\hat{\theta}_2\) 更有效。

有效性关注无偏估计量的波动程度，方差越小的估计，其“代表性”就越好，在同一参数的无偏估计量中更优。

例8-4：

以正态总体 \(N(\mu, \sigma^2)\) 为例，\(\mu\) 的两个常用无偏估计分别是：

1. 样本均值： \(\bar{X} = \frac{1}{n} \sum_{i=1}^n X_i\)；

2. 样本中位数： \(M\)。（没错，因为正态分布均值、中位数都是一个数，所以样本中位数也是\(\mu\)的无偏估计量）

那么，样本均值和样本中位数，哪个是更有效的均值估计量呢？

可以证明两者的期望均是 \(\mu\)，但方差有明显差异（计算过程略）：

• \(\operatorname{Var}(\bar{X}) = \dfrac{\sigma^2}{n}\)；

• 当\(n\) 较大时，样本中位数 \(M\) 的方差近似为
  \(\operatorname{Var}(M) \approx \dfrac{\pi \sigma^2}{2n} \approx 1.57 \dfrac{\sigma^2}{n}\)。

显然 \(\operatorname{Var}(\bar{X}) < \operatorname{Var}(M)\)，因此 \(\bar{X}\) 比 \(M\) 更有效，意味着用样本均值估计总体均值通常具有更高的精度。

接下来我们通过用R语言进行模拟验证。

# 比较样本均值和样本中位数的有效性
library(tidyverse)
library(ggplot2)

# ---数据准备 ---
set.seed(123)
mu <- 5       # 总体均值真值
sigma <- 2
n <- 30       # 样本量
simulations <- 1000

# 生成模拟数据
results <- data.frame(
  simulation = 1:simulations,
  sample_mean = numeric(simulations),
  sample_median = numeric(simulations)
)

# 循环模拟
for (i in 1:simulations) {
  sample_data <- rnorm(n, mean = mu, sd = sigma)
  results$sample_mean[i] <- mean(sample_data)
  results$sample_median[i] <- median(sample_data)
}

# 转换为长表，方便作图
tidy_results <- pivot_longer(
  results, 
  cols = c(sample_mean, sample_median), 
  names_to = "Estimator", 
  values_to = "estimate"
) %>%
  mutate(
    Estimator = factor(Estimator, 
                       levels = c("sample_mean", "sample_median"),
                       labels = c("样本均值", "样本中位数"))
  )

# 计算统计量用于标注
stats_summary <- tidy_results %>%
  group_by(Estimator) %>%
  summarise(
    Mean_Val = mean(estimate),
    Std_Val = sd(estimate),
    .groups = "drop"
  )

# --- 画图 ---

plot_data <- tidy_results

p <- ggplot(plot_data, aes(x = estimate, color = Estimator)) +
  geom_histogram(aes(y = ..density.., fill = Estimator), 
                 bins = 20, alpha = 0.3, color = NA, position = "identity") +
  geom_density(linetype = 1, size = 1.2) +
  
  # 红色实线：真值
  geom_vline(xintercept = mu, linetype = "solid", color = "red", size = 1) +
  geom_vline(data = stats_summary, 
             aes(xintercept = Mean_Val), 
             linetype = "dashed", size = 1,color = c("orange","steelblue")) +

  labs(
    subtitle = paste("总体分布：N(", mu, ",", sigma^2, ") | 样本量 n =", n, "| 模拟次数:", simulations),
    x = "估计值",
    y = "密度",
    caption = "说明：红色实线为总体真值；彩色虚线为模拟的期望位置。"
  ) +
  # 颜色使用 Manual 指定
  scale_fill_manual(values = c("样本均值" = "orange", "样本中位数" = "steelblue")) +
  scale_color_manual(values = c("样本均值" = "orange", "样本中位数" = "steelblue"),
                     guide = "none") + # 因为我们加了 scale_color_identity，这里要禁用原本图例以免冲突
  
  theme_minimal(base_size = 14) +
  theme(
    plot.subtitle = element_text(hjust =0, size = 11),
    legend.position = "top",
    panel.grid.minor = element_blank()
  )

print(p)

# --- 输出结果 ---
{
cat("=== 模拟结果汇总 ===\n")
print(as.data.frame(
  stats_summary %>%
    rename(
      `估计量` = Estimator, 
      `模拟均值` = Mean_Val, 
      `标准差` = Std_Val
    ) %>%
    mutate(`模拟均值` = round(`模拟均值`, 3), `标准差` = round(`标准差`, 3))
), row.names = FALSE) # row.names = FALSE 可以顺便把行号 1, 2 去掉

cat("\n直观结论：\n1. 两个分布的中心（模拟均值）都紧密围绕在真值", mu, "附近 -> 两者都具有无偏性。\n")
cat("2. 样本均值（绿色）的分布更窄更高，标准差更小 -> 样本均值更有效。\n")
}

## === 模拟结果汇总 ===
##      估计量 模拟均值 标准差
##    样本均值    4.988  0.355
##  样本中位数    4.982  0.436
## 
## 直观结论：
## 1. 两个分布的中心（模拟均值）都紧密围绕在真值 5 附近 -> 两者都具有无偏性。
## 2. 样本均值（绿色）的分布更窄更高，标准差更小 -> 样本均值更有效。

c). 一致性（Consistency）

无偏性和有效性评价的是估计量在固定样本量的表现，而一致性则关注样本量增大时的长期趋势。一个好的估计应当在我们获得越来越多的数据时，无限逼近于参数真值。

对估计量序列 \(\hat{\theta}_n\)，如果对任意小的 \(\epsilon > 0\)，都有
\[ \lim_{n \to \infty}P\bigl(|\hat{\theta}_n - \theta| < \epsilon \bigr) = 1, \]
则称 \(\hat{\theta}_n\) 是 \(\theta\) 的一致估计量。

一致性描述的是大样本下的可靠性：当样本容量 \(n\)足够大时，估计值落在真值附近任意小范围内的概率趋近于1。换言之，只要数据量充分，我们就几乎可以确定估计结果与目标参数非常接近。

一致性并不要求估计必须无偏。实际上，只要偏差与方差均随 \(n\) 趋于 0，估计仍然是一致的。
与有效性不同，一致性一般不用于比较两个估计的优劣，而是作为一种“基本门槛”：不一致的估计会在样本量增加时依然远离真值，这在应用中往往是不可接受的。

例8-5：

利用R语言模拟并可视化样本均值作为总体均值估计量的一致性。

# 一致性：随着n增大，估计值越来越接近真值
set.seed(123)
mu <- 10  # 总体均值真值
sigma <- 5
max_n <- 1000
trials <- 50  # 重复试验次数

# 模拟不同样本量下的估计
all_trials <- list()

for (trial in 1:trials) {
  # 生成一大组数据
  data_pool <- rnorm(max_n, mean = mu, sd = sigma)
  
  # 计算累积均值（即不同n时的估计）
  cum_means <- cumsum(data_pool) / (1:max_n)
  
  all_trials[[trial]] <- data.frame(
    n = 1:max_n,
    trial = trial,
    estimate = cum_means
  )
}

# 合并所有试验
consistency_df <- do.call(rbind, all_trials)

# 绘制一致性图
ggplot(consistency_df, aes(x = n, y = estimate, group = trial)) +
  geom_line(alpha = 0.2, color = "gray70") +
  geom_hline(yintercept = mu, color = "red", linetype = "solid", size = 1) +
  # 添加平均趋势线
  stat_summary(aes(group = 1, y = estimate), 
               fun = mean, geom = "line", 
               color = "blue", size = 1) +
  labs(title = "估计量的一致性演示",
       subtitle = paste("蓝色实线：试验平均值|红色水平线：真值 μ =", mu,
                       "|灰色线条：各次试验路径"),
       x = "样本量 n",
       y = "样本均值估计值",
       caption = "在所有情况下，随着n增大，估计值越来越接近真值且波动减小") +
  theme_minimal() +
  coord_cartesian(ylim = c(mu - 3, mu + 3))

无偏性、有效性和一致性是从三个不同维度评价和选择估计量的准则。无偏性指估计量的评估基点是否准确，即估计值的期望是否等于真值，避免系统性偏差。有效性指估计量的评估精度是否高，即在所有无偏估计中，方差是否最小，结果是否足够稳定、波动小。一致性指估计量的评估过程是否可靠，即随着样本量增大，估计值是否必然会无限接近真值，确保大样本下的可信赖性。

这三个准则在数学上并没有必然的包含关系，但如果需要在理论层面建立优先级，一致性通常处于首要地位，它是大样本（渐近）理论的基石，也是统计推断可靠性的根本保证。相比之下，无偏性和有效性更多关注估计量在有限样本下的表现，有时为了满足一致性，可以适当放松无偏性（这正是最大似然估计方法广受应用的原因之一：它通常具有良好的渐近性质，即使在小样本下可能存在一定的偏差）。在满足无偏性的前提下，有效性则成为进一步筛选估计量的次要标准，因为它衡量的实质是在所有无偏估计中寻找方差最小、精度最高的那一个。因此，从优先逻辑来看，一般先追求一致性，再保障无偏性，最后在无偏类中寻求有效性。

8.3 区间估计

8.3.1 区间估计的基本思想

点估计的基本思想是选择一个合适的统计量（如样本均值、样本方差），用它的观测值作为总体参数的估计值。这种适合是通过无偏性、有效性和一致性来衡量。然而，这种方式存在一个根本性缺陷：缺乏对不确定性的度量。一个点估计值可以告诉我们参数的”最佳猜测”，但没有回答“这个估计有多可靠？可能的误差范围有多大？如果再做一次抽样，结果会有多大变化…”等重要信息。

例如，用样本均值 \(\bar{X} = 5\) 估计总体均值 \(\mu\)，我们可能知道5是这个样本给出的最佳总体均值参数估计，但却无法知道这种最佳估计有着怎样的误差（本质是忽略了抽样误差：这次样本的均值是5，下次样本均值是6。5是正确的还是6是正确的？这个问题不解决，统计推断就会退回到定性论战的阶段：试想，一个出名的团队做了一次调查，估计参数是5。而尚未成名的你，看着很努力的调查结果，估计的参数是6，会不会深深地怀疑自己）。

为了解决点估计的局限性，统计学家提出了区间估计的方法。这种方法基于三个基本认识：

• 正视抽样误差：承认由于随机抽样的偶然性，任何统计量都有波动。

• 利用抽样分布：研究统计量在重复抽样下的分布规律。

• 构造概率保证：基于抽样分布，构建区间，以高概率包含参数真值（而非确定真值）的方式进行估计。

我们把以上思想具象化，来推导区间估计的步骤：

首先，我们需要确定，这个区间有多高的概率包含参数真值。这里引入个概念，置信水平（在第五章我们定性地了解过）。置信水平用 \(1-\alpha\) 表示，其中 \(\alpha\) 是显著性水平，通常取 0.05、0.01 或 0.10。（这里需要解释一下，置信水平是用显著性水平来定义的，而显著性水平是后面假设检验的一个术语。它们本来就是一个硬币的两面。但由于区间估计的思想出现时，假设检验思想就已经有雏形了。区间估计就直接从假设检验里面拿了 显著性水平\(\alpha\)来定义置信水平。）置信水平 \(1-\alpha\) 的含义是：在重复抽样下，构造的区间中包含参数真值的比例约为 \((1-\alpha)\times 100\%\)。

以估计正态总体均值 \(\mu\)为例：我们可能希望构造一个区间，使得它有95%的把握包含真实的 \(\mu\)，此时置信水平 \(1-\alpha = 0.95\)，\(\alpha = 0.05\)。

其次，我们需要回到抽样分布。由于我们可能要估计的参数太多，不可能为每个参数都单独研究其点估计量的精确抽样分布（因为点估计量的分布通常依赖于未知参数本身，使其难以直接用于概率推断），我们于是通过构造枢纽量（Pivotal Quantity） 来解决这一难题。

枢纽量是以一个良好的点估计量为核心，通过将其与待估参数进行恰当的数学组合（如标准化）而构造出来的。具体来说，枢轴量是一个包含样本观测值（通常体现为点估计量）和待估参数\(\theta\)的函数。其最关键的特性在于，经过这种组合变换后，它的概率分布变得完全已知且不依赖于任何未知参数。这使我们能够利用已知的分布（如标准正态分布、t分布）进行概率计算。

在上面的例子中（已知方差 \(\sigma^2\)）：我们选择样本均值 \(\bar{X}\) 作为点估计。已知其抽样分布为 \(\bar{X} \sim N(\mu, \sigma^2/n)\)。据此可以构造枢轴量： \[ Z = \frac{\bar{X} - \mu}{\sigma / \sqrt{n}} \] 该统计量服从标准正态分布，即 \(Z \sim N(0,1)\)。其分布完全已知，不依赖于未知参数 \(\mu\)。

最后，我们将第一步确定的置信概率应用于第二步中构造的枢轴量的已知分布上,就可以找出一个区间了。虽然满足该置信概率的区间有无数种选择，但为了规范性和通常情况下的优良性（如区间长度最短），我们普遍采用对称的方式，即将犯错的风险 \(\alpha\) 平均分配在分布的两侧。

继续上面的例子，对于标准正态分布，我们找到双侧分位数 \(z_{\alpha/2}\)，使得 \(P(-z_{\alpha/2} \leq Z \leq z_{\alpha/2}) = 1-\alpha\)。代入枢轴量 \(Z\) 的表达式： \[ P\left(-z_{\alpha/2} \leq \frac{\bar{X} - \mu}{\sigma / \sqrt{n}} \leq z_{\alpha/2}\right) = 1-\alpha \] 通过解这个关于 \(\mu\) 的不等式，即可得到\(\mu\)的置信区间： \[ \left[ \bar{X} - z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}},\ \bar{X} + z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}} \right] \] 这就是总体均值 \(\mu\) 在置信水平为 \(1-\alpha\) 下的置信区间。

8.3.2 几种典型的参数区间估计

a) 总体均值的区间估计（单个总体，方差已知）

前提条件：总体服从正态分布 \(X \sim N(\mu, \sigma^2)\)，且总体方差 \(\sigma^2\) 已知。

1. 确定置信水平：设定置信水平为 \(1-\alpha\)（如 0.95 或 0.99）。

2. 构造枢轴量：选择样本均值 \(\bar{X}\) 作为点估计。其枢轴量为Z统计量： \[ Z = \frac{\bar{X} - \mu}{\sigma / \sqrt{n}} \sim N(0, 1) \] 该统计量服从标准正态分布，其分布不依赖于未知参数 \(\mu\)。

3. 构建对称置信区间： 根据标准正态分布的双侧分位数 \(z_{\alpha/2}\)，有： \[ P\left(-z_{\alpha/2} \le \frac{\bar{X} - \mu}{\sigma / \sqrt{n}} \le z_{\alpha/2}\right) = 1-\alpha \] 将不等式等价变形为关于 \(\mu\) 的表达式： \[ P\left(\bar{X} - z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}} \le \mu \le \bar{X} + z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}}\right) = 1-\alpha \] 因此，总体均值 \(\mu\) 的 \(1-\alpha\) 置信区间为： \[ \left[ \bar{X} - z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}},\ \bar{X} + z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}} \right] \]

例8-6：

研究某大学学生每月学习时间（小时）。根据历年数据，学习时间服从正态分布，总体标准差σ=12小时。现随机抽取100名学生（n=100），测得样本均值为48小时。

求总体均值μ的95%置信区间。

1. 确定置信水平： 设定置信水平 \(1-\alpha = 0.95\)，则显著性水平 \(\alpha = 0.05\)。

2. 构造枢轴量： 选择样本均值 \(\bar{X}\) 作为点估计。其枢轴量为Z统计量： \[ Z = \frac{\bar{X} - \mu}{\sigma / \sqrt{n}} \sim N(0, 1) \]

3. 构建对称置信区间并代入数值： 根据标准正态分布的双侧分位数 \(z_{\alpha/2}\)，有： \[ P\left(-z_{\alpha/2} \le \frac{\bar{X} - \mu}{\sigma / \sqrt{n}} \le z_{\alpha/2}\right) = 1-\alpha \] 将不等式等价变形为关于 \(\mu\) 的表达式，并代入已知数值 \(\bar{X} = 48\)，\(\sigma = 12\)，\(n = 100\)，以及 \(\alpha = 0.05\) 对应的临界值 \(z_{0.025} = 1.96\)：

\[ \begin{aligned} P\left(\bar{X} - z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}} \le \mu \le \bar{X} + z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}}\right) &= 1-\alpha \\ \left[ 48 - 1.96 \times \frac{12}{\sqrt{100}}, \quad 48 + 1.96 \times \frac{12}{\sqrt{100}} \right] & \\ \end{aligned} \] 因此，总体均值 \(\mu\) 的 \(1-\alpha\) 置信区间为： \[ \left[ 45.648, \quad 50.352 \right] \] R语言可以帮助我们很快地完成计算：

# 已知参数
n <- 100
sample_mean <- 48
sigma <- 12
alpha <- 0.05

# 计算标准误和临界值
standard_error <- sigma / sqrt(n)
z_critical <- qnorm(1 - alpha/2)  # 约为1.96

# 计算置信区间
ci_lower <- sample_mean - z_critical * standard_error
ci_upper <- sample_mean + z_critical * standard_error

cat(sprintf("总体均值μ的95%%置信区间为: [%.3f, %.3f] 小时\n", ci_lower, ci_upper))

## 总体均值μ的95%置信区间为: [45.648, 50.352] 小时

b) 总体均值的区间估计（单个总体，方差未知）

前提条件：总体服从正态分布 \(X \sim N(\mu, \sigma^2)\)，且总体方差 \(\sigma^2\) 未知。

1. 确定置信水平：设定置信水平为 \(1-\alpha\)。
2. 构造枢轴量：由于 \(\sigma\) 未知，用样本标准差 \(S\) 代替。构造的枢轴量为： \[ T = \frac{\bar{X} - \mu}{S / \sqrt{n}} \sim t(n-1) \] 该统计量服从自由度为 \(n-1\) 的 t 分布，其分布不依赖于未知参数 \(\mu\) 和 \(\sigma\)。
3. 构建对称置信区间： *根据 t 分布的双侧分位数 \(t_{\alpha/2}(n-1)\)，有： \[ P\left(-t_{\alpha/2}(n-1) \le \frac{\bar{X} - \mu}{S / \sqrt{n}} \le t_{\alpha/2}(n-1)\right) = 1-\alpha \] 将不等式等价变形为关于 \(\mu\) 的表达式： \[ P\left(\bar{X} - t_{\alpha/2}(n-1) \cdot \frac{S}{\sqrt{n}} \le \mu \le \bar{X} + t_{\alpha/2}(n-1) \cdot \frac{S}{\sqrt{n}}\right) = 1-\alpha \] 因此，总体均值 \(\mu\) 的 \(1-\alpha\) 置信区间为： \[ \left[ \bar{X} - t_{\alpha/2}(n-1) \cdot \frac{S}{\sqrt{n}},\ \bar{X} + t_{\alpha/2}(n-1) \cdot \frac{S}{\sqrt{n}} \right] \]

例8-7：

研究城市居民通勤时间。假设通勤时间服从正态分布，但方差未知。随机调查25名居民（n=25），得到通勤时间数据如下： [25, 30, 35, 40, 45, 28, 32, 38, 42, 33, 29, 36, 41, 27, 34, 39, 31, 37, 43, 26, 44, 30, 35, 32, 38]

估计通勤时间的95%置信区间。

1. 确定置信水平：

设定置信水平 \(1-\alpha = 0.95\)，则显著性水平 \(\alpha = 0.05\)。

2. 构造枢轴量：

选择样本均值 \(\bar{X}\) 作为点估计，样本标准差 \(S\) 代替总体标准差。其枢轴量为t统计量： \[ t = \frac{\bar{X} - \mu}{S / \sqrt{n}} \sim t(n-1) \]

3. 构建对称置信区间并代入数值： 根据 \(t\) 分布的双侧分位数 \(t_{\alpha/2}(n-1)\)，有： \[ P\left(-t_{\alpha/2}(n-1) \le \frac{\bar{X} - \mu}{S / \sqrt{n}} \le t_{\alpha/2}(n-1)\right) = 1-\alpha \] 将不等式等价变形为关于 \(\mu\) 的表达式，并代入已知数值 \(\bar{X} = 34\)，\(S \approx 5.717\)，\(n = 25\)，以及 \(\alpha = 0.05\)、自由度 \(df=24\) 对应的临界值 \(t_{0.025}(24) \approx 2.064\)：

\[ \begin{aligned} P\left(\bar{X} - t_{\alpha/2}(n-1) \cdot \frac{S}{\sqrt{n}} \le \mu \le \bar{X} + t_{\alpha/2}(n-1) \cdot \frac{S}{\sqrt{n}}\right) &= 1-\alpha \\ \left[ 34 - 2.064 \times \frac{5.717}{\sqrt{25}}, \quad 34 + 2.064 \times \frac{5.717}{\sqrt{25}} \right] & \\ \end{aligned} \] 因此，总体均值 \(\mu\) 的 \(95%\) 置信区间为： \[ \left[ 32.399, \quad 37.201 \right] \]

# 模拟通勤时间数据（分钟）
set.seed(123)
commute_time <- c(25, 30, 35, 40, 45, 28, 32, 38, 42, 33, 
                  29, 36, 41, 27, 34, 39, 31, 37, 43, 26, 
                  44, 30, 35, 32, 38)

# 使用更简单的t检验计算置信区间
t_test_result <- t.test(commute_time, conf.level = 0.95)
cat(sprintf("总体均值μ的95%%置信区间为: [%.3f, %.3f] 小时\n", t_test_result$conf.int[1], t_test_result$conf.int[2]))

## 总体均值μ的95%置信区间为: [32.399, 37.201] 小时

c) 总体比例的区间估计

前提条件：大样本情况下，对总体中具有某种特征的个体所占比例 \(p\) 进行估计。

1. 确定置信水平：设定置信水平为 \(1-\alpha\)。
2. 构造枢轴量：选择样本比例 \(\hat{p} = \frac{X}{n}\)（\(X\) 为具有该特征的个体数）作为点估计。由中心极限定理，在大样本下，枢轴量近似服从标准正态分布： \[ Z = \frac{\hat{p} - p}{\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}} \overset{\text{近似}}{\sim} N(0, 1) \] （注：标准误中的 \(p\) 用其点估计 \(\hat{p}\) 代替）。
3. 构建对称置信区间： 根据标准正态分布的双侧分位数 \(z_{\alpha/2}\)，有近似概率： \[ P\left(-z_{\alpha/2} \le \frac{\hat{p} - p}{\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}} \le z_{\alpha/2}\right) \approx 1-\alpha \] 将不等式等价变形为关于 \(p\) 的表达式： \[ P\left(\hat{p} - z_{\alpha/2} \cdot \sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \le p \le \hat{p} + z_{\alpha/2} \cdot \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\right) \approx 1-\alpha \] 因此，总体比例 \(p\) 的 \(1-\alpha\) 置信区间为： \[ \left[ \hat{p} - z_{\alpha/2} \cdot \sqrt{\frac{\hat{p}(1-\hat{p})}{n}},\ \hat{p} + z_{\alpha/2} \cdot \sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \right] \]

例8-8：

调查某社区1000名选民（n=1000），其中580人支持某项政策。

求支持率p的95%置信区间。

1. 确定置信水平：

设定置信水平 \(1-\alpha = 0.95\)，则显著性水平 \(\alpha = 0.05\)。

2. 构造枢轴量： 选择样本比例 \(\hat{p}\) 作为总体支持率 \(p\) 的点估计。在样本量足够大时，其枢轴量为Z统计量： \[ Z = \frac{\hat{p} - p}{\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}} \sim N(0, 1) \] （注：由于总体比例 \(p\) 未知，计算标准误差时通常用样本比例 \(\hat{p}\) 代替）

3. 构建对称置信区间并代入数值： 根据标准正态分布的双侧分位数 \(z_{\alpha/2}\)，有： \[ P\left(-z_{\alpha/2} \le \frac{\hat{p} - p}{\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}} \le z_{\alpha/2}\right) = 1-\alpha \] 将不等式等价变形为关于 \(p\) 的表达式，并代入已知数值 \(\hat{p} = \frac{580}{1000} = 0.58\)，\(n = 1000\)，以及 \(\alpha = 0.05\) 对应的临界值 \(z_{0.025} = 1.96\)：

\[ \begin{aligned} P\left(\hat{p} - z_{\alpha/2} \cdot \sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \le p \le \hat{p} + z_{\alpha/2} \cdot \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\right) = 1-\alpha \\ \left[ 0.58 - 1.96 \times \sqrt{\frac{0.58 \times (1-0.58)}{1000}}, \quad 0.58 + 1.96 \times \sqrt{\frac{0.58 \times 0.42}{1000}} \right] \\ \end{aligned} \]

因此，总体支持率 \(p\) 的 \(1-\alpha\) 置信区间为： \[ \left[ 0.549, \quad 0.611 \right] \]

n <- 1000
x <- 580  # 支持人数
p_hat <- x / n
alpha <- 0.05

# 计算标准误和临界值
standard_error <- sqrt(p_hat * (1 - p_hat) / n)
z_critical <- qnorm(1 - alpha/2)

# 计算置信区间
ci_lower <- p_hat - z_critical * standard_error
ci_upper <- p_hat + z_critical * standard_error

{
  cat(sprintf("总体支持率p的95%%置信区间为: [%.3f, %.3f]\n", ci_lower, ci_upper))
  cat(sprintf("即支持率在%.1f%%到%.1f%%之间\n", ci_lower*100, ci_upper*100))
}

## 总体支持率p的95%置信区间为: [0.549, 0.611]
## 即支持率在54.9%到61.1%之间

d) 总体方差的区间估计

前提条件：总体服从正态分布 \(X \sim N(\mu, \sigma^2)\)。

1. 确定置信水平：设定置信水平为 \(1-\alpha\)。

2. 构造枢轴量：选择样本方差 \(S^2\) 作为点估计。构造的枢轴量为： \[ \chi^2 = \frac{(n-1)S^2}{\sigma^2} \sim \chi^2(n-1) \] 该统计量服从自由度为 \(n-1\) 的卡方分布，其分布不依赖于未知参数 \(\sigma^2\)。

3. 构建置信区间： 卡方分布非对称，但仍将 \(\alpha\) 平分到两侧，找到两个分位数 \(\chi^2_{1-\alpha/2}(n-1)\) 和 \(\chi^2_{\alpha/2}(n-1)\)，使得： \[ P\left(\chi^2_{1-\alpha/2}(n-1) \le \frac{(n-1)S^2}{\sigma^2} \le \chi^2_{\alpha/2}(n-1)\right) = 1-\alpha \] 将不等式等价变形为关于 \(\sigma^2\) 的表达式： \[ P\left(\frac{(n-1)S^2}{\chi^2_{\alpha/2}(n-1)} \le \sigma^2 \le \frac{(n-1)S^2}{\chi^2_{1-\alpha/2}(n-1)}\right) = 1-\alpha \] 因此，总体方差 \(\sigma^2\) 的 \(1-\alpha\) 置信区间为： \[ \left[ \frac{(n-1)S^2}{\chi^2_{\alpha/2}(n-1)},\ \frac{(n-1)S^2}{\chi^2_{1-\alpha/2}(n-1)} \right] \]

例8-9：

研究学院开会时长（分钟）的变异性。假设时长服从正态分布，抽取20次会议记录时长，记录如下： [85, 120, 65, 110, 30, 95, 75, 105, 45, 115, 50, 80, 100, 70, 90, 108, 55, 85, 118, 42]

求开会时长方差和标准差的95%置信区间。

1. 确定置信水平：

设定置信水平 \(1-\alpha = 0.95\)，则显著性水平 \(\alpha = 0.05\)。

2. 构造枢轴量：

选择样本方差 \(S^2\) 作为总体方差 \(\sigma^2\) 的点估计。其枢轴量为： \[ \chi^2 = \frac{(n-1)S^2}{\sigma^2} \sim \chi^2(n-1) \]

3. 构建对称置信区间并代入数值： 根据卡方分布的上侧分位数 \(\chi^2_{\alpha/2}(n-1)\) 和 \(\chi^2_{1-\alpha/2}(n-1)\)，有： \[ P\left(\chi^2_{1-\alpha/2} \le \frac{(n-1)S^2}{\sigma^2} \le \chi^2_{\alpha/2}\right) = 1-\alpha \] 将不等式等价变形为关于 \(\sigma^2\) 的表达式。首先计算样本统计量：\(n = 20\)，样本均值 \(\bar{X} = 82.15\)，样本方差 \(S^2 \approx 752.871\)，\(\chi^2_{0.025}(19) \approx 32.852\)，\(\chi^2_{0.975}(19) \approx 8.907\)，代入数值计算方差的置信区间：

\[ \begin{aligned} P\left(\frac{(n-1)S^2}{\chi^2_{\alpha/2}} \le \sigma^2 \le \frac{(n-1)S^2}{\chi^2_{1-\alpha/2}}\right) &= 1-\alpha \\ \left[ \frac{19 \times 752.871}{32.852}, \quad \frac{19 \times 752.871}{8.907} \right] & \\ \end{aligned} \] 因此，总体方差 \(\sigma^2\) 的 95% 置信区间为： \[ \left[ 435.420, \quad 1606.077 \right] \]

进而，总体标准差 \(\sigma\) 的 95% 置信区间为上述区间的平方根： \[ \left[ \sqrt{435.420}, \quad \sqrt{1606.077} \right] \approx \left[ 20.867, \quad 40.076 \right] \]

# 模拟数据
session_length <- c(85, 120, 65, 110, 30, 95, 75, 105, 45, 115,
                   50, 80, 100, 70, 90, 108, 55, 85, 118, 42)

n <- length(session_length)
alpha <- 0.05
sample_var <- var(session_length)

# 计算卡方分布临界值
chi2_lower <- qchisq(1 - alpha/2, df = n-1)
chi2_upper <- qchisq(alpha/2, df = n-1)

# 计算方差置信区间
var_ci_lower <- (n-1) * sample_var / chi2_lower
var_ci_upper <- (n-1) * sample_var / chi2_upper

# 计算标准差置信区间
sd_ci_lower <- sqrt(var_ci_lower)
sd_ci_upper <- sqrt(var_ci_upper)

{
  cat(sprintf("总体方差σ²的95%%置信区间: [%.3f, %.3f]\n", var_ci_lower, var_ci_upper))
  cat(sprintf("总体标准差σ的95%%置信区间: [%.3f, %.3f] 分钟\n", sd_ci_lower, sd_ci_upper))
}

## 总体方差σ²的95%置信区间: [435.420, 1606.077]
## 总体标准差σ的95%置信区间: [20.867, 40.076] 分钟

e) 两总体方差比的区间估计

适用条件：两总体分别服从正态分布 \(X \sim N(\mu_1, \sigma_1^2)\)， \(Y \sim N(\mu_2, \sigma_2^2)\)，且相互独立。

1. 确定置信水平：设定置信水平为 \(1-\alpha\)。
2. 构造枢轴量：选择两样本方差的比值 \(S_1^2 / S_2^2\) 作为点估计。构造的枢轴量为： \[ F = \frac{S_1^2 / \sigma_1^2}{S_2^2 / \sigma_2^2} = \frac{S_1^2 / S_2^2}{\sigma_1^2 / \sigma_2^2} \sim F(n_1-1, n_2-1) \] 该统计量服从第一自由度为 \(n_1-1\)、第二自由度为 \(n_2-1\) 的 F 分布。
3. 构建置信区间： F 分布非对称，将 \(\alpha\) 平分到两侧，利用 F 分布的性质 \(F_{1-\alpha/2}(n_1-1, n_2-1) = 1 / F_{\alpha/2}(n_2-1, n_1-1)\)，有： \[ P\left(F_{1-\alpha/2}(n_1-1, n_2-1) \le \frac{S_1^2 / S_2^2}{\sigma_1^2 / \sigma_2^2} \le F_{\alpha/2}(n_1-1, n_2-1)\right) = 1-\alpha \] 将不等式等价变形为关于方差比 \(\frac{\sigma_1^2}{\sigma_2^2}\) 的表达式： \[ P\left(\frac{S_1^2 / S_2^2}{F_{\alpha/2}(n_1-1, n_2-1)} \le \frac{\sigma_1^2}{\sigma_2^2} \le \frac{S_1^2 / S_2^2}{F_{1-\alpha/2}(n_1-1, n_2-1)}\right) = 1-\alpha \] 因此，两总体方差比 \(\frac{\sigma_1^2}{\sigma_2^2}\) 的 \(1-\alpha\) 置信区间为： \[ \left[ \frac{S_1^2 / S_2^2}{F_{\alpha/2}(n_1-1, n_2-1)},\ \frac{S_1^2 / S_2^2}{F_{1-\alpha/2}(n_1-1, n_2-1)} \right] \] 或等价地写为： \[ \left[ \left(\frac{S_2^2}{S_1^2}\right) \cdot F_{1-\alpha/2}(n_2-1, n_1-1),\ \left(\frac{S_1^2}{S_2^2}\right) \cdot F_{\alpha/2}(n_1-1, n_2-1) \right] \]

例8-9：

假设两个班级的成绩都服从正态分布: 班级1使用传统教学法，30名同学的成绩为： [71, 73, 87, 76, 76, 89, 79, 65, 70, 71, 85, 78, 78, 76, 71, 89, 79, 59, 81, 71, 66, 73, 67, 69, 70, 62, 82, 76, 66, 85] 班级2使用新教学法，25名同学成绩为： [83, 74, 89, 89, 88, 86, 85, 77, 74, 73, 70, 76, 63, 104, 92, 65, 73, 72, 87, 77, 81, 78, 77, 94, 75]

求两个总体方差比σ₁²/σ₂²的95%置信区间。

1. 确定置信水平：

设定置信水平 \(1-\alpha = 0.95\)，则显著性水平 \(\alpha = 0.05\)。

2. 构造枢轴量：

两个样本的样本方差分别为 \(S_1^2\) 和 \(S_2^2\)。其枢轴量为 F 统计量： \[ F = \frac{S_1^2 / \sigma_1^2}{S_2^2 / \sigma_2^2} = \frac{S_1^2}{S_2^2} \cdot \frac{\sigma_2^2}{\sigma_1^2} \sim F(n_1-1, n_2-1) \]

3. 构建对称置信区间并代入数值：

计算样本统计量：

班级 1：样本量 \(n_1 = 30\)，样本均值 \(\bar{X}_1 = 74.667\) ，样本方差\(S_1^2 \approx 61.402\)。

班级 2：样本量 \(n_2 = 25\)，样本均值 \(\bar{X}_2 = 79.040\)，样本方差\(S_2^2 \approx 90.909\)。

临界值分别为：\(F_{\alpha/2}(29, 24) = F_{0.025}(29, 24) \approx 2.218\)，\(F_{1-\alpha/2}(29, 24) = F_{0.975}(29, 24) = 1 / F_{0.025}(24, 29) \approx 0.457\)。

代入置信区间公式： \[ \left[ \frac{S_1^2 / S_2^2}{F_{\alpha/2}(n_1-1, n_2-1)},\ \frac{S_1^2 / S_2^2}{F_{1-\alpha/2}(n_1-1, n_2-1)} \right] \]

计算置信区间：

\[ [ 0.304, \quad 1.455] \]

# 班级A：传统教学法，成绩变异性较小
class_a_scores <- c(71, 73, 87, 76, 76, 89, 79, 65, 70, 71, 85, 78, 78, 76, 71, 89, 79, 59, 81, 71, 66, 73, 67, 69, 70, 62, 82, 76, 66, 85)
# 班级B：新教学法，成绩变异性较大  
class_b_scores <- c(83, 74, 89, 89, 88, 86, 85, 77, 74, 73, 70, 76, 63, 104, 92, 65, 73, 72, 87, 77, 81, 78, 77, 94, 75)

# 计算样本方差
s2_a <- var(class_a_scores)
s2_b <- var(class_b_scores)
n_a <- length(class_a_scores)
n_b <- length(class_b_scores)

# 方差比点估计
variance_ratio <- s2_a / s2_b

# 计算F分布的临界值
alpha <- 0.05
f_lower <- qf(alpha/2, df1=n_a-1, df2=n_b-1)    # 下侧分位数
f_upper <- qf(1-alpha/2, df1=n_a-1, df2=n_b-1)  # 上侧分位数

# 计算方差比的置信区间
ci_lower <- variance_ratio / f_upper
ci_upper <- variance_ratio / f_lower

cat(sprintf("两个总体方差比σ₁²/σ₂²的95%%置信区间: [%.3f, %.3f]\n", ci_lower, ci_upper))

## 两个总体方差比σ₁²/σ₂²的95%置信区间: [0.305, 1.455]

典型参数区间估计方法总结

待估参数	条件	枢轴量	使用分布
单个总体均值 (μ)	总体服从正态分布，方差 σ² 已知	$Z = $	标准正态分布 (Z分布)
单个总体均值 (μ)	总体服从正态分布，方差 σ² 未知	$T = $	t分布 (自由度 df = n-1)
单个总体比例 (p)	大样本情况 (如 $np > 5, n(1-p) > 5 $)	$Z = $ (近似)	标准正态分布 (Z分布) (近似)
单个总体方差 (σ²)	总体服从正态分布	$^2 = $	卡方分布 (χ²分布) (自由度 df = n-1)
两总体方差比 (σ₁²/σ₂²)	两总体独立且服从正态分布	$F = = $	F分布 (自由度 df₁ = n₁-1, df₂ = n₂-1)

8.3.3 其它区间估计方法

前两节我们介绍了经典的区间估计方法——枢轴量法。该方法理论严密、结果精确，尤其适用于正态总体等标准模型。然而，枢轴量法的应用前提是能够构造出分布已知且不依赖于未知参数的枢轴量，这一条件在实际问题中往往难以满足。为此，统计学家发展出了多种更为灵活的替代方法。本节简要介绍两种常用的区间估计方法：大样本法和自助法。

a) 大样本法

大样本法基于大样本理论，特别是中心极限定理。该定理指出，当样本容量足够大时，许多统计量的抽样分布会趋近于正态分布，而不依赖于总体分布的具体形式。因此，即使无法获得精确的枢轴量，我们仍可以借助渐近分布构造近似的置信区间。这种方法适用性广泛，尤其适用于样本量较大、总体分布未知或形式复杂的场景。尽管其结果为近似区间，但在实际应用中具有良好的实用性和可接受性。

例8-10：

调查200名员工工作满意度得分（1-10分），分布未知且非正态。具体如下（数据量较大，采样模拟程序给出）

set.seed(789) satisfaction <- c(rnorm(150, mean=7, sd=1.5), rnorm(50, mean=4, sd=2)) #大部分员工满意，部分员工不满意

请用大样本法求中位数的95%置信区间。

本例中 \(n=200\)，基于中心极限定理（CLT），当样本量足够大（通常 \(n>30\)）时，即便总体分布的具体形式未知（如本例中的混合正态分布），样本中位数 \(\tilde{X}\) 的抽样分布仍近似服从正态分布。由此，利用正态分布分位数公式构建置信区间：\([ \tilde{X} - z_{\alpha/2} \cdot SE(\tilde{X}), \tilde{X} + z_{\alpha/2} \cdot SE(\tilde{X}) ]\)。其中，在大样本条件下，中位数的标准误可通过样本标准差近似估算，公式为：\(SE(\tilde{X}) \approx 1.253 \cdot \frac{S}{\sqrt{n}}\)。

利用R语言计算如下：

set.seed(789)
satisfaction <- c(rnorm(150, mean=7, sd=1.5),  # 大部分员工满意
                  rnorm(50, mean=4, sd=2))     # 部分员工不满意

# 大样本中位数置信区间（基于正态近似）
n <- length(satisfaction)
median_val <- median(satisfaction)

# 计算中位数的标准误（大样本近似公式）
se_median <- 1.253 * sd(satisfaction) / sqrt(n)

z_critical <- qnorm(1 - 0.05/2)
ci_lower <- median_val - z_critical * se_median
ci_upper <- median_val + z_critical * se_median

cat(sprintf("中位数的95%%大样本置信区间: [%.2f, %.2f]\n", ci_lower, ci_upper))

## 中位数的95%大样本置信区间: [6.11, 6.86]

b) 自助法（Bootstrap）

自助法是一种基于计算机模拟的现代统计方法，其核心思想是通过对原始样本进行有放回的重复抽样，模拟出统计量的经验分布。具体而言，自助法从样本中生成大量“自助样本”，并计算每个样本的统计量值，从而构建出该统计量的经验分布。基于这一分布，可以直接计算置信区间，如使用百分位数法确定区间上下限。自助法的优势在于不依赖于总体分布假设和复杂的理论推导，适用于小样本、分布未知或统计量复杂的场合，具有较强的通用性和灵活性。

例8-11：

研究教育年限与年收入的关系。两组数据如下：

set.seed(321) n <- 150 education <- rnorm(n, mean=16, sd=3) income <- 20000 + 5000*education + 2000*rnorm(n) + 100*(education-16)^2

由于数据分布复杂，不适合传统的参数估计方法。请用自助法估计两者相关系数的95%置信区间。

在本例中，收入与教育年限之间存在非线性关系，且误差项并非标准的正态分布。这种复杂的分布特征使得经典的参数估计方法（如基于正态假设的 t 检验）在构建置信区间时可能失效。因此，我们采用自助法。我们将该 \(n=150\) 的观测样本视为“虚拟总体”，通过有放回的重复抽样构建大量的自助样本，计算每个样本的 Pearson 相关系数。利用这些模拟得出的相关系数的经验分布，我们可以直接选取其 2.5% 和 97.5% 分位数（百分位数法）作为置信区间的上下限，从而在不依赖总体分布具体形式的前提下，准确估计两者相关系数的 95% 置信区间。

我们利用R程序实现如下：

library(ggplot2)  
library(boot)     # 用于自助法

# 模拟教育年限和年收入数据
set.seed(321)
n <- 150
education <- rnorm(n, mean=16, sd=3)  # 教育年限
# 收入与教育年限相关但存在非线性
income <- 20000 + 5000*education + 2000*rnorm(n) + 100*(education-16)^2

# 创建数据框
data <- data.frame(education, income)

# 定义相关系数统计量函数
cor_stat <- function(data, indices) {
  sample_data <- data[indices, ]
  return(cor(sample_data$education, sample_data$income))
}

# 执行自助法
set.seed(123)
boot_results <- boot(data, cor_stat, R = 9999)  # R=重复抽样次数

# 计算95%置信区间（百分位数法）
boot_ci <- boot.ci(boot_results, conf = 0.95, type = "perc")

{
cat("教育年限与年收入的相关系数:\n")
cat(sprintf("点估计: %.3f\n", boot_results$t0))
cat("95%自助法置信区间:[",boot_ci$percent[4],",",boot_ci$percent[5],"]")
}

## 教育年限与年收入的相关系数:
## 点估计: 0.985
## 95%自助法置信区间:[ 0.9815737 , 0.9888549 ]

# 可视化自助分布
hist(boot_results$t, breaks=30, main="Bootstrap相关系数分布",
     xlab="相关系数", col="lightblue")
abline(v=boot_ci$percent[4:5], col="red", lwd=2, lty=2)
abline(v=boot_results$t0, col="blue", lwd=2)
legend("topright", legend=c("点估计", "95%置信区间"), 
       col=c("blue", "red"), lwd=2)