##########################################################################
#                 BINOMICKE ROZDELENI                                   ##
##########################################################################

p <- 0.2
n <- 10
x <- 0:10
pq <- seq(0.001, 1 - 0.001, length = 100)
px <- dbinom(x, size = n, prob = p) # pravdepodobnosti funkce
Fx <- pbinom(x, size = n, prob = p)#distribucni funkce
Finvp <- qbinom(pq, size = n, prob = p) #kvantilova funkce
ukaz <- data.frame(Hodnoty = x, px = px, Fx = Fx)
print(ukaz)

postscript("fig01.ps", width = 7, height = 10, horizontal = FALSE)
par(mfrow = c(3, 1))
plot(x, px, type = "h", xlab = "k", ylab = "P(X=k)", col = "red")
title(main = paste("Binom(", n, ", ", p, ")", sep = ""))
plot(x, Fx, type = "s", xlab = "x", ylab = "P(X<=x)", col = "blue")
plot(pq, Finvp, type = "s", xlab = "p", ylab = "Kvantilova funkce",
     col = "darkgreen")
dev.off()

# nahodne binomicke rozdeleni

rn <- rbinom(1000, size = n, prob = p)
rn[1:10]
table(rn)
prop.table(table(rn))
round(px, 3)


##################################################################################
#                   normalni rozdeleni                                          ##
##################################################################################

mu <- 100
sigma <- 15
x <- seq(mu - 3 * sigma, mu + 3 * sigma, length = 100)
pq <- seq(0.001, 1 - 0.001, length = 100)
fx <- dnorm(x, mean = mu, sd = sigma) #hustota
Fx <- pnorm(x, mean = mu, sd = sigma) #distribucni fce
Finvp <- qnorm(pq, mean = mu, sd = sigma) #kvantilova fce
ukaz <- data.frame(Hodnoty = x, fx = fx, Fx = Fx)
print(ukaz[c(1:5, 45:55, 96:100), ])

pdf("fig02.pdf", width = 7, height = 10)
par(mfrow = c(3, 1))
plot(x, fx, type = "l", xlab = "x", ylab = "f(x)", col = "red")
title(main = paste("N(", mu, ", ", sigma^2, ")", sep = ""))
plot(x, Fx, type = "l", xlab = "x", ylab = "P(X<=x)", col = "blue")
plot(pq, Finvp, type = "l", xlab = "p", ylab = "Kvantilova funkce",
     col = "darkgreen")
dev.off()


#############################################################################
#                   IRIS                                                  ###
#############################################################################

str(iris)
names(iris)
summary(iris)
var(iris$Sepal.Length) #rozptyl
hist(iris$Sepal.Length)
plot(density(iris$Sepal.Length))

cor(iris[,1:4]) #korelacni matice

with(iris, plot(Sepal.Length, Sepal.Width, col=Species, pch=as.numeric(Species)))

library(scatterplot3d)
scatterplot3d(iris$Petal.Width, iris$Sepal.Length, iris$Sepal.Width)
library(rgl)
plot3d(iris$Petal.Width, iris$Sepal.Length, iris$Sepal.Width) #pohyblivy graf

distMatrix <- as.matrix(dist(iris[,1:4]))
heatmap(distMatrix)


####################################################################################
#                       regrese                                                   ##
####################################################################################

year <- rep(2008:2010, each=4)
quarter <- rep(1:4, 3)
cpi <- c(162.2, 164.6, 166.5, 166.0,   #Consumer Price Index
         166.2, 167.0, 168.6, 169.5,
         171.0, 172.1, 173.3, 174.0)
plot(cpi, xaxt="n", ylab="CPI", xlab="")
# draw x-axis
axis(1, labels=paste(year,quarter,sep="Q"), at=1:12, las=3)

fit <- lm(cpi ~ year + quarter)
fit
plot(fit)

library(scatterplot3d)
s3d <- scatterplot3d(year, quarter, cpi, highlight.3d=T, type="h", lab=c(2,3))
s3d$plane3d(fit)

#GL
data("bodyfat", package="mboost")
myFormula <- DEXfat ~ age + waistcirc + hipcirc + elbowbreadth + kneebreadth
bodyfat.glm <- glm(myFormula, family = gaussian("log"), data = bodyfat)
summary(bodyfat.glm)
pred <- predict(bodyfat.glm, type="response")
plot(bodyfat$DEXfat, pred, xlab="Observed Values", ylab="Predicted Values")
abline(a=0, b=1)


#CLUSTERING

iris2 <- iris
iris2$Species <- NULL
(kmeans.result <- kmeans(iris2, 3))
table(iris$Species, kmeans.result$cluster)

plot(iris2[c("Sepal.Length", "Sepal.Width")], col = kmeans.result$cluster)
# plot cluster centers> 
points(kmeans.result$centers[,c("Sepal.Length", "Sepal.Width")], col = 1:3,
             pch = 8, cex=2)

idx <- sample(1:dim(iris)[1], 40)
irisSample <- iris[idx,]
irisSample$Species <- NULL
hc <- hclust(dist(irisSample), method="ave")
plot(hc, hang = -1, labels=iris$Species[idx])
# cut tree into 3 clusters
rect.hclust(hc, k=3)
groups <- cutree(hc, k=3)

