原創(chuàng)：hxj7
本文繼續(xù)介紹一些和概率統(tǒng)計(jì)相關(guān)的模擬。
前文《R-概率統(tǒng)計(jì)與模擬》介紹了一些用 R 進(jìn)行概率模擬的實(shí)驗(yàn)靠欢，本文繼續(xù)上次的工作，并在此過(guò)程中回顧一些相關(guān)的概率統(tǒng)計(jì)知識(shí)盐肃。
一共五題：
- 對(duì)pi值的估計(jì)（蒙特卡洛模擬經(jīng)典示例）
- 貝葉斯公式練習(xí)
- 多個(gè)獨(dú)立并符合同一個(gè)正態(tài)分布的變量的平方和符合卡方分布
- 多個(gè)獨(dú)立且符合同一個(gè)柯西分布的變量的平均值仍符合柯西分布
- 馬爾可夫鏈練習(xí)
題目一：對(duì)pi值的估計(jì)（蒙特卡洛模擬經(jīng)典示例）
圓周率 pi ()的值大家很早就知道了俯艰。其實(shí)我們可以用模擬的方法來(lái)自行估算。方法是：
假設(shè)有一個(gè)單位圓（半徑為1）懂更，我們?cè)谒耐饨泳匦卫镫S機(jī)產(chǎn)生一個(gè)點(diǎn)眨业，那么這個(gè)點(diǎn)落在單位圓中的理論概率是
其中和
分別是單位圓和外接矩形的面積。所以沮协，理論上：
那么龄捡，我們只需要模擬出這個(gè)值，就可以得到
的估計(jì)值慷暂。
如果按上述隨機(jī)產(chǎn)生點(diǎn)的方法聘殖，是：
其中，是落在單位圓里的點(diǎn)的數(shù)量，
是落在外接矩形里的點(diǎn)的數(shù)量奸腺。
得到
按照上述方法，筆者寫了一段隨機(jī)產(chǎn)生點(diǎn)的代碼：
# Q1: the value of pi
# 也是蒙特卡洛方法的經(jīng)典示例
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
產(chǎn)生的點(diǎn)如下：
題目二：貝葉斯公式練習(xí)
假設(shè)變量
和
互相獨(dú)立突照，且都服從
上的均勻分布坠宴，并且令事件
和
為：
那么試計(jì)算。
從理論上計(jì)算：
由題設(shè)可知绷旗，
是：
因?yàn)?img class="math-inline" src="https://math.jianshu.com/math?formula=X" alt="X" mathimg="1">和是：
所以，
很容易就知道：
由貝葉斯定理可知：
用于模擬的代碼是：
# 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
可以看出衔肢，當(dāng)模擬的重復(fù)結(jié)果達(dá)到一百萬(wàn)次時(shí)庄岖，模擬的結(jié)果和實(shí)際值很接近。
另外角骤，筆者還嘗試另一種計(jì)算方式：
令
隅忿，那么當(dāng)計(jì)算出
的
后，那么
但是沒(méi)有成功邦尊。如果有朋友找到計(jì)算的方法背桐，希望能指導(dǎo)一下。
題目三：多個(gè)獨(dú)立并符合同一個(gè)正態(tài)分布的變量的平方和符合卡方分布
正如標(biāo)題所說(shuō)蝉揍，模擬的任務(wù)就是看看多個(gè)獨(dú)立并符合同一個(gè)正態(tài)分布的變量的平方和是否符合卡方分布链峭。我們會(huì)嘗試不同的變量數(shù)目進(jìn)行模擬。
我們看一個(gè)變量是否符合某個(gè)特定的概率分布又沾，可以看這個(gè)變量的
或者
是否符合那個(gè)特定的概率分布弊仪。因?yàn)樯鲜?img class="math-inline" src="https://math.jianshu.com/math?formula=%5Ctext%7Bp.d.f.%2Fp.f.%2Fc.d.f.%2Fm.g.f.%7D" alt="\text{p.d.f./p.f./c.d.f./m.g.f.}" mathimg="1">都是一個(gè)概率分布的特征函數(shù)。
所以杖刷，這次會(huì)模擬出“多個(gè)獨(dú)立并符合同一個(gè)正態(tài)分布的變量的平方和”這個(gè)變量的曲線励饵。
用于模擬的代碼：
# Q3: 模擬多個(gè)獨(dú)立同分布（正態(tài)分布）變量的平方和（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")
模擬出的曲線：
我們可以看看卡方分布理論上的曲線：
（圖片引自網(wǎng)頁(yè)https://wiki.mbalib.com/wiki/%E5%8D%A1%E6%96%B9%E5%88%86%E5%B8%83）
可以看出，模擬曲線和理論曲線很相像滑燃。
題目四：多個(gè)獨(dú)立且符合同一個(gè)柯西分布的變量的平均值仍符合柯西分布
如同題目三役听，這次是看柯西分布的平均值。同樣是看曲線：
模擬的代碼：
# 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")
模擬的曲線：
可以看出：
- 模擬出的曲線非常像
函數(shù)的曲線表窘，而理論上一個(gè)柯西分布的
- 當(dāng)
值（
變量的數(shù)目）變化時(shí)，它們的模擬
題目五：馬爾可夫鏈練習(xí)
馬爾可夫鏈大家應(yīng)該都很熟悉了麦备。我們用一個(gè)小題目回顧一下：
假設(shè)有一個(gè)四宮格，從一個(gè)格子到另一個(gè)格子的概率有一個(gè)概率轉(zhuǎn)移矩陣。如果最初是在格子
里面凛篙，每走一步都是隨機(jī)的黍匾，那么
步后到達(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
可以看出呛梆，當(dāng)模擬的重復(fù)次數(shù)到達(dá)十萬(wàn)次時(shí)锐涯，模擬值很接近理論值了。
小結(jié)
從前文到本文填物，我們共通過(guò)八個(gè)小題目回顧了一些概率統(tǒng)計(jì)相關(guān)的知識(shí)纹腌，并嘗試用 R 去做一些模擬，希望能對(duì)朋友們有所幫助滞磺。如果文中有任何錯(cuò)誤升薯，期望大家能指正！
（公眾號(hào)：生信了）
