Визуализируйте функцию, используя двойную интеграцию в R - Wacky Result

Я не совсем уверен, что вы просите заговорить. Но я знаю, что вы хотите визуализировать свою скалярную функцию двух аргументов.

Вот некоторые подходы. Сначала мы определяем вашу функцию.

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

В пакете lattice есть несколько функций, которые помогут вам рисовать трехмерные графики, включая 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)