# Ze zwracaniem
n <- 10  
sample(1:10, n, replace = TRUE) 
##  [1] 8 2 5 2 3 7 6 2 6 4
# Bez zwracania
n <- 10  
sample(1:10, n, replace = FALSE) 
##  [1] 10  9  8  4  3  7  2  5  6  1

Zadanie 1 Estymacja błędu standardowego średniej próbkowej metodą bootstrap.

Wygeneruj przykładową populację liczności 25 o takim rozkładzie (funkcja rnorm()).

n <- 25  # Liczność populacji
mean <- 3  # Średnia
sd <- 1  # Odchylenie standardowe

populacja <- rnorm(n, mean, sd)  # Generowanie populacji o rozkładzie normalnym

print(populacja)
##  [1] 4.289807 4.172686 3.529093 3.400101 3.451881 2.726090 4.682704 4.190293
##  [9] 3.537358 2.944086 1.675225 1.782996 4.119147 4.349633 3.239754 2.248108
## [17] 3.105042 2.735123 1.842323 1.336595 2.107751 2.271808 2.821928 4.279090
## [25] 3.243105

Wygeneruj 1000 prób o liczności 25 z otrzymanych danych (funkcja replicate()).

próby <- replicate(1000, rnorm(25, 3, 1))

Oblicz średnią każdej z 1000 wylosowanych prób bootstrapowych oraz przedstaw histogram rozkładu tych średnich. Porównaj ten rozkład z rozkładem cechy (rozkładem normalnym N (3;1)). Oceń, czy przybliżenie rozkładu metodą bootstrap jest adekwatne (czy kształty obu rozkładów są podobne)

# Obliczenie średniej dla każdej kolumny/próby
srednie<-c()
for (i in 1:1000){
  srednie<-append(srednie,mean(próby[i]))
}

# Wygenerowanie histogramu rozkładu średnich
hist(srednie, col = "lightsteelblue", main = "Histogram rozkładu średnich", xlab = "Średnia")

n <- 1000  # Liczba obserwacji
mean <- 3  # Średnia
sd <- 1  # Odchylenie standardowe

cecha <- rnorm(n, mean, sd)  # Generowanie rozkładu cechy o rozkładzie normalnym

srednie_cechy<-c()
for (i in 1:1000){
  srednie_cechy<-append(srednie,mean(cecha[i]))
}  

# Wygenerowanie histogramu rozkładu średnich
hist(srednie_cechy, col = "thistle", main = "Histogram rozkładu średnich", xlab = "Średnia")

Porównanie wykresów:

par(mfrow=c(1,2))

hist(srednie, col = "lightsteelblue", main = "Histogram rozkładu średnich", xlab = "Średnia")

hist(srednie_cechy, col = "thistle", main = "Rozkład normalny średnich", xlab = "Średnia")

Porównaj średnią wygenerowanej próby wyjściowej ze średnią wektora średnich z prób bootstrapowych. Oceń ich różnicę.

library(kableExtra)

# Obliczenie średniej próby wyjściowej
średnia_próby_wyjściowej <- mean(próby)

# Obliczenie średniej wektora średnich z prób bootstrapowych
średnia_bootstrap <- mean(srednie)

różnica <- średnia_próby_wyjściowej-średnia_bootstrap

# Porównanie średnich
Wynik<-c(średnia_próby_wyjściowej, średnia_bootstrap, różnica)
Działanie<-c("Średnia próby wyjściowej", "Średnia wektora średnich z prób bootstrapowych", "Różnica")
ramka <- data.frame(Wynik, Działanie)
ramka %>%
  kbl(caption="Porównanie średniej wygenerowanej próby wyjściowej ze średnią wektora średnich z prób bootstrapowych") %>%
  kable_material(c("striped", "hover"))
Porównanie średniej wygenerowanej próby wyjściowej ze średnią wektora średnich z prób bootstrapowych
Wynik Działanie
2.986462 Średnia próby wyjściowej
2.930564 Średnia wektora średnich z prób bootstrapowych
0.055898 Różnica

Różnica między średnią próby wyjściowej a średnią wektora średnich z prób bootstrapowych wynosi 0.000826987524339829. Różnica jest bardzo mała, zatem można wnioskować, że średnia próby wyjściowej jest zbliżona do średniej wektora średnich z prób bootstrapowych. Oznacza to, że wygenerowane próby bootstrapowe dobrze odzwierciedlają próbę wyjściową i można na ich podstawie wnioskować o rozkładzie populacji.

Przy użyciu funkcji boot (pakiet boot) wyznacz błąd standardowy dla średniej.

library(boot)
## Warning: pakiet 'boot' został zbudowany w wersji R 4.2.3
funkcja <- function(próby, i){mean(próby[i])}

# Obliczenie błędu standardowego z użyciem funkcji boot()
boot_results <- boot(próby, funkcja, R = 1000)

# Wyciągnięcie błędu standardowego z wyników
se <- sd(boot_results$t)

# Wyświetlenie wyniku
cat("Błąd standardowy dla średniej:", se, "\n")
## Błąd standardowy dla średniej: 0.2305955

Zadanie 2 Estymacja przedziałowa średniej metodą bootstrap.

Wygeneruj próbkę losową 50 obserwacji wybranej zmiennej z pliku dane.sav (funkcja sample()). Jako populację potraktuj cały zbiór N-obserwacji w pliku.

# Instalacja pakietu haven (jeśli jeszcze nie jest zainstalowany)
#install.packages("haven")

# Ładowanie pakietu haven
library(haven)
## Warning: pakiet 'haven' został zbudowany w wersji R 4.2.3
# Wczytanie pliku .sav
dane <- read_sav("dane (1).sav")
# Wygenerowanie próbki losowej z całej populacji
próbka_sav <- sample(dane$Powierzchnia_mieszkalna, size = 50)
print(próbka_sav)
##  [1] 105  48  76  31 167  46  45  64  74  53  60  54  90  70  43  63 240 115  40
## [20]  53  48  90  30  70  70 130  35  53 189  76  90 110  63  61  54 120  40  80
## [39] 110  60  80  43 160  52  41  99 100  48  41 170

Oszacuj przedziałowo (klasycznie – stosując aproksymację rozkładem normalnym) oraz nieklasycznie (stosując funkcję boot.ci()) średnią wybranej zmiennej w próbce. Oblicz oraz porównaj oba błędy standardowe (względny %) estymacji przedziałowej klasycznej i nieklasycznej. Podaj wnioski.

Szacowanie klasyczne – stosując aproksymację rozkładem normalnym:

# Obliczenie średniej próbkowej i odchylenia standardowego próbkowego
x_bar <- mean(próbka_sav)
s <- sd(próbka_sav)

# Określenie poziomu istotności (np. 95%)
alpha <- 0.05

# Obliczenie granic przedziału
lower <- x_bar - qnorm(1 - alpha/2) * (s / sqrt(length(próbka_sav)))
upper <- x_bar + qnorm(1 - alpha/2) * (s / sqrt(length(próbka_sav)))

# Wyświetlenie przedziału ufności
cat("Przedział ufności dla średniej w próbce:", lower, "-", upper, "\n")
## Przedział ufności dla średniej w próbce: 66.70851 - 91.29149

Szacowanie nieklasyczne:

library(boot)

# Obliczenie przedziału ufności z użyciem funkcji boot() i boot.ci()
boot_results <- boot(data = próbka_sav, statistic = function(próbka_sav, i) {return(mean(próbka_sav[i]))}, R = 50)
boot_ci <- boot.ci(boot_results, type = "basic", conf = 0.95)

# Wyświetlenie przedziału ufności
boot_ci
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 50 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_results, conf = 0.95, type = "basic")
## 
## Intervals : 
## Level      Basic         
## 95%   (65.71, 89.65 )  
## Calculations and Intervals on Original Scale
## Some basic intervals may be unstable
# Obliczenie błędu standardowego dla estymacji przedziałowej klasycznej
klasyczne <- c(61.19856, 85.76144)
se_classic <- (klasyczne[2] - klasyczne[1]) / 2

# Obliczenie błędu standardowego dla estymacji przedziałowej nieklasycznej
nieklasyczne <- c(60.89, 85.97)
se_nonclass <- (nieklasyczne[2] - nieklasyczne[1]) / 2

# Obliczenie błędów standardowych jako względne %
se_classic_percent <- (se_classic / x_bar) * 100
se_nonclass_percent <- (se_nonclass / x_bar) * 100

# Porównanie błędów standardowych
Wynik<-c(se_classic, se_nonclass, se_classic_percent, se_nonclass_percent)
Działanie<-c("Błąd standardowy dla estymacji przedziałowej klasycznej", "Błąd standardowy dla estymacji przedziałowej nieklasycznej", "Błąd standardowy dla estymacji przedziałowej klasycznej jako względny %", "Błąd standardowy dla estymacji przedziałowej nieklasycznej jako względny%")
ramka <- data.frame(Wynik, Działanie)
ramka %>%
  kbl(caption="Porównanie błędów standardowych") %>%
  kable_material(c("striped", "hover"))
Porównanie błędów standardowych
Wynik Działanie
12.28144 Błąd standardowy dla estymacji przedziałowej klasycznej
12.54000 Błąd standardowy dla estymacji przedziałowej nieklasycznej
15.54613 Błąd standardowy dla estymacji przedziałowej klasycznej jako względny %
15.87342 Błąd standardowy dla estymacji przedziałowej nieklasycznej jako względny%

błąd standardowy dla estymacji przedziałowej nieklasycznej jest nieco większy niż dla estymacji przedziałowej klasycznej. Oznacza to, że estymacja przedziałowa klasyczna daje nieco bardziej precyzyjne oszacowanie średniej w próbie. Jednak różnica między tymi błędami standardowymi jest niewielka, więc obie metody dają podobne wyniki.