【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学习之蒙特卡罗积分

  • 23
    点赞
  • 194
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
非常感谢您的提问。根据您的描述,我会尽力回答您的问题。 针对您的问题,我了解到您想利用切尔诺夫界为计算π和定积分的数值随机算法之一建立置信区间,使得我们可以根据置信水平和置信区间估计所需随机实验的次数。这是一个比较复杂的问题,需要进行一些数学推导和分析。 首先,我们需要了解一下切尔诺夫界和置信区间的概念。切尔诺夫界是一种用于估计随机变量与其期望值之间差距的界限,也称为Chebyshev不等式。它基于一个简单的观察:对于任何一组随机变量,平均值与其标准差的比例总是小于一个特定的值。根据切尔诺夫界,对于任何一个随机变量X,对于任何k>0,有 Pr[|X-E(X)| >= k*σ] <= 1/k^2 其中E(X)表示随机变量的期望值,而σ表示随机变量的标准差。该不等式意味着,对于任意给定的k,我们可以推断出X落在以E(X)为中心,以kσ为半径的区间内的概率不小于1-1/k^2 。这个区间就是切尔诺夫置信区间。 其次,我们需要了解一下数值随机算法的基本原理。数值随机算法是一种基于随机模拟的算法,可以用于计算一些复杂的数学问题,例如π的值和定积分等。算法的基本思路是,通过随机抽样来获得原问题的近似解,然后通过估计误差来确定这个近似解的精度。 针对您的问题,我们可以考虑用Monte Carlo方法来估计π和定积分。Monte Carlo方法是一种基于随机抽样和统计推断的方法,可以用于计算复杂的数学问题。其基本思路是,通过对随机点的采样进行统计推断,从而得到原问题的近似解。对于π的计算,我们可以随机产生一批坐标落在单位圆内的点,然后通过比较圆的面积和正方形面积的比例来估计π的值。对于定积分的计算,我们可以随机采样一些点,并计算这些点在被积函数下的概率密度函数值的均值。通过用这个均值乘以积分区间的长度,我们就可以估算出定积分的值。 在使用Monte Carlo方法时,我们需要进行一些统计分析,以确定估计值的精度。根据切尔诺夫界,我们可以得到如下估计公式: Pr[|π-π_est| >= ε] <= σ^2 / (n*ε^2) 其中,π是真实值,π_est是估计值,ε是置信水平,σ是标准差,n是采样次数。通过这个公式,我们可以确定所需的采样次数,以达到所要求的置信区间。例如,如果我们要求ε=0.01,置信水平为95%,则根据公式,我们需要进行40000次采样。 希望这些信息对您有所帮助。如果您有任何其他问题,请随时向我提出,我将尽力回答。谢谢!
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值