1.1 问题背景

数据来源于 2021 年邵阳市统计年鉴发布的常住人口（修订版），2010-2020年邵阳市常住人口数（单位：万人）、城市化率（单位：%）。

一座城市的常住人口数受复杂的政治、经济因素影响，具有不确定性，可以看作一个随机变量，在相对稳定的内外部环境下，不考虑时间趋势性，每年的常住人口数近似服从正态分布。简单起见，将以上这组常住人口数据视为来自某正态总体的一组样本观察值。

import pandas as pd
shaoyang = pd.DataFrame(
    {
        "年份": [2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020],
        "常住人口": [707.17, 704.01, 702.91, 698.78, 693.63, 690.33, 688.84, 687.21, 679.04, 665.39, 656.15],
        "城市化率": [32.84, 34.16, 35.77, 37.94, 39.97, 42.44, 44.71, 46.76, 48.46, 49.74, 52.16]
    }
)

1.2 提出假设

在方差 \(\sigma^2\) 未知的情况下，检验正态分布的均值 \(\mu\) 的情况。

\[ \begin{aligned} \quad H_0: \mu \geq \mu_0 \quad vs. \quad H_1: \mu < \mu_0 \end{aligned} \]

上式中 \(H_0\) 称为零假设（或称原假设），\(H_1\) 称为备择假设。不妨设 \(\mu_0 = 700\) 。

1.3 构造检验统计量

在样本量充足的情况下，根据中心极限定理，可以推导出如下结果，此处推导过程复杂，略。

\[ \begin{aligned} & \frac{\bar{x} - \mu}{\sigma / \sqrt{n}} \sim \mathcal{N}(0,1) \quad \frac{(n-1)s^2}{\sigma^2} \sim \chi^2(n-1) \\ & \mathsf{E}\{s^2\} = \sigma^2 \quad \mathsf{Var}\{s^2\} = \frac{2\sigma^4}{n-1} \end{aligned} \]

其中，样本均值 \(\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i\) ，样本方差 \(s^2 = \frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2\) 。

为了检验上面提出的问题，在总体的方差\(\sigma^2\)未知的情况下，用样本方差来替换，构造检验统计量

\[ T=\frac{\bar{x} - \mu_0}{s / \sqrt{n}} \]

该检验统计量服从自由度为 \(n-1\) 的 \(t\) 分布，此处推导过程复杂，略。

对于上述检验问题，拒绝域的形式是根据备择假设来定的。

\[ \{T < t_{\alpha}(n-1)\} \]

其中，\(\alpha\) 称为检验的显著性水平，一般取 0.05。控制犯第一类错误（拒真错误）的概率，即原假设成立的情况下，样本观察值落入拒绝域的概率，这个概率要控制在0.05 以内。只有这样，检验的结果才有说服力。说服力来源于小概率原理，即小概率事件不可能在一次实验中发生。若发生了，就不是小概率事件。有点反证法的味道！其中这个小概率的量化，一般取 0.05，事件发生的概率小于或等于 0.05，即认为此事件为小概率事件。

1.4 计算检验统计量的值和 P 值

t 分布的数学公式是非常复杂的，常见的统计软件中内置了有关的计算函数。t 分布的概率密度函数如下：

\[ f(x; n) = \frac{\Gamma\left(\frac{n+1}{2}\right)}{\sqrt{n\pi} \Gamma\left(\frac{n}{2}\right)} \left(1 + \frac{x^2}{n}\right)^{-\frac{n+1}{2}} \]

其中，\(n\) 是自由度。t 分布的累积分布函数如下：

\[ F(x; n) = \int_{-\infty}^{x} f(w; n) dw \]

# 样本均值
shaoyang["常住人口"].mean()

## np.float64(688.4963636363636)

# 样本标准差
shaoyang["常住人口"].std()

## np.float64(16.13687409179812)

# 检验统计量
T = (shaoyang["常住人口"].mean() - 700) / (shaoyang["常住人口"].std()/len(shaoyang) ** 0.5)
T

## np.float64(-2.3643516907814464)

# 临界值
import scipy.stats as stats
# 计算 alpha 处的分位数
ppf = stats.t.ppf(0.05, df=len(shaoyang)-1)  # 5% 分位数
ppf

## np.float64(-1.8124611228116767)

注意到检验统计量 T 的值确实小于临界值，样本观察值落在拒绝域中。说明，要拒绝原假设，即原正态总体的均值是小于 700 的。

# 计算 P 值
stats.t.cdf(T, df=len(shaoyang)-1)

## np.float64(0.019828842503383064)

其实，scipy 模块的统计部分有一些假设检验的函数。下面直接调用单样本 T 检验函数stats.ttest_1samp来计算 T 值和 P 值。

t_stat, p_value = stats.ttest_1samp(a=shaoyang["常住人口"], popmean=700, alternative='less')
# T 检验统计量
print(t_stat)

## -2.364351690781446

# P 值
print(p_value)

## 0.019828842503383095

计算结果和前面手搓的一样。邵阳市常住人口数总体上大于 700 万吗？根据现有的数据，这个问题的答案是否定的。

2.1 置信区间

现在的问题是根据这组样本观察值，构造一个区间，使得正态总体的均值/真值落入该区间的可能性/概率为 95% 。换句话说，我们要找一个区间，这个区间覆盖总体的真值的可能性/概率不低于 95% 。置信区间覆盖总体的真值的概率即覆盖概率，这个覆盖概率要大于 95% ，或者接近置信水平。

res = stats.ttest_1samp(a=shaoyang["常住人口"], popmean=700, alternative='less')
# 置信区间
res.confidence_interval(confidence_level=0.95)

## ConfidenceInterval(low=np.float64(-inf), high=np.float64(697.3148037023432))

补充在 R 软件中的代码实现，如下：

x <- c(
  707.17, 704.01, 702.91, 698.78, 693.63,
  690.33, 688.84, 687.21, 679.04, 665.39, 656.15
)
sd(x)

## [1] 16.14

t.test(x, mu = 700, alternative = "less")

## 
## 	One Sample t-test
## 
## data:  x
## t = -2.4, df = 10, p-value = 0.02
## alternative hypothesis: true mean is less than 700
## 95 percent confidence interval:
##   -Inf 697.3
## sample estimates:
## mean of x 
##     688.5

计算结果与前面完全一致。做假设检验的同时，可以获得参数真值的置信区间。

2.2 功效水平

一个检验的功效函数（Power function）是样本观察值落在拒绝域的概率。功效是用来衡量检验统计量的检验效果的，功效越大，表示检验统计量越灵敏，实际表现是用更少的样本量就可以对原假设问题作出结论。在业务中，减少样本量意味着减少时间和实验成本。在 Python 模块 scipy 的统计部分中有功效计算函数 stats.power，归在重抽样和蒙特卡罗方法栏目下。

在这个单样本的 t 检验问题下，检验的对象是正态分布的均值参数 \(\mu\)，功效函数记为 \(g(\mu)\) ，它与显著性水平 \(\alpha\) （也是犯第一类错误的概率）和犯第二类错误的概率 \(\beta\) 的关系。

\[ \begin{aligned} g(\mu) = \left\{ \begin{array}{l} \alpha(\mu) & H_0 ~\mathrm{成立} \\ 1-\beta(\mu) & H_1 ~\mathrm{成立} \end{array} \right. \end{aligned} \]

\[ \begin{aligned} g(\mu) &= P_{\mu}\{\bar{x} < \mu_0\} \\ &= P_{\mu}\{\frac{\bar{x}-\mu}{s/\sqrt{n}} < \frac{ \mu_0-\mu}{s/\sqrt{n}}\} \\ &= P_{\mu}\{t < \frac{ \mu_0-\mu}{s/\sqrt{n}} \} \\ &= \int_{-\infty}^{\frac{ \mu_0-\mu}{s/\sqrt{n}}} f(x;n-1)dx \end{aligned} \] 其中，将 \(\frac{\bar{x}-\mu}{s/\sqrt{n}}\) 看作自由度为 \(n-1\) 的 \(t\) 分布的随机变量，而 \(f(x)\) 是相应的 \(t\) 分布的概率密度函数。当将样本均值 \(\bar{x}\) 代入 \(\mu\) 时，Power 值就是前面计算的 P 值。

import numpy as np
rng = np.random.default_rng(seed=20262026)
# 总体（常住人口）的分布是正态分布
# 样本均值和样本标准差来替代总体均值和总体标准差
# 将从这个分布中生成随机数
rvs = lambda size: rng.normal(loc=shaoyang["常住人口"].mean(), scale=shaoyang["常住人口"].std(), size=size)
# 单样本 t 检验
test = stats.ttest_1samp
# 样本量
n_obs = len(shaoyang["常住人口"])
# 默认样本量是 n_resamples=10000
res = stats.power(test=test, rvs=rvs, n_observations=n_obs, significance=0.05, kwargs={'popmean': 700, 'alternative': 'less'})
# 输出结果
res

## PowerResult(power=np.float64(0.7109), pvalues=array([0.01801163, 0.08560967, 0.00089626, ..., 0.03789451, 0.02300109,
##        0.00028921], shape=(10000,)))

# 打印 power 值
res.power

## np.float64(0.7109)

# P 值小于 0.05 的比例表示功效
np.mean(res.pvalues < 0.05)

## np.float64(0.7109)

from statsmodels.stats.power import TTestPower
# Cohen's d，即均值差除以标准差
effect_size = (shaoyang["常住人口"].mean()-700)/shaoyang["常住人口"].std()
effect_size

## np.float64(-0.7128788573422227)

# 显著性水平
alpha=0.05
res = TTestPower().power(
  effect_size=effect_size,
  nobs=n_obs,
  alpha=alpha,
  alternative='smaller' # 注意与 effect_size、备择假设的符号保持一致
)
print(res)

## 0.710660157004632

power.t.test(
  n = length(x), delta = abs(mean(x) - 700),
  sd = sd(x), sig.level = 0.05,
  type = "one.sample",
  alternative = "one.sided" # 这里没有方向可选，delta 必须为正，计算结果才对
)

## 
##      One-sample t test power calculation 
## 
##               n = 11
##           delta = 11.5
##              sd = 16.14
##       sig.level = 0.05
##           power = 0.7107
##     alternative = one.sided

# 前面 t 检验等价的功效计算
library(pwr)
pwr.t.test(
  d = (mean(x)-700)/sd(x), # Cohen's d
  n = length(x),
  sig.level = 0.05,
  type = "one.sample",
  alternative = "less" # 注意与备择假设、参数 d 的方向保持一致
)

## 
##      One-sample t test power calculation 
## 
##               n = 11
##               d = -0.7129
##       sig.level = 0.05
##           power = 0.7107
##     alternative = less

pwr.t.test(
  d = (700-mean(x))/sd(x), # Cohen's d
  n = length(x),
  sig.level = 0.05,
  type = "one.sample",
  alternative = "greater" # 同上
)

## 
##      One-sample t test power calculation 
## 
##               n = 11
##               d = 0.7129
##       sig.level = 0.05
##           power = 0.7107
##     alternative = greater