R-概率统计与模拟

本文记录了三个概率统计相关的小题目，以回顾一些概率统计的知识。

（公众号：生信了）

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

这几个题目都是和概率统计相关，本来都是可以推演出精确的解，但是有意思的是，笔者从一位网友处得知这类题目可以用 R 来做模拟求得一个近似解。这是笔者之前从未尝试过的，所以动手一做：

题目一：X10的期望值

假设$X_1$在$[0,1]$上服从均匀分布，$X_2$在$[0,1+X_1]$上服从均匀分布，$X_3$在$[0,1+X_2]$上服从均匀分布，依此类推，请问$X_{10}$的期望值是多少？

其实我们可以比较容易地推断出$X_i(i=1,2,\dots )$的期望值是
$$E(X_i)=1-\frac{1}{2^i}$$

所以，$X_{10}$的期望值是：
$$E(X_{10})=1-\frac{1}{2^{10}}=1-\frac{1}{1024}=0.9990234$$

假设两个随机变量$X$和$Y$，满足$Y=aX+b$，那么$E(Y)=aE(X)+b$.

这是精确解，那么如何做模拟呢？笔者没有实际动手做过模拟，但是记得“抛十万次硬币，正面朝上的次数会非常接近于五万”，所以笔者对模拟的初步认识就是用大量的随机实验去模拟，每一次随机实验会得到一个结果，这个结果要么符合我们的要求，要么不符合。所有实验的结果中符合我们要求的结果的次数除以总次数就是我们想要的概率值。

要想让模拟的结果接近真实值，模拟的总次数要足够多。为了解决这个问题，同时看看不同模拟次数的效果如何，笔者编写了一小段 R 代码：

# Q1
oxn <- function(n) {
  x <- 0
  for (i in 1:n)
    x <- runif(1, 0, 1 + x)
  x
}

sxn <- function(n, k) {
  set.seed(SEED)
  mean(replicate(n, oxn(k)))
}

SEED <- 123
repNum.xn <- 10 ^ (1:5)
xn <- 10
res.s.xn <- sapply(repNum.xn, sxn, k=xn)
res.t.xn <- 1 - 1 / 2 ^ xn   # theoretical value
plot(log10(repNum.xn), res.s.xn, type="b", main=paste0("Simulation for X", xn),
     xlab="Number of Replications (log10)", ylab="Probability", col="blue")
abline(a=res.t.xn, b=0, col="red")
legend("bottomright", legend = c("theoretical", "simulative"), 
       col=c("red", "blue"), lty=1, bty="n")

当模拟次数不同时，结果如下：

从图中可以看出，当模拟次数达到10万次时，模拟的结果已经很接近真实值了。

题目二：球投盒子

假设10个球随机投入16个盒子中，请问每个盒子的球数都小于等于1的概率是多少？

这个问题的精确解是：
$$\frac{P_{16,10}}{16^{10}}=\frac{16!}{6!16^{10}}$$

$$P_{n,k}=\frac{n!}{(n-k)!}.$$

其模拟的代码是：

# Q2
oballs <- function(k, n) {
  1 - as.numeric(any(table(sample(1:n, k, replace = T)) > 1))
}

sballs <- function(N, k, n) {
  set.seed(SEED)
  res <- replicate(N, oballs(k, n))
  mean(res)
}

SEED <- 123
nballs <- 10
nboxes <- 16
repNum.balls <- 10 ^ (1:5)
res.s.balls <- sapply(repNum.balls, sballs, k=nballs, n=nboxes)
res.t.balls <- factorial(nboxes) / factorial(nboxes - nballs) / nboxes ^ nballs
plot(log10(repNum.balls), res.s.balls, type="b", main="Simulation for Balls",
     xlab="Number of Replications (log10)", ylab="Probability", col="blue")
abline(a=res.t.balls, b=0, col="red")
legend("bottomright", legend = c("theoretical", "simulative"), 
       col=c("red", "blue"), lty=1, bty="n")

当模拟次数不同时，结果是：

从图中可以看出，当模拟次数达到1000次时，模拟的结果已经很接近真实值了。

题目三：信封问题

假设有$n$个信封，每个信封上都有地址。现在将$n$封信随机放到这$n$个信封里，请问至少有一封信放到了写着正确地址的信封里的概率是多少？

这个问题比较难，它的精确解是：
$$p_n=1-\frac{1}{2!}+\frac{1}{3!}-\frac{1}{4!}+\cdots +(-1{)}^{n+1}\frac{1}{n!}.$$

当$n$足够大时，上述结果会收敛：
$$\underset{n\to \infty }{\lim }{p}_n=1-\frac{1}{e}.$$

$$\begin{array}{rl}\mathrm{Pr}\left(\bigcup _{i=1}^n{A}_i\right)=& \sum _{i=1}^n\mathrm{Pr}(A_i)-\sum _{i<j}\mathrm{Pr}(A_i\cap A_j)+\sum _{i<j<k}\mathrm{Pr}(A_i\cap A_j\cap A_k)\\ & -\sum _{i<j<k<l}\mathrm{Pr}(A_i\cap A_j\cap A_k\cap A_l)+\cdots \\ & +(-1{)}^{n+1}\mathrm{Pr}(A_1\cap A_2\cap \cdots \cap A_n).\end{array}.$$

相对应地，其模拟代码如下：

# Q3
oletter <- function(n) {
  as.numeric(any(sample(1:n, n, replace = F) == 1:n))
}

sletter <- function(n, N) {
  set.seed(SEED)
  res <- replicate(N, oletter(n))
  mean(res)
}

rletter <- function(n) {
  sum(sapply(1:n, function(x) (-1) ^ (x %% 2 + 1) / factorial(x)))
}

SEED <- 123
nrep <- 100000
letterNum <- 2 ^ (1:8)
res.s.letter <- sapply(letterNum, sletter, N=nrep)
res.t.letter <- sapply(letterNum, rletter)
plot(log2(letterNum), res.s.letter, type="b", col="blue", 
     main="Simulation for Letters", xlab="Letter Number (log2)", 
     ylab="Probability")
lines(log2(letterNum), res.t.letter, type="b", col="red")
legend("bottomright", legend = c("theoretical", "simulative"), 
       col=c("red", "blue"), lty=1, bty="n")

不同的$n$数值对应的结果是：

从图中可以看出，当n达到8以后，概率已经趋于稳定了。