On-the-fly variance correction

stat <- function(v)
{
  list(mean = mean(v), sqsum = sum(v^2), var = var(v), sd = sd(v))
}

add <- function(l, n, x)
{
  l$mean <- (n * l$mean + x) / (n + 1)
  l$sqsum <- l$sqsum + x^2
  l$var <- (l$sqsum - (n + 1) * l$mean^2) / n
  l$sd <- sqrt(l$var)
  l
}

add_many <- function(l, n, x)
{
  k <- length(x)
  l$mean <- (n * l$mean + sum(x)) / (n + k)
  l$sqsum <- l$sqsum + sum(x^2)
  l$var <- (l$sqsum - (n + k) * l$mean^2) / (n + k - 1)
  l$sd <- sqrt(l$var)
  l
}

delete <- function(l, n, x)
{
  l$mean <- (n * l$mean - x) / (n - 1)
  l$sqsum <- l$sqsum - x^2
  l$var <- (l$sqsum - (n - 1) * l$mean^2) / (n - 2)
  l$sd <- sqrt(l$var)
  l
}

delete_many <- function(l, n, x)
{
  k <- length(x)
  if (k > n) stop()
  l$mean <- (n * l$mean - sum(x)) / (n - k)
  l$sqsum <- l$sqsum - sum(x^2)
  l$var <- (l$sqsum - (n - k) * l$mean^2) / (n - k - 1)
  l$sd <- sqrt(l$var)
  l
}

change <- function(l, n, x_out, x_in)
{
  l$mean <- l$mean + (x_in - x_out) / n
  l$sqsum <- l$sqsum + x_in^2 - x_out^2
  l$var <- (l$sqsum - n * l$mean^2) / (n - 1)
  l$sd <- sqrt(l$var)
  l
}

change_many <- function(l, n, x_out, x_in)
{
  if (length(x_in) != length(x_out)) stop()
  l$mean <- l$mean + (sum(x_in) - sum(x_out)) / n
  l$sqsum <- l$sqsum + sum(x_in^2) - sum(x_out^2)
  l$var <- (l$sqsum - n * l$mean^2) / (n - 1)
  l$sd <- sqrt(l$var)
  l
}

change_universal <- function(l, n, x_out, x_in)
{
  k_out <- length(x_out)
  k_in <- length(x_in)
  n_total <- (n + k_in - k_out)
  if (k_out > k_in + n) stop()
  l$mean <- (n * l$mean + sum(x_in) - sum(x_out)) / n_total
  l$sqsum <- l$sqsum + sum(x_in^2) - sum(x_out^2)
  l$var <- (l$sqsum - n_total * l$mean^2) / (n_total - 1)
  l$sd <- sqrt(l$var)
  l
}

v <- 1:100
l <- stat(v)
n <- length(v)
x <- 1:5 / 10

all.equal(stat(c(v, 0)), add(l, n, 0))
all.equal(stat(v[-1]), delete(l, n, 1))
all.equal(stat(replace(v, 1, 5)), change(l, n, 1, 5))
all.equal(stat(c(v, x)), add_many(l, n, x))
all.equal(stat(v[-(1:10)]), delete_many(l, n, 1:10))
all.equal(stat(replace(v, 1:10, 1:10/10)), change_many(l, n, 1:10, 1:10/10))

w <- v
w <- replace(w, 40:50, 300:310)
w[101:110] <- 201:210
w <- w[-(1:5)]
x_out <- c(40:50, 1:5)
x_in <- c(300:310, 201:210)

all.equal(stat(w), change_universal(l, n, x_out, x_in))
all.equal(stat(w), change_universal(l, n, rev(x_out), rev(x_in)))