set.seed(1) # Določimo začetno točko vzorčenje.
lower <- numeric(10000) # Ustvarimo vektor, v katerega bomo shranili meje intervalnih ocen.
upper <- numeric(10000) # Ustvarimo vektor, v katerega bomo shranili meje intervalnih ocen.

n=100 # Določimo velikost vzorca

# Ponavljamo vzorčenje 100 enot 10 000 krat; za vsak vzorec ocenimo povprečje in intervalno oceno (alfa=0,05). Pri izračunu uporabimo oceno povprečja in oceno standardnega odklona, ki ga opazimo na vzorcu.
for(i in 1:10000) {
  
  Y <- runif(n, min=0, max=10)
  lower[i] <- mean(Y) + qt(0.025, df=n-1) * sd(Y)/sqrt(n)
  upper[i] <- mean(Y) + qt(0.975, df=n-1) * sd(Y)/sqrt(n)
  
}
# Združimo spodnjo in zgornjo mejo intervala v matriko
CIs <- cbind(lower, upper)
head(CIs)
##         lower    upper
## [1,] 4.647524 5.709417
## [2,] 4.635375 5.714061
## [3,] 3.797064 4.903242
## [4,] 4.455056 5.605279
## [5,] 4.445534 5.652940
## [6,] 4.782810 5.891401
# Izmed 10 000 vzorcev preverimo delež vzorcev, ki vsebuje pravo vrednost povprečja (Mu=5).
mean(CIs[, 1] <= 5 & 5 <= CIs[, 2])
## [1] 0.9479
# Poiščemo vzorce, ki med prvimi stotimi vzorci ne vsebujejo prave vrednosti povprečja
ID <- which(!(CIs[1:100, 1] <= 5 & 5 <= CIs[1:100, 2]))
ID
## [1]  3 27 37 69 96
# Določimo obliko grafa
plot(0, 
     xlim = c(0, 10), 
     ylim = c(1, 100), 
     ylab = "Vzorec", 
     xlab = expression(mu), 
     main = "Intervalne ocene za aritmetično sredino")

# S sivo označimo dobre intervalne ocene; tiste, ki zgrešijo, označimo z rdečo barvo
colors <- rep(gray(0.7), 100)
colors[ID] <- "red"

# Narišemo navpično črto pri pravi vrednosti parametra (Mu=5)
abline(v = 5, lty = 2)

# V graf vrišemo prvih 100 intervalnih ocen.
for(j in 1:100) {
  
    lines(c(CIs[j, 1], CIs[j, 2]), 
        c(j, j), 
        col = colors[j], 
        lwd = 2)
 
}