2013-12-15 4 views
3

몬테카를로 메서드를 사용하여 일식 영역 (a = 6b = 3)을 계산해야합니다. 또한 내부 점이 빨간색이고 바깥 점이 검은 색 인 결과의 플롯 (다이어그램)을 만들어야합니다. 마지막에 내가 "정기적 인 결과"와 "몬테카를로 결과"를 비교해야몬테카를로 메서드가있는 R - 타원 영역

방정식은 방법 100000 개 응답이 있어야합니다 (x^2)/36+(y^2)/9=1

입니다.

이것은 내가하는 일입니다. 분명히 작동하지 않습니다.

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 

답변

3

, 그러나 여기에 우리가 간다 :

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) 

enter image description here

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

정말 고마워요! 정말! – fiore

1
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) 

enter image description here