r_populacji <- 1000
populacja <- runif(r_populacji, 5, 17)
?apply
Następnie wygenerujemy 200 prób o rozmiarze 15:
r_proby <- 15
l_prob <- 200
losowe_indeksy <- sample(1:r_populacji, size=r_proby*l_prob, replace=TRUE)
proby <- matrix(populacja[losowe_indeksy], r_proby, l_prob)
colnames(proby) <- paste(rep("próba",ncol(proby)),c(1:ncol(proby)))
proby
Oblicz statystyki opisowe dla populacji:
pop_var <- function(data)
{
mean((data - mean(data))^2)
}
opis_popul <- c("średnia"=mean(populacja), "wariancja"=pop_var(populacja), "odchylenie standardowe"=sqrt(pop_var(populacja)))
opis_popul
Dla każdej próby wyznacz estymator średniej, estymator wariancji oraz obciążoną wersję estymatora wariancji. Możesz użyć funkcji apply.
estymator_srednia <- apply(proby, 2, mean)
estymator_wariancja <- apply(proby, 2, var)
obciazony_estymator_wariancja <- apply(proby, 2, pop_var)
Następnie dla każdego z estymatorów wyznacz wartość oczekiwaną (uśredniając po próbach) oraz obciążenie:
wiersze <- c("Wartość oczekiwana", "Obciążenie")
kolumny <- c("Estymator średnia", "Estymator wariancja", "Estymator obciążony wariancja")
wartosci_oczekiwane <- c(mean(estymator_srednia), mean(estymator_wariancja), mean(obciazony_estymator_wariancja))
obciazenia <- c(wartosci_oczekiwane[1]-opis_popul[1], wartosci_oczekiwane[2]-opis_popul[2], wartosci_oczekiwane[3]-opis_popul[2])
matrix(c(rbind(wartosci_oczekiwane, obciazenia)), length(wiersze), length(kolumny), dimnames=list(wiersze, kolumny))
Dla estymatora wartości średniej policz dodatkowo odchylenie standardowe i teoretyczne odchylenie standardowe:
c("Odchylenie standardowe"=sd(estymator_srednia),"Teoretyczne"=opis_popul[3]/sqrt(r_proby))
Utwórz histogram wartości estymatora średniej:
hist(estymator_srednia, main="Estymator średnia", xlab="Wartość estymatora")
Ćwiczenie 2: Badanie estymacji przedziałowej.¶
Na początek wygenerujemy 200 prób o rozmiarze 20 z rozkładu normalnego.
r_proby <- 15
l_prob <- 200
m <- 0
sigma <- 1
proby <- matrix(rnorm(r_proby*l_prob, m, sigma), r_proby, l_prob)
colnames(proby) <- paste(rep("próba",ncol(proby)),c(1:ncol(proby)))
proby
Wyznacz kolejno dla każdej próby estymator średniej (możesz użyć funkcji apply), lewy i prawy koniec przedziału, szerokość przedziału oraz wypisz czy prawdziwa średnia znalazła się w przedziale.
poziom_istotnosci <- 0.05
kwantyl <- qnorm(1-poziom_istotnosci/2, m, sigma)
estymator_srednia <- apply(proby, 2, mean)
l_koniec_przedzialu <- estymator_srednia - kwantyl*sigma/sqrt(r_proby)
p_koniec_przedzialu <- estymator_srednia + kwantyl*sigma/sqrt(r_proby)
szerokosc_przedzialu <- p_koniec_przedzialu - l_koniec_przedzialu
czy_srednia_w_przedziale <- (m > l_koniec_przedzialu) & (m < p_koniec_przedzialu)
round(matrix(c(rbind(estymator_srednia, l_koniec_przedzialu, p_koniec_przedzialu, szerokosc_przedzialu, czy_srednia_w_przedziale)), 5, l_prob, dimnames = list(c("Estymator średniej", "L koniec przedziału", "P koniec przedziału", "Szerokość przedziału", "Średnia w przedziale?"),colnames(proby))),3)
Następnie wyznacz prawdopodobieństwo uśredniając po 200 próbach. Czy zmiana m i sigma wpływa na to prawdopodobieństwo?
c("Pr(średnia w przedziale)"=sum(czy_srednia_w_przedziale)/l_prob)
Ćwiczenie 3: Badanie zgodności estymatorów średniej i wariancji¶
Zacznij od wygenerowania n losowych liczb z rozkładu normalnego:
n <- 1000
m <- 0
sigma <- 1
losowe <- rnorm(n, m, sigma)
Wyznacz wartości estymatora średniej oraz wariancji i porównaj ją z prawdziwymi wartościami dla rosnącej ilości zmiennych (możesz użyć funkcji sample). Błąd przyjmij w sensie kwadratu różnicy między wartością prawdziwą a estymowaną.
estymatory <- matrix(rep(0, (n-1)*4), n-1, 4, dimnames=list(1:(n-1),c("Średnia", "Wariancja", "Błąd średnia", "Błąd wariancja")))
for (i in (2:n))
{
estymatory[i-1,"Åšrednia"] <- mean(losowe[sample(1:n, i, replace=FALSE)])
estymatory[i-1,"Wariancja"] <- var(losowe[sample(1:n, i, replace=FALSE)])
}
estymatory[,"Błąd średnia"] <- (estymatory[,"Średnia"] - m)^2
estymatory[,"BÅ‚Ä…d wariancja"] <- (estymatory[,"Wariancja"] - sigma^2)^2
estymatory
plot_estimate <- function(data, estimate_name, error_name, plot_name, true_value)
{
y_limits <- c(min(0,min(data[,estimate_name])),max(max(data[,estimate_name]),max(data[,error_name])))
plot(data[,estimate_name], t='l',col='blue', ylab="wartość", xlab="n", main=plot_name, ylim=(y_limits))
lines(data[,error_name], col='red')
abline(h=true_value, col='green')
e_name <- tolower(estimate_name)
legend('topright', c(paste("Estymator",e_name), error_name, paste("Prawdziwa wartość",e_name)), lty=1, col=c('blue', 'red', 'green'))
}
Utwórz wykresy estymatora średniej oraz wariancji i ich błędu dla rosnącej liczby zmiennych. Skorzystaj z funkcji plot_estimate. Wykres dla estymatora średniej:
plot_estimate(estymatory, "Średnia", "Błąd średnia", "Estymator średniej", m)
Wykres dla estymatora wariancji:
plot_estimate(estymatory, "Wariancja", "BÅ‚Ä…d wariancja", "Estymator wariancji", sigma^2)
