Integrating/Integral in R: Find the catch

Question

Ok, so this has me stumped for over more then 3 days and after still not being one step closer to the solution, I'm gonna try my luck here.

In the past, I've written some code for one sorted dataset in specific and it goes like this:

n <- length(data)
maxobs <- max(data)
minobs <- min(data)
FG <- function(m=NULL, h = NULL){
    n<- length(data) #Number of observations
    if (m<minobs){FG = 0} else { 
        if (m >maxobs){FG = 1} else {
                    FG = sum(pnorm((m - data)/h)-pnorm((minobs-data)/h))/sum(pnorm((maxobs - data)/h)-pnorm((minobs-data)/h))   
    }}
    return(FG)
}

f<- function(m,h){
    f<- FG(m,h)^n
    return(f)
}

##Integration
max <- NULL
delta<- function(h,max=maxobs){
    delta <- integrate(Vectorize(f), minobs, max, h)$value
    return (delta)
}

which works PERFECTLY. For example, if one selects data := c(1,2,3,4,5), one gets the correct result for

> delta(0.1, maxobs)
[1] 0.6300001

However, now I'm trying to generalize it for each sorted dataset, so what I did was (to be clear: the dataset x is sorted before exercising all these functions)

FG <- function(x, m=NULL, h = NULL){
  n<- length(x) #Number of observations
  maxobs <- max(x)
  minobs <- min(x)
  if (m<minobs){FG = 0} else { 
    if (m >maxobs){FG = 1} else {
      FG = sum(pnorm((m - x)/h)-pnorm((minobs-x)/h))/sum(pnorm((maxobs - x)/h)-pnorm((minobs-x)/h)) 
    }}
  return(FG)
}

f<- function(x,m,h){
  n <- length(x)
  f<- FG(x,m,h)^n
  return(f)
}

##Integration
delta<- function(x,h,maxu= max(x)){
minobs <- min(x)
  delta <- integrate(Vectorize(f), minobs, maxu, h)$value
  return (delta)
}

But now, delta(data,0.1) gives

delta(data,0.1)
[1] 0.

Which doesn't make any sense to me. Same function, same dataset, but now with a wrong value. What am I doing wrong?

Any help would be dearly appreciated.

EDIT: After taking a closer look to the Vectorize function and integrate function, I've now edited my delta function to:

delta<- function(x,h,maxu= max(x)){
minobs <- min(x)
  delta <- integrate(Vectorize(f, vectorize.args= c("m","h")), minobs, maxu, h)$value
  return (delta)
}

but this now just results in another error:

Error in integrate(Vectorize(f, vectorize.args = c("m", "h")), lower = minobs, : evaluation of function gave a result of wrong length

I thought Vectorize was supposed to prevent such errors?

Answer 1

The main problem here is that integrate expects you to have the variable you are integrating over as the first argument. In the first set of code, you are integrating over m. In the second set, you are attempting to integrate over x.

The shortest edit is to make a helper function to place the arguments in the necessary order for integrate:

delta<- function(x,h,maxu= max(x)){
  minobs <- min(x)
  g <- function(m) f(x,m,h)
  return( integrate(Vectorize(g), minobs, maxu)$value )
}

Now you'll get your desired results

delta(data,0.1)
# [1] 0.6300001

I believe the source of your second error was due to attempting to vectorize over h, whereas you really only wanted to vectorize over m. The helper function approach above eliminates this problem as well by only exposing the variable you wish to integrate over.

Note that I cannot tell what you are really trying to do here, but I'll also offer this rewrite that should be equivalent to your implementation, but is perhaps a bit easier to follow:

FG <- function(m, x, h) {
  n <- length(x)
  d <- function(t) pnorm((t-x)/h)

  if(m < x[1]) return(0)
  if(m > x[n]) return(1)

  return( sum(d(m)-d(x[1]))/sum(d(x[n])-d(x[1])) )
}

f<- function(m, x, h){
  n <- length(x)
  mapply(function(m) FG(m,x,h)^n, m)
}

delta<- function(x, h, lb=x[1], ub=x[length(x)]) {
  return( integrate(f, lb, ub, x, h)$value )
}