蒙特卡罗模拟计算定积分(R)

最新推荐文章于 2024-06-06 20:00:18 发布

一大盘洋芋

最新推荐文章于 2024-06-06 20:00:18 发布

阅读量1.5k

点赞数 1

文章标签：概率论

本文链接：https://blog.csdn.net/weixin_51362945/article/details/127713214

版权

参考资料：概率论与数理统计教程第二版（茆诗松）4.3

定积分线性变换（换元法）

对于一般区间[a,b]上的定积分：

$J=\int_{a}^{b}g(x)dx$

可以作线性变换y=(x-a)/(b-a)，转化为[0,1]区间上的积分：

若 $c=g(a)\leq g(x)\leq g(b)= d$ ，

令 $f(y)=\frac{1}{d-c}[g(a+(b-a)y)-c]$

则 $0\leq f(y) \leq 1$ ，此时：

$J=\int_{a}^{b}g(x)dx=S_{0}J'+c(b-a)$

其中， $S_{0}=(b-a)(d-c)$ ， $J'=\int_{0}^{1}f(y)dy$

蒙特卡罗模拟

随机投点法（伯努利大数定律）

设二维随机变量(X,Y)服从 $\left\{0\leq x\leq 1,0\leq y\leq 1\right\}$ 上的均匀分布且独立。

记事件 $A=\left\{Y\leq f(X)\right\}$ ，其概率为：

$P(Y\leq f(X))=\int_{0}^{1}\int_{0}^{f(x)}dydx=\int_{0}^{1}f(x)dx=J'$

用蒙特卡罗方法随机投点，将(X,Y)看成正方形 $\left\{0\leq x\leq 1,0\leq y\leq 1\right\}$ 内以均匀分布随机投的点，计算随机点落在区域A中的频率（即 $y_{i}\leq f(x_{i})i=1,..n$ 发生的次数占比），当n很大时，该频率作为 $J'$ 的近似概率值（伯努利大数定律）。

例如，计算 $\int_{2}^{3}x^{2}dx=\frac{19}{3}=6.333$

a=2
b=3
g=function(x){
  t=x**2
  return(t)
}
c=g(a)
d=g(b)
f=function(y){
  t=1/(d-c)*(g(a+(b-a)*y)-c)
  return(t)
}
n=100000
m=function(n){
  m1=runif(n)
  m2=runif(n)
  s=0
  for(i in 1:n){
    if(m2[i]<=f(m1[i])){
      s=s+1
    }
  }
  p=s/n
  s0=(b-a)*(d-c)
  J=s0*p+c*(b-a)
  return(J)
}
m(n)

结果为：6.34315 6.3347 6.33895 6.3346等

平均值法（辛钦大数定理）

设随机变量X服从(0,1)上的均匀分布，则f(x)的数学期望为：

$E(f(X))=\int_{0}^{1}f(x)dx=J'$

由辛钦大数定理，可以用f(x)观察值的平均值去估计数学期望。使用蒙特卡罗模拟，先生成n个(0,1)上的均匀分布随机数（n足够大），再由这些随机数确定f(x)观察值的平均值：

$J'=\frac{1}{n}\sum_{i=1}^{n}f(x_{i})$

例如，计算 $\int_{2}^{3}x^{2}dx=\frac{19}{3}=6.333$

a=2
b=3
g=function(x){
  t=x**2
  return(t)
}
c=g(a)
d=g(b)
f=function(y){
  t=1/(d-c)*(g(a+(b-a)*y)-c)
  return(t)
}
n=100000
m=function(n){
  m1=runif(n)
  p=f(m1[1])
  for(i in 2:n){
    p=(p+f(m1[i]))/2
  }
  s0=(b-a)*(d-c)
  J=s0*p+c*(b-a)
  return(J)
}
m(n)

结果为：6.350545 5.169826 7.269659 6.478476等

注：平均值容易受极端值影响，误差可能较大。