title: Integration inconsistency with adaptIntegrate
body: <p>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 <a href="http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2011.00117.x/full" rel="nofollow noreferrer">this (open access) paper</a>. 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.</p>
  
  <p>I've implemented this with <code>adaptIntegrate</code> from the <code>cubature</code> 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.</p>
  
  <p><img src="https://i.stack.imgur.com/z1o4w.png" alt="enter image description here"></p>
  
  <pre><code>library(cubature)
  f &lt;- function(xmin, xmax, ymin, ymax) {
    adaptIntegrate(function(x) {
      r &lt;- sqrt((x[3] - x[1])^2 + (x[4] - x[2])^2)
      ifelse(r &gt; 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
  </code></pre>
  
  <p>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.</p>
  
  <p><img src="https://i.stack.imgur.com/zRtYw.png" alt="enter image description here"></p>
  
  <pre><code>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
  </code></pre>
  
  <p>Why are these integrals so different (<code>0.01949567</code> vs <code>0.01016105</code>)? 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?</p>
  
  <p>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.</p>
  
  <p><hr />
  <strong>EDIT:</strong>
  For the <code>A -&gt; B</code> case, nested <code>integrate</code> returns a solution similar to the first <code>adaptIntegrate</code> solution. For the <code>A -&gt; C</code> case, it returns <code>Error in integrate(function(ky) { : the integral is probably divergent</code>.</p>
  
  <pre><code>g &lt;- function(Bx, By, Ax, Ay) {
    r &lt;- sqrt((Ax - Bx)^2 + (Ay - By)^2)
    ifelse(r &gt; 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
  </code></pre>