R-概率统计与模拟（二）

一只羊

发布于 2019-10-11 17:38:40

7710

发布于 2019-10-11 17:38:40

文章被收录于专栏：生信了生信了

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

前文《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

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