R-概率統計與模擬

2019 年 10 月 8 日
筆記

本文記錄了三個概率統計相關的小題目，以回顧一些概率統計的知識。

正如筆者在前文《公眾號一歲啦》中所說，近期在複習概率統計相關的知識。機緣巧合，筆者遇到了幾個比較有意思的題目，和朋友們分享一下：

這幾個題目都是和概率統計相關，本來都是可以推演出精確的解，但是有意思的是，筆者從一位網友處得知這類題目可以用 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以後，概率已經趨於穩定了。