Nieklasyczne metody statystyki

Lab 6

Wygeneruj ciąg losowo wybranych liczb z zakresu od 1 do 10. Zastosuj losowanie ze zwracaniem i bez zwracania.

sample(1:10) #bez zwracania

##  [1]  4 10  9  7  8  5  2  1  6  3

sample(1:10, replace=TRUE) #ze zwracaniem==

##  [1] 9 5 4 5 8 2 4 6 8 3

Zadanie 1. Estymacja błędu standardowego średniej próbkowej metodą bootstrap. Załóżmy, że chcemy badać populację o rozkładzie normalnym ze średnią 3 i odchyleniem standardowym 1.

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

rn <- rnorm(25, mean=3, sd=1) 

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

proby_1000 <- replicate(1000, rnorm(25, mean=3, sd=1))  

Oblicz średnią każdej z 1000 wylosowanych prób bootstrapowych

srednia <- apply(proby_1000, 2, mean) 

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)

hist(srednia) 

hist(rn)

hist(rnorm(25,3,1))

Powyzej przedstawiono 3 histogramy. Histogram srednia nie jest symetryczny. Histogram rn (rozklad normlany) pokazuje jak ksztaltuja sie wyniki rn. Histogram rnorm (25,3,1) prezentuje wyniki dla średniej 3 i odchylenia 1.

Kształt 1 i 2 histogramu jest podobny.

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

mean(rn)

## [1] 2.94824

mean(srednia)

## [1] 2.996395

Średnie różnią się od siebie.

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

sd.boot=function(srednia, idx){
  ans=sd(srednia[idx])
}

boot(srednia, statistic = sd.boot, R=1000)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = srednia, statistic = sd.boot, R = 1000)
## 
## 
## Bootstrap Statistics :
##      original        bias    std. error
## t1* 0.1988978 -0.0001617877 0.004512747

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.

setwd("C://Users//kinit//Downloads")
dane <- read.spss("dane.sav", to.data.frame = TRUE, use.value.labels = TRUE)

pm <- sample(dane$Powierzchnia_mieszkalna, 50)
pm

##  [1]  60  74  59  50  50 190  35  54  63  70  42 100 130  96  54  48  80 168 120
## [20] 220  25 120  65  58  51  84  87  63  60  57  30  48  44  25  45  43 250 100
## [39]  74  32  70  70  38  60  63 180  62 114  39 200

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.

klasycznie <- sd(pm)/sqrt(pm)
klasycznie

##  [1]  6.722234  6.053038  6.778963  7.363839  7.363839  3.777570  8.801471
##  [8]  7.085857  6.560229  6.223580  8.034607  5.207020  4.566859  5.314393
## [15]  7.085857  7.515687  5.821626  4.017303  4.753338  3.510572 10.414041
## [22]  4.753338  6.458514  6.837153  7.291287  5.681325  5.582512  6.560229
## [29]  6.722234  6.896867  9.506675  7.515687  7.849879 10.414041  7.762168
## [36]  7.940632  3.293209  5.207020  6.053038  9.204799  6.223580  6.223580
## [43]  8.446902  6.722234  6.560229  3.881084  6.612922  4.876821  8.337906
## [50]  3.681919

dolny <- mean(mean(pm)-1.96*klasycznie)
gorny <- mean(mean(pm)+1.96*klasycznie)

dolny

## [1] 67.61807

gorny

## [1] 93.18193

Odchylenie standardowe

sd(pm)

## [1] 52.0702

Nieklasycznie

mean.boot=function(pm, idx){
  ans=mean(pm[idx])
}

srednia_boot_pm=boot(pm, statistic = mean.boot, R=50)
class(srednia_boot_pm)

## [1] "boot"

names(srednia_boot_pm)

##  [1] "t0"        "t"         "R"         "data"      "seed"      "statistic"
##  [7] "sim"       "call"      "stype"     "strata"    "weights"

plot(srednia_boot_pm)

boot.ci(srednia_boot_pm, conf = 0.95, type=c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 50 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = srednia_boot_pm, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (66.88, 94.94 )   (67.04, 94.70 )  
## Calculations and Intervals on Original Scale
## Some percentile intervals may be unstable

Średnia oszacowana w sposób nieklasyczny zawiera się w przedziale od 71,60 do 98,02, średnia określona klasycznie mieści się w tymże przedziale.

Nieklasycznie - odchylenie standardowe

sd.boot1=function(x, idx){
  ans=sd(x[idx])
}

odch_nieklas_pm=boot(pm,statistic = sd.boot1, R=50)
class(odch_nieklas_pm)

## [1] "boot"

names(odch_nieklas_pm)

##  [1] "t0"        "t"         "R"         "data"      "seed"      "statistic"
##  [7] "sim"       "call"      "stype"     "strata"    "weights"

plot(odch_nieklas_pm)

Odchylenie standardowe nieklasyczne

boot.ci(odch_nieklas_pm, conf = 0.95, type=c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 50 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = odch_nieklas_pm, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (38.58, 66.16 )   (36.84, 64.23 )  
## Calculations and Intervals on Original Scale
## Some percentile intervals may be unstable

Odchylenie mieści się między 30,44 i 54,67. Odchylenie standardowe klasyczne również mieści się w tym zakresie.

Można zatem stwierdzić, że zarówno średnie i odchylenia standardowe obliczone metodą klasyczną i nieklasyczną są poprawne i pokrywają sie.