Plots and tables summarizing results of Poisson simulations

Analysis settings

We will extract the results from these methods:

methods <- c("ash","BMSM","haarfisz_R")

Specify the row and column names for the tables:

table.row.names <- c("SMASH","BMSM","Haar-Fisz")
table.col.names <- c("(1/100,3)","(1/8,8)","(1/128,128)")

These are settings used in plotting the test functions:

n <- 1024
t <- 1:n/n

Set up environment

Add text here.

library(ggplot2)
library(cowplot)
library(xtable)

Some of the test functions are defined in signals.R:

source("../code/signals.R")

Load results

Load the results of the simulation experiments.

load("../output/pois.RData")

Spikes data

This is the function used to simulate the “Spikes” data sets at different ranges of intensities:

mu.s <- spike.f(t)
plot(t,mu.s,xlab = "",ylab = "",type = "l")

This table summarizes the results from the Spikes simulations:

mise.s.table   <- cbind(mise.s.1[methods],
                        mise.s.8[methods],
                        mise.s.128[methods])
rownames(mise.s.table) <- table.row.names
colnames(mise.s.table) <- table.col.names
print(xtable(mise.s.table,caption = " "),type = "html",
      html.table.attributes = "border=0")

Each column shows results at a different range of intensities. The individual table entries give the average error (MISE) in the estimates, in which the average is taken over the 100 data sets simulated at the given range of intensities.

The combined violin-boxplots provide a visualization of the same results:

m                  <- length(mise.ash.s.1)
method.labels      <- c("Haar-Fisz","BMSM","SMASH")
mise.hf.ti.r.s.1   <- colMeans(rbind(mise.hf.ti.r.4.s.1,
                                     mise.hf.ti.r.5.s.1,
                                     mise.hf.ti.r.6.s.1,
                                     mise.hf.ti.r.7.s.1))
mise.hf.ti.r.s.8   <- colMeans(rbind(mise.hf.ti.r.4.s.8,
mise.hf.ti.r.5.s.8,
                                     mise.hf.ti.r.6.s.8,
                                     mise.hf.ti.r.7.s.8))
mise.hf.ti.r.s.128 <- colMeans(rbind(mise.hf.ti.r.4.s.128,
                                     mise.hf.ti.r.5.s.128,
                                     mise.hf.ti.r.6.s.128,
                                     mise.hf.ti.r.7.s.128))
pdat1 <- data.frame(method = factor(rep(method.labels,each = m),method.labels),
                    mise   = c(mise.hf.ti.r.s.1,mise.BMSM.s.1,mise.ash.s.1))
pdat8 <- data.frame(method = factor(rep(method.labels,each = m),method.labels),
                    mise   = c(mise.hf.ti.r.s.8,mise.BMSM.s.8,mise.ash.s.8))
pdat128 <-
  data.frame(method = factor(rep(method.labels,each = m),method.labels),
             mise   = c(mise.hf.ti.r.s.128,mise.BMSM.s.128,mise.ash.s.128))
create.violin.plots <- function (pdat1, pdat8, pdat128) {
  p1 <- ggplot(pdat1,aes(x = method,y = mise)) +
        geom_violin(fill = "skyblue",color = "skyblue") +
        geom_boxplot(width = 0.15,outlier.shape = NA) +
        coord_flip() +
        labs(x = "",y = "MISE",title = "intensities (1/100,3)") +
        theme(axis.line = element_blank())
  p8 <- ggplot(pdat8,aes(x = method,y = mise)) +
        geom_violin(fill = "skyblue",color = "skyblue") +
        geom_boxplot(width = 0.15,outlier.shape = NA) +
        scale_x_discrete(breaks = NULL) +
        coord_flip() +
        labs(x = "",y = "MISE",title = "intensities (1/8,8)") +
        theme(axis.line = element_blank())
  p128 <- ggplot(pdat128,aes(x = method,y = mise)) +
        geom_violin(fill = "skyblue",color = "skyblue") +
        geom_boxplot(width = 0.15,outlier.shape = NA) +
        scale_x_discrete(breaks = NULL) +
        coord_flip() +
        labs(x = "",y = "MISE",title = "intensities (1/128,128)") +
        theme(axis.line = element_blank())
  return(plot_grid(p1,p8,p128,nrow = 1,ncol = 3))
}
create.violin.plots(pdat1,pdat8,pdat128)

Combine the results of the simulation experiments into several larger tables.

mise.ang.table <- cbind(mise.ang.1[methods],
                        mise.ang.8[methods],
                        mise.ang.128[methods])
mise.bur.table <- cbind(mise.bur.1[methods],
                        mise.bur.8[methods],
                        mise.bur.128[methods])
mise.cb.table  <- cbind(mise.cb.1[methods],
                        mise.cb.8[methods],
                        mise.cb.128[methods])
mise.b.table   <- cbind(mise.b.1[methods],
                        mise.b.8[methods],
                        mise.b.128[methods])
mise.hs.table  <- cbind(mise.hs.1[methods],
                        mise.hs.8[methods],
                        mise.hs.128[methods])
rownames(mise.ang.table) <- table.row.names
rownames(mise.b.table)   <- table.row.names
rownames(mise.cb.table)  <- table.row.names
rownames(mise.hs.table)  <- table.row.names
rownames(mise.bur.table) <- table.row.names
colnames(mise.ang.table) <- table.col.names
colnames(mise.b.table)   <- table.col.names
colnames(mise.cb.table)  <- table.col.names
colnames(mise.hs.table)  <- table.col.names
colnames(mise.bur.table) <- table.col.names

Angles data

This is the function used to simulate the “Angles” data sets at different ranges of intensities:

mu.ang <- dop.f(t)
sig <- ((2 * t + 0.5) * (t <= 0.15)) +
  ((-12 * (t - 0.15) + 0.8) * (t > 0.15 & t <= 0.2)) +
  0.2 * (t > 0.2 & t <= 0.5) + 
  ((6 * (t - 0.5) + 0.2) * (t > 0.5 & t <= 0.6)) +
  ((-10 * (t - 0.6) + 0.8) * (t > 0.6 & t <= 0.65)) +
  ((-0.5 * (t - 0.65) + 0.3) * (t > 0.65 & t <= 0.85)) +
  ((2 * (t - 0.85) + 0.2) * (t > 0.85))
mu.ang <- 3/5 * ((5/(max(sig) - min(sig))) * sig - 1.6) - 0.0419569
plot(t,mu.ang,xlab = "",ylab = "",type = "l")

This table summarizes the results from the Angles simulations:

mise.ang.table   <- cbind(mise.ang.1[methods],
                          mise.ang.8[methods],
                          mise.ang.128[methods])
rownames(mise.ang.table) <- table.row.names
colnames(mise.ang.table) <- table.col.names
print(xtable(mise.ang.table,caption = " "),type = "html",
      html.table.attributes = "border=0")

Each column shows results at a different range of intensities. The individual table entries give the average error (MISE) in the estimates, in which the average is taken over the 100 data sets simulated at the given range of intensities.

The combined violin-boxplots provide a visualization of the same results:

mise.hf.ti.r.ang.1   <- colMeans(rbind(mise.hf.ti.r.4.ang.1,
                                       mise.hf.ti.r.5.ang.1,
                                       mise.hf.ti.r.6.ang.1,
                                       mise.hf.ti.r.7.ang.1))
mise.hf.ti.r.ang.8   <- colMeans(rbind(mise.hf.ti.r.4.ang.8,
                                       mise.hf.ti.r.5.ang.8,
                                       mise.hf.ti.r.6.ang.8,
                                       mise.hf.ti.r.7.ang.8))
mise.hf.ti.r.ang.128 <- colMeans(rbind(mise.hf.ti.r.4.ang.128,
                                       mise.hf.ti.r.5.ang.128,
                                       mise.hf.ti.r.6.ang.128,
                                       mise.hf.ti.r.7.ang.128))
pdat1 <-
  data.frame(method = factor(rep(method.labels,each = m),method.labels),
             mise   = c(mise.hf.ti.r.ang.1,mise.BMSM.ang.1,mise.ash.ang.1))
pdat8 <-
  data.frame(method = factor(rep(method.labels,each = m),method.labels),
             mise   = c(mise.hf.ti.r.ang.8,mise.BMSM.ang.8,mise.ash.ang.8))
pdat128 <-
  data.frame(method=factor(rep(method.labels,each = m),method.labels),
             mise  =c(mise.hf.ti.r.ang.128,mise.BMSM.ang.128,mise.ash.ang.128))
create.violin.plots(pdat1,pdat8,pdat128)

Heavisine data

This is the function used to simulate the “Heavisine” data sets at different ranges of intensities:

heavi <- 4 * sin(4 * pi * t) - sign(t - 0.3) - sign(0.72 - t)
mu.hs <- heavi/sqrt(var(heavi)) * 1 * 2.99/3.366185
mu.hs <- mu.hs - min(mu.hs)
plot(t,mu.hs,xlab = "",ylab = "",type = "l")

This table summarizes the results from the Angles simulations:

mise.ang.table   <- cbind(mise.hs.1[methods],
                          mise.hs.8[methods],
                          mise.hs.128[methods])
rownames(mise.hs.table) <- table.row.names
colnames(mise.hs.table) <- table.col.names
print(xtable(mise.hs.table,caption = " "),type = "html",
      html.table.attributes = "border=0")

Each column shows results at a different range of intensities. The individual table entries give the average error (MISE) in the estimates, in which the average is taken over the 100 data sets simulated at the given range of intensities.

print(xtable(mise.hs.table,caption="Comparison of methods for denoising Poisson data for the ``Heavisine'' test function for 3 different (min,max) intensities ((0.01,3), (1/8,8), (1/128,128)). Performance is measured using MISE over 100 independent datasets, with smaller values indicating better performance. Values colored in red indicates the smallest MISE amongst all methods for a given (min, max) intensity.",label="table:pois_hs",digits=2),type = "html")

Bursts data

This is the function used to simulate the “Bursts” data sets at different ranges of intensities:

I_1 <- exp(-(abs(t - 0.2)/0.01)^1.2) * (t <= 0.2) +
       exp(-(abs(t - 0.2)/0.03)^1.2) * (t > 0.2)
I_2 <- exp(-(abs(t - 0.3)/0.01)^1.2) * (t <= 0.3) +
       exp(-(abs(t - 0.3)/0.03)^1.2) * (t > 0.3)
I_3 <- exp(-(abs(t - 0.4)/0.01)^1.2) * (t <= 0.4) +
       exp(-(abs(t - 0.4)/0.03)^1.2) * (t > 0.4)
mu.bur <- 2.99/4.51804 * (4 * I_1 + 3 * I_2 + 4.5 * I_3)
plot(t,mu.bur,xlab = "",ylab = "",type = "l")

print(xtable(mise.bur.table,caption="Comparison of methods for denoising Poisson data for the ``Bursts'' test function for 3 different (min,max) intensities ((0.01,3), (1/8,8), (1/128,128)). Performance is measured using MISE over 100 independent datasets, with smaller values indicating better performance. Values colored in red indicates the smallest MISE amongst all methods for a given (min, max) intensity.",label="table:pois_bur",digits=2),type = "html")

Clipped Blocks data

This is the function used to simulate the “Clipped Blocks” data sets at different ranges of intensities:

pos   <- c(0.1,0.13,0.15,0.23,0.25,0.4,0.44,0.65,0.76,0.78,0.81)
hgt   <- 2.88/5 * c(4,-5,3,-4,5,-4.2,2.1,4.3,-3.1,2.1,-4.2)
mu.cb <- rep(0,n)
for (j in 1:length(pos))
  mu.cb <- mu.cb + (1 + sign(t - pos[j])) * (hgt[j]/2)
mu.cb[mu.cb < 0] <- 0
plot(t,mu.cb,xlab = "",ylab = "",type = "l")

print(xtable(mise.cb.table,caption="Comparison of methods for denoising Poisson data for the ``Clipped Blocks'' test function for 3 different (min,max) intensities ((0.01,3), (1/8,8), (1/128,128)). Performance is measured using MISE over 100 independent datasets, with smaller values indicating better performance. Values colored in red indicates the smallest MISE amongst all methods for a given (min, max) intensity.",label="table:pois_cb",digits=2),type = "html")

Bumps data

This is the function used to simulate the “Bumps” data sets at different ranges of intensities:

pos  <- c(0.1,0.13,0.15,0.23,0.25,0.4,0.44,0.65,0.76,0.78,0.81)
hgt  <- 2.97/5 * c(4,5,3,4,5,4.2,2.1,4.3,3.1,5.1,4.2)
wth  <- c(0.005,0.005,0.006,0.01,0.01,0.03,0.01,0.01,0.005,0.008,0.005)
mu.b <- rep(0, n)
for (j in 1:length(pos))
  mu.b <- mu.b + hgt[j]/((1 + (abs(t - pos[j])/wth[j]))^4)
plot(t,mu.b,xlab = "",ylab = "",type = "l")

print(xtable(mise.b.table,caption="Comparison of methods for denoising Poisson data for the ``Bumps'' test function for 3 different (min,max) intensities ((0.01,3), (1/8,8), (1/128,128)). Performance is measured using MISE over 100 independent datasets, with smaller values indicating better performance. Values colored in red indicates the smallest MISE amongst all methods for a given (min, max) intensity.",label="table:pois_b",digits=2),type = "html")

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] xtable_1.8-2  cowplot_0.9.3 ggplot2_3.1.0
# 
# loaded via a namespace (and not attached):
#  [1] Rcpp_1.0.0        compiler_3.4.3    pillar_1.2.1     
#  [4] git2r_0.23.0      plyr_1.8.4        workflowr_1.1.1  
#  [7] bindr_0.1.1       R.methodsS3_1.7.1 R.utils_2.6.0    
# [10] tools_3.4.3       digest_0.6.17     evaluate_0.11    
# [13] tibble_1.4.2      gtable_0.2.0      pkgconfig_2.0.2  
# [16] rlang_0.2.2       yaml_2.2.0        bindrcpp_0.2.2   
# [19] withr_2.1.2       stringr_1.3.1     dplyr_0.7.6      
# [22] knitr_1.20        rprojroot_1.3-2   grid_3.4.3       
# [25] tidyselect_0.2.4  glue_1.3.0        R6_2.2.2         
# [28] rmarkdown_1.10    purrr_0.2.5       magrittr_1.5     
# [31] whisker_0.3-2     backports_1.1.2   scales_0.5.0     
# [34] htmltools_0.3.6   assertthat_0.2.0  colorspace_1.4-0 
# [37] labeling_0.3      stringi_1.2.4     lazyeval_0.2.1   
# [40] munsell_0.4.3     R.oo_1.21.0