R: Scope in nested functions, part 2: nested do.call's

Question

Here I have a real question and a greatly simplified question that may or may not encapsulate my confusion about the first question.

The function multi_quad below is supposed to integrate dfun, here dgb2 from the GB2 package, over the limits from lower bound lb to upper bound ub. The base-R integrate function successfully does this integration for the parameter values and limits given here, but fails for ub=4.55e8 or more. I have a question posted about this problem, here: R: Convergence problems with numerical integration

I an effort to solve the convergence problem, I have the function multi_quad, below, which is supposed to attempt integration with various alternative quadrature methods from the pracma package. I am using safely( ) from the purrr package in an attempt to try all the methods, return values if they converge, and error messages otherwise.

For each iteration of map, I get an error that argument "lb" is missing, with no default. Here the both the enclosing and the calling environment of the outer do.call is the execution environment of m_quad, and both the enclosing and the calling environment of m_quad is the execution environment of multi_quad. Since I supply multi_quad with a value for lb, I expect m_quad and do.call to be able to find it, regardless of how lb in the outer do.call is scoped. This is incorrect, and I do not know why.

library(GB2)
library(pracma)
library(purr)

multi_quad <- function(dfun, params, lb, ub){
  m_quad <- function(q_method, dfun, params, lb, ub){
    do.call(what=q_method, 
            args=list(f=function(X){do.call(what=dfun, args=params)}, a=lb, b=ub))
    }
  safe_m_quad <- safely(m_quad)
  q_methods<- list("quadgk",  "quadcc",  "simpadpt", "renshaw", "quadgr") %>% 
    map(safe_m_quad)  -> out
  out
}

GB2_params <- list(shape1 = 3.652, scale = 65797, shape2 = 0.3, shape3 = 0.8356)

results <- multi_quad(dfun=dgb2, params=GB2_params, lb=1, ub=3e8)

results

Below is a vastly simplified example where I think my expectations of when and how an inner function can find a value supplied to a calling function are frustrated in a similar way. Here the variables with names starting in "extra" are to prevent positional matching.

In this example, the execution environment of fa is both the calling and the enclosing environment of fb. Why isn't the value of arg1, under that name, from fa in scope for arg1 of fb? Whether fb's arg1 is looking for arg1 in its calling environment or it's enclosing environment, or in fb's execution environment and up from there, in every case it should arrive in fa's execution environment. But it doesn't. Why is this? Where is it looking instead?

fa <- function(extra0=0, arg1, extra2=2){
  extra0
  extra2
  fb(arg1)
  fb <- function(extra3=3, arg1){
    extra3
    print(arg1)}  
}

fa(arg1=1)

Answer 1

(1) For me, integrate works up to 7.0e08. Think of it what huge integration domain this is. You're applying functions from pracma that are more meant as method demonstrations, not for pretentious tests. quadgkis similar to integrate and works up to 1.0e09, the integral function in pracma up to 1.0e10 !

If you want to integrate over the unbounded interval [0. Inf], then integrate will not work again. But integral recognizes the unboundedness and reverts to some other approach:

fun = function(x) dgb2(x, shape1 = 3.652, scale = 65797,
                          shape2 = 0.3, shape3 = 0.8356)
pracma::integral(fun, 0, Inf)
## For infinite domains Gauss integration is applied!
## [1] 1

It actually employs pracma::quadinf that transforms the infinite domain into a finite domain before applying a Gaussian Quadrature rule.

(2) In your simplified text example, the first bug is that you call fb(arg1) before having defined it. If you put fb(arg1) after the function definition, arg1 will be assigned to extra3 and the parameter arg1 is still undefined.

A possible correction could be:

fa <- function(extra0=0, arg1, extra2=2){
    fb <- function(extra3=3, arg1){
        print(arg1)
    }
    fb(arg1 = arg1)
}

fa(arg1=1)
## [1] 1

Or change fb to fb <- function(arg1, extra3=3) {...}. What the lines extra0, extra1, etc. are good for remains unclear.

(3) Added to provide some reasonable integral value:

Intgral functions integral and quadgk are able to return reasonable values for intervals up 5e08 and larger by selecting appropriate parameters.

I1 = integrate(fun, 0, 5.0e08, rel.tol=1e-14, subdivisions=1000)$value; I1 - 1
## [1] -4.792833e-13

I2 = integral(fun, 0, 5.0e08, reltol=1e-14); I2 -1
## [1] -4.775069e-13

I3 = quadgk(fun, 0, 5.0e08, tol=1e-14); I3 - 1
## [1] -4.773959e-13

So that we get the following integral values:

| Function          | Value             |
|-------------------|-------------------|
| stats::integrate  | 0.999999999999521 |
| pracma::integral  | 0.999999999999522 |
| pracma::quadgk    | 0.999999999999523 |

These functions use similar "adaptive Gaussian quadrature" rules, so no wonder they return similar results. Other competely different integration routines may come up with (slightly ?) different results.