R语言笔记9: 模拟——随机数、抽样、线性模型
R语言基础系列：
Simulation
今天学习的内容是模拟 —— Simulation，在统计和一些其他的应用中很重要，所以在这里介绍一些R语言中可以做模拟的函数。
用于模拟已知概率分布的数字和变量
有些函数可以直接生成符合某种概率分布的随机数字或变量，例如：
rnorm()：指定一个均值和标准差，即可生成符合正态分布的随机数字变量rpois()：从已知平均发生率（rate）的泊松分布中生成泊松随机变量
一共有四类基本的函数和概率分布函数相关，它们的前缀分别是d, r, p以及q：
d：用来估计密度（density）r：用来产生随机数字（random）p：估计累计分布（cumulative distribution）q：估计分位数（quantile）
每种分布都有以上这四种前缀构成的函数，比如rpois(), dpois(), ppois(), qpois()等。
每种函数都有不同的参数，那正态分布的四个函数举例：
dnorm(x, mean = 0, sd = 1, log = FALSE)
## 可以求密度的对数值
pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
## 可以对概率求对数；计算该分布的左尾，如果填FALSE就是计算右尾
qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
## 可以对概率求对数；计算该分布的左尾，如果填FALSE就是计算右尾
rnorm(n, mean = 0, sd = 1)
## n是你想生成的随机变量的个数所有函数都需要我们给定均值和标准差，以此确定实际的概率分布，如果不指定，那么默认属于标准正态分布（均值为0，标准差为1）
生成最简单的服从标准正态分布的随机数：
> ## Simulate standard Normal random numbers
> x <- rnorm(10)
> x
[1] 0.01874617 -0.18425254 -1.37133055 -0.59916772 0.29454513
[6] 0.38979430 -1.20807618 -0.36367602 -1.62667268 -0.25647839修改参数：
> x <- rnorm(10, 20, 2)
> x
[1] 22.20356 21.51156 19.52353 21.97489 21.48278 20.17869 18.09011
[8] 19.60970 21.85104 20.96596
> summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
18.09 19.75 21.22 20.74 21.77 22.20设置随机数字生成种子
set.seed()函数可以用来设置随机数字生成种子（seed）。
这种方法可产生相同的随机数，也叫“伪随机数”。它怎么用？
举例
我们可以设置种子为任意整数，然后生成随机数字如下：
> set.seed(1)
> rnorm(5)
[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078
> set.seed(2)
> rnorm(5)
[1] -0.89691455 0.18484918 1.58784533 -1.13037567 -0.08025176
> set.seed(5)
> rnorm(5)
[1] -0.84085548 1.38435934 -1.25549186 0.07014277 1.71144087我们看到种子设定值不同时随机数字也不同，但是再设定相同种子就会出现完全一样的随机数字：
> set.seed(1)
> rnorm(5)
[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078这个操作能让你重复得到之前生成的随机数字，可以让我们的模拟过程可重复，所以在生成随机数字之前最好设定一个种子。
再拿泊松分布举例：
生成不同发生率的10个随机变量：
> rpois(10, 1) ## Counts with a mean of 1
[1] 0 0 1 1 2 1 1 4 1 2
> rpois(10, 2) ## Counts with a mean of 2
[1] 4 1 2 0 1 1 0 1 4 1
> rpois(10, 20) ## Counts with a mean of 20
[1] 19 19 24 23 22 24 23 20 11 22估计泊松分布的累积分布函数：
## 在平均发生率为2的泊松分布中，出现小于等于2的随机变量的概率是多少
> ppois(2,2)
[1] 0.6766764
## 在平均发生率为2的泊松分布中，出现小于等于4的随机变量的概率是多少
> ppois(4,2)
[1] 0.947347
## 在平均发生率为2的泊松分布中，出现小于等于6的随机变量的概率是多少
> ppois(6,2)
[1] 0.9954662Simulating a Linear Model
我们再来看怎样从简单的线性模型中提取随机值
如果x是单一自变量
举例：
## 首先记得先设定一个种子
> set.seed(20)
## 设置自变量x取值，属于标准正态分布
> x <- rnorm(100)
## 随机噪声e是属于标准差为2的正态分布
> e <- rnorm(100, 0, 2)
## 设置线性方程的回归系数和截距
> y <- 0.5 + 2 * x + e
## 输出概要结果
> summary(y)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-6.4084 -1.5402 0.6789 0.6893 2.9303 6.5052
## 画出图表
> plot(x,y)这就是通过回归模型模拟出的x和y的图表，可以看出明显的线性关系。
如果x是二元变量
x如果是二元变量，可以代表两种性别、两种实验处理（可以是实验组和对照组）这一类情况
使用二项分布函数rbinom()
## 重新设定种子
> set.seed(10)
## 随机取100个二元变量数据，并设置参数
> x <- rbinom(100, 1, 0.5)
## 取符合正态分布的随机变量
> e <- rnorm(100, 0, 2)
## 生成线性模型
> y <- 0.5 + 2 *x + e
## 查看结果和作图
> summary(y)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.4936 -0.1409 1.5767 1.4322 2.8397 6.9410
> plot(x, y)可以看到，x是二元变量，y依旧是连续的，符合正态分布
从复杂模型中模拟取值
假设结果y服从于一个平均值为μ的泊松分布，log(μ)服从一个截距为β0，斜率β1的线性函数，x是其中的自变量
> set.seed(1)
> x <- rnorm(100)
## 模拟线型变量log(μ)
> log.mu <- 0.5 + 0.3 * x
## 使用rpois函数来计算y，取幂
> y <- rpois(100, exp(log.mu))
## 查看结果
> summary(y)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 1.00 1.00 1.55 2.00 6.00
> plot(x, y)可以看到x越大y越大，呈线性关系
Random Sampling
这里使用的函数是sample()，可以从你给定的一组对象中，随机抽取样本
假如我们提供一个数值向量，sample()可以从中随机抽取样本。因此我们可以人为确定一个分布，指定给一个向量并从中取样
举例：
从整数1到10中取样，取出的数字不放回，如果重复去取样两次，得到的数字是不重复的，对字母取样也是同理：
> set.seed(1)
> sample(1:10, 4)
[1] 3 4 5 7
> sample(1:10, 4)
[1] 3 9 8 5
> ## Doesn't have to be numbers
> sample(letters, 5)
[1] "q" "b" "e" "x" "p"如果只提供向量，不指定任何条件，返回结果就是这些整数重新排列，数字内部无重复：
> ## Do a random permutation
> sample(1:10)
[1] 4 7 10 6 9 2 8 3 1 5
> sample(1:10)
[1] 2 3 4 1 9 5 10 8 6 7有放回的取样，设置参数replace = TRUE，因为有放回，所以可以看到数字内部有重复：
> ## Sample w/replacement
> sample(1:10, replace = TRUE)
[1] 2 9 7 8 2 8 5 9 7 8Simulation 小结
下面敲黑板
- 使用R中的函数从指定的概率分布中取样
使用rnorm(), rpois(), rbinom()等。
常用分布包括正态分布（normal）、泊松分布（poission）、二项分布（binomial）、指数分布（exponential）、伽马分布（gamma）等。
- sample函数可用来从指定向量中随机取样
- 无论何时记得用seed生成种子
参考资料：
视频课程 R Programming by Johns Hopkins University：https://www.coursera.org/learn/r-programming/home/welcome
讲义 Programming for Data Science ：https://bookdown.org/rdpeng/rprogdatascience/R
猜你喜欢
10000+：肠道细菌 人体上的生命 宝宝与猫狗 梅毒狂想曲 提DNA发Nature 实验分析谁对结果影响大 Cell微生物专刊
文献阅读 热心肠 SemanticScholar Geenmedical
16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外150+ PI，1500+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”