以下是使用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, "\n")

最后是分层抽样法,它将取样区间划分为多个子区间,并在每个子区间中进行采样。以下是将(-1, 0)和(0, 1)两个区间各自作为一个层级的示例代码:

set.seed(123)

n <- 100000

x1 <- runif(n/2, -1, 0)
y1 <- exp(x1)

x2 <- runif(n/2, 0, 1)
y2 <- exp(x2)

integral <- (mean(y1)*(1/2) + mean(y2)*(1/2)) * 2

cat("The estimated integral using stratified sampling is", integral, "\n")

由于我们无法知道真实的积分值,因此无法确定哪种方法的误差最小。但是,我们可以比较估计值与蒙特卡罗方法所生成的不同随机样本数量的真实积分值之间的距离来评估误差。通常,误差会随着样本数量的增加而减小。

以下是使用10个不同的样本大小(从1到10万)计算每种估计方法的误差并绘制误差图的示例代码:

set.seed(123)

f <- function(x) exp(x)

exact_integral <- integrate(f, -1, 1)$value

n_samples <- 10^((1:5)*2)
exact_integrals <- rep(exact_integral, length(n_samples))

random_sampling_errors <- c()
mean_estimation_errors <- c()
importance_sampling_errors <- c()
stratified_sampling_errors <- c()

for (n in n_samples) {
  x <- runif(n, -1, 1)
  y <- runif(n, 0, exp(1))
  
  count <- sum(y <= f(x))
  integral <- count/n * (2*exp(1))
  random_sampling_errors <- c(random_sampling_errors, abs(integral - exact_integral))
  
  integral <- mean(f(x))
  mean_estimation_errors <- c(mean_estimation_errors, abs(integral - exact_integral))
  
  x <- rnorm(n, 0, 1)
  y <- f(x)
  
  integral <- mean(y/exp(0))
  importance_sampling_errors <- c(importance_sampling_errors, abs(integral - exact_integral))
  
  x1 <- runif(n/2, -1, 0)
  y1 <- f(x1)
  x2 <- runif(n/2, 0, 1)
  y2 <- f(x2)
  
  integral <- (mean(y1)*(1/2) + mean(y2)*(1/2)) * 2
  stratified_sampling_errors <- c(stratified_sampling_errors, abs(integral - exact_integral))
}

plot(log10(n_samples), log10(random_sampling_errors), type="b", col="red", ylim=range(log10(c(random_sampling_errors, mean_estimation_errors, importance_sampling_errors, stratified_sampling_errors))))
points(log10(n_samples), log10(mean_estimation_errors), type="b", col="blue")
points(log10(n_samples), log10(importance_sampling_errors), type="b", col="green")
points(log10(n_samples), log10(stratified_sampling_errors), type="b", col="purple")
legend("topright", legend=c("random sampling", "mean estimation", "importance sampling", "stratified sampling"), col=c("red", "blue", "green", "purple"), lty=1)

该代码会生成一个图形,其中包含每种方法的误差随着样本数量的变化而变化的情况。根据该图形,我们可以看到重要抽样法和分层抽样法通常比随机投点法和平均值估计法具有更小的误差。

请注意,这里提供的所有示例代码都是基于特定函数和参数的。实际上,您需要根据具体问题调整这些代码以获得准确的结果。此外,使用蒙特卡罗方法时,您还需要考虑到关于采样分布,步长等方面的其他细节。

Logo

魔乐社区(Modelers.cn) 是一个中立、公益的人工智能社区,提供人工智能工具、模型、数据的托管、展示与应用协同服务,为人工智能开发及爱好者搭建开放的学习交流平台。社区通过理事会方式运作,由全产业链共同建设、共同运营、共同享有,推动国产AI生态繁荣发展。

更多推荐