欢迎

Book in early development. Planned release in 202X.

上帝在掷骰子吗?

图 0.1: 上帝在掷骰子吗?

本书风格

可以说,点估计、区间估计、假设检验、统计功效是每一个学数理统计的学生都绕不过去的坎,离开学校从事数据相关的工作,它们仍然是必备的工具。所以,本书会覆盖相关内容,但是和高校的教材最大的区别是更加注重它们之间的区别和联系,毕竟每一个统计概念都是经过了千锤百炼,而我们的主流教材始终如一地遵循的一个基本套路,就是突然给出一大堆定义、命题或定理,紧接着冗长的证明过程,然后给出一些难以找到实际应用背景的例子。三板斧抡完后就是给学生布置大量的习题,这种教学方式无论对于立志从事理论工作的还是将来投身于工业界的学生都是不合适的。

极大似然估计最初由德国数据学家 Gauss 于 1821 年提出,但未得到重视,后来, R. A. Fisher 在 1922 年再次提出极大似然的思想,探讨了它的性质,使它得到广泛的研究和应用。(茆诗松, 程依明, and 濮晓龙 2006)

这是国内某著名数理统计教材在极大似然估计开篇第一段的内容,后面是各种定义、定理、公式推导。教材简短一句话,这里面有很多信息值得发散,一个数学家提出了统计学领域极其重要的一个核心思想,他是在研究什么的时候提出了这个想法,为什么后来没有得到重视,整整 100 年以后,Fisher 又是怎么提出这一思想的呢?他做了什么使得这个思想被广泛接受和应用?虽然这可能有点离题,但是读者可以获得很多别的启迪,要知道统计领域核心概念的形成绝不是一蹴而就的,这一点也绝不局限于统计科学,任何一门科学都是这样的,比如物理学之于光的波粒二象性。历史上,各门各派的学者历经多年的思想碰撞才最终沉淀出现在的结晶。笔者认为,学校要想培养出有原创理论创新的人才,在对待前辈的成果上,我们要不吝笔墨和口水,传道不等于满堂灌和刷分机,用寥寥数节课或者数页纸来梳理学者们几十年乃至上百年的智慧结晶是非常值得的,我们甚至可以从当时的社会、人文去剖析。

下面以区间估计为例,希望能为传道做一点事情。区间估计的意义是解决点估计可靠性问题,它用置信系数解决了对估计结果的信心问题,弥补了点估计的不足。置信系数是最大的置信水平。非常欣赏有人在收集关于统计学历史的材料,读者不妨去看看 https://github.com/sctyner/history_of_statistics

1934 年 C. J. Clopper 和 E. S. Pearson 给出二项分布 \(B(n, p)\) 参数 \(p\) 的置信带 (Clopper and Pearson 1934),图 0.2 提炼了文章的主要结果。

给定置信系数 \(1- \alpha = 0.95\) 和样本量 \(n = 10\) 的情况下,二项分布参数 \(p\) 的置信带。样本量为 10,正面朝上的次数为 2,置信水平为 0.95 的情况下,参数 \(p\) 的精确区间估计为 \((p_1, p_2) = (0.03, 0.55)\)。

图 0.2: 给定置信系数 \(1- \alpha = 0.95\) 和样本量 \(n = 10\) 的情况下,二项分布参数 \(p\) 的置信带。样本量为 10,正面朝上的次数为 2,置信水平为 0.95 的情况下,参数 \(p\) 的精确区间估计为 \((p_1, p_2) = (0.03, 0.55)\)

区间半径这么长,区间估计的意义何在?增加样本量可以使得半径更短,那么至少应该有多少样本量才可以让估计变得有意义呢?就是说用估计比不用估计更好呢?答案是 39 个,留给读者思考一下为什么?读者可能已经注意到,置信带是关于点 \((5, 0.5)\) 中心对称的,这又是为什么,并且两头窄中间胖,像个酒桶?

Base R 提供的 uniroot() 函数只能求取一元非线性方程的一个根,而 rootSolve 包提供的 uniroot.all() 函数可以求取所有的根。在给定分位点下,我们需要满足方程的最小的概率值。

Base R 提供的 binom.test() 函数可以精确计算置信区间,而 prop.test() 函数可近似计算置信区间。

# 近似计算 Wilson 区间
prop.test(x = 2, n = 10, p = 0.95, conf.level = 0.95, correct = TRUE)
## Warning in prop.test(x = 2, n = 10, p = 0.95, conf.level = 0.95,
## correct = TRUE): Chi-squared approximation may be incorrect
## 
## 	1-sample proportions test with continuity correction
## 
## data:  2 out of 10, null probability 0.95
## X-squared = 103, df = 1, p-value <2e-16
## alternative hypothesis: true p is not equal to 0.95
## 95 percent confidence interval:
##  0.03543 0.55782
## sample estimates:
##   p 
## 0.2
# 精确计算
binom.test(x = 2, n = 10, p = 0.95, conf.level = 0.95)
## 
## 	Exact binomial test
## 
## data:  2 and 10
## number of successes = 2, number of trials = 10, p-value =
## 2e-09
## alternative hypothesis: true probability of success is not equal to 0.95
## 95 percent confidence interval:
##  0.02521 0.55610
## sample estimates:
## probability of success 
##                    0.2

二项分布 \(B(n,\theta)\) 成功概率 \(\theta\)

library(ggplot2)
ggplot(data.frame(x = c(0, 1)), aes(x)) +
  stat_function(
    fun = pbinom, geom = "path", 
    args = list(size = 10, q = 2),
    color = "gray70", alpha = .8 # 颜色最浅
  ) +
  stat_function(
    fun = pbinom, geom = "path",
    args = list(size = 10, q = 4),
    color = "gray50", alpha = .8
  ) +
  stat_function(
    fun = pbinom, geom = "path",
    args = list(size = 10, q = 6),
    color = "gray30", alpha = .8
  ) +
  labs(x = expression(theta), y = expression(p[theta]), 
       title = "pbinom() with fixed sample size = 10") +
  annotate("text", label = "q = 2", x = 0.32, y = 0.50, colour = "gray70") +
  annotate("text", label = "q = 4", x = 0.50, y = 0.50, colour = "gray50") +
  annotate("text", label = "q = 6", x = 0.70, y = 0.50, colour = "gray30") +
  theme_minimal(base_size = 16)

实际达到的置信度水平随真实的未知参数值和样本量的变化而剧烈波动,这意味着这种参数估计方法在实际应用中不可靠、真实场景中参数真值是永远未知的,样本量是可控的,并且是可以变化的。根本原因在于这类分布是离散的,比如这里的二项分布。当数据 \(x\) 是离散的情况,置信区间的端点\(\ell(x)\)\(u(x)\) 也是离散的。这种缺陷是无法避免的,清晰的置信区间和离散的数据之间存在无法调和的冲突

覆盖概率 \(P_{\theta}(X = x)\) 和参数真值 \(\theta\) 的关系 (Brown, Cai, and DasGupta 2001; Geyer and Meeden 2005)

prop.test(x = 2, n = 10, p = 0.95, conf.level = 0.95, correct = TRUE)
## Warning in prop.test(x = 2, n = 10, p = 0.95, conf.level = 0.95,
## correct = TRUE): Chi-squared approximation may be incorrect
## 
## 	1-sample proportions test with continuity correction
## 
## data:  2 out of 10, null probability 0.95
## X-squared = 103, df = 1, p-value <2e-16
## alternative hypothesis: true p is not equal to 0.95
## 95 percent confidence interval:
##  0.03543 0.55782
## sample estimates:
##   p 
## 0.2
set.seed(2020)
rbinom(1, size = 30, prob = 0.2) # 得到观测值 7 
## [1] 7
7 + qnorm(1-0.95/2)*sqrt(0.2*0.8/30)
## [1] 7.005
prop.test(x = 7, n = 30, p = 0.2, conf.level = 0.95, correct = TRUE) # 得到观测值 7 对应的区间估计
## 
## 	1-sample proportions test with continuity correction
## 
## data:  7 out of 30, null probability 0.2
## X-squared = 0.052, df = 1, p-value = 0.8
## alternative hypothesis: true p is not equal to 0.2
## 95 percent confidence interval:
##  0.1064 0.4270
## sample estimates:
##      p 
## 0.2333
pbinom(7, size = 30, prob = 0.2, lower.tail = TRUE)
## [1] 0.7608
# 计算分位点
qbinom(p = 0.95, size = 30, prob = 0.2, lower.tail = TRUE)
## [1] 10
binom.test(x = 10, n = 30, p = 0.2, conf.level = 0.95)
## 
## 	Exact binomial test
## 
## data:  10 and 30
## number of successes = 10, number of trials = 30, p-value =
## 0.1
## alternative hypothesis: true probability of success is not equal to 0.2
## 95 percent confidence interval:
##  0.1729 0.5281
## sample estimates:
## probability of success 
##                 0.3333
prop.test(x = 5, n = 10, p = 0.2, conf.level = 0.95, correct = TRUE) # 得到观测值 7 对应的区间估计
## Warning in prop.test(x = 5, n = 10, p = 0.2, conf.level = 0.95,
## correct = TRUE): Chi-squared approximation may be incorrect
## 
## 	1-sample proportions test with continuity correction
## 
## data:  5 out of 10, null probability 0.2
## X-squared = 3.9, df = 1, p-value = 0.05
## alternative hypothesis: true p is not equal to 0.2
## 95 percent confidence interval:
##  0.2014 0.7986
## sample estimates:
##   p 
## 0.5

\[ \sum_{x = C_1 + 1}^{C_2 -1} x \binom{n}{x} p^{x}(1-p)^{n-x} = np \sum_{x = C_1 + 1}^{C_2 -1} \binom{n -1}{x -1} p^{x -1}(1-p)^{(n -1)-(x-1)} \]

n = 30 
c2 = 20 
c1 = 10
p = 0.2
n * p * (pbinom(c2 - 2, n - 1, p) - pbinom(c1 - 1, n - 1, p))
## [1] 0.2956

多重比较与检验

多重比较 p.adjust() 函数 Adjust P-values for Multiple Comparisons 单因素多重比较 oneway.test()

set.seed(123)
x <- rnorm(50, mean = c(rep(0, 25), rep(3, 25)))
p <- 2 * pnorm(sort(-abs(x)))
# ?p.adjust
round(p, 3)
##  [1] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## [11] 0.001 0.002 0.003 0.004 0.005 0.007 0.007 0.009 0.009 0.011
## [21] 0.021 0.049 0.061 0.063 0.074 0.083 0.086 0.119 0.189 0.206
## [31] 0.221 0.286 0.305 0.466 0.483 0.492 0.532 0.575 0.578 0.619
## [41] 0.636 0.645 0.656 0.689 0.719 0.818 0.827 0.897 0.912 0.944
# round(p.adjust(p), 3)
# round(p.adjust(p, "BH"), 3)

混合正态分布的参数估计

两个二元正态分布的碰撞,点的密度估计值代表概率密度值,

散点图:faithful 数据集

图 0.3: 散点图:faithful 数据集

统计检验,决策风险,显著性水平

Charles J. Geyer 的文章 Fuzzy and Randomized Confidence Intervals and P-Values (Geyer and Meeden 2005) 文章中的图 1 名义覆盖概率的计算见 (Blyth and Hutchinson 1960)

统计计算与软件实现

binom.test(x = 2, n = 10, p = 0.95)
## 
## 	Exact binomial test
## 
## data:  2 and 10
## number of successes = 2, number of trials = 10, p-value =
## 2e-09
## alternative hypothesis: true probability of success is not equal to 0.95
## 95 percent confidence interval:
##  0.02521 0.55610
## sample estimates:
## probability of success 
##                    0.2
binom.test(x = 1, n = 10, p = 0.95)
## 
## 	Exact binomial test
## 
## data:  1 and 10
## number of successes = 1, number of trials = 10, p-value =
## 2e-11
## alternative hypothesis: true probability of success is not equal to 0.95
## 95 percent confidence interval:
##  0.002529 0.445016
## sample estimates:
## probability of success 
##                    0.1

置信区间

# 置信系数 p 0.95
binom.test(x = 1, n = 10, p = 0.95)$conf.int
## [1] 0.002529 0.445016
## attr(,"conf.level")
## [1] 0.95
binom.test(x = 1, n = 10, p = 0.95)$p.value
## [1] 1.865e-11

如果此时计算的 P 值是 0.05 反推真实的成功概率 \(\theta\) 为多少

q 分位点

# 计算覆盖概率
pbinom(q = 1, size = 10, prob = 0.95, lower.tail = TRUE)
## [1] 1.865e-11
pbinom(q = 0:10/10, size = 10, prob = 0.95, lower.tail = TRUE)
##  [1] 9.766e-14 9.766e-14 9.766e-14 9.766e-14 9.766e-14 9.766e-14
##  [7] 9.766e-14 9.766e-14 9.766e-14 9.766e-14 1.865e-11
round(pbinom(0:4, 10, 1 / 2), 5)
## [1] 0.00098 0.01074 0.05469 0.17188 0.37695
# power.prop.test() # 比例检验的功效

现代统计建模的三重境界

现代统计建模的三重境界:修改自 2019 年冬季 Bradley Efron 的课程笔记(第一部分) http://statweb.stanford.edu/~ckirby/brad/STATS305B_Part-1_corrected-2.pdf

图 0.4: 现代统计建模的三重境界:修改自 2019 年冬季 Bradley Efron 的课程笔记(第一部分) http://statweb.stanford.edu/~ckirby/brad/STATS305B_Part-1_corrected-2.pdf

本书定位

学习本书需要读者具备基本的概率、统计知识,比如上过一学期的概率论和数理统计学,也需要读者接触过编程知识,比如至少上过一学期的 C 语言、Python 语言或 Matlab 语言。了解基本的线性代数,比如矩阵的加、减、乘、逆四则运算、线性子空间、矩阵的 LU、SVD、Eigen 等分解。

内容概要

6 章介绍数据导入导出 data.table 用于 csv 文件,而 openxlsx 用于 xlsx 文件,第 5 章介绍数据操作 Base R 和 data.table,第 7 章介绍数据可视化 ggplot2

致谢

特别感谢 XX,还有很多人通过提交 PR 或 Issues 的方式参与了本书的创作过程,没有这一点一滴的持续改进,本书不会达到现在的样子,所以我将他们列在致谢名单中,详见表 0.1,人名按照提交量(commit 的个数)倒序排列。

表 0.1: 致谢名单
贡献者 提交量
Yihui Xie 1

黄湘云
于北京

授权

本书采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 许可,请君自重,别没事儿拿去传个什么新浪爱问、百度文库以及 XX 经济论坛,项目中代码使用 MIT 协议 开源

参考文献

Blyth, Colin R., and David W. Hutchinson. 1960. “Table of Neyman-Shortest Unbiased Confidence Intervals for the Binomial Parameter.” Biometrika 47 (3/4): 381–91. http://www.jstor.org/stable/2333308.

Brown, Lawrence D., T. Tony Cai, and Anirban DasGupta. 2001. “Interval Estimation for a Binomial Proportion.” Statistical Science, no. 2: 101–33. https://projecteuclid.org/euclid.ss/1009213286.

Clopper, C. J., and E. S. Pearson. 1934. “The Use of Confidence or Fiducial Limits Illustrated in the Case of the Binomial.” Biometrika 26 (4): 404–13. https://doi.org/10.1093/biomet/26.4.404.

Geyer, Charles J., and Glen D. Meeden. 2005. “Fuzzy and Randomized Confidence Intervals and P -Values.” Statistical Science 20 (4): 358–66. https://doi.org/10.1214/088342305000000340.

茆诗松, 程依明, and 濮晓龙. 2006. 高等数理统计. 2nd ed. 北京: 高等教育出版社.