前言

Book in early development. Planned release in 202X.

上帝在掷骰子吗?

图 0.1: 上帝在掷骰子吗?

本书风格

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

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

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

下面以区间估计为例,希望能为传道做一点事情。区间估计的意义是解决点估计可靠性问题,它用置信系数解决了对估计结果的信心问题,弥补了点估计的不足。置信系数是最大的置信水平。

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 等分解。

内容概要

语言抉择

行业内可以做统计分析和建模的软件汗牛充栋,比较顶级的收费产品有 SAS 和 SPSS,在科学计算领域的 Matlab 和 Mathematica 也有相当强的统计功能,而用户基数最大的是微软 Excel,抛开微软公司的商业手段不说,Excel 的市场份额却是既成事实。 Brian D. Ripley 20 多年前的一句话很有意思,放在当下也是适用的。

Let’s not kid ourselves: the most widely used piece of software for statistics is Excel.

— Brian D. Ripley (Ripley 2002)

有鉴于 Excel 在人文、社会、经济和管理等领域的影响力,熟悉 R 语言的人把它看作超级收费版的 Excel,这实在是一点也不过分。事实上,我司就是一个很好的明证,一个在线教育类的互联网公司,各大业务部门都在使用 Excel 作为主要的数据分析工具。然而,Excel 的不足也十分突出,工作过程无法保存和重复利用,Excel 也不是数据库,数据集稍大,操作起来愈发困难,对于复杂的展示,需要借助内嵌的 VBA,由于缺乏版本控制,随着时间的推移,几乎不可维护。所以,我们还是放弃 Excel 吧,Jenny Bryan 更在 2016 年国际 R 语言大会上的直截了当地喊出了这句话1。Nathan Stephens 对 Excel 的缺陷不足做了全面的总结2

Some people familiar with R describe it as a supercharged version of Microsoft’s Excel spreadsheet software.

— Ashlee Vance3

另一方面,我们谈谈开源领域的佼佼者 — R (https://cran.r-project.org/),Python (https://www.python.org/) 和 Octave (http://www.gnu.org/software/octave/)。Python 号称万能的胶水语言,从系统运维到深度学习都有它的广泛存在,它被各大主流 Linux 系统内置,语言风格上更接近于基数庞大的开发人员,形成了强大的生态平台。 Octave 号称是可以替代 Matlab 的科学计算软件,在兼容 Matlab 的方面确实做的很不错,根据 Julia 官网给出的各大编程语言的测试 https://julialang.org/benchmarks/,性能上不能相提并论。

R 语言扩展包生态系统

图 0.5: R 语言扩展包生态系统

R 提供了丰富的图形接口,包括 Tcl/Tk , Gtk, Shiny 等,以及基于它们的衍生品 rattle(RGtk2)、Rcmdr(tcl/tk)、radiant(shiny)。更多底层介绍,见 John Chamber 的著作《Extending R》。

TikZ 在绘制示意图方面有很大优势,特别是示意图里包含数学公式,这更是 LaTeX 所擅长的方面

JASP https://jasp-stats.org 是一款免费的统计软件,源代码托管在 Github 上 https://github.com/jasp-stats/jasp-desktop,主要由阿姆斯特丹大学 E. J. Wagenmakers 教授 https://www.ejwagenmakers.com/ 领导的团队维护开发,实现了很多贝叶斯和频率统计方法,相似的图形用户界面使得 JASP 可以作为 SPSS 的替代,目前实现的功能见 https://jasp-stats.org/current-functionality/,统计方法见博客 https://www.bayesianspectacles.org/

Patrick Burns 收集整理了 R 语言中奇葩的现象,写成 The R Inferno 直译过来就是《R 之炼狱》。这些奇葩的怪现象可以看做是 R 风格的一部分,对于编程人员来说就是一些建议和技巧,参考之可以避开某些坑。 Paul E. Johnson 整理了一份真正的 R 语言建议,记录了他自己从 SAS 转换到 R 的过程中遇到的各种问题 http://pj.freefaculty.org/R/Rtips.html。Michail Tsagris 和 Manos Papadakis 也收集了 70 多条 R 编程的技巧和建议,力求以更加 R 范地将语言特性发挥到极致 (Tsagris and Papadakis 2018)。 Python 社区广泛流传着 Tim Peters 的 《Python 之禅》,它已经整合进每一版 Python 软件中,只需在 Python 控制台里执行 import this 可以获得。

  1. Beautiful is better than ugly.
  2. Explicit is better than implicit.
  3. Simple is better than complex.
  4. Complex is better than complicated.
  5. Flat is better than nested.
  6. Sparse is better than dense.
  7. Readability counts.
  8. Special cases aren’t special enough to break the rules.
  9. Although practicality beats purity.
  10. Errors should never pass silently.
  11. Unless explicitly silenced.
  12. In the face of ambiguity, refuse the temptation to guess.
  13. There should be one– and preferably only one –obvious way to do it.
  14. Although that way may not be obvious at first unless you’re Dutch.
  15. Now is better than never.
  16. Although never is often better than right now.
  17. If the implementation is hard to explain, it’s a bad idea.
  18. If the implementation is easy to explain, it may be a good idea.
  19. Namespaces are one honking great idea – let’s do more of those!

— The Zen of Python

总之,编程语言到一定境界都是殊途同归的,对美的认识也是趋同的,道理更是相通的,Python 社区的 Pandas https://github.com/pandas-dev/pandas 和 Matplotlib https://github.com/matplotlib/matplotlib 也有数据框和图形语法的影子。Pandas https://github.com/pandas-dev/pandas 明确说了要提供与 data.frame 类似的数据结构和对应统计函数等,而 Matplotlib 偷了 ggplot2 绘图样式 https://matplotlib.org/3.2.1/gallery/style_sheets/ggplot.html

获取帮助

R 社区提供了丰富的帮助资源,可以在 R 官网搜集的高频问题 https://cran.r-project.org/faqs.html 中查找,也可在线搜索 https://cran.r-project.org/search.htmlhttps://rseek.org/https://stackoverflow.com/questions/tagged/r,更多获取帮助方式见 https://www.r-project.org/help.html。在当下信息爆炸的时代,唯一不缺的就是各种学习资源:

适合入门的书籍:

与数据可视化相关:

与 shiny 相关:

与统计推断相关:

以及课程资源:

还有各类食谱:

以及中外博客:

除了 R 语言,我们还需要掌握一点和命令行相关的东西,比如 Bash 和 Makefile 等。

写作环境

书籍项目架构图

图 0.6: 书籍项目架构图

本书 R Markdown 源文件托管在 Github 仓库里,本地使用 RStudio IDE 编辑,bookdown 组织各个章节的 Rmd 文件和输出格式,使用 Git 进行版本控制。每次提交修改到 Github 上都会触发 Travis 自动编译书籍,将一系列 Rmd 文件经 knitr 调用 R 解释器执行里面的代码块,并将输出结果返回,Pandoc 将 Rmd 文件转化为 md 、 html 或者 tex 文件。若想输出 pdf 文件,还需要准备 TeX 排版环境,最后使用 Netlify 托管书籍网站,和 Travis 一起实现连续部署,使得每次修改都会同步到网站。最近一次编译时间 2020年05月24日02时24分16秒,本书用 R version 3.6.3 (2020-02-29) 编译,完整运行环境如下:

xfun::session_info(packages = c(
  "knitr", "rmarkdown", "bookdown"
), dependencies = FALSE)
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 8 (Core)
## 
## Locale:
##   LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##   LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##   LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##   LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##   LC_ADDRESS=C               LC_TELEPHONE=C            
##   LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## Package version:
##   bookdown_0.19   knitr_1.28.5    rmarkdown_2.1.3
## 
## Pandoc version: 2.7.3

借助 bookdown (Xie 2016) 可以将 Rmd 文件组织起来, rmarkdown (Allaire et al. 2020)knitr (Xie 2015) 将源文件编译成 Markdown 文件, Pandoc 将 Markdown 文件转化成 HTML 和 TeX 文件, TinyTeX (Xie 2019) 可以将 TeX 文件进一步编译成 PDF 文档,书中大量的图形在用 ggplot2 包制作 (Wickham 2016),而统计理论相关的示意图用 Base R 创作。

本书在后续章节中陆续用到新的 R 包,其安装过程不会在正文中呈现,下面以在 CentOS 8 上安装 sf 包为例介绍。首先需要安装一些系统依赖,具体安装哪些依赖参见 sf 包开发站点 https://github.com/r-spatial/sf

sudo dnf copr enable simc/stable # gdal-devel
sudo dnf install -y sqlite-devel gdal-devel \
  proj-devel geos-devel udunits2-devel

然后,在 R 命令行窗口中,执行安装命令:

install.packages('sf')

至此,安装完成。如遇本地未安装的新 R 包,可从其官方文档中找寻安装方式。如果你完全不知道自己应该安装哪些,考虑把下面的依赖都安装上

sudo dnf install -y \
  # magick
  ImageMagick-c++-devel \ 
  # pdftools
  poppler-cpp-devel \ 
  # gifski
  cargo 

最后,本书在三个位置提供网页版, 网站 Github Pages 发布最近一次在 Travis 构建成功的版本 https://xiangyunhuang.github.io/masr/,网站 Bookdown 发布本地手动创建的版本 https://bookdown.org/xiangyun/masr/ ,网站 Netlify 发布最新的开发版 https://masr.netlify.app/

记号约定

正文中的代码、函数、参数及参数值以等宽正体表示,如 data(list = c('iris', 'BOD')), 其中函数名称 data(),参数及参数值 list = c('iris', 'BOD') ,R 程序包用粗体表示,如 graphics

ruler()
----+----1----+----2----+----3----+----4----+----5----+----6----+----
123456789012345678901234567890123456789012345678901234567890123456789

Winston Chang 整理了一份 LaTeX 常用命令速查小抄 https://wch.github.io/latexsheet/latexsheet.pdf

参考文献

Allaire, JJ, Yihui Xie, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, Hadley Wickham, Joe Cheng, Winston Chang, and Richard Iannone. 2020. Rmarkdown: Dynamic Documents for R. https://github.com/rstudio/rmarkdown.

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.

Hornik, Kurt. 2020. “R FAQ: Frequently Asked Questions on R.” https://CRAN.R-project.org/doc/FAQ/R-FAQ.html.

Ripley, Brian D. 2002. “Statistical Methods Need Software: A View of Statistical Computing.” Opening Lecture Royal Statistical Society. Plymouth. https://www.stats.ox.ac.uk/~ripley/RSS2002.pdf.

Tsagris, Michail, and Manos Papadakis. 2018. “Taking R to Its Limits: 70+ Tips.” PeerJ Preprints 6: e26605v1. https://doi.org/10.7287/peerj.preprints.26605v1.

Wickham, Hadley. 2016. ggplot2: Elegant Graphics for Data Analysis. 2nd ed. New York: Springer-Verlag. https://ggplot2-book.org/.

Xie, Yihui. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

Xie, Yihui. 2016. Bookdown: Authoring Books and Technical Documents with R Markdown. Boca Raton, Florida: Chapman; Hall/CRC. https://github.com/rstudio/bookdown.

Xie, Yihui. 2019. “TinyTeX: A Lightweight, Cross-Platform, and Easy-to-Maintain Latex Distribution Based on TeX Live.” TUGboat, no. 1: 30–32. https://tug.org/TUGboat/Contents/contents40-1.html.

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