In shabbychef/fromo: Fast Robust Moments

# set the knitr options ... for everyone!
# if you unset this, then vignette build bonks. oh, joy.
#opts_knit$set(progress=TRUE)
opts_knit$set(eval.after='fig.cap')
# for a package vignette, you do want to echo.
# opts_chunk$set(echo=FALSE,warning=FALSE,message=FALSE)
opts_chunk$set(warning=FALSE,message=FALSE)
#opts_chunk$set(results="asis")
opts_chunk$set(cache=TRUE,cache.path="cache/ref_timings_")

#opts_chunk$set(fig.path="figure/",dev=c("pdf","cairo_ps"))
opts_chunk$set(fig.path="figure/ref_timings_",dev=c("png"))
opts_chunk$set(fig.width=7,fig.height=6,dpi=100,out.width='700px',out.height='600px')

# doing this means that png files are made of figures;
# the savings is small, and it looks like shit:
#opts_chunk$set(fig.path="figure/",dev=c("png","pdf","cairo_ps"))
#opts_chunk$set(fig.width=4,fig.height=4)
# for figures? this is sweave-specific?
#opts_knit$set(eps=TRUE)

# this would be for figures:
#opts_chunk$set(out.width='.8\\textwidth')
# for text wrapping:
options(width=160,digits=2)
opts_chunk$set(size="small")
opts_chunk$set(tidy=TRUE,tidy.opts=list(width.cutoff=50,keep.blank.line=TRUE))
library(ggplot2)
library(fromo)
library(dplyr)
library(moments)
library(microbenchmark)

fromo timings against reference implementations

Let us check timings of simple sums:

library(fromo)
library(microbenchmark)

print(Sys.info())
print(sessionInfo())


set.seed(12345)
x    <- rnorm(1e5)
wins <- 250

ref_running_sum <- function(x,wins) {
    cx <- cumsum(x)
    cx - c(rep(0,wins),cx[1:(length(cx)-wins)])
}

ref_running_sum2 <- function(x,wins) {
    cx <- cumsum(x)
    c(cx[1:wins],cx[(wins+1):length(cx)] - cx[1:(length(cx)-wins)])
}


# check first
blah <- running_sum(x,wins) - ref_running_sum(x,wins)
print(summary(blah[4:length(blah)]))

# timings
checkit <- microbenchmark(sum(x),
                                                    mean(x),
                                                    gc(),
                                                    running_sum(x,wins),
                                                    running_sum(x,wins,na_rm=FALSE,restart_period=50000L),
                                                    running_mean(x,wins),
                                                    ref_running_sum(x,wins),
                                                    ref_running_sum2(x,wins),
                                                    cumsum(x))

print(checkit)

checkit %>%
    group_by(expr) %>%
        dplyr::summarize(meant=mean(time,na.rm=TRUE)) %>%
    ungroup() %>%
    dplyr::filter(grepl('running_sum',expr)) %>%
    mutate(timeover=meant / min(meant,na.rm=TRUE)) %>%
    kable()

Welford standard deviation is easy to compute quickly:

library(fromo)
library(microbenchmark)

print(Sys.info())
print(sessionInfo())


set.seed(12345)
x    <- rnorm(1e5)
wins <- 250

# check first
blah <- running_sd(x,wins) - ref_running_sd(x,wins)
print(summary(blah[4:length(blah)]))

# timings
checkit <- microbenchmark(sd(x),
                                                    sd3(x), 
                                                    ref_sd(x),
                                                    ref_sd_objecty(x),
                                                    running_sd(x,wins),
                                                    running_sd(x,wins,na_rm=FALSE,restart_period=50000L),
                                                    gc(),
                                                    ref_running_sd(x,wins),
                                                    ref_running_sd_narm(x,wins),
                                                    ref_running_sd_intnel(x,wins),
                                                    ref_running_sd_objecty(x,wins),
                                                    ref_running_sd_onecheck(x,wins),
                                                    ref_running_sd_fooz(x,wins),
                                                    ref_running_sd_barz(x,wins))


print(checkit)

checkit %>%
    group_by(expr) %>%
        dplyr::summarize(meant=mean(time,na.rm=TRUE)) %>%
    ungroup() %>%
    dplyr::filter(grepl('running_sd',expr)) %>%
    mutate(timeover=meant / min(meant,na.rm=TRUE)) %>%
    kable()