Несогласованность интеграции с адаптациейIntegrate

Я оцениваю вероятность расселения вида по ландшафту с координатной сеткой при заданном ядре расселения (функции расстояния) с максимальным расстоянием расселения. Я пытаюсь рассчитать вероятность рассредоточения от области к области, как описано в уравнении. 8 из этой статьи (в открытом доступе) . Это включает в себя четырехкратное интегрирование, оценивая значение функции рассеивания для каждой возможной комбинации исходной и целевой точек в исходной и целевой ячейках соответственно.

Я реализовал это с помощью adaptIntegrate из пакета cubature следующим образом для исходной ячейки A, целевой ячейки B и упрощенного ядра рассеивания, где рассеивание равно 1, когда межточечное расстояние > 1,25, и 0 в противном случае. Это показано графически ниже, где красная область ячейки B недостижима, поскольку ни одна точка в ячейке A не находится на расстоянии 1,25.

введите здесь описание изображения

library(cubature)
f <- function(xmin, xmax, ymin, ymax) {
  adaptIntegrate(function(x) {
    r <- sqrt((x[3] - x[1])^2 + (x[4] - x[2])^2)
    ifelse(r > 1.25, 0, 1)
  }, 
  lowerLimit=c(-0.5, -0.5, xmin, ymin), 
  upperLimit=c(0.5, 0.5, xmax, ymax), 
  maxEval=1e5)
}

f(xmin=1.5, xmax=2.5, ymin=-0.5, ymax=0.5)

# $integral
# [1] 0.01949567
# 
# $error
# [1] 0.001225998
# 
# $functionEvaluations
# [1] 100035
# 
# $returnCode
# [1] 0

Я получаю другой интеграл, когда рассматриваю целевую ячейку C, расположенную на том же расстоянии, но выше, а не справа от ячейки A.

введите здесь описание изображения

f(xmin=-0.5, xmax=0.5, ymin=1.5, ymax=2.5)

# $integral
# [1] 0.01016105
# 
# $error
# [1] 0.0241325
# 
# $functionEvaluations
# [1] 100035
# 
# $returnCode
# [1] 0

Почему эти интегралы такие разные (0.01949567 против 0.01016105)? Я неправильно закодировал? Изменение допуска и максимального количества оценок, по-видимому, не имеет большого значения. В качестве альтернативы, есть ли лучший подход к кодированию решения этой проблемы?

Я понимаю, что вопросы об общем подходе, вероятно, лучше подходят для stats.stackexchange.com, но я разместил здесь, так как подозреваю, что может быть что-то, что я упускаю из виду с самим кодированием.

EDIT: для случая A -> B вложенный integrate возвращает решение, аналогичное первому решению adaptIntegrate. Для случая A -> C возвращается Error in integrate(function(ky) { : the integral is probably divergent.

g <- function(Bx, By, Ax, Ay) {
  r <- sqrt((Ax - Bx)^2 + (Ay - By)^2)
  ifelse(r > 1.25, 0, 1)
}

integrate(function(Ay) {
  sapply(Ay, function(Ay) {
    integrate(function(Ax) {
      sapply(Ax, function(Ax) {
        integrate(function(By) {
          sapply(By, function(By) {
            integrate(function(Bx) g(Bx, By, Ax, Ay), 1.5, 2.5)$value # Bx
          })
        }, -0.5, 0.5)$value # By
      })
    }, -0.5, 0.5)$value # Ax
  })
}, -0.5, 0.5)$value # Ay

# [1] 0.019593

r numerical-integration
person jbaums    schedule 13.07.2014    source источник

Ответы (3)


arrow_upward
3
arrow_downward

Причина этого, по-видимому, заключается в том, как работает adaptIntegrate, поскольку, очевидно, единственное, что вы меняете, — это порядок интегрирования. Неидентичные результаты вероятны только из-за приблизительного интегрирования (см. первый ответ здесь), но это больше похоже на ошибку.

Вот значения r при вычислении f(xmin=1.5, xmax=2.5, ymin=-0.5, ymax=0.5)

введите здесь описание изображения

и f(xmin=-0.5, xmax=0.5, ymin=1.5, ymax=2.5)

введите здесь описание изображения

поэтому внутри функции должно что-то происходить, поскольку диапазон значений сильно различается.

Одной из альтернатив для этого является интеграция Монте-Карло, которая хороша в этом случае, поскольку ваши очки распределяются равномерно.

MCI <- function(Ax, Ay, Bx, By, N, r) {
  d <- sapply(list(Ax, Ay, Bx, By), function(l) runif(N, l[1], l[2]))
  sum(sqrt((d[, 1] - d[, 3])^2 + (d[, 2] - d[, 4])^2) <= r) / N
}

set.seed(123)
MCI(c(-0.5, 0.5), c(-0.5, 0.5), c(1.5, 2.5), c(-0.5, 0.5), 100000, 1.25)
# [1] 0.0194
MCI(c(-0.5, 0.5), c(-0.5, 0.5), c(-0.5, 0.5), c(1.5, 2.5), 100000, 1.25)
# [1] 0.01929
person Julius Vainora    schedule 15.07.2014
comment
Спасибо @Julius - похоже, интеграция Монте-Карло может быть подходящим вариантом. Не могли бы вы рассказать мне, как вы сохранили значения r? (Я могу cat записать их в файл, но есть ли более быстрый способ?) Я оставлю вопрос активным еще немного, на случай, если у других появится больше информации о причинах adaptIntegrate несоответствия. - person jbaums; 16.07.2014
comment
@jbaums, конечно, я поставил rs <- numeric(100035); cnt <- 1 перед функцией и rs[cnt] <<- r; cnt <<- cnt + 1 внутри функции. - person Julius Vainora; 16.07.2014

arrow_upward
2
arrow_downward

Обычно меры расстояния (x1-x2)^2+(y1-y2)^2. Можете ли вы объяснить, почему вы вычитаете x из y при построении r? Рассмотрим альтернативную версию:

f <- function(xmin, xmax, ymin, ymax) {
   adaptIntegrate(function(x) {
     r <- sqrt((x[4] - x[3])^2 + (x[2] - x[1])^2)
     ifelse(r > 1.25, 0, 1)
   }, 
   lowerLimit=c(-0.5, -0.5, xmin, ymin), 
   upperLimit=c(0.5, 0.5, xmax, ymax), 
   maxEval=1e5)
 }

 f(xmin=1.5, xmax=2.5, ymin=-0.5, ymax=0.5)
#-------------
$integral
[1] 0.01016105

$error
[1] 0.0241325

$functionEvaluations
[1] 100035

$returnCode
[1] 0
#---------
 f(xmin=-0.5, xmax=0.5, ymin=1.5, ymax=2.5)
#---------
$integral
[1] 0.01016105

$error
[1] 0.0241325

$functionEvaluations
[1] 100035

$returnCode
[1] 0
person IRTFM    schedule 13.07.2014
comment
Я считаю, что правильно рассчитываю расстояние, поскольку порядок элементов x определяется порядком элементов двух векторов *Limit. Они относятся к положению вдоль каждого из четырех измерений, т. е. когда мы движемся от -0,5 до 0,5 в первом и втором измерениях (представляющих границы координат x и y ячейки A соответственно), от xmin до xmax в третьем (x-границы ячейки B) и от ymin до ymax в четвертом (y-границы ячейки B). Моя функция-оболочка f просто передает xmin, xmax, ymin и ymax в adaptIntegrate. - person jbaums; 19.07.2014

arrow_upward
1
arrow_downward

Сопровождающий пакета R cubature (Нарас) сообщил мне, что библиотека Cubature C дает те же результаты, о которых я сообщаю в вопросе выше, и что это вряд ли является ошибкой; скорее, h-адаптивная кубатурная подпрограмма (для которой пакет R является интерфейсом) в некоторых случаях менее точна, чем p-адаптивная подпрограмма Cubature, которая удваивает количество точек выборки в соответствующих регионах.

Нарас также предоставил следующий код julia, который демонстрирует согласованные pcubature решения для двух случаев, представленных в моем вопросе (элементы возвращаемого значения - расчетный интеграл, за которым следует расчетная абсолютная ошибка).

using Cubature

# integrand
f = x -> ifelse(sqrt((x[3] - x[1])^2 + (x[4] - x[2])^2) > 1.25, 0, 1)

# A to B case
pcubature(f, [-0.5, -0.5, 1.5, -0.5], [0.5, 0.5, 2.5, 0.5], abstol=1e-5)    
# (0.019593408732917292,3.5592555263398717e-6)

# A to C case
pcubature(f, [-0.5, -0.5, -0.5, 1.5], [0.5, 0.5, 0.5, 2.5], abstol=1e-5)   
# (0.019593408732918302,3.559255527241928e-6)
person jbaums    schedule 19.07.2014