# 05 Odhady a hypoteza

# nahraj datový soubor Stulong
attach(Stulong)
summary(Stulong)
dim(Stulong)

# jake je rozlozeni vysek v zakladním souboru (populaci 954 muzu)?
hist(vyska)  # rozložení výšek je souměrné kolem hodnoty cca 175 cm
mean(vyska)     # 175.8 cm
abline(v=175.8, col=2)
sd(vyska)  # směrodatná odchylka v datech je sd(X)=6.26 cm

# konfidencni interval, kdyz znam variabilitu zakladnich dat:
# levá mez: prumer - kvantil_N(0,1) * sigma/odmocnina(n)
mean(vyber)-qnorm(0.975,mean=0,sd=1)*sd(vyska)/sqrt(30)
# pravá mez: prumer + kvantil_N(0,1) * sigma/odmocnina(n)
mean(vyber)+qnorm(0.975,mean=0,sd=1)*sd(vyska)/sqrt(30)
# zapamatuj si, že pro souměrná rozdělení (normální a t-rozdělení) platí:
qnorm(0.025,0,1)  # 2.5% kvantil je stejné číslo jako 97.5% kvantil, až na znaménko
qnorm(0.975,0,1)
# podobně pro t-rozdeleni
qt(0.025,df=29)
qt(0.975, df=29)

# konfidencní interval, když neznám variabilitu zakladnich dat
# sigma odhadnu pomocí sd(vyber) a kvantil beru z t-rozdeleni
# levá mez: prumer - kvantil_t(n-1) * sd(vyber)/odmocnina(n)
mean(vyber)-qt(0.975,df=29)*sd(vyber)/sqrt(30)
# pravá mez: prumer + kvantil_t(n-1) * sd(vyber)/odmocnina(n)
mean(vyber)+qt(0.975,df=29)*sd(vyber)/sqrt(30)


# vyber vzorek 30 vysek
set.seed(21)   # mean=178.1 cm
vyber=sample(vyska, 30)
mean(vyber)
hist(vyber, freq=F, xlab="výška [cm]", ylab="odhad hustoty", main="Výška mužů ve výběru (rel. četnosti)",xlim=c(150,200))

# máme reprezentativni vzorek? Odpovida vyberovy prumer vysek prumeru populacnimu?
# zformuluj hypotezu:
# H0:
# H1:

# grafické znázornění:
abline(v=mean(vyber), col="blue")
lines(x=seq(150,200, by=0.5), y=dnorm(seq(150,200, by=0.5),mean=175.8, sd=6.3), col=2)
abline(v=175.8, col="red",lty=2)
text(x=190,y=0.06,labels = "hustota očekávaného rozdělení", col="red")
text(x=190,y=0.055,labels =expression(N(175.8, 6.26)), col="red")
text(190, 0.05, labels = "průměr výběru", col="blue")
# toto byl histogram namerenych vysek
# pro rozhodovani o prumeru potrebujeme hustotni krivku pro prumer, tedy N(175.8, 6.26/odmoc(30))
6.26/sqrt(30)  # 1.14
# novy graf:
x<-seq(170,182,by=0.1)
y<-dnorm(x,mean=175.8, sd=1.14)
plot(x,y,type="l",xlab="výška [cm]", ylab="odhad hustoty", lwd=2)
text(172,0.3, labels=expression(N(175.8, 6.26/sqrt(30))))
abline(v=175.8,lty=2)
abline(v=mean(vyber), col="blue", lty=2)
abline(h=0,lty=3)
# pokud je vyberovy prumer "dostatecne" vzdaleny od mi0=175.8 cm, hypotezu zamitneme
# co je to "dostatecne" vzdaleny? Zalezi, jakou zvolime pravdepodobnost alfa chyby 1. druhu 
# pravdepodobnosti alfa se take rika hladina testu
# v biologii volime typicky alfa = 0.05 ~ 5 %
# a protoze dopredu nevim, na kterou stranu bude vyberovy prumer vychyleny, 
# rozdelim tuto pravdepodobnost 2.5 % vlevo a 2.5 % pravo
# plochu pravdepodobnosti pod hustotni krivkou urcuje kvantil(alfa/2):
qnorm(0.025, mean=175.8, sd=1.14)  # kvantil vlevo
abline(v=173.57, lty=2, col="red")
qnorm(0.975, mean=175.8, sd=1.14)  # kvantil vpravo
abline(v=178.03, lty=2, col="red")
# kde leží vyberovy prumer? Je "dostatecne" vzdaleny? Potom zamitame hypotezu, 
# ze nas vyber ma populacni stredni vysku muzu 175.8 cm (plati pro vyber ze set.seed(21))


# kdyz neznam variabilitu dat (prumeru)?
# -> pouziju vyberovy odhad S a t-rozdeleni
sd(vyber)  # pro set.seed(21) je S = 7.1
# vypocet testove statistiky t:
(178.1-175.8)/7.1*sqrt(30)
# graf zobrazujici situaci v t-rozdeleni
x<-seq(-3,3,by=0.1)
y<-dt(x, df=29)
plot(x,y,type="l",xlab="hodnoty t-statistiky", ylab="hustota t-rozdeleni", main="Hustota t-rozdělení s 29 stupni volnosti", lwd=2)
text(-2,0.35,labels = "t-rozdělení, df = 29")
abline(v=0, lty=2)
abline(v=(178.1-175.8)/7.1*sqrt(30), lty=2,col="blue")
text(-2,0.32, labels = "testová statistika: t = 1.77", col="blue")
# pravy kvantil t-rozdeleni:
qt(0.975,df=29)
abline(v=qt(0.975,df=29), lty=2, col="red")
text(-2,0.29, labels = "97.5% kvantil = 2.04", col="red")
# testova statistika t není větší než kvantil(0.975) => hypotezu H0 nezamitam
# tento vysledek odrazi fakt, ze jsme neznali variabilitu hypotetickych dat a museli jsme ji odhadovat
#
# dopocitame p-hodnotu pro ziskanou testovou statistiku:
pt(q = 1.77, df = 29)    # toto odpovida P(X<= testova stat.)
1-pt(q = 1.77, df = 29)  # p-hodnota je ta "doplnkova pravdepodobnost, tj. P(X> test.stat.)
2*(1-pt(q = 1.77, df = 29))  # a protoze mame oboustranny test, musíme dopocitat i plochu na druhe strane hustotni krivky