#######################################################################################
################ Statistická laboratoř - 7. workshop ##################################
############### Adéla Vrtková, Martina Litschmannová ##################################
#######################################################################################
library(moments)
library(signmedian.test)

## načtení dat
data=read.csv("http://am-nas.vsb.cz/lit40/DATA/aku2.csv",sep=";", dec=".", header=TRUE)
data=na.omit(data)

## Zadání - určete intervalové odhady střední hodnoty a směrodatné odchylky pro kapacitu po 5 cyklech 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

a5=data$kap5[data$Výrobce=="A"]

summary(a5)

boxplot(a5)
hist(a5)
qqnorm(a5)
qqline(a5)

skewness(a5)
kurtosis(a5)-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(a5,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(a5,conf.level=0.95)  # pro střední hodnotu
test$conf.int

# R nemá pro intervalový odhad sm. odchylky funkci -> klikátko, nebo "ručně"
n=length(a5)
s=sd(a5)
alpha=0.05
sqrt(((n-1)*s^2)/qchisq(1-alpha/2,n-1) )    # dolní mez IO
sqrt(((n-1)*s^2)/qchisq(alpha/2,n-1) )      # horní mez IO


## Co když normalita splněna nebude? -> podívejme se na kapacitu po 5 cyklech od všech výrobců
summary(data$kap5)

boxplot(data$kap5)
hist(data$kap5)
qqnorm(data$kap5)
qqline(data$kap5)

skewness(data$kap5)
kurtosis(data$kap5)-3

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

wilcox.test(data$kap5, conf.int = TRUE)
signmedian.test(data$kap5, conf.level=0.95)


## Intervalový odhad pravděpodobnosti, že změna v kapacitě bude o více než 10% -> byla definována nová proměnná ID
table(data$ID)
n=length(data$ID)
x=259
p=x/n
alfa=0.05

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

#nebo

p-qnorm(1-alfa/2)*sqrt((p*(1-p))/n)   # dolní mez IO - Waldův 
p+qnorm(1-alfa/2)*sqrt((p*(1-p))/n)   # horní mez IO - Waldův

#nebo

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

## interpretace - pravděpodobnost, že pokles v kapacitě bude o více než 10% je s 95% spolehlivostí z vypočteného intervalu

##############################################################################
## podívejme se ještě na pokles v kapacitě akumulátorů
summary(data$pokles)

boxplot(data$pokles)
hist(data$pokles)
qqnorm(data$pokles)
qqline(data$pokles)

skewness(data$pokles)
kurtosis(data$pokles)-3

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

wilcox.test(data$pokles, conf.int = TRUE)
signmedian.test(data$pokles, conf.level=0.95)