所以我使用蒙特卡罗方法来评估一堆函数的定积分。 首先，
y = x ^ (-0.5) ; for x in [0.01,1]
为此，我在 R 中的代码如下所示
#s <- NULL
m<- 100
a<- 0.01
b<- 1
set.seed(5)
x<-runif(m,a,b)
y<-runif(m,0,1)
for (i in 1:m){
if(y[i]<(x[i]^(-0.5))){
s[i] <- 1
}
else{
s[i] <-0
}
}
nn<-sum(s==1)*(b-a)/m
print(nn)
答案 (nn)：0.99
实际答案：1.8
我无法弄清楚我哪里出了问题。我是不是做错了什么？
【问题讨论】：
标签： r statistics simulation
一个小于 1 的负数次方的数总是大于任何小于 1 的数，所以当你得到一个全为 1 的向量时，你不应该感到惊讶。
您使用的矩形太短（高度为 1）。实际上，它应该是 10 高（因为 0.01^-0.5=10）是最大值。
然后你把矩形的总面积乘以s的平均值，所以修改后的代码是这样的：
s <- NULL
m<- 100
a<- 0.01
b<- 1
set.seed(5)
x<-runif(m,a,b)
y<-10*runif(m,0,1)
for (i in 1:m){
if(y[i]<(x[i]^(-0.5))){
s[i] <- 1
}
else{
s[i] <-0
}
}
nn<-sum(s)*(b-a)/m*10#note that the addition of the area of the rectangle
print(nn)
我得到了 1.683 的结果，这更接近真实答案。
编辑：做了一个多余的乘法，答案稍作修改
【讨论】：
m<-100000 a<-0.01 b<-1 set.seed(7063257) x<-runif(m,a,b) y<-10*runif(m,0,1) s<-ifelse(y<(x^(-0.05)),1,0) s<-sum(s)*10*(b-a)/m print(s) 我得到大约 1.038
正如 user1362215 指出的那样，您的函数应该包含在矩形中。如果增加 n，您将更接近解决方案。这是一个矢量化的解决方案。结果在范围内。
# Hit and miss
f <- function(x) x ^ (-0.5)
n <- 1000000
a <- 0.01
b <- 1
#ceiling(max(f((seq(0.01,1,by=0.001)))))
#[1] 10
set.seed(5)
x <- runif(n,a,b)
y <- 10*runif(n,0,1)
R <- sum(y < f(x))/n
(b-a)*10*R
#[1] 1.805701
# Repeat a few times to look at the distribution
set.seed(5)
n <- 100000
r <- replicate(1000,sum(10*runif(n,0,1) < f(runif(n,a,b)))/n *(b-a)*10)
hist(r)
summary(r)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1.755 1.792 1.800 1.800 1.809 1.845
# Sample mean method for comparison
set.seed(5)
r <- replicate(1000, mean(f(runif(n, a,b)))*(b-a))
hist(r)
summary(r)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1.788 1.798 1.800 1.800 1.803 1.813
重新编辑：我假设 x*2 + y^2, [-1,1] 您指的是一个圆而不是函数 f(z)。所以真的要通过模拟来估计单位圆的面积/Pi。
f2 <- function(x) sqrt(1-x^2)
s <- seq(-1 , 1 ,by=0.001)
plot(s,f2(s))
# Get the max value of function within the range
c <- ceiling(max(f2(s)))
# [1] 1
n <- 1000000
a <- -1
b <- 1
set.seed(5)
x <- runif(n,a,b)
y <- c*runif(n,0,1)
R <- sum(y < f2(x))/n
(b-a)*c*R
#[1] 1.57063 # multiply it by 2 to get full area
pi/2
#[1] 1.570796
【讨论】：
接受/拒绝的 Monte Carlo 替代方法是统一生成 x 值，平均得到的 y = f(x) 值以估计平均高度，然后将其乘以区间长度以获得估计面积。我对 R 不太了解，所以这里用 Ruby 来说明算法：
def f(x)
x ** -0.5
end
sum = 0.0
10000.times { sum += f(0.01 + 0.99 * rand) }
print (1.0 - 0.01) * (sum / 10000)
我得到的结果在 1.8 +/- 0.02 范围内
您还可以通过使用对立随机变量来提高估算器的精度 - 对于您生成的每个 x，还可以使用与 x 的中位数镜像的对称 x 值。
使用@user20650 的代码指导如何在 R 中执行此操作，您可以估计 Pi / 2 如下：
f <- function(x) sqrt(1-x^2)
n <- 100000
a <- -1
b <- 1
range <- b-a
set.seed(5)
r <- replicate(1000, mean(f(runif(n,a,b))) * range)
hist(r)
summary(r)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1.566 1.570 1.571 1.571 1.572 1.575
这种方法不需要边界函数，通常它比接受/拒绝方法产生更高的精度。
【讨论】：