R语言系列第二期（番外篇）：R先生教你统计概率与分布

本文首发于 ”百味科研芝士“ 微信公众号，转载请注明：百味科研芝士，Focus科研人的百味需求

点击蓝字关注我们

还记得我们在系列2开始的时候为大家介绍的几个特别的函数吗，rnorm(),dnorm()...？如果你忘记了，详情点击：

R语言系列第二期：②R编程、函数、数据输入等功能

在这个部分，我们会给大家介绍一下概率与分布的统计知识以及R中包含的关于随机抽样和处理理论分布的函数，这个部分的内容同时也是下一个系列描述性统计和图表的基础。

1

随机抽样

> sample(1:6,3)

[1] 2 4 6

#Tips：sample()函数第一个参数是被抽取的值向量，第二个参数是抽样次数，就是样本量大小。每次使用相同的参数，结果也是不同的，因为每次抽样都是随机的。这里还需要注意一点的是，sample()函数默认是无放回的抽样方式，也就是说，样本不会包含同一个数字两次，并且样本量不能超过被抽取的值向量长度。如果想有放回的抽样，那么需要加上参数replace=TRUE。

> sample(1:6,7)

Error in sample.int(length(x), size, replace, prob) :

'replace = FALSE'，因此不能取比总体要大的样本

> sample(1:6,7,replace=T)

[1] 5 4 3 5 2 3 4

同时，R可以实现对称抽样和不对称抽样，比如，多次投硬币是有放回的抽样，实际中我们通常认为正反两面的概率是一致的，都是50%。我们可以模拟10次投掷硬币：

> sample(c("Z","F"),10,replace=T)

 [1] "Z" "F" "F" "F" "F" "F" "F" "Z" "Z" "F"

公平地掷硬币的话，出现正反面的概率应该是一样的，但是随机事件的思想告诉我们实际结果不一定就像我们预想的对称情形。它同样适用于其他情形，比如说一种外科手术的成功率是90%，我们可以通过sample()函数中prob参数来模拟这种不相等概率的数据：

>sample(c("success","failure"),10,replace=T,prob=c(0.9,0.1))

 [1] "success" "failure" "failure" "success" "success" "success" "failure"

 [8] "success" "success" "success"

当然，这个函数不仅仅局限与两种结局的事情，也可以是多分类的事件，比如：

> sample(c("优","良","及格","不及格"),30,replace=T,prob=c(0.4,0.4,0.15,0.05))

[1] "良" "良" "不及格" "良" "及格""良" "及格" 

 [8] "优" "良" "优" "不及格" "优" "良" "良"   

[15] "优" "优" "优" "良" "良" "优" "优"   

[22] "良" "良" "良" "良" "优" "及格" "良"   

[29] "良" "良"   

#Tips：需要注意的是，事件数目要和概率数目保持一致，不然会提示错误，另外，你可以让概率和不为1，只是剩余概率指示的事件不会出现在结果里。

不过，这不是产生样本集合的最好方法，因为我们在实际中往往不会关注每个个体的情况，而是关注总体汇总的情况。后文会给大家介绍。

2

概率计算以及排列组合

先前我们提到的掷色子的例子 > sample(1:6,3) ，得到第一个样本的概率为1/6，第二个是1/5，以此类推。那么给定一个样本的整体概率就是1/（6*5*4）。在R中，使用prod()函数，可以用于计算数字向量的乘积，即排列A63。那么概率为：

> 1/prod(6:4)

[1] 0.008333333

但是有的时候给定结果的数字不一定要求按照特定顺序排列，也许可以打乱顺序，就像乐透彩票一样，我们只是想知道抽取的每个数字我们是否含有，而不关注它是第几个，这种情况下，我们就需要在上面概率的基础上乘上出现期望组合的次数，第一个数字有3种可能，第二个有2种，最后一个1种，从而可能数字是3*2*1，也就是3！，那么赢得乐透大奖的概率就是：

> prod(3:1)/prod(6:4)

[1] 0.05

有其他的方法可以得到同样的结果，比如说我们能够计算出所有可能的数量是：组合C63,概率就是它的倒数。在R里，可以使用choose()函数来解决组合问题，这个概率就可以写成：

> 1/choose(6,3)

[1] 0.05

3

统计知识：离散分布和连续分布

当观察一个独立重复的二项试验时，通常对每次试验的成功或失败并不感兴趣，更感兴趣的是成功或者失败的总数(概率)，显然这个数字是随机的，因为它依赖每一次随机结果，因此被称为随机变量。这个整体数据的分布就是二项分布（图示见下文）。而如果实验结果为多种不连续的可能，我们可以认为整体数据的分布为离散分布。

有些数据来自于对实质连续尺度的测量，比如温度、浓度等。实际中，它们只能被记录成有限精度的值。这种随机波动会遵循某种模式，通常会集中在某个中心值附近，这里我们不能像离散分布那样去定义每个点的概率，因为在连续分布中，任何特定值的概率为零。就像连续的函数中，每个点的积分都是零。所以这里不能用点概率，而采用的是概率密度的概念。

正态分布也是连续分布的一种理想化状态，它有标志性的钟形曲线，也是连续变量统计中常用的一种概率分布。（图示见后文）

4

R内置的分布

与建模和统计检验有关的常用标准分布都已经嵌入在R中，因此可以完全取代传统的统计表格。这里我们着重了解一下正态分布和二项分布，二者比较常见，且其他分布都遵循相同的模式。

对一个统计分布可以计算4项基本内容：

1. 密度或点概率
2. 累计概率分布函数
3. 分位数
4. 随机数

在R的所有的分布，关于上面列出的4项都对应一个相应的函数。比如对于正态分布，它们分别为dnorm，pnorm，qnorm，rnorm（分别对应密度、概率、分位数和随机数）。

Part1.密度

连续分布的密度是指得到一个接近x的值的相对可能性的度量。在一个特定区间得到一个值得概率是在相应曲线下的面积。对于离散分布，密度用点概率描述，也就是得到x值的概率。

密度函数是4中函数中最少在实际中应用的。如果要绘制正态分布钟形曲线，可以这样：

> x<-seq(-4,4,0.1)

> plot(x,dnorm(x),type="l")

#Tips：这里是字母“L”的小写，不是数字“1”。Type=有几个选项，详细情况可以参考 > ?plot 。Seq()产生从-4到4，步长为0.1的等距数值。Type=“l”使得函数在点与点之间画线而不是画出点本身来。 

创建图形的另一个方法是使用curve()函数：

> curve(dnorm(x),from=-4,to=4)

这个方法通常更加方便，但是需要满足y的值能够通过x的简单函数表达式表示出来的条件。

对于离散分布，变量只能取那些明确的值，更倾向于画竖线图。下面是n=50,p=0.4的二项分布图形：

> x<-0:50

> plot(x,dbinom(x,size=50,prob=0.4),type="h") 

在这个dbinom()函数里有3个参数。除了x（出现阳性事件的次数），还需要具体说明实验次数n和概率参数p。比如，画出投掷一枚枚硬50次出现正面的数量。其实，dnorm还有其他参数，即均值和标准差，他们分别默认0和1，因为通常我们默认的是标准正态分布。

Part2.累积分布函数

累积分布函数描述的是对一个给定分布小于或等于x的累积概率。相应的R函数按惯例以“p”(probability第一个字母)开头。

正如可以对密度作图，也可以对累积分布函数作图，但是我们更需要的是实际的数字，即我们计算的概率到底是多少。比如健康人群某项生化指标可以很好地用均值为132、标准差为13的正态分布来描述。那么，如果一个患者的检测值为160，那么：

> 1-pnorm(160,132,13)

[1] 0.01562612

> 1-pnorm(160,mean=132,sd=13)

这就表示，在这个健康人群中，只有0.0156的概率能够找到检测值是160或者更高的结果。pnorm()返回一个在给定分布下取得小于第一个参数事件的概率。

对于二项分布，同样可以计算尾部概率。20个病人每人进行2种治疗，问治疗A还是治疗B更好，结果16个病人觉得A好。问题是这是否能作为A治疗好于B的充分证据呢，还是说只是偶然发生的，如果我们设定喜欢治疗A的人数服从n=20，p=0.5的二项分布。那么累积概率是：

> pbinom(16,size=20,p=0.5)

[1] 0.9987116

然后需要用1减去这部分概率，但是这样就剔掉了数量为16的概率，我们需要的是观测到16个或者更极端的概率和，那么我们可以用15来替代：

> 1-pbinom(15,size=20,p=0.5)

[1] 0.005908966

这个是典型的单侧检验，就是说A比B好，当然我们也可以计算一下B比A好的概率（尽管根据我们收集的样本来看可能性很小），这就意味着4个或者更少的人喜欢A，那么两种情况的概率合计：

> 1-pbinom(15,size=20,p=0.5)+pbinom(4,20,0.5)

[1] 0.01181793

#Tips：你会发现这个数值正好是上面的两倍，因为这两种情况在密度函数是是对称的。另外有的时候并不一定需要写参数的关键字，只要参数值顺序正确即可。

当然这些都是复杂的计算，之后的系列我们会给大家介绍直接计算实际概率的函数比如说t.test，binom.test等，他们可以直接给出我们想要的结果。

Part3.分位数

分位数函数是累积分布函数的反函数。P-分位数是具有这样性质的一个值：得到小于等于它的概率为P。

#Tips：统计分布表几乎都是根据分位数函数结果给出的。

分位数常用于置信区间的计算，置信区间的公式是：

> xbar<-83

> sigma<-12

> n<-5

> stderr<-sigma/sqrt(n)

> stderr

[1] 5.366563

> xbar+stderr*qnorm(0.025)

[1] 72.48173

> xbar+stderr*qnorm(0.975)

[1] 93.51827

我们就可以得到μ的95%置信区间，从72.48173到93.51827。

#Tips：这里是已知总体标准差的计算方式，如果总体标准差未知，需要引入自由度。这个部分我们后续讨论。另外，正态分布是对称的，所以N0.025 = - N0.975.

> qnorm(0.025)

[1] -1.959964

> qnorm(0.975)

[1] 1.959964

Part.4 随机数字

使用随机数字产生函数非常简单，第一个参数指定产生随机数字的数量，后续参数类似于其他与相同分布有关的函数中相应位置的参数。

> rnorm(10)

[1] -0.38433961  0.07077883  1.69446851  2.13576539  0.21176494

[6]  1.24767280 -0.44994990  0.11796112 -0.98289727  0.96196080

> rnorm(10,mean=7,sd=5)

[1] -1.6260339 17.1604979 10.7228844 10.3688051  7.0181209

[6] -1.3064159  0.5518916 -3.5541634  7.1747866  3.7209794

> rbinom(10,size=20,prob=0.5)

[1] 10 11  9  9 11 12 11 10  8 12

到这里，大家应该大致明白了统计里的一些基本统计概念以及如何使用R计算一些简单的概率和位点。不过给你一个大型的样本使用这样的方法似乎很难计算，好在统计学家已经为我们设计好了相应统计方法，R中也纳入了这部分的内容，因此之后的系列会给大家介绍如何使用R语言直接计算我们需要的统计量和P值，敬请期待。

参考资料：
1.《Introductory Statistics with R Second Edition》Peter Dalgaard著
2.《R语言初学者指南》 人民邮电出版社  Brian Dennis著
3.Vicky的小笔记本《blooming for you》by Vicky