Pretpostavimo da je uzorak uzet iz neke populacije. Bit bootsrapinga je da se iz orginalnog uzorka uzorkuje, sa vaćanjem, novi uzorak iste duljine, pa pomoću novog uzorka računamo neku statistiku. Proceduru uzorkovanja i računanja statistike ponavljamo mnogo puta tako da dobijemo niz rezultata primjene statistike.
VAŽNO, bootstraping polazi od principa:
Bootstraped uzorak se prema orginalnom uzorku odnosi kao
orginalni uzorak prema populaciji.
To znači da je populacija za bootstraped uzorak jednaka orginalnom uzorku.
Na temelju ova važne činjenice podaci o uzorkovanim rezultatima mogu se primjeniti na populaciju.
Odredite onaj prinos žita za koji vrijedi da je 1% vjerojatno da ćemo dobiti manji prinos.
quantile(mhw$grain, p = 0.01, type = 5)
## 1%
## 2.945
Budući mhw sadrži 500 mjerenja to je taj kvantil podatak između
5. i 6. mjerenja
mean(sort(mhw$grain)[5:6])
## [1] 2.945
Zadatak Napišite funkciju koju ćete predati funkciji boot, a
koja će vraćati 1% kvantil uzorka kojeg generira bootstraping u
funkciji boot. Funkcija mora imati dva parametra, jedan je ime uzorka
(koji je data frame) a drugi vektor indeksa uzorka. Taj vektor slučajno
generira funkcija boot (uzorkovanje sa vraćanjem).
boot.q01 = function (data, indices) {
obs = data[indices, ] # bootstrap iz podataka data
return(quantile(obs$grain, p = 0.01, type = 5))
}
Zadatak Procjenite 1% kvantil populacije žita pomoću neparametarskog bootstrapa.
Funkcije boot() ima, najmanje tri agumenta. Prvi je
data.frame podataka, drugi je funkcije FUN koja će
iz uzorkovanog data.framea izračunati neku statistiku a
treći je broj uzorkovanja R.
Rezultat je složen (vidi ?boot) ali najbitniji dio je sadržan
t0 koja sadrži primjenu funkcije FUN na orginalni uzorak it (varijabli t) koji sadrže R rezultata primjene
funkcije FUN na uzorkovane uzorke
(R rezultata primjene statistike). Taj vektor se dobije naredbomrezultat$t
Provedimo bootstraping sa 1024 uzorkovanja, a na svakom uzorku će se računati 1% kvantil
require (boot)
## Loading required package: boot
b.q01 = boot (mhw, boot.q01, R = 1024)
print (b.q01)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = mhw, statistic = boot.q01, R = 1024)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 2.945 -0.008242 0.07224
Grafički prikažimo rezultat b.q01 (funkcija plot
automatski prepoznaje da je argument b.q01 klase boot).
Biti će nacrtana dva grafa, histogram od b.q01$t i
normal qq plot varijable b.q01$t.
plot(b.q01)
Na slici histograma je isprekidanim pravcem prikazan položaj 1% kvantila orginalnog (non-bootstrapped) skupa podataka. I histogram i qq-plot su isprekidani što je posljedica činjenice da se, prorodno, resampling generira iz konačnog (diskretnog) skupa podataka.
Zadatak Izračunajte bootstrap procjenu 1% kvantila prinosa zrna. Izračunajte bootstrap procjenu normalnog i bazičnog intervala povjerenja uz 0,95% pouzdanosti.
Orginalna procjena, u uzorku, 1% kvantila je sadržana u
rezultatu b.q01 u varijabli t0:
b.q01$t0
## 1%
## 2.945
1% kvantili resamplanih uzoraka su sadržani u varijabli
t pa ćemo, za točkovnu procjenu 1% kvantila uzeti aritemtičku
sredinu od t:
mean(b.q01$t)
## [1] 2.937
Funkcija, iz biblioteke boot, koja računa intervale
povjerenja zove se boot.ci() i zahtjeva tri
parametra
boot, u našem slučaju b.q01conf koja je pouzdanost, u našem slučaju 0,95.typeb.q01.ci = boot.ci(b.q01, conf = 0.95, type = c("norm", "basic"))
b.q01.ci
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1024 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = b.q01, conf = 0.95, type = c("norm", "basic"))
##
## Intervals :
## Level Normal Basic
## 95% ( 2.812, 3.095 ) ( 2.840, 3.080 )
## Calculations and Intervals on Original Scale
Pokušajmo razumijeti relaciju između prinosa slame i zrna. Pretpostavimo da ta relacija ima zajednički uzrok tj. proces rasta pšenice utiče i na prinos slame i na prinos zrna na sistematski način. U strukturnoj analizi pokušavamo otkriti tu povezanost i postavljamo hipotezu da je povezanost u obliku \( y=\alpha + \beta x \) gdje su i \( x \) i \( y \) matematičke varijable. No, ako sa \( X \) i \( Y \) označimo prave vrijednosti tada vrijedi \[ \begin{aligned} X &= x +\xi \\ Y &= y +\eta \end{aligned} \]
gdje su \( \xi \) i \( \eta \) slučajne greške, tj. umjesto da smo izmjerili pravu vrijednost \( X \) mi smo izmjerili \( x \) i napravili slučajnu grešku \( \xi \). Slično za \( \eta \). Greške nam nisu poznate kao ni njihove varijance \( \sigma_{\xi}^2 \) i \( \sigma_{\eta}^2 \).
Označimo sa \( \lambda \) omjer \[ \lambda = \frac{\sigma_{\eta}^2}{\sigma_{\xi}^2}. \] Procjenitelj koeficijenta smjera \( \beta \), ako je \( X \) ulaz (prediktor) a \( Y \) izlaz (prediktant), po metodi maksimalne vjerodostojnosti, je \[ {\hat{\beta}}_{X.Y} = \frac{1}{2 s_{XY}}\left\{\left(s_Y^2 - \lambda s_X^2\right) + \sqrt{\left(s_Y^2 - \lambda s_X^2\right)^2 + 4\lambda s_{XY}^2} \right\} \]
Veličina \( \lambda \) nije poznata pa se nadamo ćemo je dobro aproksimirati omjerom varijanci uzorka \[ \lambda \approx \frac{s_{y}^2}{s_{x}^2}. \]
Napišimo funkciju koja će računati procjenitelj od \( \beta \).
struct.beta <- function(y, x, lambda) {
a <- var(y)-lambda*var(x);
c <- var(x,y);
return((a + sqrt(a^2 + 4 * lambda * c^2))/(2*c))
}
Sada su procjene za \( \beta \) i \( \alpha \) dane sa
beta = struct.beta (mhw$straw, mhw$grain, var(mhw$straw)/var(mhw$grain))
alpha = mean (mhw$straw) - beta * mean (mhw$grain)
Zadatak Odredite 95% interval povjerenja za koeficijent smjera \( \beta \) i odsječak \( \alpha \) pomoću neparametarskog bootstrapa.
Prvo moramo napisati funkciju koja će iz uzorkovanog niza računati
statistiku. Zapravo moramo malo modificirati funkciju struct.beta.
Ta funkcija će vraćati listu dvije vrijednosti; funkcija boot() će
bilježiti te liste i sve njih vratiti u varijabli t rezultata.
boot.sr = function (data, indices) {
obs = data[indices, ]
beta = struct.beta (obs$straw, obs$grain, var(obs$straw)/var(obs$grain))
alpha = mean (obs$straw) - beta * mean (obs$grain)
return (c(beta, alpha))
}
Provedimo bootstrap
b.sr = boot(mhw, boot.sr, R = 1024)
print(b.sr)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = mhw, statistic = boot.sr, R = 1024)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 1.960 0.001782 0.05511
## t2* -1.225 -0.006737 0.21615
Vizualno prikažimo rezultat bootstrapa; prvo za koeficijent smjera \( \beta \)
(vektor t rezultata b.sr ima dva stupca, pa funkciji boot.ci treba, parametrom index reći koji stupac treba koristiti za grafički prikaz)
plot(b.sr, index = 1)
a zatim prikažimo rezultate za \( \alpha \)
plot(b.sr, index = 2)
Sa slika se može zaključiti da su dobivene distribucije približno normalne.
Odredimo osnovi (basic) i normal interval povjerenja za \( \beta \) i
\( \alpha \). Budući su prema predhodnim slikama distribucije približno normalne
očekujemo da neće biti velike razlike između basic i normal intervala
povjerenja.
Prvo točkovna procjena za \( \beta \)
mean(b.sr$t[, 1])
## [1] 1.962
zatim interval povjerenja (vektor t rezultata b.sr ima dva stupca,
pa funkciji boot.ci treba, parametrom index reći koji stupac treba
uzeti u obzir pri računu)
b.sr.ci = boot.ci(b.sr, conf = 0.95, type = c("norm", "basic"), index = 1)
b.sr.ci
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1024 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = b.sr, conf = 0.95, type = c("norm", "basic"),
## index = 1)
##
## Intervals :
## Level Normal Basic
## 95% ( 1.850, 2.066 ) ( 1.852, 2.064 )
## Calculations and Intervals on Original Scale
Točkovna procjena za \( \alpha \)
mean(b.sr$t[, 2])
## [1] -1.232
95% interval povjeranje za \( \alpha \)
a.sr.ci = boot.ci(b.sr, conf = 0.95, type = c("norm", "basic"), index = 2)
a.sr.ci
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1024 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = b.sr, conf = 0.95, type = c("norm", "basic"),
## index = 2)
##
## Intervals :
## Level Normal Basic
## 95% (-1.642, -0.795 ) (-1.625, -0.790 )
## Calculations and Intervals on Original Scale
Ako želim samo podatke o intervalima povjeranje, umjesto da ispišemo
cijeli sadržao od a.sr.ci možemo izvršiti naredbe
a.sr.ci$normal
## conf
## [1,] 0.95 -1.642 -0.7948
a.sr.ci$normal[2:3]
## [1] -1.6421 -0.7948
a.sr.ci$basic
## conf
## [1,] 0.95 999.4 25.63 -1.625 -0.7899
a.sr.ci$basic[4:5]
## [1] -1.6246 -0.7899
Zadatak Odredite točkovnu procjenu i interval povjerenja koeficijenta korelacije između prinosa zrna i slame.
Početi ćemo tako da definiramo funkciju koja računa koeficijent korelacije.
boot.cc = function(data, indices) {
obs = data[indices, ]
cc = cor(obs$grain, obs$straw)
return(cc)
}
Primjenimo bootstraping:
b.cc = boot(mhw, boot.cc, R = 1024)
b.cc
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = mhw, statistic = boot.cc, R = 1024)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.7298 -0.0007324 0.02645
Grafički prikaz bootstrapinga za koeficjent korelacije
plot(b.cc)
Slike pokazuju da postoje manja odstupanje od normalne razdiobe na rubovima područja definicije kao i asimetrija histograma.
Koeficijent korelacije originalnih podataka (bez bootstrapinga)
cor(mhw$grain, mhw$straw)
## [1] 0.7298
Predhodni podatak bi trebao biti isti kao i vrijednost varijabel t0
b.cc$t0
## [1] 0.7298
Točkovna procjena koeficijenta korelacije na temelj bootstapniga
mean(b.cc$t)
## [1] 0.729
Provedimo test normalnosti rezultata bootstrapinga
shapiro.test(b.cc$t)
##
## Shapiro-Wilk normality test
##
## data: b.cc$t
## W = 0.9924, p-value = 4.283e-05
Dobivana vrijednost p-value pokazuje da hipotezu o normalnosti rezultata bootstrapinga treba odbaciti, što znači da se možemo više pouzdati u interval povjeranja tipa “basic” nego “normal”.
Intervali povjerena, na temelju bootstrapinga
b.cc.ci = boot.ci(b.cc, conf = 0.95, type = c("norm", "basic"))
b.cc.ci
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1024 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = b.cc, conf = 0.95, type = c("norm", "basic"))
##
## Intervals :
## Level Normal Basic
## 95% ( 0.6787, 0.7824 ) ( 0.6830, 0.7889 )
## Calculations and Intervals on Original Scale
Numerički podaci o intervalima povjerenja
b.cc.ci$normal[2:3]
## [1] 0.6787 0.7824
b.cc.ci$basic[4:5]
## [1] 0.6830 0.7889