I am estimating the probability of a species dispersing across a gridded landscape, given a dispersal kernel (a function of distance) with a maximum dispersal distance. I'm attempting to calculate area-to-area dispersal probabilities as described in eqn. 8 of [this (open access) paper](http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2011.00117.x/full). This involves quadruple integration, evaluating the value of the dispersal function for every possible combination of source and target point in the source and target cells, respectively.

I've implemented this with `adaptIntegrate` from the `cubature` package, as follows, for source cell A, target cell B, and a simplified dispersal kernel where dispersal is 1 when the inter-point distance > 1.25 and 0 otherwise. This is shown graphically below, where the red region of cell B is unreachable since no point in cell A is within a distance of 1.25.

![enter image description here](https://i.stack.imgur.com/z1o4w.png)

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

```

I get a different integral when considering a target cell, C, that is placed the same distance away, but above rather than to the right of cell A.

![enter image description here](https://i.stack.imgur.com/zRtYw.png)

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

```

Why are these integrals so different (`0.01949567` vs `0.01016105`)? Have I coded it incorrectly? Changing the tolerance and maximum number of evaluations appears to make no great difference. Alternatively, is there a better approach to coding a solution to this type of problem?

I realise that questions about the general approach are probably better suited to stats.stackexchange.com, but I've posted here since I suspect there may be something that I'm overlooking with the coding itself.

---

**EDIT:**
For the `A -> B` case, nested `integrate` returns a solution similar to the first `adaptIntegrate` solution. For the `A -> C` case, it returns `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

```