####################################################################################
###### POPISNA STATISTIKA ##########################################################
####################################################################################


library(openxlsx)
pisky <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "Pisky", startRow = 1, colNames = TRUE, rowNames = TRUE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
summary(pisky)
pisky = pisky[,-1]
colnames(pisky)
x =  pisky[,"Sr"]

########################################################################################################################################################

### arithmetic mean

mean(x, na.rm=TRUE)
apply(pisky,2,mean)

# trimming
mean(x,trim=0.1)
apply(pisky,2,function(xx) mean(xx,trim=0.1))
library(asbio)
trim.me(x, trim = 0.1) # trimmovany soubor


### rozptyl (variance)
# populace n / vyber n-1 viz Excel
var(x, na.rm=TRUE)
library(asbio)
var(trim.me(x, trim = 0.1)) # var trim. souboru

#### SD
sd(x, na.rm=TRUE)
sqrt(var(x, na.rm=TRUE))
sd(trim.me(x, trim = 0.1)) # sd trim souboru


### geometric mean
prod(x)^(1/length(x)) #compute the geometric mean
library(FSA)
geomean(x, na.rm = FALSE, zneg.rm = FALSE)

# lognormal mean
library(fuel)
lognormalmean(x, estimator="ml", base = exp(1), n = length(x), m = n - 1, d = 1/n)


#### geometric SD
library(FSA)
sg = geosd(x, na.rm = FALSE, zneg.rm = FALSE)
10^(sg)

#### lognormal SD
library(fuel)
lognormalsd(x, estimator="ml", base = exp(1),n = length(x),m = n - 1, d = 1/n)

### variacni rozpeti
range(x)
max(x)
min(x)
max(x)-min(x)

max(x)/min(x) 
# ratio exceeds 2 orders of magnitude, logarithmic plots will be informative
# if the ratio is between 1.5 and 2 orders of magnitude, logarithmically scaled plots will likely be informative;


#### Variacni koeficient
# CV > 100 % = plots on a logarithmic scale
# CV = 70% - 100%, the inspection of logarithmically scaled plots will likely be informative. 

sd(x)/mean(x)
library(DescTools)
CoefVar(x, weights = NULL, unbiased = FALSE, conf.level = 0,95, na.rm = FALSE)


### sikmost
library(moments)
skewness(x)# sikmost
skewness(x)/(sqrt(6/length(x))) # sikmost delena odhadem str. chyby
#hodnoty pod -2 a nad 2 indikují nenormalitu (Havranek)


### spicatost
library(moments)
kurtosis(x)# spicatost
kurtosis(x)-3
(kurtosis(x)-3)/(sqrt(24/length(x))) # sikmost delena odhadem str. chyby
#hodnoty pod -2 a nad 2 indikují nenormalitu (Havranek)


#### kvantily ####################################################################

### kvartily
Q1 = quantile(x,0.25)
Q2 = quantile(x,0.5)
Q3 = quantile(x,0.75)

quantile(x,c(0.25,0.5,0.75))


### median
median(x)
quantile(x, 0.5)


### interkvartilove rozpeti
IQR(x)
Q3-Q1


### absolutni medianova odchylka (MAD)
mad(x)



### SUMMARY ###########

# pro dany sloupec
summary(x)

library(pastecs)
stat.desc(x)

library(fields)
stats(x)

library(s20x)
summaryStats(x)

library(epiR)
epi.descriptives(x, conf.level = 0.95)


# pro celou matici
summary(pisky)

library(psych)
describe(pisky)

library(Hmisc)
describe(pisky)

library(pastecs)
stat.desc(pisky)

library(fields)
stats(pisky)


### konfidencni intervaly
t.test(x,conf.level=0.95)$conf.int

library(DescTools)
MeanCI(x,conf.level=0.95)

library(Rmisc)
CI(x,ci=0.95)


### bootstrap CI
library(DescTools)
BootCI(x, y = NULL, mean, bci.method = "norm", conf.level = 0.95, sides = "two.sided", R = 999)
# bci.method = c("norm", "basic", "stud", "perc", "bca")
# sides = c("two.sided", "left", "right")

library(boot)
Mboot = boot(x,function(x,i) mean(x[i]), R=10000)
mean(Mboot$t[,1])
boot.ci(Mboot,conf = 0.95,type = "norm")


#### Predikcni intervaly
library(predictionInterval)
pi.m(M=mean(x),SD=sd(x),n=length(x),rep.n=length(x))


#### Tolerancni intervaly
library(tolerance)
nptol.int(y,alpha=0.05,P=0.95,side=2,method="WALD",upper=NULL,lower=NULL)
# method = c("WILKS", "WALD", "HM","YM")


####################################################################################
######  GRAFY EDA ###################################################################
####################################################################################

library(MASS)
data(chem)
x = chem
# A numeric vector of 24 determinations of copper in wholemeal flour, in parts per million.

pisky <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "Pisky", startRow = 1, colNames = TRUE, rowNames = TRUE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
summary(pisky)
pisky = pisky[,-1]
colnames(pisky)
x =  pisky[,"Ba"]

########################################################################################################################################################

# Cullen Frey plot
library(fitdistrplus)
descdist(x, discrete=FALSE)
descdist(x, discrete=FALSE, boot=900)
descdist(trim.me(x, trim = 0.1), discrete=FALSE, boot=900) # asbio


### stripchart ############################################################################################################

stripchart(x,method="overplot",pch=1,col=1,cex=1.5)         # vertical=FALSE
stripchart(x,method="overplot",pch=1,col=1,cex=1.5,log="x")

stripchart(x,method="jitter",jitter=0.1, pch=1,col=1,cex=1.5)
stripchart(x,method="jitter",jitter=0.1, pch=1,col=1,cex=1.5,log="x")

stripchart(x,method="stack",pch=1,col=1,cex=1.5)
stripchart(x,method="stack",pch=1,col=1,cex=1.5,log="x")


library(beeswarm)
beeswarm(x, log = FALSE,pch = 1, col = 1, cex = 1.2, vertical=FALSE)
beeswarm(x, log = TRUE,pch = 1, col = 1, cex = 1.2, vertical=FALSE)



### boxplot #####################################################################

boxplot(x, horizontal=TRUE, notch = FALSE, xlab = "x")

boxplot(x, horizontal=TRUE, notch = TRUE, xlab = "x")

library(beeswarm)
boxplot(x, horizontal=TRUE, notch = FALSE, xlab = "x", log = "x", cex = 1.2)
beeswarm(x, log = TRUE,pch = 1, col = 1, cex = 1.2, vertical=FALSE, add=TRUE)

boxplot(x,xlab="x",las=2, horizontal=TRUE, cex=1.2,log="x") 
stripchart(x,pch=1, col=1, vert=FALSE, add=TRUE, cex=1.2,log="x")

boxplot.stats(x)
bs=boxplot.stats(x)
bs$out



### histogram #####################################################################

hist(x, freq=TRUE, breaks = 20, plot = TRUE,main = NULL) # pro counts

hist(x, freq=FALSE, breaks = 10, plot = TRUE,main = NULL)
rug(x)


#### Jadrovy odhad hustoty pravdepodobnosti (Kernel density estimator, KDE) ####################################################

den1 = density(x)
den1$bw
den1$x
den1$y
plot(den1$x, den1$y, type="l", ylab = "density", xlab = "x")
rug(x, lwd=2)

plot(density(x), main="")  #A bit bumpy
rug(x, lwd=2)

density(x,adjust = 10)
plot(density(x,adjust = 10), main="") #Very smooth
rug(x, lwd=2)

density(x,adjust = 0.1)
plot(density(x,adjust = 0.1), main="") #crazy bumpy
rug(x, lwd=2)


#### houslovy graf ##########################################################################################

library(vioplot)

vioplot(x, range=1.5, horizontal=TRUE, col=0)
stripchart(x, pch=1, col=2, vert=FALSE, add=TRUE, cex=1.2)

vioplot(x, range=1.5, horizontal=TRUE, col=0)
stripchart(x, method="stack", pch=1, col=2, vert=FALSE, add=TRUE, cex=1.2)

vioplot(x, range=1.5, horizontal=TRUE, col=0)
beeswarm(x, log = FALSE,pch = 1, col = 2, cex = 1.2, vertical=FALSE, add=TRUE)


#### QQ graf ############################################################################################################################

my.qq <- qqnorm(sort(x), type = "p", xlab = "x", ylab = "Quantiles of the Normal Distribution",main=NULL)
qqline(x, col = "red")

library(car)
qqPlot(x, dist="norm",line="robust",envelope=.95)
# line=c("quartiles", "robust", "none")


#### EDA plots ################################################################################################

library(UsingR)
simple.eda(x)

####################################################################################################################################################################


### Testy na odlehlé body ###########################################################################

library(outliers)
dixon.test(x, type=22, opposite=FALSE, two.sided=TRUE)
# 10 pro 3-7, 11 pro 8-10, 21 pro 11-13, 22 pro 14 a vice hodnot
grubbs.test(x, type=10, opposite=FALSE, two.sided=TRUE)
grubbs.test(x, type=11, opposite=FALSE, two.sided=TRUE)
grubbs.test(x, type=20, opposite=FALSE, two.sided=TRUE)

# 10 pro 1 odl. bod, 11 pro 1 odl. bod na kazdem konci, 20 pro 2 odl. body na 1 konci

library(climtrends)
FindOutliersTietjenMooreTest(x,k=1,alpha=0.05)

library(EnvStats)
rosnerTest(x, k = 2, alpha = 0.05, warn = TRUE)

#### TESTY NORMALITY #####################################################################################

library(s20x)
normcheck(x, xlab = c("Theoretical Quantiles", ""),
          ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue",
          bootstrap = TRUE, B = 5, bpch = 3, bcol = "lightgrey",
          shapiro.wilk = FALSE, whichPlot = 1:2, usePar = TRUE)


# Shapiro Francia test
library(nortest)
ShapiroFranciaTest(x)

# Shapiro Wilk test
shapiro.test(x)

# Anderson Darling test
library(nortest)
ad.test(x)

# DAgostinuv test sikmosti
library(moments)
agostino.test(x, alternative = "two.sided")
# alternative = c("two.sided", "less", "greater")

# Anscombe-Glynnuv test spicatosti
library(moments)
anscombe.test(x, alternative = "two.sided")
# alternative = c("two.sided", "less", "greater")


#### TESTY GOODNESS OF FIT #####################################################################################

# Anderson Darlinguv test
library(DescTools)
AndersonDarlingTest(x, "pnorm", mean=mean(x), sd=sd(x))

# Cramer von Mises
library(goftest)
cvm.test(x, "pnorm", mean=mean(x), sd=sd(x), estimated=TRUE)

# Kolmogorov Smirnov test
library(KScorrect)
LcKS(x, cdf="pnorm", nreps = 4999, G = 1:9, varModel = "E", parallel = FALSE, cores = NULL)
# varModel = c("E", "V")
# "pnorm" for normal,"pmixnorm" for (univariate) normal mixture,"plnorm" for lognormal (log-normal, log normal), "punif" for uniform, "plunif" for loguniform (log-uniform, log uniform),
library(dgof)
ks.test(x, y="pnorm",alternative = "two.sided", exact = NULL, tol=1e-8, simulate.p.value=FALSE, B=2000)
# alternative = c("two.sided", "less", "greater")

### Fitting distribution / Goodness of fit #################################################################

### normal
library(MASS)
dfp <- fitdistr(x, densfun="normal")
(dffl <- dfp$loglik)
dff <- dfp$estimate
mean <- as.numeric(dff[1]); mean
sd <- as.numeric(dff[2]); sd
mean(x)
sd(x)

ncd = seq(min(x), max(x), 0.05)
npd = dnorm(ncd, mean = mean, sd = sd)
hist(x,freq = FALSE)
lines(ncd,npd, type="l",col=2)
rug(x)

library(EnvStats)
norm.est <- enorm(x, method = "mvue", ci = TRUE, ci.type = "two-sided",ci.method = "exact", conf.level = 0.95, ci.param = "mean")
norm.m <- as.vector(norm.est[[3]])[1]
norm.sd <- as.vector(norm.est[[3]])[2]
###
sw.norm <- gofTest(x, test="sw", dist = "norm"); sw.norm
plot(sw.norm)
norm.W <- sw.norm[[6]]
norm.p <- sw.norm[[10]]


### lognormal
library(MASS)
dfp <- fitdistr(x, densfun="lognormal")
(dffl <- dfp$loglik)
dff <- dfp$estimate
meanlog <- as.numeric(dff[1])
sdlog <- as.numeric(dff[2])

library(EnvStats)
lnorm.est <- elnorm(x, method = "mvue", ci = TRUE, ci.type = "two-sided",ci.method = "exact", conf.level = 0.95)
lnorm.mlog <- as.vector(lnorm.est[[3]])[1]
lnorm.sdlog <- as.vector(lnorm.est[[3]])[2]
#
sw.lnorm <- gofTest(x, test="sw", dist = "lnorm",param.list = list(meanlog=lnorm.mlog, sdlog = lnorm.sdlog)); sw.lnorm
plot(sw.lnorm)
lnorm.W <- sw.lnorm[[6]]
lnorm.p <- sw.lnorm[[10]]


#### MULTIMODALITA #####################################################################################

##### modelovani kontaminovanych distribuci

# multimodalita 1
xfit<-seq(-20, 20, length=10000)
yfit1<-dnorm(xfit, mean=-5, sd=3)
yfit2<-dnorm(xfit, mean=5, sd=3)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

# multimodalita 2
xfit<-seq(-10, 20, length=10000)
yfit1<-dnorm(xfit, mean=-2, sd=1)
yfit2<-dnorm(xfit, mean=5, sd=4)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

# zmena sikmosti  1
xfit<-seq(-10, 20, length=10000)
yfit1<-dnorm(xfit, mean=0, sd=2)
yfit2<-dnorm(xfit, mean=4, sd=4)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

# zmena sikmosti  2
xfit<-seq(-35,50, length=10000)
yfit1<-dnorm(xfit, mean=-6, sd=10)
yfit2<-dnorm(xfit, mean=12, sd=15)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

# zmena spicatosti  1
xfit<-seq(-20, 20, length=10000)
yfit1<-dnorm(xfit, mean=-2.5, sd=4)
yfit2<-dnorm(xfit, mean=5, sd=4)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

# zmena spicatosti 2
xfit<-seq(-8, 8, length=10000)
yfit1<-dnorm(xfit, mean=0, sd=0.8)
yfit2<-dnorm(xfit, mean=0, sd=2.4)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

########################################################################################################

kyjov<- read.delim("c:\\Users\\lubop\\Documents\\kurs\\kyjov.txt", header=T)
x = kyjov[which(kyjov[,3]=="1034"),]
x = x[,6]

####################################################################################

hist(x, main="")
rug(x)

library(diptest)
dtb = dip(x, full.result = TRUE, min.is.0 = TRUE, debug = FALSE); dtb
dip.test(x, simulate.p.value = FALSE, B = 2000)
plot(dtb)

library(Rfolding)
folding.test(as.matrix(x))


#### Kernel density estimate ########################################################################################

library(multimode)
locmodes(x,mod0=1,display=TRUE,xlab="x")
locmodes(x,mod0=2,display=TRUE,xlab="x")
locmodes(x,mod0=3,display=TRUE,xlab="x")
locmodes(x,mod0=4,display=TRUE,xlab="x")

modetree(x,bws=NULL,gridsize=NULL,cbw1=NULL,cbw2=NULL,display=TRUE,logbw=FALSE,logbw.regulargrid=NULL,xlab="x")
modetree(x,bws=c(0.1,1),gridsize=NULL,cbw1=NULL,cbw2=NULL,display=TRUE,logbw=FALSE,logbw.regulargrid=NULL,xlab="x")

bw1 = bw.crit(x,mod0=1,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),full.result=F)
bw2 = bw.crit(x,mod0=2,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),full.result=F)
bw3 = bw.crit(x,mod0=3,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),full.result=F)
bw4 = bw.crit(x,mod0=4,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),full.result=F)


library(silvermantest) # https://www.mathematik.uni-marburg.de/~stochastik/R_packages/
modes=1:4
densities.plot(x,modes=modes,mark_modes=TRUE,in_one=TRUE*(length(modes)<4))
#  null hypothesis: underlying density has at most k modes
silverman.test(x, k=1, M = 999, adjust = FALSE, digits = 6)
silverman.test(x, k=2, M = 999, adjust = FALSE, digits = 6)
silverman.test(x, k=3, M = 999, adjust = FALSE, digits = 6)
silverman.plot(x, kmin=1, kmax=5, alpha=0.05, adjust=FALSE)


library(mixtools)
fit <- normalmixEM(x, k = 3) #try to fit k Gaussians
summary(fit)


###################################################################################################################
######### Testy shody ###############################################################################################
##########################################################################################################

b314<- read.table("c:\\Users\\lubop\\Documents\\kurs\\Data1\\ascii03\\B\\B314.txt", header=T)
# dekl obsah 2.4 mg/l
is.list(b314)
x = unlist(b314)

#### t-test ############################################################################################
# one sample t-test: t.test(x,mu=3, alternative = "two.sided") 

t.test(x, y = NULL, alternative = "two.sided", mu = 2.4, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

library(pwr)
p.t.one <- pwr.t.test(d=(abs(mean(x)-2.4)), power=0.9, type="one.sample", alternative="two.sided")
plot(p.t.one)


#### medianovy test ############################################################################################

library(BSDA)
SIGN.test(x, y = NULL, md = 2.4, alternative = "two.sided", conf.level = 0.95)


#### Wilcoxonuv (Mann-Whitney) test ############################################################################################

wilcox.test(x, y = NULL, alternative = "two.sided", mu = 2.4, paired = FALSE, exact = NULL, correct = TRUE, conf.int = TRUE, conf.level = 0.95) 

########################################################################################################################################################################


##### Dvouvyberove testy  ###########################################################################################################################################


###### Scale #######################################################################################

### F test
var.test(x, y, ratio = 1, alternative = "two.sided", conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

### Ansari Bradley two sample test
ansari.test(x, y, alternative = "two.sided", exact = FALSE, conf.int = FALSE, conf.level = 0.95)


###### Location ###################################################################################

b325ab<- read.table("c:\\Users\\lubop\\Documents\\kurs\\Data1\\ascii03\\B\\B325ab.txt", header=T)
# dekl obsah 2.4 mg/l
is.list(b325ab)
X = as.data.frame(b325ab)
# x = unlist(b325ab)
x= unlist(X[,1])
y= unlist(X[,2])


#### t-test ############################################################################################

# t.test(x,mu=3, alternative = "two.sided") 

sd(x)
sd(y)
# srovnani pomoci F-testu nebo jineho testu shody rozptylu

boxplot(x,y)

t.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95)
t.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) # Welchuv test
# alternative = c("two.sided", "less", "greater")

library(pwr)
p.t.two <- pwr.t.test(d=(abs(mean(x)-mean(y))), power=0.9, type="two.sample", alternative="two.sided")
plot(p.t.two)
plot(p.t.two, xlab="sample size per group")

p.t.two2 <- pwr.t2n.test(n1 = length(x), n2= length(y), d = (abs(mean(x)-mean(y))), sig.level = 0.05, power = NULL,
                         alternative ="two.sided")
# alternative = c("two.sided","less","greater")
plot(p.t.two2)
plot(p.t.two2, xlab="sample size per group")


#### 2 sample medianovy test ############################################################################################

# one sample sign test 
library(BSDA)
SIGN.test(x, y, md = 0, alternative = "two.sided", conf.level = 0.95)


#### 2 sample Wilcoxonuv (Mann-Whitney) test ############################################################################################


wilcox.test(x, y,alternative = "two.sided", mu = 0, paired = FALSE, exact = FALSE, correct = TRUE,conf.int = FALSE, conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

library(exactRankTests)
wilcox.exact(x, y, mu = 0,alternative = "two.sided",paired = FALSE, exact = FALSE,conf.int = TRUE, conf.level = 0.95)

######################################################################################################################################################

####################################################################################################################################################

### Srovnani 2 distribuci ##########################################################################################


### Kolmogorov Smirnov test
ks.test(x,y, alternative = "two.sided", exact = FALSE)
# alternative = c("two.sided", "less", "greater")

library(KScorrect)
ks_test_stat(x, y)

library(dgof)
ks.test(x, y,alternative = "two.sided", exact = NULL, tol=1e-8, simulate.p.value=TRUE, B=2000)
# alternative = c("two.sided", "less", "greater")

library(twosamples)
ks_test(x, y, nboots = 2000, p = 0.025)


### Cramer von Mises
library(CDFt)
CramerVonMisesTwoSamples(x, y)

library(cramer)
cramer.test(x,y,conf.level=0.95,replicates=1000,sim="ordinary",just.statistic=FALSE,kernel="phiCramer", maxM=2^14, K=160)

library(twosamples) # nefunguje
cvm_test(x, y, nboots = 2000, p = 0.05)

#######################################################################################################################


