--
title: R语言中dnorm, pnorm, qnorm与rnorm以及随机数
date: 2018-09-07 12:02:00
type: "tags"
tags:
- 随机数
categories: - 生物统计
mathjax: true
前言
在R语言中，与正态分布（或者说其它分布）有关的函数有四个，分别为dnorm,pnorm,qnorm和rnorm，其中，dnorm表示密度函数，pnorm表示分布函数，qnorm表示分位数函数，rnorm表示生成随机数的函数。在R中与之类似的函数还有很多，具体的可以通过help(Distributions)命令去查看，对于分位数或百分位数的一些介绍可以看这篇笔记《分位数及其应用》，关于正态分布的知识可以看这篇笔记《正态分布笔记》。
现在这篇笔记就介绍一下这些函数的区别。
R中的随机数背景
R提供了多种随机数生成器(random number generators, RNG)，默认采用的是Mersenne twister方法产生的随机数，该方法是由Makoto Matsumoto和Takuji Nishimura于1997年提出来的，其循环周期是。R里面还提供了了Wichmann-Hill、Marsaglia-Multicarry、Super-Duper、Knuth-TAOCP-2002、Knuth-TAOCP和L'Ecuyer-CMRG等几种随机数生成方法，可以通过
RNGkind()函数进行更改，例如，如果要改为WIchmann-Hill方法，就使用如下语句：
RNGkind(kind="Wich")
set.seed()随机数种子
在R中使用随机数函数，例如rnorm()函数来生成的随机数是不一样的，有时我们在做模拟时，为了比较不同的方法，就需要生成的随机数都一样，即重复生成相同的随机数，此时就可以使用set.seed()来设置随机数种子，其参数为整数，如下所示：
> set.seed(1)
> runif(5)
[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
> set.seed(2)
> runif(5)
[1] 0.1848823 0.7023740 0.5733263 0.1680519 0.9438393
> set.seed(1)
> runif(5)
[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
dnorm
dnorm中的d表示density，norm表示正态贫，这个函数是正态分布的概率密度(probability density)函数。
正态分布的公式如下所示：
给定x，μ和σ后，dnorm()这个函数返回的就是会返回上面的这个公式的值，这个值就是Z-score，如果是标准正态分布，那么上述的公式就变成了这个样子，如下所示：
现在看一个案例，如下所示：
dnorm(0,mean=0,sd=1)
# 这个是标准正态分布函数
# [1] 0.3989423
dnorm(0,mean=0,sd=1)由于是标准正态分布函数的概率密度，这个命令其实可以直接写为dnorm(0)即可，如下所示：
dnorm(0)
# [1] 0.3989423
再看一个非标准正态分布的案例，如下所示：
dnorm(2, mean = 5, sd = 3)
# [1] 0.08065691
虽然在dnorm()中，x是一个概率密度函数(PDF,Probability Density Function)的独立变量，但它也能看作是一组经过Z转换后的一组变量，现在我们看一下使用dnorm来绘制一个正态分布的概率密度函数曲线，如下所示：
z_scores <- seq(-3,3,0.1)
# 生成一个Z-score的向量
z_scores
# [1] -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5
# [17] -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1
# [33] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7
# [49] 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0
现在使用dnorm()函数计算一下Z_scores的概率密度，如下所示：
dvalues <- dnorm(z_scores)
dvalues
# [1] 0.004431848 0.005952532 0.007915452 0.010420935 0.013582969 0.017528300
# [7] 0.022394530 0.028327038 0.035474593 0.043983596 0.053990967 0.065615815
# [13] 0.078950158 0.094049077 0.110920835 0.129517596 0.149727466 0.171368592
# [19] 0.194186055 0.217852177 0.241970725 0.266085250 0.289691553 0.312253933
# [25] 0.333224603 0.352065327 0.368270140 0.381387815 0.391042694 0.396952547
# [31] 0.398942280 0.396952547 0.391042694 0.381387815 0.368270140 0.352065327
# [37] 0.333224603 0.312253933 0.289691553 0.266085250 0.241970725 0.217852177
# [43] 0.194186055 0.171368592 0.149727466 0.129517596 0.110920835 0.094049077
# [49] 0.078950158 0.065615815 0.053990967 0.043983596 0.035474593 0.028327038
# [55] 0.022394530 0.017528300 0.013582969 0.010420935 0.007915452 0.005952532
# [61] 0.004431848
现在绘图，如下所示：
plot(z_scores,dvalues, # Plot where y = values and x = index of the value in the vector
type = "l", # Make it a line plot
main = "pdf of the Standard Normal",
xlab= "Z-score")
从上面的结果可以看出，在每个Z-score处，dnorm可以绘制出这个Z-score对应的正态分布的pdf的高度。
pnorm
pnorm函数中的p表示Probability，它的功能是，在正态分布的PDF曲线上，返回从负无穷到q的积分，其中这个q指的是一个Z-score。现在我们大概就可以猜测出pnorm(0)的值是0.5，因为在标准正态分布曲线上，当Z-score等于0时，这个点正好在标准正态分布曲线的正中间，那么从负无穷到0之间的曲线面积就是整个标准正态分布曲线下面积的一半，如下所示：
pnorm(0)
# pnorm()默认的参数与dnorm()一样，都是标准正态分布，即平均数为0，标准差为1的正态分布
# [1] 0.5
pnorm函数还能使用lower.tail参数，如果lower.tail设置为FALSE，那么pnorm()函数返回的积分就是从q到正无穷区间的PDF下的曲线面积，因此我们就知道了，pnorm(q)与1-pnorm(q,lower.tail=FALSE)的结果是一样的，如下所示：
pnorm(2)
# [1] 0.9772499
pnorm(2, mean = 5, sd = 3)
# [1] 0.1586553
pnorm(2, mean = 5, sd = 3, lower.tail = FALSE)
# [1] 0.8413447
1 - pnorm(2, mean = 5, sd = 3, lower.tail = FALSE)
# [1] 0.1586553
在计算机出现之前的时代里，统计学家们使用正态分布进行统计时，通常是要查正态分布表的，但是，在计算机时代，通常都不使用正态分布表了，在R中，pnorm()这个函数完全可以取代正态分布表了，现在我们使用一个Z-scores的向量来计算一下相应的累积概率，如下所示：
# 此处还使用前面生成的z_scores
z_scores
# [1] -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5
# [17] -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1
# [33] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7
# [49] 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0
pvalues <- pnorm(z_scores)
pvalues
# [1] 0.001349898 0.001865813 0.002555130 0.003466974 0.004661188 0.006209665
# [7] 0.008197536 0.010724110 0.013903448 0.017864421 0.022750132 0.028716560
# [13] 0.035930319 0.044565463 0.054799292 0.066807201 0.080756659 0.096800485
# [19] 0.115069670 0.135666061 0.158655254 0.184060125 0.211855399 0.241963652
# [25] 0.274253118 0.308537539 0.344578258 0.382088578 0.420740291 0.460172163
# [31] 0.500000000 0.539827837 0.579259709 0.617911422 0.655421742 0.691462461
# [37] 0.725746882 0.758036348 0.788144601 0.815939875 0.841344746 0.864333939
# [43] 0.884930330 0.903199515 0.919243341 0.933192799 0.945200708 0.955434537
# [49] 0.964069681 0.971283440 0.977249868 0.982135579 0.986096552 0.989275890
# [55] 0.991802464 0.993790335 0.995338812 0.996533026 0.997444870 0.998134187
# [61] 0.998650102
plot(pvalues,
xaxt = "n", # 不绘制x轴的刻度
type = "l", # 使用曲线将各点连接起来
main = "标准正态分布的CDF曲线",
xlab= "分位数(Quantiles)",
ylab="概率密度(Probability Density)")
# 以下是添加x轴的刻度
# These commands label the x-axis
axis(1, at=which(pvalues == pnorm(-2)), labels=round(pnorm(-2), 2))
axis(1, at=which(pvalues == pnorm(-1)), labels=round(pnorm(-1), 2))
axis(1, at=which(pvalues == pnorm(0)), labels=c(.5))
axis(1, at=which(pvalues == pnorm(1)), labels=round(pnorm(1), 2))
axis(1, at=which(pvalues == pnorm(2)), labels=round(pnorm(2), 2))
以上就是标准正态分布的累积分布函数(CDF,Cumulative Distribution Function)曲线。
qnorm
简单来说，qnorm是正态分布累积分布函数(CDF,Cumulative Distribution Function）的反函数，也就是说它可以视为pnorm的反函数，这里的q指的是quantile，即分位数。
使用qnorm这个函数可以回答这个问题：正态分布中的第p个分位数的Z-score是多少？
现在我们来计算一下，在正态分布分布中，第50百分位数的Z-score是多少，如下所示：
qnorm(0.5)
# [1] 0
# 这里计算的就是在标准正态分布中，第50百分位数的Z-score是多少
pnorm(0)
# [1] 0.5
# 这里计算的就是在标准正态分布中，当Z-score是0时，它对应的曲线下面积是多少，也就是对应的是哪个百分位数
再来看一个案例：在正态分布中，第96个百分位的Z-score是多少，如下所示：
qnorm(.96)
# [1] 1.750686
再来看一个案例：在正态分布中，第99个百分位的Z-score是多少，如下所示：
qnorm(0.99)
# [1] 2.326348
再来看一下pnorm()这个函数，如下所示：
pnorm(qnorm(0))
# [1] 0
pnorm(2.326348)
# [1] 0.99
从上面我们可以看到，pnorm这个函数的功能是，我们知道某个Z-score是多少，它位于哪个分位数上。
接着我们进一步举例来说明一下qnorm和pnorm的具体功能，如下所示：
oldpar <- par()
par(mfrow=c(1,2))
# 设置绘图页面，页面布局是一行两列
quantiles <- seq(0, 1, by = .05)
# 以5%为步长，生成0到1的百分数
quantiles
# [1] 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75
# [17] 0.80 0.85 0.90 0.95 1.00
qvalues <- qnorm(quantiles)
# 计算每个百分位数对应的Z-score
qvalues
# [1] -Inf -1.6448536 -1.2815516 -1.0364334 -0.8416212 -0.6744898 -0.5244005
# [8] -0.3853205 -0.2533471 -0.1256613 0.0000000 0.1256613 0.2533471 0.3853205
# [15] 0.5244005 0.6744898 0.8416212 1.0364334 1.2815516 1.6448536 Inf
现在进行绘图，如下所示：
plot(qvalues,
type = "l", # We want a line graph
xaxt = "n", # No x-axis
xlab="概率密度(Probability Density)",
ylab="Z-scores")
# Same pnorm plot from before
plot(pvalues, # Plot where y = values and x = index of the value in the vector
xaxt = "n", # Don't label the x-axis
type = "l", # Make it a line plot
main = "标准正态分布的CDF曲线",
xlab= "分位数(Quantiles)",
ylab="概率密度(Probability Density)")
# 绘制x轴的刻度
axis(1, at=which(pvalues == pnorm(-2)), labels=round(pnorm(-2), 2))
axis(1, at=which(pvalues == pnorm(-1)), labels=round(pnorm(-1), 2))
axis(1, at=which(pvalues == pnorm(0)), labels=c(.5))
axis(1, at=which(pvalues == pnorm(1)), labels=round(pnorm(1), 2))
axis(1, at=which(pvalues == pnorm(2)), labels=round(pnorm(2), 2))
rnorm
rnomr()函数的功能用于生成一组符合正态分布的随机数，在学习各种统计学方法时，rnorm这个函数应该是最常用的，它的参数有n,mean，sd，其中n表示生成的随机数，mean与sd分别表示正态分布的均值与标准差，现在举个例子，如下所示：
set.seed(1000)
# 设定随身数种子
rnorm(5)
# 生成5个服从标准正态分布的随机数
# [1] -0.44577826 -1.20585657 0.04112631 0.63938841 -0.78655436
n10 <- rnorm(10, mean = 70, sd = 5);n10
# 生成10个，服从均值为70，标准差为5的正态分布的随机数
# [1] 71.80351 60.47042 73.38245 74.64034 73.20814 75.98877 78.17864 68.84888 72.23209
# [10] 81.19844
n100 <- rnorm(100, mean = 70, sd = 5);n100[1:10]
# 生成100个，服从均值为70，标准差为5的正态分布的随机数
# [1] 62.09470 72.75453 65.42062 68.65545 76.84043 74.17830 75.40543 65.32326 70.71987
# [10] 70.81957
n10000 <- rnorm(10000, mean = 70, sd = 5);n10000[1:10]
# 生成1000个，服从均值为70，标准差为5的正态分布的随机数
# [1] 66.14695 70.16345 69.54554 76.15484 74.54789 71.78985 67.85345 73.85163 67.58083
# [10] 73.98425
现在我们绘制一下上面的几个向量的直方图，看一下它们的均值是否在70附近，如下所示：
oldpar <- par()
par(mfrow=c(1,3))
# breaks为设置直方图的宽度
hist(n10, breaks = 5)
hist(n100, breaks = 20)
hist(n10000, breaks = 100)
在R语言中，生成不同分布的各种类型的函数都是以d,p,q,r开头的，使用原理跟上面的正态分布都一样。
runif-生成均匀分布的随机数
runif(5) # 生成 5 个介于 0 和 1 之间的均匀分布的随机数
runif(5, 1,10) # 生成 5 个介于 0 和 10 之间的均匀分布的随机数
rnorm(5) # 生成 5 个正态分布的随机数，它们的中位数为 0，标准差为 1
rnorm(5, 3, 7) # 生成 5 个正态分布的随机数，它们的中位数为 3，标准差为 7
R语言中的随机数
sample()函数是一个用于生成随机数的重要的核心函数，如果仅传递一个数值n给它，就会返回一个从1到n的自然数的排列，如果传递是n:m就是生成从n到m的随机数，如是是7,5，则会生成5个小于7的随机数，如下所示：
> sample(7)
[1] 5 1 2 6 4 3 7
> sample(7:5)
[1] 5 6 7
> sample(7, 5)
[1] 6 4 1 3 2
从上面的结果可以看出来，这些数字都是不同的，也就是说，sample函数默认情况下是不重复抽样，每个值只出现一次，如果允许有重复抽样，需要添加参数replace = TRUE，如下所示：
> sample(7, 10, replace = TRUE)
[1] 3 7 7 2 5 7 1 4 5 5
sample函数通常会从某些向量中随机挑一些参数，如下所示：
> sample(colors(), 5)
[1] "honeydew4" "wheat4" "deepskyblue3" "lightgray" "royalblue1"
也可以挑日期，如下所示：
> sample(.leap.seconds, 4)
[1] "1994-07-01 08:00:00 CST" "1977-01-01 08:00:00 CST" "1982-07-01 08:00:00 CST" "1972-07-01 08:00:00 CST"
R语言中其它分布的随机数
| 分布 | 中文名称 | R中的表达式 | 参数 |
|---|---|---|---|
| Beta | 贝塔分布 | beta(a,b) |
shape1、shape2 |
| Binomial | 二项分布 | binom(n,p) |
size、prob |
| Cauchy | 柯西分布 | cauchy() |
location、scale |
| Chi-square | 卡方分布 | chisq(df) |
df |
| Exponential | 指数分布 | exp(lambda) |
rate |
| F | F分布 | f(df1,df2) |
df1、df2 |
| Gamma | 伽玛分布 | gamma() |
shape、rate |
| Geometric | 几何分布 | geom() |
prob |
| Hypergeometric | 超几何分布 | hyper() |
m、n、k |
| Logistic | 逻辑分布 | logis() |
location、scale |
| Negative binomial | 负二项分布 | nbinom() |
size、prob |
| Normal | 正态分布 | norm() |
mean、sd |
| Multivariate normal | 多元正态分布 | mvnorm() |
mean、cov |
| Poisson | 泊松分布 | pois() |
lambda |
| T | t分布 | t() |
df |
| Uniform | 均匀分布 | unif() |
min、max |
| Weibull | 威布尔分布 | weibull() |
shape、scale |
| Wilcoxon | 威尔考可森分布 | wilcox() |
m、n |
上述分布函数前面加上r，p、q、d就可以表示相应的目的：
-
r：生成相应分布的随机数； -
d：生成相应分布的密度函数； -
p：生成相应分布的累积概率密度函数； -
q：生成相应分布的分位数函数。
参考资料
- Introduction to dnorm, pnorm, qnorm, and rnorm for new biostatisticians
- R数据分析——方法与案例详解.电子工业出版社出版.方匡南 朱建平 姜叶飞.2015
