calcL <- function(x) {
  n <- length(x)
  (x - mean(x))^2 - (n-1)/n*var(x)
}
approxL <- function(x, eps=1e-6) {
  n <- length(x)
  theta <- (n-1)/n*var(x)
  L <- numeric(n)
  for (i in 1:n) {
    p <- rep(1/n,n)*(1-eps)
    p[i] <- p[i] + eps
    mean.i <- sum(p*x)
    L[i] <- (sum(p*(x-mean.i)^2)-theta)/eps    
  }
  L
}

x <- rnorm(100)
L <- calcL(x)
L.approx <- approxL(x, eps=1e-6)
mean(abs(L-L.approx))