# Materialy do zajec z estymatorow
# http://www.biecek.pl/semestr/stats
#
# ---------------------------------------------------------------
# estymacja sredniej i wariancji z proby
#
ocen.mu.var <- function(x) {
  mu <- mean(x)
  S1 <- sum((x - mu)^2)/length(x)
  S2 <- sum((x - mu)^2)/(length(x)-1)
  list(mu = mu, S1 = S1, S2 = S2)
}

# zobaczmy czy S1 jest obciazonym estymatorem variancji
#
n <- 20
m <- 10^5
mu  <- 0
var <- 1
obs <- matrix(rnorm(n*m,mu,sqrt(var)),n,m)
S1 <- NULL
for (i in 1:m) {
    S1[i] <- (ocen.mu.var(obs[,i]))$S1
}

hist(S1,100)
lines(c(S1,S1),c(0,m),lwd=3)
# policz ocene obciazenia
mean(S1 - var)

# zobaczmy czy S2 jest obciazonym estymatorem variancji
#
S2 <- NULL
for (i in 1:m) {
    S2[i] <- (ocen.mu.var(obs[,i]))$S2
}

hist(S2,100)
lines(c(S2,S2),c(0,m),lwd=3)
# policz ocene obciazenia
mean(S2 - var)


# zobaczmy ktory estymator ma mniejszy blad sredniokwadratowy
#
mean((S1 - var)^2)
mean((S2 - var)^2)

#
# ktory estymator jest lepszy?
#