R语言笔记9: 模拟——随机数、抽样、线性模型

宏基因组2018-10-29 20:18:05

R语言基础系列:

  • 你知道R中的赋值符号箭头(<-)和等号(=)的区别吗?

  • 1数据类型(向量、数组、矩阵、 列表和数据框)

  • 2读写数据所需的主要函数、与外部环境交互

  • 3数据筛选——提取对象的子集

  • 4向量、矩阵的数学运算

  • 5控制结构

  • 6函数及作用域

  • 7认识循环函数lapply和sapply

  • 8分解数据框split和查看对象str

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.9954662

Simulating 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 8

Simulation 小结

下面敲黑板

- 使用R中的函数从指定的概率分布中取样

使用rnorm(), rpois(), rbinom()等。

常用分布包括正态分布(normal)、泊松分布(poission)、二项分布(binomial)、指数分布(exponential)、伽马分布(gamma)等。

- sample函数可用来从指定向量中随机取样
- 无论何时记得用seed生成种子

参考资料:

  1. 视频课程 R Programming by Johns Hopkins University:https://www.coursera.org/learn/r-programming/home/welcome

  2. 讲义 Programming for Data Science :https://bookdown.org/rdpeng/rprogdatascience/R

猜你喜欢

10000+:肠道细菌 人体上的生命 宝宝与猫狗 梅毒狂想曲 提DNA发Nature 实验分析谁对结果影响大  Cell微生物专刊

系列教程:微生物组入门 Biostar 微生物组  宏基因组

专业技能:生信宝典 学术图表 高分文章 不可或缺的人

一文读懂:宏基因组 寄生虫益处 进化树

必备技能:提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析:图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具:16S预测培养基 生信绘图

科研经验:云笔记  云协作 公众号

编程模板 Shell  R Perl

生物科普  生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难,我们建立了“宏基因组”专业讨论群,目前己有国内外150+ PI,1500+ 一线科研人员加入。参与讨论,获得专业解答,欢迎分享此文至朋友圈,并扫码加主编好友带你入群,务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助,首先阅读《如何优雅的提问》学习解决问题思路,仍末解决群内讨论,问题不私聊,帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战,关注“宏基因组”

Copyright © 温县电话机虚拟社区@2017