蒙特卡罗模拟 - 怎么了？

Question

我想用函数的蒙特卡罗模拟计算曲线下面积

。我想在区间 [-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）。我做错了什么？当然，经过百万次采样后，算法应该更接近 9.3333

Answer 1

我们可以像这样尝试蒙特卡洛模拟

> 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

Answer 2

这是一个显示您正在使用的内容的情节。我希望它能帮助您更好地理解我在评论中写的内容：

您似乎忽略了 y=1 下方的矩形。它的面积 (=4) 是缺失的数量。所以代码对于计算非偏移表达式 x^2 是正确的。更改为 y_n

评论是答案的一半，即您没有模拟 y_n 的 0 和 1 之间的点。那些需要在蒙特卡洛模型中整合一个区域。另一个修改是将正确的总面积 [-2

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 点图：

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