Gaussian mean estimation in simulated data sets

Set up environment

Load the ggplot2 and cowplot packages, and the functions definining the mean and variances used to simulate the data.

library(ggplot2)
library(cowplot)
source("../code/signals.R")

Load results

Load the results of the simulation experiments.

load("../output/gaus-dscr.RData")

Simulated data with constant variances

This plot reproduces Fig. 2 of the manuscript comparing the accuracy of estimated mean curves in the data sets simulated from the “Spikes” mean function with constant variance.

First, extract the results used to generate this plot.

homo.data.smash <-
  res[res$.id    == "sp.3.v1" &
      res$method == "smash.s8",]
homo.data.smash.homo <-
  res[res$.id    == "sp.3.v1" &
      res$method == "smash.homo.s8",]
homo.data.tithresh <-
  res[res$.id == "sp.3.v1" &
      res$method == "tithresh.homo.s8",]
homo.data.ebayes <-
  res[res$.id    == "sp.3.v1" &
      res$method == "ebayesthresh",]
homo.data.smash.true <-
  res[res$.id == "sp.3.v1" &
  res$method  == "smash.true.s8",]
homo.data <-
  res[res$.id == "sp.3.v1" &
  (res$method == "smash.s8" |
   res$method == "ebayesthresh" |
   res$method == "tithresh.homo.s8"),]

Transform these results into a data frame suitable for ggplot2.

pdat <-
  rbind(data.frame(method      = "smash",
                   method.type = "est",
                   mise        = homo.data.smash$mise),
        data.frame(method      = "smash.homo",
                   method.type = "homo",
                   mise        = homo.data.smash.homo$mise),
        data.frame(method      = "tithresh",
                   method.type = "homo",
                   mise        = homo.data.tithresh$mise),
        data.frame(method      = "ebayesthresh",
                   method.type = "homo",
                   mise        = homo.data.ebayes$mise),
        data.frame(method      = "smash.true",
                   method.type = "true",
                   mise        = homo.data.smash.true$mise))
pdat <-
  transform(pdat,
            method = factor(method,
                            names(sort(tapply(pdat$mise,pdat$method,mean),
                                       decreasing = TRUE))))

Create the combined boxplot and violin plot using ggplot2.

p <- ggplot(pdat,aes(x = method,y = mise,fill = method.type)) +
     geom_violin(fill = "skyblue",color = "skyblue") +
     geom_boxplot(width = 0.15,outlier.shape = NA) +
     scale_y_continuous(breaks = seq(6,16,2)) +
     scale_fill_manual(values = c("darkorange","dodgerblue","gold"),
                       guide = FALSE) +
     coord_flip() +
     labs(x = "",y = "MISE") +
     theme(axis.line = element_blank(),
           axis.ticks.y = element_blank())
print(p)

From this plot, we see that three versions of SMASH outperformed EbayesThresh and TI thresholding.

Next, we compare the same methods in simulated data sets with heteroskedastic errors.

Simulated data with heteroskedastic errors: “Spikes” mean signal and “Clipped Blocks” variance

In this scenario, data sets were simulated using the “Spikes” mean function and the “Clipped Blocks” variance function. The next couple plots reproduce part of Fig. 3 in the manuscript.

This plot shows the mean function as a block line, and the +/- 2 standard deviations as orange lines:

t         <- (1:1024)/1024
mu        <- spikes.fn(t,"mean")
sigma.ini <- sqrt(cblocks.fn(t,"var"))
sd.fn     <- sigma.ini/mean(sigma.ini) * sd(mu)/3
par(cex.axis = 1,cex.lab = 1.25)
plot(mu,type = "l", ylim = c(-0.05,1),xlab = "position",ylab = "",
     lwd = 1.75,xaxp = c(0,1024,4),yaxp = c(0,1,4))
lines(mu + 2*sd.fn,col = "darkorange",lty = 5,lwd = 1.75)
lines(mu - 2*sd.fn,col = "darkorange",lty = 5,lwd = 1.75)

Extract the results from running the simulations.

hetero.data.smash <-
  res[res$.id == "sp.3.v5" & res$method == "smash.s8",]
hetero.data.smash.homo <-
  res[res$.id == "sp.3.v5" & res$method == "smash.homo.s8",]
hetero.data.tithresh.homo <-
  res[res$.id == "sp.3.v5" & res$method == "tithresh.homo.s8",]
hetero.data.tithresh.rmad <-
  res[res$.id == "sp.3.v5" & res$method == "tithresh.rmad.s8",]
hetero.data.tithresh.smash <-
  res[res$.id == "sp.3.v5" & res$method == "tithresh.smash.s8",]
hetero.data.tithresh.true <-
  res[res$.id == "sp.3.v5" & res$method == "tithresh.true.s8",]
hetero.data.ebayes <-
  res[res$.id == "sp.3.v5" & res$method == "ebayesthresh",]
hetero.data.smash.true <-
  res[res$.id == "sp.3.v5" & res$method == "smash.true.s8",]

Transform these results into a data frame suitable for ggplot2.

pdat <-
  rbind(data.frame(method      = "smash",
                   method.type = "est",
                   mise        = hetero.data.smash$mise),
        data.frame(method      = "smash.homo",
                   method.type = "homo",
                   mise        = hetero.data.smash.homo$mise),
        data.frame(method      = "tithresh.rmad",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.rmad$mise),
        data.frame(method      = "tithresh.smash",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.smash$mise),
        data.frame(method      = "tithresh.true",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.true$mise),
        data.frame(method      = "ebayesthresh",
                   method.type = "homo",
                   mise        = hetero.data.ebayes$mise),
        data.frame(method      = "smash.true",
                   method.type = "true",
                   mise        = hetero.data.smash.true$mise))
pdat <-
  transform(pdat,
            method = factor(method,
                            names(sort(tapply(pdat$mise,pdat$method,mean),
                                       decreasing = TRUE))))

Create the combined boxplot and violin plot using ggplot2.

p <- ggplot(pdat,aes(x = method,y = mise,fill = method.type)) +
     geom_violin(fill = "skyblue",color = "skyblue") +
     geom_boxplot(width = 0.15,outlier.shape = NA) +
     scale_fill_manual(values=c("darkorange","dodgerblue","limegreen","gold"),
                       guide = FALSE) +
     coord_flip() +
     scale_y_continuous(breaks = seq(10,70,10)) +
     labs(x = "",y = "MISE") +
     theme(axis.line = element_blank(),
           axis.ticks.y = element_blank())
print(p)

In the “Spikes” scenario, we see that SMASH, when allowing for heteroskedastic errors, outperforms EbayesThresh and all variants of TI thresholding (including TI thresholding with the true variance). Further, SMASH performs almost as well when estimating the variance compared to when provided with the true variance.

Simulated data with heteroskedastic errors: “Corner” mean signal and “Doppler” variance

In this next scenario, the data sets were simulated using the “Corner” mean function and the “Doppler” variance function. These plots were also used for Fig. 3 of the manuscript.

This plot shows the mean function as a block line, and the +/- 2 standard deviations as orange lines:

mu        <- cor.fn(t,"mean") 
sigma.ini <- sqrt(doppler.fn(t,"var"))
sd.fn     <- sigma.ini/mean(sigma.ini) * sd(mu)/3
plot(mu,type = "l", ylim = c(-0.05,1),xlab = "position",ylab = "",
     lwd = 1.75,xaxp = c(0,1024,4),yaxp = c(0,1,4))
lines(mu + 2*sd.fn,col = "darkorange",lty = 5,lwd = 1.75)
lines(mu - 2*sd.fn,col = "darkorange",lty = 5,lwd = 1.75)

Extract the results from running these simulations.

hetero.data.smash.2 <-
  res[res$.id == "cor.3.v3" & res$method == "smash.s8",]
hetero.data.smash.homo.2 <-
  res[res$.id == "cor.3.v3" & res$method == "smash.homo.s8",]
hetero.data.tithresh.homo.2 <-
  res[res$.id == "cor.3.v3" & res$method == "tithresh.homo.s8",]
hetero.data.tithresh.rmad.2 <-
  res[res$.id == "cor.3.v3" & res$method == "tithresh.rmad.s8",]
hetero.data.tithresh.smash.2 <-
  res[res$.id == "cor.3.v3" & res$method == "tithresh.smash.s8",]
hetero.data.tithresh.true.2 <-
  res[res$.id == "cor.3.v3" & res$method == "tithresh.true.s8",]
hetero.data.ebayes.2 <-
  res[res$.id == "cor.3.v3" & res$method == "ebayesthresh",]
hetero.data.smash.true.2 <-
  res[res$.id == "cor.3.v3" & res$method == "smash.true.s8",]

Transform these results into a data frame suitable for ggplot2.

pdat <-
  rbind(data.frame(method      = "smash",
                   method.type = "est",
                   mise        = hetero.data.smash.2$mise),
        data.frame(method      = "smash.homo",
                   method.type = "homo",
                   mise        = hetero.data.smash.homo.2$mise),
        data.frame(method      = "tithresh.rmad",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.rmad.2$mise),
        data.frame(method      = "tithresh.smash",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.smash.2$mise),
        data.frame(method      = "tithresh.true",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.true.2$mise),
        data.frame(method      = "ebayesthresh",
                   method.type = "homo",
                   mise        = hetero.data.ebayes.2$mise),
        data.frame(method      = "smash.true",
                   method.type = "true",
                   mise        = hetero.data.smash.true.2$mise))
pdat <-
  transform(pdat,
            method = factor(method,
                            names(sort(tapply(pdat$mise,pdat$method,mean),
                                       decreasing = TRUE))))

Create the combined boxplot and violin plot using ggplot2.

p <- ggplot(pdat,aes(x = method,y = mise,fill = method.type)) +
     geom_violin(fill = "skyblue",color = "skyblue") +
     geom_boxplot(width = 0.15,outlier.shape = NA) +
     scale_fill_manual(values=c("darkorange","dodgerblue","limegreen","gold"),
                       guide = FALSE) +
     coord_flip() +
     scale_y_continuous(breaks = seq(1,5)) +
     labs(x = "",y = "MISE") +
     theme(axis.line = element_blank(),
           axis.ticks.y = element_blank())
print(p)

Similar to the “Spikes” scenario, we see that the SMASH method, when allowing for heteroskedastic variances, outperforms both the TI thresholding and EbayesThresh approaches.

Session information

sessionInfo()
# R version 3.4.3 (2017-11-30)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS High Sierra 10.13.6
# 
# Matrix products: default
# BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
# 
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] cowplot_0.9.3 ggplot2_3.0.0
# 
# loaded via a namespace (and not attached):
#  [1] Rcpp_0.12.18      later_0.7.2       dscr_0.1-7       
#  [4] compiler_3.4.3    pillar_1.2.1      git2r_0.21.0     
#  [7] plyr_1.8.4        workflowr_1.1.1   bindr_0.1.1      
# [10] R.methodsS3_1.7.1 R.utils_2.6.0     tools_3.4.3      
# [13] digest_0.6.16     evaluate_0.10.1   tibble_1.4.2     
# [16] gtable_0.2.0      pkgconfig_2.0.1   rlang_0.2.1      
# [19] shiny_1.1.0       yaml_2.2.0        bindrcpp_0.2.2   
# [22] withr_2.1.2       stringr_1.3.0     dplyr_0.7.5      
# [25] knitr_1.20        rprojroot_1.3-2   grid_3.4.3       
# [28] tidyselect_0.2.4  glue_1.2.0        R6_2.2.2         
# [31] rmarkdown_1.9     purrr_0.2.5       magrittr_1.5     
# [34] whisker_0.3-2     promises_1.0.1    backports_1.1.2  
# [37] scales_0.5.0      htmltools_0.3.6   assertthat_0.2.0 
# [40] xtable_1.8-2      mime_0.5          colorspace_1.4-0 
# [43] httpuv_1.4.3      stringi_1.1.7     lazyeval_0.2.1   
# [46] munsell_0.4.3     R.oo_1.21.0