蒙特卡洛集成在R：得到錯誤的答案[使用命中或小姐]

Question

所以我使用蒙特卡洛方法來評估定義積分的一堆功能。 首先，

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

我不知道我要去哪裏錯了。我做錯了什麼？

Answer 1

小於1的負數小於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的結果，更接近很多真正的答案。

編輯：由多餘乘法，答案略微

訂正

Answer 2

蒙特卡羅替代接受/拒絕是均勻地產生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，可以估算π/ 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 

無邊界功能，一般它產生比更大的精度接受/拒絕方法。

Answer 3

正如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