7 随机模拟积分
7.1 习题
习题1
- 利用恒等式\(EK=\sum_{n=0}^{\infty} P\{K>n\}\)。
不难得到，\(P\{K>n\}=P\{\sum_{i=1}^n U_i<1\}\)，\(U_i\)iid标准均匀分布，这涉及到标准均匀分布随机数之和的分布，学名为Irwin-Hall分布，我们不妨设\(X_n=\sum_{i=1}^n U_i\)，其密度函数为\(f_n(x)\)，由Irwin-Hall分布可以知道， \[ P\{K>n\}=P\{\sum_{i=1}^n U_i<1\}=\int_0^1 f_n(x)\mathrm dx=\frac1{n!} \]
再利用恒等式\(EK=\sum_{n=0}^{\infty} P\{K>n\}=\sum_{n=0}^{\infty} \frac1{n!}=e\)
- 定义函数
sampling_K <- function(n)
{
k <- 10 # 期望多少个均匀分布随机数之和能大于1
U <- runif(k*n)
U <- matrix(U,n,k)
S <- t(apply(U,1,cumsum))
x <- vector("integer",n)
# 按行循环
for(i in seq_along(x)){
j <- match(TRUE,S[i,]>1,nomatch = -1)
# 匹配失败
if(j==-1){
jj <- 0
while(j!=-1){
u <- cumsum(runif(k))
j <- match(TRUE,u>1,nomatch = -1)
jj <- jj+1
}
x[[i]] <- jj*k+j
}
# 匹配成功
else{
x[[i]] <- j
}
}
list(x=x,S=S)
}测试一下：
sampling_K(10)
#> $x
#> [1] 4 2 5 3 2 3 2 2 2 2
#>
#> $S
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> [1,] 0.21808916 0.2357843 0.2658886 1.2296056 1.592953 1.887922 2.140974
#> [2,] 0.98763424 1.9464288 2.3330987 2.9809780 3.850568 4.545099 4.878368
#> [3,] 0.34846189 0.6999084 0.8075462 0.9006072 1.690359 2.342677 2.983584
#> [4,] 0.38104699 0.8358012 1.3227212 2.1230489 3.029562 3.403099 3.458583
#> [5,] 0.02098596 1.0019890 1.8790498 2.2702906 3.216077 4.026435 4.870650
#> [6,] 0.74972687 0.9305819 1.6982691 2.4187291 2.497471 2.874808 3.553516
#> [7,] 0.16295804 1.0107804 1.5122450 1.7406904 2.666825 2.871587 3.278336
#> [8,] 0.31904129 1.1742926 1.4711369 2.3982398 2.572193 3.375264 3.975827
#> [9,] 0.59789026 1.2797888 1.7417666 1.8840650 2.333020 3.101203 4.013269
#> [10,] 0.58799732 1.0799603 1.6595000 2.0488883 3.037629 3.134606 3.559394
#> [,8] [,9] [,10]
#> [1,] 2.535080 2.973997 3.701060
#> [2,] 5.344494 5.515443 5.869804
#> [3,] 3.829691 4.718623 5.003419
#> [4,] 3.758749 4.060928 4.226667
#> [5,] 5.260501 5.512067 5.698450
#> [6,] 4.420283 4.758983 5.289225
#> [7,] 3.546575 3.659780 3.665562
#> [8,] 4.769557 5.550969 5.796674
#> [9,] 4.322063 5.169168 5.891357
#> [10,] 3.672081 4.365124 4.733642估计e：
n <- 100000
x <- sampling_K(n)$x
e <- mean(x)
e
#> [1] 2.71649- 直接用R求标准差：
sd(x)
#> [1] 0.874894195%置信区间为\(\hat\theta \pm 2 \frac{S_N}{ \sqrt{N}}\)：
d <- sd(x)/sqrt(n)
cat("95%置信区间为(",e-d,",",e+d,")",sep = "")
#> 95%置信区间为(2.713723,2.719257)习题2
- 容易知道： \[ P\{M>n\} = P\{U_1\le U_2\le\cdots\le U_n\} \]
由于\(U_i\overset{iid}{\sim}U(0,1)\)，记事件集合\(A=\{U_{i_1}\le U_{i_2}\le\cdots\le U_{i_n},i_1,\cdots,i_n\text{是}1,\cdots,n\text{的一个排列}\}\)，则\((U_1,\cdots,U_n)\)在集合\(A\)上均匀分布，从而 \[ P\{M>n\} = P\{U_1\le U_2\le\cdots\le U_n\}=\frac1{n!} \]
显然、易证。
定义函数：
sampling_M <- function(n)
{
k <- 10 # 期望值
U <- runif(k*n)
U <- matrix(U,n,k)
S <- U[,1:(k-1)]>U[,2:k]
x <- vector("integer",n)
# 按行循环
for(i in seq_along(x)){
j <- match(TRUE,S[i,],nomatch = -1)
# 匹配失败
if(j==-1){
jj <- 0
while(j!=-1){
u <- runif(k)
u <- u[1:(k-1)]>u[2:k]
j <- match(TRUE,u,nomatch = -1)
jj <- jj+1
}
x[[i]] <- jj*k+j+1
}
# 匹配成功
else{
x[[i]] <- j+1
}
}
list(x=x,S=S)
}测试一下：
sampling_M(10)
#> $x
#> [1] 3 2 2 3 3 2 3 2 3 2
#>
#> $S
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] FALSE TRUE TRUE FALSE TRUE FALSE FALSE FALSE TRUE
#> [2,] TRUE FALSE FALSE TRUE TRUE FALSE TRUE TRUE FALSE
#> [3,] TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE TRUE
#> [4,] FALSE TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE
#> [5,] FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE TRUE
#> [6,] TRUE FALSE FALSE TRUE FALSE TRUE TRUE FALSE TRUE
#> [7,] FALSE TRUE TRUE FALSE TRUE FALSE TRUE FALSE FALSE
#> [8,] TRUE TRUE FALSE TRUE FALSE TRUE FALSE TRUE FALSE
#> [9,] FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE FALSE
#> [10,] TRUE TRUE TRUE FALSE TRUE FALSE TRUE FALSE TRUE估计e：
n <- 100000
x <- sampling_M(n)$x
e <- mean(x)
e
#> [1] 2.71777d <- sd(x)/sqrt(n)
cat("估计的标准差：",
sd(x),
"，近似95%的置信区间：(",
e-d,",",
e+d,")",
sep="")
#> 估计的标准差：0.8767234，近似95%的置信区间：(2.714998,2.720542)