在R中计算每个周期的方差

Question

我正在处理一组如下所示的数据:

team runs_scored       date
LAN           3        2014-03-22
ARI           1        2014-03-22
LAN           7        2014-03-23
ARI           5        2014-03-23
LAN           1        2014-03-30
SDN           3        2014-03-30

我试图在这个集合上测试一个预测模型,其中一个输入参数是runs_scoredt-1 的方差.换句话说,为了预测第四次观察的结果变量,我需要LAN基于数据集中先前观察的方差.

我可以计算累积平均值和总和,但我无法弄清楚如何计算数据集中的累积方差.我正在进行大部分数据操作dplyr,但是如果它能得到我需要的东西,我不反对使用替代解决方案

Answer 1

将方差公式写为,(sum(x^2)-length(x)*mean(x)^2)/(length(x)-1)您可以看到它可以很容易地推广到累积方差,只需将其中的每个函数替换为其累积版本(cummean来自dplyr).所以,

library(dplyr)
cum_var <- function(x){
    n <- 1:length(x)
    (cumsum(x^2)-n*cummean(x)^2)/(n-1)
}

与@ MrFlick的速度比较cumvar似乎令人鼓舞.

x <- rnorm(1e6)
all.equal(cum_var(x), cumvar(x))
#[1] TRUE
system.time(cumvar(x))[3]
elapsed 
   5.52 
system.time(cum_var(x))[3]
elapsed 
   0.04 

Answer 2

如果您想要累积方差,可以实现方差的在线算法.主要的好处是它可以线性扩展而不是指数扩展,就像迭代所有可能的子集一样.

如果你有

x<-c(3,1,7,5,1,3)

你可以做

cumvar<-function(x) {
   tail(Reduce(local({mm<-0; nn<-0; function(a,b) 
        {nn<<-nn+1; d<-b-mm; mm<<-mm+d/nn; a+d*(b-mm)}}), 
        x, 0, accumulate=TRUE), -1)/(seq_along(x)-1)
}
cumvar(x)
# [1]       NaN 24.500000 14.333333 10.000000  7.700000  6.166667  5.333333   4.696429  4.111111  3.777778

返回相同的结果

cumvar2 <- function(x)  {
    sapply(seq_along(x), function(i) var(x[1:i]))
}
cumvar2(x)
# [1]        NA 24.500000 14.333333 10.000000  7.700000  6.166667  5.333333  4.696429  4.111111  3.777778

我们可以将效率与效率进行比较

set.seed(15)
x<-rpois(100, 5)
microbenchmark:::microbenchmark(cumvar(x), cumvar2(x))

# Unit: microseconds
#        expr      min        lq      mean   median       uq      max neval cld
#   cumvar(x)  272.502  297.2425  335.2058  315.490  339.625  957.728   100  a 
#  cumvar2(x) 1672.323 1793.0960 2089.8104 1865.838 1956.208 6386.863   100   b

但是如果你想使用这个算法,我建议你阅读维基页面,如果你只计算方差一,那么双通法更健壮.

你可以用用它dplyr与

dd<-read.table(text="team runs_scored       date
LAN           3        2014-03-22
ARI           1        2014-03-22
LAN           7        2014-03-23
ARI           5        2014-03-23
LAN           1        2014-03-30
SDN           3        2014-03-30", header=T)

dd %>% mutate(cvar=lag(cumvar(runs_scored)))

#   team runs_scored       date     cvar
# 1  LAN           3 2014-03-22       NA
# 2  ARI           1 2014-03-22      NaN
# 3  LAN           7 2014-03-23 2.000000
# 4  ARI           5 2014-03-23 9.333333
# 5  LAN           1 2014-03-30 6.666667
# 6  SDN           3 2014-03-30 6.800000