--
title: R語言中dnorm, pnorm, qnorm與rnorm以及隨機數(shù)
date: 2018-09-07 12:02:00
type: "tags"
tags:
- 隨機數(shù)
categories: - 生物統(tǒng)計
mathjax: true
前言
在R語言中，與正態(tài)分布（或者說其它分布）有關(guān)的函數(shù)有四個碟绑，分別為dnorm,pnorm,qnorm和rnorm呕臂，其中，dnorm表示密度函數(shù)，pnorm表示分布函數(shù)，qnorm表示分位數(shù)函數(shù)，rnorm表示生成隨機數(shù)的函數(shù)狸驳。在R中與之類似的函數(shù)還有很多，具體的可以通過help(Distributions)命令去查看缩赛，對于分位數(shù)或百分位數(shù)的一些介紹可以看這篇筆記《分位數(shù)及其應(yīng)用》耙箍，關(guān)于正態(tài)分布的知識可以看這篇筆記《正態(tài)分布筆記》。
現(xiàn)在這篇筆記就介紹一下這些函數(shù)的區(qū)別酥馍。
R中的隨機數(shù)背景
R提供了多種隨機數(shù)生成器(random number generators, RNG)辩昆，默認采用的是Mersenne twister方法產(chǎn)生的隨機數(shù)，該方法是由Makoto Matsumoto和Takuji Nishimura于1997年提出來的旨袒，其循環(huán)周期是汁针。R里面還提供了了Wichmann-Hill、Marsaglia-Multicarry砚尽、Super-Duper施无、Knuth-TAOCP-2002、Knuth-TAOCP和L'Ecuyer-CMRG等幾種隨機數(shù)生成方法必孤，可以通過
RNGkind()函數(shù)進行更改猾骡，例如，如果要改為WIchmann-Hill方法，就使用如下語句：
RNGkind(kind="Wich")
set.seed()隨機數(shù)種子
在R中使用隨機數(shù)函數(shù)卓练，例如rnorm()函數(shù)來生成的隨機數(shù)是不一樣的，有時我們在做模擬時购啄，為了比較不同的方法襟企，就需要生成的隨機數(shù)都一樣，即重復(fù)生成相同的隨機數(shù)狮含，此時就可以使用set.seed()來設(shè)置隨機數(shù)種子，其參數(shù)為整數(shù)，如下所示：
> 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表示正態(tài)貧青柄，這個函數(shù)是正態(tài)分布的概率密度(probability density)函數(shù)。
正態(tài)分布的公式如下所示：
給定x映胁，μ和σ后木羹，dnorm()這個函數(shù)返回的就是會返回上面的這個公式的值，這個值就是Z-score解孙，如果是標準正態(tài)分布坑填，那么上述的公式就變成了這個樣子，如下所示：
現(xiàn)在看一個案例弛姜，如下所示：
dnorm(0,mean=0,sd=1)
# 這個是標準正態(tài)分布函數(shù)
# [1] 0.3989423
dnorm(0,mean=0,sd=1)由于是標準正態(tài)分布函數(shù)的概率密度脐瑰，這個命令其實可以直接寫為dnorm(0)即可，如下所示：
dnorm(0)
# [1] 0.3989423
再看一個非標準正態(tài)分布的案例廷臼，如下所示：
dnorm(2, mean = 5, sd = 3)
# [1] 0.08065691
雖然在dnorm()中苍在，x是一個概率密度函數(shù)(PDF,Probability Density Function)的獨立變量，但它也能看作是一組經(jīng)過Z轉(zhuǎn)換后的一組變量荠商，現(xiàn)在我們看一下使用dnorm來繪制一個正態(tài)分布的概率密度函數(shù)曲線寂恬，如下所示：
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
現(xiàn)在使用dnorm()函數(shù)計算一下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
現(xiàn)在繪圖结啼，如下所示：
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")
從上面的結(jié)果可以看出掠剑，在每個Z-score處，dnorm可以繪制出這個Z-score對應(yīng)的正態(tài)分布的pdf的高度郊愧。
pnorm
pnorm函數(shù)中的p表示Probability朴译，它的功能是，在正態(tài)分布的PDF曲線上属铁，返回從負無窮到q的積分眠寿，其中這個q指的是一個Z-score。現(xiàn)在我們大概就可以猜測出pnorm(0)的值是0.5焦蘑，因為在標準正態(tài)分布曲線上盯拱，當(dāng)Z-score等于0時，這個點正好在標準正態(tài)分布曲線的正中間，那么從負無窮到0之間的曲線面積就是整個標準正態(tài)分布曲線下面積的一半狡逢，如下所示：
pnorm(0)
# pnorm()默認的參數(shù)與dnorm()一樣宁舰，都是標準正態(tài)分布，即平均數(shù)為0奢浑，標準差為1的正態(tài)分布
# [1] 0.5
pnorm函數(shù)還能使用lower.tail參數(shù)蛮艰，如果lower.tail設(shè)置為FALSE，那么pnorm()函數(shù)返回的積分就是從q到正無窮區(qū)間的PDF下的曲線面積雀彼，因此我們就知道了壤蚜，pnorm(q)與1-pnorm(q,lower.tail=FALSE)的結(jié)果是一樣的，如下所示：
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
在計算機出現(xiàn)之前的時代里徊哑，統(tǒng)計學(xué)家們使用正態(tài)分布進行統(tǒng)計時袜刷，通常是要查正態(tài)分布表的，但是莺丑，在計算機時代著蟹，通常都不使用正態(tài)分布表了，在R中窒盐，pnorm()這個函數(shù)完全可以取代正態(tài)分布表了草则，現(xiàn)在我們使用一個Z-scores的向量來計算一下相應(yīng)的累積概率，如下所示：
# 此處還使用前面生成的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 = "標準正態(tài)分布的CDF曲線",
xlab= "分位數(shù)(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))
以上就是標準正態(tài)分布的累積分布函數(shù)(CDF,Cumulative Distribution Function)曲線蟹漓。
qnorm
簡單來說炕横，qnorm是正態(tài)分布累積分布函數(shù)(CDF,Cumulative Distribution Function）的反函數(shù)，也就是說它可以視為pnorm的反函數(shù)葡粒，這里的q指的是quantile份殿，即分位數(shù)。
使用qnorm這個函數(shù)可以回答這個問題：正態(tài)分布中的第p個分位數(shù)的Z-score是多少嗽交？
現(xiàn)在我們來計算一下卿嘲，在正態(tài)分布分布中，第50百分位數(shù)的Z-score是多少夫壁，如下所示：
qnorm(0.5)
# [1] 0
# 這里計算的就是在標準正態(tài)分布中拾枣，第50百分位數(shù)的Z-score是多少
pnorm(0)
# [1] 0.5
# 這里計算的就是在標準正態(tài)分布中，當(dāng)Z-score是0時盒让，它對應(yīng)的曲線下面積是多少梅肤，也就是對應(yīng)的是哪個百分位數(shù)
再來看一個案例：在正態(tài)分布中，第96個百分位的Z-score是多少邑茄，如下所示：
qnorm(.96)
# [1] 1.750686
再來看一個案例：在正態(tài)分布中姨蝴，第99個百分位的Z-score是多少，如下所示：
qnorm(0.99)
# [1] 2.326348
再來看一下pnorm()這個函數(shù)肺缕，如下所示：
pnorm(qnorm(0))
# [1] 0
pnorm(2.326348)
# [1] 0.99
從上面我們可以看到左医，pnorm這個函數(shù)的功能是授帕，我們知道某個Z-score是多少，它位于哪個分位數(shù)上浮梢。
接著我們進一步舉例來說明一下qnorm和pnorm的具體功能跛十，如下所示：
oldpar <- par()
par(mfrow=c(1,2))
# 設(shè)置繪圖頁面，頁面布局是一行兩列
quantiles <- seq(0, 1, by = .05)
# 以5%為步長秕硝，生成0到1的百分數(shù)
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)
# 計算每個百分位數(shù)對應(yīng)的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
現(xiàn)在進行繪圖偶器，如下所示：
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 = "標準正態(tài)分布的CDF曲線",
xlab= "分位數(shù)(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()函數(shù)的功能用于生成一組符合正態(tài)分布的隨機數(shù)，在學(xué)習(xí)各種統(tǒng)計學(xué)方法時缝裤，rnorm這個函數(shù)應(yīng)該是最常用的，它的參數(shù)有n,mean颊郎，sd憋飞，其中n表示生成的隨機數(shù)，mean與sd分別表示正態(tài)分布的均值與標準差姆吭，現(xiàn)在舉個例子榛做，如下所示：
set.seed(1000)
# 設(shè)定隨身數(shù)種子
rnorm(5)
# 生成5個服從標準正態(tài)分布的隨機數(shù)
# [1] -0.44577826 -1.20585657 0.04112631 0.63938841 -0.78655436
n10 <- rnorm(10, mean = 70, sd = 5);n10
# 生成10個，服從均值為70内狸，標準差為5的正態(tài)分布的隨機數(shù)
# [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的正態(tài)分布的隨機數(shù)
# [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的正態(tài)分布的隨機數(shù)
# [1] 66.14695 70.16345 69.54554 76.15484 74.54789 71.78985 67.85345 73.85163 67.58083
# [10] 73.98425
現(xiàn)在我們繪制一下上面的幾個向量的直方圖，看一下它們的均值是否在70附近昂灵，如下所示：
oldpar <- par()
par(mfrow=c(1,3))
# breaks為設(shè)置直方圖的寬度
hist(n10, breaks = 5)
hist(n100, breaks = 20)
hist(n10000, breaks = 100)
在R語言中避凝，生成不同分布的各種類型的函數(shù)都是以d,p,q,r開頭的，使用原理跟上面的正態(tài)分布都一樣眨补。
runif-生成均勻分布的隨機數(shù)
runif(5) # 生成 5 個介于 0 和 1 之間的均勻分布的隨機數(shù)
runif(5, 1,10) # 生成 5 個介于 0 和 10 之間的均勻分布的隨機數(shù)
rnorm(5) # 生成 5 個正態(tài)分布的隨機數(shù)管削，它們的中位數(shù)為 0，標準差為 1
rnorm(5, 3, 7) # 生成 5 個正態(tài)分布的隨機數(shù)撑螺，它們的中位數(shù)為 3含思，標準差為 7
R語言中的隨機數(shù)
sample()函數(shù)是一個用于生成隨機數(shù)的重要的核心函數(shù)，如果僅傳遞一個數(shù)值n給它甘晤，就會返回一個從1到n的自然數(shù)的排列含潘，如果傳遞是n:m就是生成從n到m的隨機數(shù)，如是是7,5安皱，則會生成5個小于7的隨機數(shù)调鬓，如下所示：
> 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
從上面的結(jié)果可以看出來，這些數(shù)字都是不同的酌伊，也就是說腾窝，sample函數(shù)默認情況下是不重復(fù)抽樣缀踪，每個值只出現(xiàn)一次，如果允許有重復(fù)抽樣虹脯，需要添加參數(shù)replace = TRUE驴娃，如下所示：
> sample(7, 10, replace = TRUE)
[1] 3 7 7 2 5 7 1 4 5 5
sample函數(shù)通常會從某些向量中隨機挑一些參數(shù)，如下所示：
> 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語言中其它分布的隨機數(shù)
| 分布 | 中文名稱 | R中的表達式 | 參數(shù) |
|---|---|---|---|
| Beta | 貝塔分布 | beta(a,b) |
shape1唇敞、shape2 |
| Binomial | 二項分布 | binom(n,p) |
size、prob |
| Cauchy | 柯西分布 | cauchy() |
location咒彤、scale |
| Chi-square | 卡方分布 | chisq(df) |
df |
| Exponential | 指數(shù)分布 | 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 | 正態(tài)分布 | norm() |
mean鞋屈、sd |
| Multivariate normal | 多元正態(tài)分布 | mvnorm() |
mean、cov |
| Poisson | 泊松分布 | pois() |
lambda |
| T | t分布 | t() |
df |
| Uniform | 均勻分布 | unif() |
min故觅、max |
| Weibull | 威布爾分布 | weibull() |
shape厂庇、scale |
| Wilcoxon | 威爾考可森分布 | wilcox() |
m、n |
上述分布函數(shù)前面加上r输吏，p权旷、q、d就可以表示相應(yīng)的目的：
-
r：生成相應(yīng)分布的隨機數(shù)贯溅； -
d：生成相應(yīng)分布的密度函數(shù)炼杖； -
p：生成相應(yīng)分布的累積概率密度函數(shù)； -
q：生成相應(yīng)分布的分位數(shù)函數(shù)盗迟。
參考資料
- Introduction to dnorm, pnorm, qnorm, and rnorm for new biostatisticians
- R數(shù)據(jù)分析——方法與案例詳解.電子工業(yè)出版社出版.方匡南 朱建平 姜葉飛.2015
