Fast Robust Moments in R with Rcpp
Switch branches/tags
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Permalink
Failed to load latest commit information.
R
docker
m4
man-roxygen
man
nodist
rpkg_make @ 628c25c
src
tests
tools/figure
.Rbuildignore
.gitignore
.gitmodules
.travis.yml
ChangeLog
DESCRIPTION
Makefile
NAMESPACE
README.Rmd
README.md

README.md

fromo

Build Status codecov.io CRAN Downloads Total RCpp is true

Fast Robust Moments -- Pick Three!

Fast, numerically robust, higher order moments in R, computed via Rcpp, mostly as an exercise to learn Rcpp. Supports computation on vectors and matrices, and Monoidal append (and unappend) of moments. Computations are via the Welford-Terriberry algorithm, as described by Bennett et al.

-- Steven E. Pav, shabbychef@gmail.com

Installation

This package can be installed from CRAN, via drat, or from github:

# via CRAN:
install.packages("fromo")
# via drat:
if (require(drat)) {
    drat:::add("shabbychef")
    install.packages("fromo")
}
# get snapshot from github (may be buggy)
if (require(devtools)) {
    install_github("shabbychef/fromo")
}

Basic Usage

Currently the package functionality can be divided into the following:

  • Functions which reduce a vector to an array of moments.
  • Functions which take a vector to a matrix of the running moments.
  • Functions which transform a vector to some normalized form, like a centered, rescaled, z-scored sample, or a summarized form, like the running Sharpe or t-stat.
  • Functions for computing the covariance of a vector robustly.
  • Object representations of moments with join and unjoin methods.

Summary moments

A function which computes, say, the kurtosis, typically also computes the mean and standard deviation, and has performed enough computation to easily return the skew. However, the default functions in R for higher order moments discard these lower order moments. So, for example, if you wish to compute Merten's form for the standard error of the Sharpe ratio, you have to call separate functions to compute the kurtosis, skew, standard deviation, and mean.

The summary functions in fromo return all the moments up to some order, namely the functions sd3, skew4, and kurt5. The latter of these, kurt5 returns an array of length 5 containing the excess kurtosis, the skewness, the standard deviation, the mean, and the observation count. (The number in the function name denotes the length of the output.) Along the same lines, there are summarizing functions that compute centered moments, standardized moments, and 'raw' cumulants:

  • cent_moments: return a k+1-vector of the kth centered moment, the k-1th, all the way down to the 2nd (the variance), then the mean and the observation count.
  • std_moments: return a k+1-vector of the kth standardized moment, the k-1th, all the way down to the 3rd, then the standard deviation, the mean, and the observation count.
  • cent_cumulants: computes the centered cumulants (yes, this is redundant, but they are not standardized). return a k+1-vector of the kth raw cumulant, the k-1th, all the way down to the second, then the mean, and the observation count.
  • std_cumulants: computes the standardized (and, of course, centered) cumulants. return a k+1-vector of the kth standardized cumulant, all the way down to the third, then the variance, the mean, and the observation count.
library(fromo)
set.seed(12345)
x <- rnorm(1000, mean = 10, sd = 2)
show(cent_moments(x, max_order = 4, na_rm = TRUE))
## [1]   47.276   -0.047    3.986   10.092 1000.000
show(std_moments(x, max_order = 4, na_rm = TRUE))
## [1]  3.0e+00 -5.9e-03  2.0e+00  1.0e+01  1.0e+03
show(cent_cumulants(x, max_order = 4, na_rm = TRUE))
## [1]   -0.388   -0.047    3.986   10.092 1000.000
show(std_cumulants(x, max_order = 4, na_rm = TRUE))
## [1] -2.4e-02 -5.9e-03  4.0e+00  1.0e+01  1.0e+03

Speed

In theory these operations should be just as fast as the default functions, but faster than calling multiple default functions. Here is a speed comparison of the basic moment computations:

library(fromo)
library(moments)
library(microbenchmark)

set.seed(123)
x <- rnorm(1000)

dumbk <- function(x) {
    c(kurtosis(x) - 3, skewness(x), sd(x), mean(x), 
        length(x))
}

microbenchmark(kurt5(x), skew4(x), sd3(x), dumbk(x), 
    dumbk(x), kurtosis(x), skewness(x), sd(x), mean(x))
## Unit: microseconds
##         expr   min    lq  mean median    uq max neval    cld
##     kurt5(x) 141.0 143.9 151.7    146 148.4 305   100     e 
##     skew4(x)  81.5  83.8  86.0     85  86.2 110   100   c   
##       sd3(x)  11.2  12.6  14.4     13  13.8  48   100  b    
##     dumbk(x) 193.2 204.4 216.4    209 216.3 449   200      f
##  kurtosis(x)  85.5  89.0  92.7     90  92.5 174   100   cd  
##  skewness(x)  86.4  90.1  95.4     91  93.3 205   100    d  
##        sd(x)  14.0  17.3  20.4     19  20.1  82   100  b    
##      mean(x)   3.8   4.6   6.4      5   5.5 109   100 a
x <- rnorm(1e+07, mean = 1e+12)

microbenchmark(kurt5(x), skew4(x), sd3(x), dumbk(x), 
    kurtosis(x), skewness(x), sd(x), mean(x), times = 10L)
## Unit: milliseconds
##         expr  min   lq mean median   uq  max neval  cld
##     kurt5(x) 1429 1456 1509   1482 1545 1694    10   c 
##     skew4(x)  806  809  856    819  917  938    10  b  
##       sd3(x)   87   87   95     91  103  109    10 a   
##     dumbk(x) 1710 1737 1797   1763 1872 1920    10    d
##  kurtosis(x)  828  832  894    878  935  998    10  b  
##  skewness(x)  799  801  860    808  914 1057    10  b  
##        sd(x)   48   49   52     50   53   61    10 a   
##      mean(x)   17   18   19     19   22   22    10 a

Weight! Weight!

Many of the methods now support the computation of weighted moments. There are a few options around weights: whether to check them for negative values, whether to normalize them to unit mean.

library(fromo)
library(moments)
library(microbenchmark)

set.seed(987)
x <- rnorm(1000)
w <- runif(length(x))

# no weights:
show(cent_moments(x, max_order = 4, na_rm = TRUE))
## [1] 2.9e+00 1.2e-02 1.0e+00 1.0e-02 1.0e+03
# with weights:
show(cent_moments(x, max_order = 4, wts = w, na_rm = TRUE))
## [1] 3.1e+00 4.1e-02 1.0e+00 1.3e-02 1.0e+03
# if you turn off weight normalization, the last
# element is sum(wts):
show(cent_moments(x, max_order = 4, wts = w, na_rm = TRUE, 
    normalize_wts = FALSE))
## [1]   3.072   0.041   1.001   0.013 493.941
# let's compare for speed!
x <- rnorm(1e+07)
w <- runif(length(x))

slow_sd <- function(x, w) {
    n0 <- length(x)
    mu <- weighted.mean(x, w = w)
    sg <- sqrt(sum(w * (x - mu)^2)/(n0 - 1))
    c(sg, mu, n0)
}
microbenchmark(sd3(x, wts = w), slow_sd(x, w))
## Unit: milliseconds
##             expr min  lq mean median  uq max neval cld
##  sd3(x, wts = w) 100 102  108    105 111 138   100  a 
##    slow_sd(x, w) 177 184  208    195 210 450   100   b

Monoid mumbo-jumbo

The as.centsums object performs the summary (centralized) moment computation, and stores the centralized sums. There is a print method that shows raw, centralized, and standardized moments of the ingested data. This object supports concatenation and unconcatenation. These should satisfy 'monoidal homomorphism', meaning that concatenation and taking moments commute with each other. So if you have two vectors, x1 and x2, the following should be equal: c(as.centsums(x1,4),as.centsums(x2,4)) and as.centsums(c(x1,x2),4). Moreover, the following should also be equal: as.centsums(c(x1,x2),4) %-% as.centsums(x2,4)) and as.centsums(x1,4). This is a small step of the way towards fast machine learning methods (along the lines of Mike Izbicki's Hlearn library).

Some demo code:

set.seed(12345)
x1 <- runif(100)
x2 <- rnorm(100, mean = 1)
max_ord <- 6L

obj1 <- as.centsums(x1, max_ord)
# display:
show(obj1)
##           class: centsums 
##     raw moments: 100 0.0051 0.09 -0.00092 0.014 -0.00043 0.0027 
## central moments: 0 0.09 -0.0023 0.014 -0.00079 0.0027 
##     std moments: 0 1 -0.086 1.8 -0.33 3.8
# join them together
obj1 <- as.centsums(x1, max_ord)
obj2 <- as.centsums(x2, max_ord)
obj3 <- as.centsums(c(x1, x2), max_ord)
alt3 <- c(obj1, obj2)
# it commutes!
stopifnot(max(abs(sums(obj3) - sums(alt3))) < 1e-07)
# unjoin them, with this one weird operator:
alt2 <- obj3 %-% obj1
alt1 <- obj3 %-% obj2
stopifnot(max(abs(sums(obj2) - sums(alt2))) < 1e-07)
stopifnot(max(abs(sums(obj1) - sums(alt1))) < 1e-07)

We also have 'raw' join and unjoin methods, not nicely wrapped:

set.seed(123)
x1 <- rnorm(1000, mean = 1)
x2 <- rnorm(1000, mean = 1)
max_ord <- 6L
rs1 <- cent_sums(x1, max_ord)
rs2 <- cent_sums(x2, max_ord)
rs3 <- cent_sums(c(x1, x2), max_ord)
rs3alt <- join_cent_sums(rs1, rs2)
stopifnot(max(abs(rs3 - rs3alt)) < 1e-07)

rs1alt <- unjoin_cent_sums(rs3, rs2)
rs2alt <- unjoin_cent_sums(rs3, rs1)
stopifnot(max(abs(rs1 - rs1alt)) < 1e-07)
stopifnot(max(abs(rs2 - rs2alt)) < 1e-07)

For multivariate input

There is also code for computing co-sums and co-moments, though as of this writing only up to order 2. Some demo code for the monoidal stuff here:

set.seed(54321)
x1 <- matrix(rnorm(100 * 4), ncol = 4)
x2 <- matrix(rnorm(100 * 4), ncol = 4)

max_ord <- 2L
obj1 <- as.centcosums(x1, max_ord, na.omit = TRUE)
# display:
show(obj1)
## An object of class "centcosums"
## Slot "cosums":
##          [,1]    [,2]   [,3]     [,4]    [,5]
## [1,] 100.0000  -0.093  0.045  -0.0046   0.046
## [2,]  -0.0934 111.012  4.941 -16.4822   6.660
## [3,]   0.0450   4.941 71.230   0.8505   5.501
## [4,]  -0.0046 -16.482  0.850 117.3456  13.738
## [5,]   0.0463   6.660  5.501  13.7379 100.781
## 
## Slot "order":
## [1] 2
# join them together
obj1 <- as.centcosums(x1, max_ord)
obj2 <- as.centcosums(x2, max_ord)
obj3 <- as.centcosums(rbind(x1, x2), max_ord)
alt3 <- c(obj1, obj2)
# it commutes!
stopifnot(max(abs(cosums(obj3) - cosums(alt3))) < 1e-07)
# unjoin them, with this one weird operator:
alt2 <- obj3 %-% obj1
alt1 <- obj3 %-% obj2
stopifnot(max(abs(cosums(obj2) - cosums(alt2))) < 1e-07)
stopifnot(max(abs(cosums(obj1) - cosums(alt1))) < 1e-07)

Running moments

Since an online algorithm is used, we can compute cumulative running moments. Moreover, we can remove observations, and thus compute moments over a fixed length lookback window. The code checks for negative even moments caused by roundoff, and restarts the computation to correct; periodic recomputation can be forced by an input parameter.

A demonstration:

library(fromo)
library(moments)
library(microbenchmark)

set.seed(1234)
x <- rnorm(20)

k5 <- running_kurt5(x, window = 10L)
colnames(k5) <- c("excess_kurtosis", "skew", "stdev", 
    "mean", "nobs")
k5
##       excess_kurtosis  skew stdev   mean nobs
##  [1,]             NaN   NaN   NaN -1.207    1
##  [2,]             NaN   NaN  1.05 -0.465    2
##  [3,]             NaN -0.34  1.16  0.052    3
##  [4,]          -1.520 -0.13  1.53 -0.548    4
##  [5,]          -1.254 -0.50  1.39 -0.352    5
##  [6,]          -0.860 -0.79  1.30 -0.209    6
##  [7,]          -0.714 -0.70  1.19 -0.261    7
##  [8,]          -0.525 -0.64  1.11 -0.297    8
##  [9,]          -0.331 -0.58  1.04 -0.327    9
## [10,]          -0.331 -0.42  1.00 -0.383   10
## [11,]           0.262 -0.65  0.95 -0.310   10
## [12,]           0.017 -0.30  0.95 -0.438   10
## [13,]           0.699 -0.61  0.79 -0.624   10
## [14,]          -0.939  0.69  0.53 -0.383   10
## [15,]          -0.296  0.99  0.64 -0.330   10
## [16,]           1.078  1.33  0.57 -0.391   10
## [17,]           1.069  1.32  0.57 -0.385   10
## [18,]           0.868  1.29  0.60 -0.421   10
## [19,]           0.799  1.31  0.61 -0.449   10
## [20,]           1.193  1.50  1.07 -0.118   10
# trust but verify
alt5 <- sapply(seq_along(x), function(iii) {
    rowi <- max(1, iii - 10 + 1)
    kurtosis(x[rowi:iii]) - 3
}, simplify = TRUE)

cbind(alt5, k5[, 1])
##         alt5       
##  [1,]    NaN    NaN
##  [2,] -2.000    NaN
##  [3,] -1.500    NaN
##  [4,] -1.520 -1.520
##  [5,] -1.254 -1.254
##  [6,] -0.860 -0.860
##  [7,] -0.714 -0.714
##  [8,] -0.525 -0.525
##  [9,] -0.331 -0.331
## [10,] -0.331 -0.331
## [11,]  0.262  0.262
## [12,]  0.017  0.017
## [13,]  0.699  0.699
## [14,] -0.939 -0.939
## [15,] -0.296 -0.296
## [16,]  1.078  1.078
## [17,]  1.069  1.069
## [18,]  0.868  0.868
## [19,]  0.799  0.799
## [20,]  1.193  1.193

See also

If you like rolling computations, do also check out the following packages (I believe they are all on CRAN):

Running 'scale' operations

Through template magic, the same code was modified to perform running centering, scaling, z-scoring and so on:

library(fromo)
library(moments)
library(microbenchmark)

set.seed(1234)
x <- rnorm(20)

xz <- running_zscored(x, window = 10L)

# trust but verify
altz <- sapply(seq_along(x), function(iii) {
    rowi <- max(1, iii - 10 + 1)
    (x[iii] - mean(x[rowi:iii]))/sd(x[rowi:iii])
}, simplify = TRUE)

cbind(xz, altz)
##              altz
##  [1,]   NaN    NA
##  [2,]  0.71  0.71
##  [3,]  0.89  0.89
##  [4,] -1.18 -1.18
##  [5,]  0.56  0.56
##  [6,]  0.55  0.55
##  [7,] -0.26 -0.26
##  [8,] -0.23 -0.23
##  [9,] -0.23 -0.23
## [10,] -0.51 -0.51
## [11,] -0.17 -0.17
## [12,] -0.59 -0.59
## [13,] -0.19 -0.19
## [14,]  0.84  0.84
## [15,]  2.02  2.02
## [16,]  0.49  0.49
## [17,] -0.22 -0.22
## [18,] -0.82 -0.82
## [19,] -0.64 -0.64
## [20,]  2.37  2.37

A list of the available running functions:

  • running_centered : from the current value, subtract the mean over the trailing window.
  • running_scaled: divide the current value by the standard deviation over the trailing window.
  • running_zscored: from the current value, subtract the mean then divide by the standard deviation over the trailing window.
  • running_sharpe: divide the mean by the standard deviation over the trailing window. There is a boolean flag to also compute and return the Mertens' form of the standard error of the Sharpe ratio over the trailing window in the second column.
  • running_tstat: compute the t-stat over the trailing window.
  • running_cumulants: computes cumulants over the trailing window.
  • running_apx_quantiles: computes approximate quantiles over the trailing window based on the cumulants and the Cornish-Fisher approximation.
  • running_apx_median: uses running_apx_quantiles to give the approximate median over the trailing window.

Lookahead

The functions running_centered, running_scaled and running_zscored take an optional lookahead parameter that allows you to peek ahead (or behind if negative) to the computed moments for comparing against the current value. These are not supported for running_sharpe or running_tstat because they do not have an idea of the 'current value'.

Here is an example of using the lookahead to z-score some data, compared to a purely time-safe lookback. Around a timestamp of 1000, you can see the difference in outcomes from the two methods:

set.seed(1235)
z <- rnorm(1500, mean = 0, sd = 0.09)
x <- exp(cumsum(z)) - 1

xz_look <- running_zscored(x, window = 301, lookahead = 150)
xz_safe <- running_zscored(x, window = 301, lookahead = 0)
df <- data.frame(timestamp = seq_along(x), raw = x, 
    lookahead = xz_look, lookback = xz_safe)

library(tidyr)
gdf <- gather(df, key = "smoothing", value = "x", -timestamp)

library(ggplot2)
ph <- ggplot(gdf, aes(x = timestamp, y = x, group = smoothing, 
    colour = smoothing)) + geom_line()
print(ph)

plot of chunk toy_zscore

Efficiency

We make every attempt to balance numerical robustness, computational efficiency and memory usage. As a bit of strawman-bashing, here we microbenchmark the running Z-score computation against the naive algorithm:

library(fromo)
library(moments)
library(microbenchmark)

set.seed(4422)
x <- rnorm(10000)

dumb_zscore <- function(x, window) {
    altz <- sapply(seq_along(x), function(iii) {
        rowi <- max(1, iii - window + 1)
        xrang <- x[rowi:iii]
        (x[iii] - mean(xrang))/sd(xrang)
    }, simplify = TRUE)
}

val1 <- running_zscored(x, 250)
val2 <- dumb_zscore(x, 250)
stopifnot(max(abs(val1 - val2), na.rm = TRUE) <= 1e-14)

microbenchmark(running_zscored(x, 250), dumb_zscore(x, 
    250))
## Unit: microseconds
##                     expr    min     lq   mean median     uq    max neval cld
##  running_zscored(x, 250)    664    696    740    725    768    960   100  a 
##      dumb_zscore(x, 250) 224253 245629 257002 252287 258771 351249   100   b

Timing against the roll package

More seriously, here we compare the running_sd3 function, which computes the standard deviation, mean and number of elements with the roll_sd and roll_mean functions from the roll package.

# dare I?
library(fromo)
library(microbenchmark)
library(roll)

set.seed(4422)
x <- rnorm(1e+05)
xm <- matrix(x)

v1 <- running_sd3(xm, 250)
rsd <- roll::roll_sd(xm, 250)
rmu <- roll::roll_mean(xm, 250)
# compute error on the 1000th row:
stopifnot(max(abs(v1[1000, ] - c(rsd[1000], rmu[1000], 
    250))) < 1e-14)
# now timings:
microbenchmark(running_sd3(xm, 250), roll::roll_mean(xm, 
    250), roll::roll_sd(xm, 250))
## Unit: milliseconds
##                      expr  min   lq mean median uq  max neval cld
##      running_sd3(xm, 250)  7.4  7.7  7.8    7.7  8  9.2   100 a  
##  roll::roll_mean(xm, 250) 13.0 13.2 13.5   13.4 13 16.1   100  b 
##    roll::roll_sd(xm, 250) 35.2 35.6 36.0   35.9 36 40.3   100   c

OK, that's not a fair comparison: roll_mean is optimized to work columwise on a matrix. Let's unbash this strawman. I create a function using fromo::running_sd3 to compute a running mean or running standard deviation columnwise on a matrix, then compare that to roll_mean and roll_sd:

library(fromo)
library(microbenchmark)
library(roll)

set.seed(4422)
xm <- matrix(rnorm(2e+05), ncol = 100)
fromo_sd <- function(x, wins) {
    apply(x, 2, function(xc) {
        running_sd3(xc, wins)[, 1]
    })
}
fromo_mu <- function(x, wins) {
    apply(x, 2, function(xc) {
        running_sd3(xc, wins)[, 2]
    })
}
wins <- 500
v1 <- fromo_sd(xm, wins)
rsd <- roll::roll_sd(xm, wins, min_obs = 3)

v2 <- fromo_mu(xm, wins)
rmu <- roll::roll_mean(xm, wins)
# compute error on the 1000th row:
stopifnot(max(abs(v1[1000, ] - rsd[1000, ])) < 1e-14)
stopifnot(max(abs(v2[1000, ] - rmu[1000, ])) < 1e-14)

# now timings: note fromo_mu and fromo_sd do
# exactly the same work, so only time one of them
microbenchmark(fromo_sd(xm, wins), roll::roll_mean(xm, 
    wins), roll::roll_mean(xm, wins, parallel_for = "cols"), 
    roll::roll_mean(xm, wins, parallel_for = "rows"), 
    roll::roll_sd(xm, wins), times = 50L)
## Unit: milliseconds
##                                              expr min  lq mean median  uq max neval  cld
##                                fromo_sd(xm, wins)  20  21   22     22  22  26    50 a   
##                         roll::roll_mean(xm, wins)  45  46   46     46  46  49    50  b  
##  roll::roll_mean(xm, wins, parallel_for = "cols")  46  47   47     47  47  50    50   c 
##  roll::roll_mean(xm, wins, parallel_for = "rows")  45  46   46     46  46  52    50  b  
##                           roll::roll_sd(xm, wins) 123 123  124    124 124 129    50    d

This was somewhat unexpected. I did not think the fromo functions would be faster than roll_sd, much less roll_mean. From these benchmarks, I suspected that my computer was not taking advantage of parallelization, since the different calls to roll_mean had the same timings. To check, I timed roll_mean for the same data size, but for larger rolling windows. During the following timing, I confirmed that all eight cores on my laptop were at 100% utilization. I believe the problem is that roll_mean is literally recomputing the moments over the entire window for every cell of the output, instead of reusing computations, which fromo mostly does:

library(roll)
library(microbenchmark)
set.seed(91823)
xm <- matrix(rnorm(2e+05), ncol = 10)
fromo_mu <- function(x, wins, ...) {
    apply(x, 2, function(xc) {
        running_sd3(xc, wins, ...)[, 2]
    })
}

microbenchmark(roll::roll_mean(xm, 10, min_obs = 3), 
    roll::roll_mean(xm, 100, min_obs = 3), roll::roll_mean(xm, 
        1000, min_obs = 3), roll::roll_mean(xm, 10000, 
        min_obs = 3), fromo_mu(xm, 10, min_df = 3), 
    fromo_mu(xm, 100, min_df = 3), fromo_mu(xm, 1000, 
        min_df = 3), fromo_mu(xm, 10000, min_df = 3), 
    times = 100L)
## Unit: milliseconds
##                                     expr   min    lq  mean median    uq   max neval     cld
##     roll::roll_mean(xm, 10, min_obs = 3)   1.5   1.6   1.8    1.8   1.9   3.4   100 a      
##    roll::roll_mean(xm, 100, min_obs = 3)  10.6  10.9  11.3   11.2  11.4  14.4   100  b     
##   roll::roll_mean(xm, 1000, min_obs = 3) 100.6 101.5 103.5  102.1 102.9 128.2   100     e  
##  roll::roll_mean(xm, 10000, min_obs = 3) 772.3 776.0 785.9  778.7 784.7 865.9   100       g
##             fromo_mu(xm, 10, min_df = 3)  11.8  13.0  13.9   14.2  14.5  16.5   100  bc    
##            fromo_mu(xm, 100, min_df = 3)  13.3  14.6  15.8   16.1  16.4  26.1   100   c    
##           fromo_mu(xm, 1000, min_df = 3)  28.2  30.7  32.6   33.5  34.1  36.5   100    d   
##          fromo_mu(xm, 10000, min_df = 3)  97.7 102.9 108.9  109.3 112.6 129.2   100      f

The runtime for operations from roll grow with the window size. The equivalent operations from fromo appear to also consume more time. In theory they would be invariant with respect to window size, but I coded them to 'restart' the computation periodically for improved accuracy. The user has control over how often this happens, in order to balance speed and accuracy. Here I set that parameter very large to show that runtimes need not grow with window size:

library(fromo)
library(microbenchmark)
set.seed(91823)
xm <- matrix(rnorm(2e+05), ncol = 10)
fromo_mu <- function(x, wins, ...) {
    apply(x, 2, function(xc) {
        running_sd3(xc, wins, ...)[, 2]
    })
}
rp <- 1L + nrow(xm)

microbenchmark(fromo_mu(xm, 10, min_df = 3, restart_period = rp), 
    fromo_mu(xm, 100, min_df = 3, restart_period = rp), 
    fromo_mu(xm, 1000, min_df = 3, restart_period = rp), 
    fromo_mu(xm, 10000, min_df = 3, restart_period = rp), 
    times = 100L)
## Unit: milliseconds
##                                                  expr min lq mean median uq max neval cld
##     fromo_mu(xm, 10, min_df = 3, restart_period = rp)  11 11   12     12 12  14   100   b
##    fromo_mu(xm, 100, min_df = 3, restart_period = rp)  11 11   12     11 12  15   100   b
##   fromo_mu(xm, 1000, min_df = 3, restart_period = rp)  11 11   12     12 12  15   100   b
##  fromo_mu(xm, 10000, min_df = 3, restart_period = rp)  10 10   11     11 11  13   100  a

Here are some more benchmarks, also against the rollingWindow package, for running sums:

library(microbenchmark)
library(fromo)
library(RollingWindow)
library(roll)

set.seed(12345)
x <- rnorm(10000)
xm <- matrix(x)
wins <- 1000

# run fun on each wins sized window...
silly_fun <- function(x, wins, fun, ...) {
    xout <- rep(NA, length(x))
    for (iii in seq_along(x)) {
        xout[iii] <- fun(x[max(1, iii - wins + 1):iii], 
            ...)
    }
    xout
}
vals <- list(running_sum(x, wins, na_rm = FALSE), RollingWindow::RollingSum(x, 
    wins, na_method = "ignore"), roll::roll_sum(xm, 
    wins), silly_fun(x, wins, sum, na.rm = FALSE))

# check all equal?
stopifnot(max(unlist(lapply(vals[2:length(vals)], function(av) {
    err <- vals[[1]] - av
    max(abs(err[wins:length(err)]), na.rm = TRUE)
}))) < 1e-12)

# benchmark it
microbenchmark(running_sum(x, wins, na_rm = FALSE), 
    RollingWindow::RollingSum(x, wins), running_sum(x, 
        wins, na_rm = TRUE), RollingWindow::RollingSum(x, 
        wins, na_method = "ignore"), roll::roll_sum(xm, 
        wins))
## Unit: microseconds
##                                                      expr  min   lq mean median   uq  max neval  cld
##                       running_sum(x, wins, na_rm = FALSE)   44   48   56     51   60  125   100 a   
##                        RollingWindow::RollingSum(x, wins)  111  117  203    126  146 1391   100  b  
##                        running_sum(x, wins, na_rm = TRUE)  107  118  129    120  137  202   100 ab  
##  RollingWindow::RollingSum(x, wins, na_method = "ignore")  359  393  474    404  432 1536   100   c 
##                                  roll::roll_sum(xm, wins) 3887 3907 4017   3941 4029 6515   100    d

And running means:

library(microbenchmark)
library(fromo)
library(RollingWindow)
library(roll)

set.seed(12345)
x <- rnorm(10000)
xm <- matrix(x)
wins <- 1000

vals <- list(running_mean(x, wins, na_rm = FALSE), 
    RollingWindow::RollingMean(x, wins, na_method = "ignore"), 
    roll::roll_mean(xm, wins), silly_fun(x, wins, mean, 
        na.rm = FALSE))

# check all equal?
stopifnot(max(unlist(lapply(vals[2:length(vals)], function(av) {
    err <- vals[[1]] - av
    max(abs(err[wins:length(err)]), na.rm = TRUE)
}))) < 1e-12)

# benchmark it:
microbenchmark(running_mean(x, wins, na_rm = FALSE, 
    restart_period = 1e+05), RollingWindow::RollingMean(x, 
    wins), running_mean(x, wins, na_rm = TRUE, restart_period = 1e+05), 
    RollingWindow::RollingMean(x, wins, na_method = "ignore"), 
    roll::roll_mean(xm, wins))
## Unit: microseconds
##                                                          expr  min   lq mean median   uq  max neval  cld
##  running_mean(x, wins, na_rm = FALSE, restart_period = 1e+05)   53   57   95     61   73 1584   100 a   
##                           RollingWindow::RollingMean(x, wins)  142  160  261    179  213 2202   100  b  
##   running_mean(x, wins, na_rm = TRUE, restart_period = 1e+05)  103  121  170    128  139 2016   100 ab  
##     RollingWindow::RollingMean(x, wins, na_method = "ignore")  375  443  571    463  536 2990   100   c 
##                                     roll::roll_mean(xm, wins) 4621 4654 5044   4797 5156 7729   100    d