我需要使用蒙特卡洛方法来计算日蚀的面积(a = 6 b = 3)。
我还必须绘制结果图(图表)，内部点为红色，外部点为黑色。最后，我必须将“蒙特卡洛结果”与“常规结果”进行比较
等式是(x^2)/36+(y^2)/9=1
该方法必须有100000个答复。
这就是我的工作。显然，它不起作用。
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
最佳答案
我怀疑这是家庭作业，但是我们开始:
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