【问题标题】：How to approximate the solution to a double integral in R using hit and miss and Halton?如何使用命中和未命中和 Halton 逼近 R 中双积分的解？
【发布时间】：2019-03-05 08:47:40
【问题描述】：

我有以下问题：

我做的是：

set.seed(1)
N = 10000
f = function(x,y) x^y * y^x

那我不知道该怎么办了。有人可以向我解释如何做到这一点吗？谢谢！

【问题讨论】：

标签： r integral

【解决方案1】：

积分的蒙特卡罗近似得到如下：

set.seed(1)
f <- function(x,y) x^y * y^x

N <- 10000
mean(f(runif(N), runif(N)))
# 0.4293375

这大概就是题中所谓的“统一方法”吧。

与数值近似比较（更好）：

library(cubature)
f <- function(xy){
  x <- xy[1]; y <- xy[2]
  x^y * y^x
} 
integral <- pcubature(f, lowerLimit = c(0,0), upperLimit = c(1,1))
integral$integral
# 0.4280186

Wolfram|Alpha 给出0.42802。

现在，通过采样 Halton 序列获得的近似值：

library(randtoolbox)
setSeed(1)
xy <- halton(n = N, dim = 2)
f <- function(x,y) x^y * y^x
mean(f(xy[,1], xy[,2]))
# 0.4277797

请注意，我不确定这是不是想要的答案，因为我还不知道什么是“命中或未命中”方法。

编辑

我看了一下Hit-or-miss 方法。它包括在 3D 空间中采样点，并获取这些点在f 定义的表面以下的比例。这里函数f取值在[0,1]，所以近似得到如下：

library(randtoolbox)
setSeed(1)
xyz <- halton(n = N, dim = 3)
f <- function(x,y) x^y * y^x
mean(f(xyz[,1],xyz[,2]) > xyz[,3])
# 0.4287

【讨论】：