R - Wacky Result에서 이중 적분을 사용하여 함수를 시각화합니다.

Question

수분 분포 곡선을 시각화하려고합니다. 나는 R에 매우 익숙하지 않으므로 제 어리 석음으로 화를 내지 말라.

llim <- 0 
ulim <- 6.29 

f <- function(x,y) {(.156812/((2*pi)*(.000005^2)*(gamma(2/.156812)))*exp(-((sqrt(x^2+y^2))/.000005)^.156812))} 

integrate(function(y) { 
    sapply(y, function(y) { 
      integrate(function(x) f(x,y), llim, ulim)$value 
      }) 
    }, llim, ulim) 

fv <- Vectorize(f) 

curve(fv, from=0, to=1000) 

는 내가 얻을 :

Error in y^2 : 'y' is missing 

Answer 1

lattice 패키지에는 wireframe() 및 persp()을 포함하여 3 차원 그래프를 그릴 때 도움이 될 수있는 몇 가지 기능이 있습니다. 3D 플롯을 사용하지 않으려면 contour()을 사용하여 등고선 플롯을 생성 할 수 있습니다.

참고 : 의도적인지 여부는 알 수 없지만 데이터가 플롯의 한쪽 구석에 매우 큰 스파이크를 생성합니다. 이는 한 모서리에 거의 눈에 띄지 않는 스파이크가있는 모든 평면에 대한 플롯을 생성합니다. 이것은 아래 윤곽 플롯에서 특히 문제가됩니다.

library(lattice) 

x <- seq(0, 1000, length.out = 50) 
y <- seq(0, 1000, length.out = 50) 

먼저 와이어 프레임 플롯 :

df <- expand.grid(x=x, y=y) 
df$z <- with(df, f(x, y)) 
wireframe(z ~ x * y, data = df) 

다음 관점 플롯 :

dm <- outer(x, y, FUN=f) 
persp(x, y, dm) 

윤곽 플롯 :

contour(x, y, dm) 

Answer 2

난 당신이 음모에 요구하는지 확실히 모르겠어요. 하지만 당신이 두 가지 인자의 스칼라 함수를 시각화하기를 원한다는 것을 알고 있습니다.

다음은 몇 가지 접근법입니다. 먼저 함수를 정의합니다.

llim <- 0 
ulim <- 6.29 

f <- function(x,y) { 
    (.156812/((2*pi)*(.000005^2)*(gamma(2/.156812)))*exp(-((sqrt(x^2+y^2))/.000005)^.156812)) 
} 

제목에서 나는 다음을 생각했습니다. intf 아래 정의 된 함수는 사각형 [0, ul] x [0, ul]에 함수를 통합하고 값을 반환합니다. 우리는 squareise에 대한 적분을 vectorise하고 plot하여 사각형의 변의 길이로 함수로 나타냅니다.

f는 분포 인 경우

intf <- function(ul) { 
    integrate(function(y) { 
    sapply(y, function(y) { 
     integrate(function(x) f(x,y), 0, ul)$value 
     }) 
    }, 0, ul)$value 
} 
fv <- Vectorize(intf) 
curve(fv, from=0, to=1000) 

, 난 당신이이 곡선의 당신의 (다소) 좋은 확률 해석을 ​​할 수 있습니다 같아요. (200 미터 x 200 미터의 수분의 20 % 확률 (?).

logf <- function(x, y) log(f(x, y)) 
x <- y <- seq(llim, ulim, length.out = 100) 
contour(x, y, outer(x, y, logf), lwd = 2, drawlabels = FALSE) 

당신은 또한 몇 가지 프로파일을 그릴 수 있습니다)

그러나, 당신은 또한 우리가 위의 통합 된 기능을 설명하는) 로그 변환 된 값의 윤곽 플롯을 (할 수 표면 :

plot(1, xlim = c(llim, ulim), ylim = c(0, 0.005), xlab = "x", ylab = "f") 
y <- seq(llim, ulim, length.out = 6) 
for (i in seq_along(y)) { 
    tmp <- function(x) f(x, y = y[i]) 
    curve(tmp, llim, ulim, add = TRUE, col = i) 
} 
legend("topright", lty = 1, col = seq_along(y), 
     legend = as.expression(paste("y = ",y))) 

그들은 출판을 가치있게 만들기 위해 조금 수정해야하지만 아이디어를 얻으실 수 있습니다. 마지막으로, 다른 사람들이 제시 한 것처럼 3 차원 도표를 만들 수 있습니다.

편집 귀하의 의견에 따라, 당신은 같은 것을 할 수 있습니다

# Define the function times radius (this time with general a and b) 
# The default of a and b is as before 
g <- function(z, a = 5e-6, b = .156812) { 
    z * (b/(2*pi*a^2*gamma(2/b)))*exp(-(z/a)^b) 
} 

# A function that integrates g from 0 to Z and rotates 
# As g is not dependent on the angle we just multiply by 2pi 
intg <- function(Z, ...) { 
    2*pi*integrate(g, 0, Z, ...)$value 
} 

# Vectorize the Z argument of intg 
gv <- Vectorize(intg, "Z") 

# Plot 
Z <- seq(0, 1000, length.out = 100) 
plot(Z, gv(Z), type = "l", lwd = 2) 
lines(Z, gv(Z, a = 5e-5),   col = "blue", lwd = 2) 
lines(Z, gv(Z, b = .150),   col = "red", lwd = 2) 
lines(Z, gv(Z, a = 1e-4, b = .2), col = "orange", lwd = 2) 

당신은 다음 원하는 a 및 b에 대한 곡선을 그릴 수 있습니다. 둘 중 하나를 지정하지 않으면 기본값이 사용됩니다. 면책 조항 : 나의 미적분은 녹슬었고 방금 내 머리 꼭대기에서 벗어났습니다. 축 주위로 함수의 회전을 제대로했는지 확인해야합니다.