n <- 200;
n0 <- 150;
x <- rnorm(n0);

# Standard variance estimate
var(x)/n0;

M <- 100000;
X <- matrix(x, nrow=n0, ncol=M)

sigmahat <- sd(x);
muhat <- mean(x);

# Xb <- matrix(sample(x, size=M*(n-n0), replace=TRUE), n-n0, M);
# Xb <- matrix(muhat+sigmahat*rnorm(M*(n-n0)), n-n0, M);
muhat.imp <-  t(matrix(muhat + rnorm(M)*sigmahat/sqrt(n0), M, n-n0));
Xb <- muhat.imp+ matrix(sigmahat*rnorm(M*(n-n0)), n-n0, M);

X.imp <- rbind(X, Xb)

Xbar.imp <- colSums(X.imp)/n;
V <- (1/n*colSums(X.imp^2) - Xbar.imp^2)*(n)/(n-1)/n
(Vbar <- mean(V));
# teoreticky 
var(x)/n


# Prakticky 
(B <- var(Xbar.imp));
# Teoreticky 
1/n^2*(n-n0)*mean((x-mean(x))^2)


# # 
Vbar + B
var(x)/n0

# Podil 
(Vbar + B*(1+1/M))/(var(x)/n0)