专栏首页生信了R-概率统计与模拟(二)

R-概率统计与模拟(二)

本文继续介绍一些和概率统计相关的模拟。

前文《R-概率统计与模拟》介绍了一些用 R 进行概率模拟的实验,本文继续上次的工作,并在此过程中回顾一些相关的概率统计知识。

一共五题:

  • 对pi值的估计(蒙特卡洛模拟经典示例)
  • 贝叶斯公式练习
  • 多个独立并符合同一个正态分布的变量的平方和符合卡方分布
  • 多个独立且符合同一个柯西分布的变量的平均值仍符合柯西分布
  • 马尔可夫链练习

题目一:对pi值的估计(蒙特卡洛模拟经典示例)

圆周率 pi 的值大家很早就知道了。其实我们可以用模拟的方法来自行估算。方法是:

按照上述方法,笔者写了一段随机产生点的代码:

产生的点如下:

题目二:贝叶斯公式练习

从理论上计算:

用于模拟的代码是:

# Q2: 贝叶斯定理
obayes <- function() {
  x <- runif(1, min=0, max=1/3)
  y <- runif(1, min=0, max=1)
  ifelse(y < sin(pi * x), 1, 0)
}

sbayes <- function(n) {
  set.seed(SEED)
  res <- replicate(n, obayes())
  mean(res)
}

SEED <- 123
N.bayes <- 1000000
sbayes(N.bayes)   # simulative value: 0.477231

1.5 / pi  # the theoretical value: 0.4774648

可以看出,当模拟的重复结果达到一百万次时,模拟的结果和实际值很接近。

另外,笔者还尝试另一种计算方式:

但是没有成功。如果有朋友找到计算 G(z) 的方法,希望能指导一下。

题目三:多个独立并符合同一个正态分布的变量的平方和符合卡方分布

正如标题所说,模拟的任务就是看看多个独立并符合同一个正态分布的变量的平方和是否符合卡方分布。我们会尝试不同的变量数目进行模拟。

所以,这次会模拟出“多个独立并符合同一个正态分布的变量的平方和”这个变量的 c.d.f. 曲线。

用于模拟的代码:

# Q3: 模拟多个独立同分布(正态分布)变量的平方和(cdf)
ochi <- function(n) {
  sum(rnorm(n) ^ 2)
}

schi <- function(n) {
  set.seed(SEED)
  res <- replicate(N.chi, ochi(n))
  ecdf(res)
}

SEED <- 123
N.chi <- 10000
varNum.chi <- c(2, 3, 5, 10)
lineCol <- c("red", "blue", "green", "black")
res.chi <- sapply(varNum.chi, schi)
plot(res.chi[[1]], do.points=F, verticals=T, col=lineCol[1],
     main="Simulation for quadratic sum of n i.i.d. variables\n(Normal distribution)")
for (i in 2:length(res.chi)) {
  lines(res.chi[[i]], do.points=F, verticals=T, col=lineCol[i])
}
legend("bottomright", legend=paste0("n=", varNum.chi), col=lineCol, lty=1, bty="n")

模拟出的 c.d.f. 曲线:

我们可以看看卡方分布理论上的 c.d.f. 曲线:

(图片引自网页https://wiki.mbalib.com/wiki/%E5%8D%A1%E6%96%B9%E5%88%86%E5%B8%83)

可以看出,模拟曲线和理论曲线很相像。

题目四:多个独立且符合同一个柯西分布的变量的平均值仍符合柯西分布

如同题目三,这次是看柯西分布的平均值。同样是看 c.d.f. 曲线:

模拟的代码:

# Q4: 模拟柯西分布的均值(cdf)
ocauchy <- function(n, g) {
  mean(rcauchy(n, 0, g))
}

scauchy <- function(n, g) {
  set.seed(SEED)
  res <- replicate(N.cauchy, ocauchy(n, g))
  ecdf(res)
}

SEED <- 123
N.cauchy <- 10000
g <- 1    # Let Xi ~ C(g, 0).
x.max <- 5
varNum.cauchy <- c(2, 3, 5, 10)
lineCol <- c("red", "blue", "green", "black")
res.cauchy <- sapply(varNum.cauchy, scauchy, g=g)
plot(res.cauchy[[1]], do.points=F, verticals=T, col=lineCol[1], 
     xlim=c(-x.max, x.max), 
     main="Simulation for mean of n i.i.d. variables\n(Cauchy distribution)")
for (i in 2:length(varNum.cauchy))
  lines(res.cauchy[[i]], do.points=F, 
        verticals=T, col=lineCol[i], xlim=c(-x.max, x.max))
legend("bottomright", legend=paste0("n=", varNum.cauchy), col=lineCol, lty=1, bty="n")

模拟的曲线:

可以看出:

  • 模拟出的曲线非常像 arctan 函数的曲线,而理论上一个柯西分布的 c.d.f. 就是一个 arctan 函数。
  • 当 n 值(i.i.d. 变量的数目)变化时,它们的模拟曲线是重合的,说明它们不仅都符合柯西分布,而且是同分布的。

题目五:马尔可夫链练习

马尔可夫链大家应该都很熟悉了。我们用一个小题目回顾一下:

用于模拟的代码是:

# Q5: 马尔科夫链(用四宫格模拟即可)
MT <- matrix(c(
  0.5, 0.2, 0.2, 0.1,
  0.2, 0.5, 0.18, 0.12,
  0.15, 0.25, 0.5, 0.1,
  0.22, 0.18, 0.1, 0.5
), byrow = T, nrow=4)

powerOfMat <- function(M, n) {
  R <- diag(1, nrow=nrow(M))
  for (i in 1:n)
    R <- R %*% M
  R
}

omarkov <- function(mt, k, start, end, nstep) {
  res <- start
  for (i in 1:nstep)
    res <- sample(1:k, 1, replace = T, prob=mt[res, ])
  ifelse(res == end, 1, 0)
}

smarkov <- function(mt, k, start, end, nstep, N) {
  set.seed(SEED)
  res <- replicate(N, omarkov(mt, k, start, end, nstep))
  mean(res)
}

SEED <- 123
startState <- 1
endState <- 3
nstate <- 4
nstep <- 10
N.markov <- 100000
smarkov(MT, nstate, startState, endState, nstep, N.markov)  # simulative value:0.2512
powerOfMat(MT, nstep)[startState, endState]  # theoretical value:0.2519698

可以看出,当模拟的重复次数到达十万次时,模拟值很接近理论值了。

小结

从前文到本文,我们共通过八个小题目回顾了一些概率统计相关的知识,并尝试用 R 去做一些模拟,希望能对朋友们有所帮助。如果文中有任何错误,期望大家能指正!

本文分享自微信公众号 - 生信了(gh_ed36a29a9a9d),作者:hxj7

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-10-10

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R-概率统计与模拟

    正如笔者在前文《公众号一岁啦》中所说,近期在复习概率统计相关的知识。机缘巧合,笔者遇到了几个比较有意思的题目,和朋友们分享一下:

    一只羊
  • 算法(三)列举所有k-mer组合

    比如,“ATGC”的所有1-mer是:’A’, ‘T’, ‘G’, ‘C’。共4^1=4种组合。

    一只羊
  • 生信(一)对BED文件进行排序

    在处理NGS数据时,经常要对BED文件进行排序。假设BED文件长这样,分隔符是’\t’:

    一只羊
  • 通过python为ZABBIX发告警邮件

    把 以上代码保存为  send.py  放到 /etc/zabbix/alertscripts 目录下

    py3study
  • zabbix通过python脚本发告警邮件

    python脚本为敏捷开发脚本,在zabbix监控也起到重要作用,以下是使用python脚本发送告警邮件配置方法。

    菲宇
  • useState使用和原理

    但是函数组件没有实例,也没有状态。函数组件使用状态需要使用 useState 钩子。

    Qiang
  • 了解的CAP和BASE等理论

    CAP,BASE和最终一致性是NoSQL数据库存在的三大基石。而五分钟法则是内存数据存储的理论依据。这个是一切的源头。

    sunsky
  • 《深入浅出Node.js》:Node异步编程解决方案 之 async函数

    关于async函数,需要明确它是generator函数的语法糖,即将生成器函数的*换成async关键字,将yield关键字换成await关键字。使用async函...

    前端_AWhile
  • React -- 组件间通信

    分为三种类型的通信关系: 1、父组件向子组件通信 2、子组件向父组件通信 3、没有嵌套关系的组件之间的通信 父组件向子组件通信 父组件通过子组件的prop...

    前朝楚水
  • 小程序在父组件执行子组件方法,可适用于下拉刷新上拉加载之后执行子组件方法

    当父组件引用了子组件的时候,会遇到父组件执行子组件的方法,比如下拉刷新上拉加载等事件只有在页面中才能检测到,但是获取数据的方法在子组件,这时就可以执行子组件方法...

    蓓蕾心晴

扫码关注云+社区

领取腾讯云代金券