# Time-stamp: <2013-03-17 18:44 Petri Koistinen>

# Simuloidaan binomikokeen onnistumistodennäköisyyden estimaattorin
# otantajakaumaa


# Todelliset, eli simuloinnissa käytettävät parametrit
n <- 10
p.true <- 0.2013

# Riippumattomien toistojen lkm:
rep <- 1000

p.est <- numeric(rep)
for (r in 1:rep) {
  # yksi otos y_1, ..., y_n
  y <- rbinom(n, size = 1, prob = p.true)
  # lasketaan onnistumisten suhteellinen osuus
  p.est[r] <- mean(y)
}

# Kymmenen ensimmäistä estimaattorin otantajakaumasta simuloitua arvoa:

p.est[1:10]

# Estimaattorin frekvenssijakauma:
table(p.est) / rep

# Otoksesta arvioitu otantajakauman odotusarvo ja varianssi

mean(p.est)
var(p.est)

# Tarkat arvot otantajakauman odotusarvolle ja varianssille ovat
#
# p.true = 0.2013
# 1/n * p.true * (1 - p.true) = 0.01607783
#
# Yllä simuloiduista otoksista lasketut suureet suppenevat kohti näitä 
# arvoja, jos simulointien lukumäärää rep kasvatetaan rajatta.

# Tätä voitaisiin toistaa erilaisilla todellisilla parametrinarvoilla.