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

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

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

一共五题：

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

题目一：对pi值的估计（蒙特卡洛模拟经典示例）

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

假设有一个单位圆（半径为1），我们在它的外接矩形里随机产生一个点，那么这个点落在单位圆中的理论概率是
$$p=\frac{S_{circle}}{S_{rect}}=\frac{\pi r^2}{(2r{)}^2}=\frac{\pi }{4}.$$
其中$S_{circle}$和$S_{rect}$分别是单位圆和外接矩形的面积。所以，理论上：
$$\pi =4p.$$
那么，我们只需要模拟出这个$p$值，就可以得到$\pi$的估计值。
如果按上述随机产生点的方法，$p$的估计值$\hat{p}$是：
$$\hat{p}=\frac{npoint{s}_{circle}}{npoint{s}_{rect}}.$$
其中，$npoint{s}_{circle}$是落在单位圆里的点的数量，$npoint{s}_{rect}$是落在外接矩形里的点的数量。
得到$p$的估计值后，$\pi$的估计值是：
$$\hat{\pi }=4\hat{p}=4\frac{npoint{s}_{circle}}{npoint{s}_{rect}}.$$

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

# Q1: the value of pi
# 也是蒙特卡洛方法的经典示例
opi <- function(i) {
  pt <- runif(2, min=-1, max=1)
  c(x=pt[1], y=pt[2], isInCircle=as.numeric(sum(pt ^ 2) < 1))
}

spi <- function(n) {
  set.seed(SEED)
  res <- sapply(1:n, opi)
  t(res)
}

SEED <- 123
N.pi <- 100000
res.pi <- as.data.frame(spi(N.pi))
ratio.pi <- mean(res.pi$isInCircle)
color.pi <- ifelse(res.pi$isInCircle == 1, "green", "yellow")
plot(x=res.pi$x, y=res.pi$y, col=color.pi, 
     main=paste0("Simulation for pi with ", N.pi, " points"), xlab="x", ylab="y")
ratio.pi * (2 ^ 2)    # the value of pi

产生的点如下：

题目二：贝叶斯公式练习

假设变量$X$和$Y$互相独立，且都服从$[0,1]$上的均匀分布，并且令事件$A$和$B$为：
A = { X < 1 3 } , B = { Y < sin ⁡ π X } . A=\{X < \frac{1}{3} \}, \quad B=\{Y < \sin \pi X \}. A={X<31​},B={Y<sinπX}.
那么试计算$\mathrm{Pr}(B|A)$。

从理论上计算：

由题设可知，$X$和$Y$的$\text{ p.d.f.}$是：
$$f_1(x)=1,f_2(y)=1.$$
因为$X$和$Y$互相独立，所以二者的$joint\text{ p.d.f.}$是：
$$f(x,y)=f_1(x)f_2(y)=1$$
所以，
$$\mathrm{Pr}(A\cap B)={\int }_0^{\frac{1}{3}}{\int }_0^{\sin \pi x}f(x,y)dxdy=\frac{1}{2\pi }$$
很容易就知道：
$$\mathrm{Pr}(A)=\frac{1}{3}$$
由贝叶斯定理可知：
$$\mathrm{Pr}(B|A)=\frac{\mathrm{Pr}(A\cap B)}{\mathrm{Pr}(A)}=\frac{3}{2\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

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

另外，笔者还尝试另一种计算方式：

令$Z=Y-\sin \pi X$，那么当计算出$Z$的$\text{c.d.f.}$ $G(z)$后，那么$\mathrm{Pr}(B)=G(0).$

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

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

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

我们看一个变量是否符合某个特定的概率分布，可以看这个变量的$\text{p.d.f./p.f.}$或者$\text{c.d.f.}$或者$\text{m.g.f.}$是否符合那个特定的概率分布。因为上述$\text{p.d.f./p.f./c.d.f./m.g.f.}$都是一个概率分布的特征函数。

所以，这次会模拟出“多个独立并符合同一个正态分布的变量的平方和”这个变量的$\text{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")

模拟出的$\text{c.d.f.}$曲线：

我们可以看看卡方分布理论上的$\text{c.d.f.}$曲线：

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

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

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

如同题目三，这次是看柯西分布的平均值。同样是看$\text{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$函数的曲线，而理论上一个柯西分布的$\text{c.d.f.}$就是一个$\arctan$函数。
• 当$n$值（$\text{i.i.d.}$变量的数目）变化时，它们的模拟$\text{c.d.f.}$曲线是重合的，说明它们不仅都符合柯西分布，而且是同分布的。

题目五：马尔可夫链练习

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

假设有一个四宫格，从一个格子到另一个格子的概率有一个概率转移矩阵$\mathbf{P}$。如果最初是在格子$1$里面，每走一步都是随机的，那么$10$步后到达格子$3$的概率是多少？

它的理论值很容易推断：
$${\mathbf{P}}_{1,3}^{10}$$

用于模拟的代码是：

# 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 去做一些模拟，希望能对朋友们有所帮助。如果文中有任何错误，期望大家能指正！