【R】随机模拟计算定积分

题目：分别用随机投点法与平均值法计算定积分：

1 随机投点法

▲分析

▲代码

library(ggplot2) 
x <- seq(0,1,0.001) #生成[0,1]序列，步长0.001
h=function(x){exp(x^2)}  #被积函数h(x)
data <- data.frame(x,h(x))
ggplot(data,aes(x,h(x)))+geom_area(fill="green")

inte <- integrate(h,0,1)  #计准确积分
in_va <- inte$value  #将integrate()所得结果转换为数值形式

#给定积分区间[a,b],h(x)在该区间的最大值为c，随机数个数为n
f=function(h,a,b,c,n){
  x <- runif(n,a,b)  #产生n个随机数在区间[a,b]
  y <- c*runif(n,0,1)  #产生区间[0,c]之间的随机数
  phat <- sum(h(x)>=y)/n  #计算落入积分区域内的频率（估计概率）
  ihat <- phat*c*(b-a)  #计算估计的积分值
}

a=0
b=1
c=exp(1)
p <- inte$value/c  

NN <- 10^(2:6)  
xmat <- matrix(0,length(NN),4) 
xmat[,1] <- NN  
colnames(xmat) <- c("N","积分估计值","随机模拟误差","95%误差界")
for (i in seq(nrow(xmat))) {
  xmat[i,"积分估计值"] <- f(h,a,b,c,NN[i])  
  xmat[i,"随机模拟误差"] <- f(h,a,b,c,NN[i])-in_va
  xmat[i,"95%误差界"] <- c*(b-a)*2*sqrt(p*(1-p)/NN[i])
}
xmat

▲结果

reference R语言中随机模拟的例子（问题八 计算积分问题）

2 平均值法

▲分析

▲积分区间为[a,b]

set.seed(123) #设定随机数种子
h=function(x){exp(x^2)}  #被积分的函数h(x)
i <- integrate(h,0,1) 
i.true <- i$value  #将integrate()所得结果转换为数值形式
f=function(h,a,b,n){
  x=runif(n)  #产生n个随机数
  mean((b-a)*h(a+(b-a)*x)) #将通用函数化成[0,1]区间上的表达式
}
c <- exp(1)
p <- i.true/c  #由作业2中推导可知模拟误差幅度大约在±2√(varI)
N <- 10^(2:6)  #随机数的个数分别取100、1000、10000、100000、1000000个
xmat <- matrix(0,length(N),4)  #建立矩阵展示结果
xmat[,1] <- N  #矩阵第一列为随机数个数N
colnames(xmat) <- c("N","积分估计值","随机模拟误差","95%误差界")
for (i in seq(nrow(xmat))) {
  xmat[i,"积分估计值"] <- f(h,0,1,N[i])  
  xmat[i,"随机模拟误差"] <- f(h,0,1,N[i])-i.true
  xmat[i,"95%误差界"] <- c*2*sqrt(p*(1-p)/N[i]) #varI=e*2*sqrt(p*(1-p)/n)
}
xmat

结果：

> xmat
         N 积分估计值  随机模拟误差   95%误差界
[1,] 1e+02   1.454480 -0.0041276735 0.271038708
[2,] 1e+03   1.449414  0.0044823251 0.085709965
[3,] 1e+04   1.458422 -0.0070815831 0.027103871
[4,] 1e+05   1.463130  0.0006053228 0.008570997
[5,] 1e+06   1.462153 -0.0002039007 0.002710387

▲积分区间为[0,1]

set.seed(111) #设定随机数种子
h=function(x){exp(x^2)}  #被积分的函数h(x)
ii <- integrate(h,0,1) 
i.true <- ii$value  #将integrate()所得结果转换为数值形式
f=function(n){y <- exp((runif(n))^2)}

N <- 10^(2:6)  #随机数的个数分别取100、1000、10000、100000、1000000个
xmat <- matrix(0,length(N),7)  #建立矩阵，以展示结果
xmat[,1] <- N  #矩阵第一列为随机数个数NN
colnames(xmat) <- c("N","积分估计值","随机模拟误差","RMSE","MAE","MRE","95%误差界")
for (i in seq(nrow(xmat))) {
  xmat[i,"积分估计值"] <- mean(f(N[i]))  #平均值法求定积分估计值  
  xmat[i,"随机模拟误差"] <- mean(f(N[i]))-i.true
  xmat[i,"RMSE"] <- sd(f(N[i]))/sqrt(N[i])
  xmat[i,"MAE"] <- 0.80*sd(f(N[i]))/sqrt(N[i])
  xmat[i,"MRE"] <- 0.80*sd(f(N[i]))/(sqrt(N[i])*mean(f(N[i])))
  xmat[i,"95%误差界"] <- 2*sd(f(N[i]))
}
xmat

结果：

> xmat
         N 积分估计值  随机模拟误差         RMSE          MAE          MRE 95%误差界
[1,] 1e+02   1.428934  4.086051e-02 0.0501545105 0.0312008578 0.0236432426 0.8720956
[2,] 1e+03   1.445884 -6.023022e-03 0.0144318733 0.0121586553 0.0082699013 0.9367443
[3,] 1e+04   1.465593  1.196663e-04 0.0047539823 0.0038052839 0.0025772210 0.9483788
[4,] 1e+05   1.463676 -7.219394e-05 0.0015087275 0.0012015586 0.0008211674 0.9468650
[5,] 1e+06   1.462819 -9.094580e-04 0.0004740806 0.0003794884 0.0002596062 0.9482749

reference 蒙特卡洛方法与定积分计算-邓一硕 、R学习之蒙特卡罗积分