Monte carlo 模拟 - 怎么了?

Monte carlo simulation - what's wrong?

我想使用monte carlo函数模拟计算曲线下面积

。我想在区间 [-2, 2]

上计算它

我目前的工作

# Define function f
f <- function(x) x^2 + 1

# I want to close my function in rectangle (-2, 2) - x axis and (1, 5) y -axis
n <- 10^6
# Randomize from x axis of rectangle
x_n <- runif(n, min = -2, max = 2)
# Randomize from y axis of rectangle
y_n <- runif(n, min = 1, max = 5)
# Calculate function values of randomized points
values <- f(x_n)

# Formula for are under the curve for monte carlo simulation is 
# Area of rectangle * (Points below curve) / (Total number of points)

所以我的结果是:

> sum(y_n < values) / n * (4 * 4)
[1] 5.329888

这是错误的结果(正确的结果是 9.33333)。我做错了什么?确保算法在 milion 采样后应该更接近 9.3333

这里有一个图表可以显示您正在使用的内容。我希望它能帮助您更好地理解我在评论中写的内容:

You seem to be ignoring the rectangle below y=1. It's area (=4) is the missing quantity. So the code is correct for calculating the non-offset expression x^2. Change to y_n <- runif(n, min = 0, max = 5) and re-run the calculations

该评论是答案的一半,即您没有为 y_n 模拟 0 和 1 之间的点。那些需要在 Monte Carlo model 的区域整合中。另一个 mod化是将正确的总面积 [-2 < x <2]x[0

f <- function(x) x^2 + 1

# I want to close my function in rectangle (-2, 2) - x axis and (1, 5) y -axis
n <- 10^6
# Randomize from x axis of rectangle
x_n <- runif(n, min = -2, max = 2)
# Randomize from y axis of rectangle
y_n <- runif(n, min = 0, max = 5)
# Calculate function values of randomized points
values <- f(x_n)

# Formula for are under the curve for monte carlo simulation is 
# Area of rectangle * (Points below curve) / (Total number of points)
 sum(y_n < values) / n * (5 * 4)
#[1] 9.3429  inaccurate to 1 or 2 parts in 933

显示第二种情况的 100 点图:

您可能会考虑的另一个 mod 是使用 set.seed 使您的计算可重现。

我们可以像这样尝试Monte Carlo模拟

> n <- 1e6

> x <- runif(n, -2, 2)

> y <- runif(n, 0, 5)

> mean(x^2 + 1 - y >= 0) * 4 * 5
[1] 9.33014

其中面积可以计算为位于曲线下方的平均点数 x^2 + 1 -y >=0