R语言 几何概型计算积分 蒙特卡罗方法

设计几何概型实验，估算积分的值

代码

详细解释在代码注释中

# 定义一个函数，该函数为所求的积分函数
# 因此在计算上述定积分和计算sqrt(2)同理
# 在此用例为sqrt(2)
f = function(x)
  sqrt(2)
# 生成x
x = seq(1.5, 3.5, length=100)
y = rep(0, length(x))
j = 1
# 计算每个x对应的y值
for (i in x) {
  y[j] = f(i)
  j = j + 1
}
# 根据函数划出积分曲线
plot(x, y, type='o')
# 确定积分边界
abline(v = 2)
abline(v = 3)
# 简单设定随机点 (x,y)|2<=x<=3,0<=y<=100
# 即随机点的分布面积为100
s = 100
a = 0
cnt <- 1
# 有兴趣的朋友可以适当的多循环几次，这样结果更为精确，但是运行时间太长，并没有跑太久
while (cnt <= 100000) {
    # 随机点 xx，yy
    xx = 2 + runif(1)
    yy = 100 * runif(1)
    # 点在积分面积内则点为红色
    if( f(xx) > yy ) {
        a = a + 1
        points(x = xx, y = yy, pch = 20, cex = 1, col = "red")
    }
    # 否则为绿色
    else
        points(x = xx, y = yy, pch = 20, cex = 1, col = "green")
  cnt = cnt + 1
}
# 计算积分
print((a / 10000) * s)