R-概率统计与模拟
本文记录了三个概率统计相关的小题目，以回顾一些概率统计的知识。
正如笔者在前文《公众号一岁啦》中所说，近期在复习概率统计相关的知识。机缘巧合，笔者遇到了几个比较有意思的题目，和朋友们分享一下：
这几个题目都是和概率统计相关，本来都是可以推演出精确的解，但是有意思的是，笔者从一位网友处得知这类题目可以用 R 来做模拟求得一个近似解。这是笔者之前从未尝试过的，所以动手一做：
题目一：X10的期望值
这是精确解，那么如何做模拟呢？笔者没有实际动手做过模拟，但是记得“抛十万次硬币，正面朝上的次数会非常接近于五万”，所以笔者对模拟的初步认识就是用大量的随机实验去模拟，每一次随机实验会得到一个结果，这个结果要么符合我们的要求，要么不符合。所有实验的结果中符合我们要求的结果的次数除以总次数就是我们想要的概率值。
要想让模拟的结果接近真实值，模拟的总次数要足够多。为了解决这个问题，同时看看不同模拟次数的效果如何，笔者编写了一小段 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的概率是多少？
这个问题的精确解是：
其模拟的代码是：
# 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次时，模拟的结果已经很接近真实值了。
题目三：信封问题
相对应地，其模拟代码如下：
# 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以后，概率已经趋于稳定了。
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2019-09-26，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录