Визуализируйте функцию, используя двойную интеграцию в 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 имеет несколько функций, которые могут помочь вам сделать 3-мерные участки, в том числе и wireframe()persp(). Если вы предпочитаете не использовать 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] и возвращает значение. Затем мы прорисовываем и строим интеграл по квадрату как функцию длины стороны квадрата.

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) 

Если F- распределение, я думаю, вы можете сделать свой (несколько) Хорошая интерпретация вероятности этой кривой. (I. ~ ~ 20% вероятности опыления (?) В 200 на 200 метров площади.)

Однако, вы также можете сделать контур участок (из логарифмических значений), которые иллюстрируют функцию мы интегрирующие выше:

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

Их нужно немного изменить, чтобы сделать их публикации достойными, но вы получите эту идею. Наконец, вы можете сделать некоторые 3d-графики, как предложили другие.

EDIT В соответствии с вашими комментариями, вы также можете сделать что-то вроде этого:

# 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 вы хотите. Если либо не указано, используется значение по умолчанию. Отказ от ответственности: мое исчисление ржавое, и я просто сделал это с моей головы. Вы должны убедиться, что я правильно вращал функцию вокруг оси.