R - 蒙特卡羅方法的橢圓區域

Question

我需要用蒙特卡羅方法計算日食的面積（a = 6 b = 3）。 此外，我必須做一個內部點紅色和黑色的結果圖（圖）。最後我有「常規結果」比較「蒙特卡洛結果」

的公式爲(x^2)/36+(y^2)/9=1

的方法必須有10萬個回覆。

這就是我所做的。顯然它不起作用。

set.seed(157619) 

n <- 100000 

xmin <- (-6) 
xmax <- (+6) 

ymin <- (-3) 
ymax <- (+3) 

rx <- (xmax-xmin)/2 
ry <- (ymax-ymin)/2 


outa <- runif(n,min=xmin,max=xmax) 
outb <- runif(n,min=ymin,max=ymax) 

dx <- outa*2 
dy <- outb*2 

ly <- dy<=(ry^2); my <- dy>(ry^2) 
lx <- dx<=(ry^2); mx <- dx>(rx^2) 

這是爲圓工作的範例代碼：

n <- 200 
xmin <- -1; xmax <- 1 
r <- (xmax-xmin)/2 
out <- runif(n,min=xmin,max=xmax) 
x <- matrix(out,ncol=2) 
d <- x[,1]^2 + x[,2]^2 
l <- d<=(r^2); m <- d>(r^2) 
win.graph(7,7.8) # così è quadrato 
plot(c(xmin,xmax),c(xmin,xmax),type="n") 
plot(x[l,1],x[l,2]) 
points(x[m,1],x[m,2],col="red",pch=19) 
(p <- sum(l)/length(l)) 
p*4 

Answer 1

我懷疑這是功課，但在這裏我們去：

set.seed(42) 
n <- 1e5 
xmax <- 6 
ymax <- 3 

x <- runif(n, 0, xmax) 
y <- runif(n, 0, ymax) 

inside <- (x^2)/36+(y^2)/9 <= 1 

plot(x, y, pch=16, cex=0.5, col=inside+1) 

mean(inside) * (xmax*ymax) *4 
#[1] 56.54376 
pi*6*3 
#[1] 56.54867 

Answer 2

set.seed(1) 
n = 1000 
a = 6 
b = 3 
x.samp = runif(n, -a, a) 
y.samp = runif(n, -b, b) 

p.in = (x.samp/a)^2 + (y.samp/b)^2 <= 1 

S = 4*a*b*sum(p.in)/n 
print(S) 

plot(x.samp, y.samp, col = p.in + 1)