Gains post-filter: simulations

First, we create the functions used in the simulations: generators (+ their means & variances), gains, filters, etc.

library(tidyverse)
library(rpart)
generator <- function(x, type, c){        # Generator functions: all except poly
    if(type == "linear"){return(c*x)}
    if(type == "quad"){return((x-c)^2)}
    if(type == "exp"){return(exp(-x*c))}
    if(type == "log"){return(log(x+c))}
    if(type == "sin"){return(sin(x*c)/c)}
}
m_gen <- function(type, c){               # Average value of generator
    if(type == "linear") return(c/2) 
    if(type == "quad") return(1/3-c+c^2)
    if(type == "exp") return((1-exp(-c))/c)
    if(type == "log") return(-1 - c*log(c)+(1+c)*log(1+c))
    if(type == "sin") return((1-cos(c))/c^2)
    if(type == "poly") {return(0)}        # Not computed 
}
var_gen <- function(type, c){             # Variance of generator
    if(type == "linear") return(c^2/12) 
    if(type == "quad") return((15*c^2-15*c+4)/45)
    if(type == "exp") return(exp(-2*c)*(exp(c)-1)*(2+c+exp(c)*(c-2))/2/(c^2))
    if(type == "log") return(1 - c*(1+c)*(log(c)-log(1+c))^2)
    if(type == "sin") return((2*c^2-c*sin(2*c)+8*cos(c)-2*cos(2*c)-6)/4/c^4)
}
fin_gen <- function(x, type,c){           # Generator with mean removed
    if(type == "poly"){
        return(gen_poly(x, c))
    } else
    return(generator(x, type, c) - m_gen(type, c))
}
sim_gen <- function(x, type, c, target_sd){ # Simulation of one point with noise
    return(fin_gen(x, type, c) + target_sd * rnorm(length(x))) 
}
gain <- function(fit){                      # Gain extracted from tree
    out_fit <- fit$frame
    var20 <- out_fit$dev[1]
    var21 <- out_fit$dev[2]
    var22 <- out_fit$dev[3]
    return(var20-(var21+var22))
}
filter_q <- function(data, q, y){           # Filtering function
    return(filter(data, y > quantile(y, 1-q) | y < quantile(y, q)))
}
gen_poly <- function(x, cc){                # Random function for the polynomial generator
    pars <- rnorm(cc) 
    while(sum(sign(pars))!=0){pars <- rnorm(cc) } # To make sure that the behavior is not too montonous
    m <- sum(pars/(2:(cc+1)))
    x <- sapply(x, `^`, 1:cc)
    p <- matrix(rep(pars, ncol(x)), nrow = cc, byrow = F)
    out <- colSums(p*x) - m
    return(out / sd(out) / runif(1, min = 12, max = 21))
}
ggplot(data.frame(x = c(0,1)), aes(x = x)) + stat_function(fun = gen_poly, args = list(cc = 6)) # Example of polynomial generator

Then, we create the simulation canvas & loop (via pmap to speed up the process). The simulations can take a lot of time. We chose modest parameters below.

nb_sim <- 50                                    # Number of simulations: 1000 in the paper
type <- c("linear", "quad", "sin", "exp", "log", "poly")
cc <- seq(0.1,0.4, length.out = nb_sim)         # For linear only, changes below
target_sd <- c(0.05)                             # Several values of noise can be tested
q <- c(0.15)                                    # Several values for q can also be tested
pars <- expand.grid(type, cc, target_sd,q)      # Initial parameters
colnames(pars) <- c("type", "cc", "target_sd", "q")
pars <- pars %>%                                # We alter parameter ranges below to get homogeneous gains
    group_by(target_sd, q) %>%
    mutate(cc = ifelse(type == "linear", seq(0.1,0.3, length.out = nb_sim), cc),
           cc = ifelse(type == "exp", seq(0.1,0.32, length.out = nb_sim), cc),
           cc = ifelse(type == "sin", seq(7,14, length.out = nb_sim), cc),
           cc = ifelse(type == "quad", seq(0.43,0.59, length.out = nb_sim), cc),
           cc = ifelse(type == "log", seq(3,13, length.out = nb_sim), cc),
           cc = ifelse(type == "poly", 6, cc))
type <- as.character(pars$type)                 # Parameter: type
cc <- matrix(pars$cc)                           # Parameter: parameter!
target_sd <- matrix(pars$target_sd)             # Parameter: level of noise
q <- matrix(pars$q)                             # Parameter: level of filter
gainz <- function(nb_points, target_sd, q, type, cc){ # Function that wraps the simulations
    x <- runif(nb_points)                             # Simulation of x
    y <- sim_gen(x, type, cc, target_sd)              # Simulation of y
    data <- data.frame(y,x)                           # Data (general case)
    fit <- rpart(y ~ x, data = data, cp = 0.0001, maxdepth = 1)                 # Tree fitting (general case)
    gain_res <- gain(fit)                             # Gain (benchmark)
    fit <- rpart(y ~ x, data = filter_q(data, q, y), cp = 0.0001, maxdepth = 1) # Tree on filtered data
    gain_res_filtered <- gain(fit)                    # Gain (post-filter)
    if(type == "poly"){                               # Below, we keep the level of generated variance  
        var_g <- var(y) - target_sd^2
    } else {
        var_g <- var_gen(type, cc)
    }
    return(c(gain_res, gain_res_filtered, var_g))
}
grd <- pmap(list(target_sd, q, type, cc),             # Run the simulations!
            gainz,                                    # On this function
            nb_points = 10^5)                         # Number of points simulated: 10^6 in the paper
res <- grd %>% unlist %>% matrix(ncol = 3, byrow = TRUE)                   # Organize results (1)
res <- data.frame(type = type, cc = cc, q = q, target_sd, gain = res[,1],  # Organize results (2)
                  gain_filtered = res[,2], var_g = res[,3])
# save(res, file = "res_agg.RData")

Finally, the plot.

res %>% 
    ggplot(aes(x = gain, y = gain_filtered, color = type)) + geom_point(size = 0.5) +
    scale_colour_manual("Legend", values = c("#ff0000","#ff0000","#ff0000","#787878","#00ccff","#9900cc")) 

The important feature is the red enveloppe: for a given initial gain, the best improvement is obtained by generators that are strictly increasing.