#######################################################################################
################ Statistická laboratoř - Zimní škola ##################################
################   Neurčitost v datech - bodové a intervalové odhady  #################
############### Adéla Vrtková, Martina Litschmannová ##################################
#######################################################################################
install.packages("moments")
library(moments)
install.packages("signmedian.test")
library(signmedian.test)

## načtení dat a převod do datové matice
data=read.csv2(file="http://am-nas.vsb.cz/lit40/DATA/aku_new.csv")
dataS=stack(data)	
colnames(dataS)=c("kap","vyrobce") 
dataS=na.omit(dataS)

a = dataS$kap[dataS$vyrobce == "A"]
b = dataS$kap[dataS$vyrobce == "B"]
c = dataS$kap[dataS$vyrobce == "C"]
d = dataS$kap[dataS$vyrobce == "D"]

## Zadání - určete intervalové odhady střední hodnoty a směrodatné odchylky pro kapacitu výrobce A
# Intervalovým odhadům musí vždy předcházet explorační analýza!
# Jak vypadají data? Mají odlehlá pozorování? Jsou symetrická? Jaká je špičatost/šikmost?
# Můžeme předpokládat normalitu? -> boxplot, histogram, qq-grafy, výpočet charakteristik

boxplot(a)

hist(a,
     freq = F)
lines(density(a))

qqnorm(a)
qqline(a)

skewness(a)
kurtosis(a)-3

# z explorační analýzy lze předpokládat normalitu -> můžeme použít funkci t.test pro int.odhad střední hodnoty
t.test(a,conf.level=0.95)  # pro střední hodnotu

# výsledek si lze i uložit a později se k němu vrátit
test=t.test(a,conf.level=0.95)  # pro střední hodnotu
test$conf.int
ls(test) #funkce ls vypíše seznam položek daného objektu

# intervalový odhad sm. odchylky
install.packages("EnvStats")
library(EnvStats)
varTest(a) # intervalový odhad rozptylu -> potřebujeme odmocninu -> ručně nebo...

t=varTest(a, conf.level = 0.95)
sqrt(t$conf.int)

##############################################################################
##############################################################################
## Co když normalita splněna nebude? -> podívejme se na hodnoty náhodně generovaných dat
x = runif(1000,1,200) # náhodný výběr z rovnoměrného rozdělení

summary(x)

boxplot(x)

hist(x,
     freq = F)
lines(density(x))

qqnorm(x)
qqline(x)

skewness(x,na.rm = T)
kurtosis(x,na.rm = T)-3

# zde nelze předpokládat normalitu - funkci t.test nelze použít -> ověření symetrie dat -> intervalový odhad mediánu

wilcox.test(x, conf.int = TRUE,conf.level = 0.95)
signmedian.test(x, conf.level=0.95)

##############################################################################
##############################################################################
## Intervalový odhad rozdílu středních kapacit akumulátorů výrobců A a B
par(mfcol = c(3,2)) # rozdělení graf. okna na 3x2, výpis grafů po sloupcích
# výpis po řádcích pomocí parametru mfrow

# Vizuální kontrola normality (nutno kontrolovat každý výběr zvlášť)
boxplot(a,
        horizontal = T,
        ylim = c(min(a),max(a)))
hist(a,
     xlim = c(min(a),max(a)))
qqnorm(a)
qqline(a)

boxplot(b,
        horizontal = T,
        ylim = c(min(b),max(b)))
hist(b,
     xlim = c(min(b),max(b)))
qqnorm(b)
qqline(b)

# grafy pro prezentaci srovnání výrobců A a B
par(mfrow = c(1,1))
boxplot(list(a,b),
        names = c("A","B"))  # pokud není nastaven tento parametr, popisek je 1,2

# ověření shody rozptylů dle výběrových charakteristik
sd(a)
sd(b)

max(var(a),var(b))/min(var(a),var(b))
# ověřujeme, zda je poměr rozptylů větší než 2
max(var(a),var(b))/min(var(a),var(b))>2 # ukázka použití logické proměnné, vhodné pro cykly

# normalita ok, populační rozptyly nepovažujeme za shodné 
# --> odhad rozdílu středních hodnot - Aspinové-Welchův test (t.test),

mean(a)-mean(b)
t.test(a,b,
       var.equal = F) # nastavením parametru var.equal volíme mezi t testem a Aspinové-Welchovým testem

##############################################################################
##############################################################################
## Intervalový odhad pravděpodobnosti, že kapacita bude větší než 1770 mAh -> potřebujeme novou proměnnou ID
dataS$ID=dataS$kap>1770 # definování binární proměnné definující, zda akumulátory splňují danou podmínku
head(dataS)

# výpočet bodového odhadu pravděpodobnosti, že kapacita bude větší než 1770 mAh
# tj. odpovídající relativní četnosti

table(dataS$ID)
n=length(dataS$ID)
x=31
p=x/n
p

binom.test(x,n,alternative="two.sided",conf.level=0.95) # Clopper-Pearson

#nebo

install.packages("binom")
library(binom)
binom.confint(x,n)

## Interpretace - Pravděpodobnost, že kapacita bude větší než 1770 mAh je s 95% spolehlivostí 6,4 % až 13,0 % (Clopperův-Pearsonův odhad).