以下是使用R语言实现蒙特卡罗方法计算e^x在(-1,1)上的定积分,以及比较不同算法的误差:
首先定义被积函数为f(x)=exp(x),然后用以下代码生成100000个随机点,并计算它们是否位于(-1,1)之间。最后,将所有落在该区间内的点的f(x)值求和并除以总的点数得到定积分的近似值。
set.seed(123)
n <- 100000
x <- runif(n, -1, 1)
y <- runif(n, 0, exp(1))
count <- sum(y <= exp(x))
integral <- count/n * (2*exp(1))
cat("The estimated integral using random sampling is", integral, "\n")
接下来是平均值估计法,它利用取样的平均值来近似估计目标函数的期望值。代码如下:
set.seed(123)
n <- 100000
x <- runif(n, -1, 1)
integral <- mean(exp(x))
cat("The estimated integral using mean estimation is", integral, "\n")
然后是重要抽样法,它使用一个更适合于函数形状的抽样分布来提高采样效率。以下是使用正态分布作为重要分布的示例代码:
set.seed(123)
n <- 100000
x <- rnorm(n, 0, 1)
y <- exp(x)
integral <- mean(y/exp(0))
cat("The estimated integral using importance sampling is", integral,