setwd('')

# 1

data <- read.table('newborns.txt', header=TRUE)

x<-na.omit(data$weight.C)

# 1a)
mu<-mean(x) # stredni hodnotu odhadneme pomoci prumeru
sigma_squared<-var(x) # rozptyl odhadneme pomoci vyberoveho rozptylu
sigma<-sqrt(sigma_squared) # a tedy smerodatnou odchylku  pomoci vyberove smerodatne odchylky

# 1b)
hist(x, prob = T, density = 30,
col = 'blue', xlab = 'porodni hmotnost',ylab = 'hustota', main = '')

lines(density(x), col = 'blue', lwd = 2)
curve(dnorm(x,mu,sigma),add=T,col='red',lwd=2)
legend('topleft',c('normalni rozdeleni','jadrovy odhad'),col=c('red','blue'),lwd=2)

plot(ecdf(x),col = 'blue', xlab = 'porodni hmotnost',ylab = 'distribucni funkce', main = '',verticals = TRUE,do.points = FALSE,lwd=2)
curve(pnorm(x,mu,sigma),add=T,col='red',lwd=2)
legend('topleft',c('normalni rozdeleni','empiricka distribucni funkce'),col=c('red','blue'),lwd=2)

# 1c)
pnorm(3800,mu,sigma)

# 1d)
1-pnorm(4000,mu,sigma)

# 1e)
pnorm(4200,mu,sigma)-pnorm(2500,mu,sigma)

# 1f)
0

# 1g)
qnorm(0.05,mu,sigma)

# 1h)
qnorm(0.85,mu,sigma)

###########################################################
# 2
###########################################################
n<-5
# 2a)
pnorm(3000,mu,sigma/sqrt(n))

# 2b)
1-pnorm(2950,mu,sigma/sqrt(n))

# 2c)
pnorm(3900,mu,sigma/sqrt(n))-pnorm(2800,mu,sigma/sqrt(n))

# 2d)
0

###########################################################
# 3
###########################################################
n<-10
# 3a)
pnorm(30000,mu*n,sigma*sqrt(n))

# 3b)
1-pnorm(33000,mu*n,sigma*sqrt(n))

# 3c)
pnorm(33500,mu*n,sigma*sqrt(n))-pnorm(27000,mu*n,sigma*sqrt(n))

# 3d)
0

###########################################################
# 4
###########################################################

data <- read.table('head.txt', header=TRUE)
x<-data$head.L
y<-data$head.W

plot(x,y,xlab='delka hlavy',ylab='sirka hlavy',asp=1)

library(mvtnorm)

mu1<-mean(x)
mu2<-mean(y)

sigma1_squared<-var(x)
sigma2_squared<-var(y)
rho<-cor(x,y)

sigma12<-rho*sqrt(sigma1_squared*sigma2_squared)
cov(x,y)

mu<-c(mu1,mu2) # odhad vektoru strednich hodnot
sigma<-matrix(c(sigma1_squared,sigma12,sigma12,sigma2_squared),nrow=2) # odhad kovariancni matice

# vytvorime si pomocnou mrizku pro vypocet hustoty
xx<-seq(165,215,by=0.5)
yy<-seq(130,170,by=0.5)

# a v kazdem jejim bode spocitame hodnotu odhadnute hustoty
f<-matrix(0, nrow=length(xx), ncol=length(yy))
for(i in 1:length(xx)){
for(j in 1:length(yy)){
f[i,j] <- dmvnorm(c(xx[i],yy[j]),mu,sigma)
}}

contour(xx, yy, f, asp = T, drawlabels = F, nlevels = 11,
xlab='delka hlavy',ylab='sirka hlavy',col='red')
points(x,y)

# na zaver jeden 3D graf
library(GA)

persp3D(xx, yy, f, xlab = 'delka hlavy',ylab = 'sirka hlavy', zlab = 'f(x, y)')

# muzeme jej trochu vylepsit
persp3D(xx, yy, f, phi = 30, theta = 5, ticktype = 'simple',xlab = 'delka hlavy',
ylab = 'sirka hlavy', zlab = 'f(x, y)', col.palette = terrain.colors, border = 'black')


###########################################################
# 5
###########################################################
hist(x, prob = T, density = 30,
col = 'blue', xlab = 'delka hlavy',ylab = 'hustota', main = '')

lines(density(x), col = 'blue', lwd = 2)
curve(dnorm(x,mu1,sqrt(sigma1_squared)),add=T,col='red',lwd=2)
legend('topright',c('normalni rozdeleni','jadrovy odhad'),col=c('red','blue'),lwd=2)

hist(y, prob = T, density = 30,
col = 'blue', xlab = 'sirka hlavy',ylab = 'hustota', main = '')

lines(density(y), col = 'blue', lwd = 2)
curve(dnorm(x,mu2,sqrt(sigma2_squared)),add=T,col='red',lwd=2)
legend('topright',c('normalni rozdeleni','jadrovy odhad'),col=c('red','blue'),lwd=2)

# 5cd)
pnorm(0,mu1,sqrt(sigma1_squared))
pnorm(0,mu2,sqrt(sigma2_squared))

###########################################################
# 6
###########################################################

curve(dchisq(x,df=1),from=0,to=30,col=1,lwd=2,ylim=c(0,0.4),ylab='hustota',main='Chi-kvadrat rozdeleni')
curve(dchisq(x,df=5),from=0,to=30,col=2,lwd=2,add=TRUE)
curve(dchisq(x,df=10),from=0,to=30,col=3,lwd=2,add=TRUE)
curve(dchisq(x,df=20),from=0,to=30,col=4,lwd=2,add=TRUE)
legend('topright',c('1','5','10','20'),col=1:4,lwd=2,title='Stupne volnosti')

curve(dt(x,df=1),from=-5,to=5,col=1,lwd=2,ylim=c(0,0.4),ylab='hustota',main='Studentovo rozdeleni')
curve(dt(x,df=5),from=-5,to=5,col=2,lwd=2,add=TRUE)
curve(dt(x,df=10),from=-5,to=5,col=3,lwd=2,add=TRUE)
curve(dt(x,df=20),from=-5,to=5,col=4,lwd=2,add=TRUE)
legend('topright',c('1','5','10','20'),col=1:4,lwd=2,title='Stupne volnosti')


curve(df(x,df1=1,df2=1),from=0,to=5,col=1,lwd=2,ylab='hustota',main='Fisherovo-Snedecorovo rozdeleni')
curve(df(x,df1=10,df2=15),from=0,to=5,col=2,lwd=2,add=TRUE)
curve(df(x,df1=15,df2=10),from=0,to=5,col=3,lwd=2,add=TRUE)
curve(df(x,df1=20,df2=20),from=0,to=5,col=4,lwd=2,add=TRUE)
legend('topright',c('1,1','10,15','15,10','20,20'),col=1:4,lwd=2,title='Stupne volnosti')

# ukazka vypoctu 10% kvantilu chi-kvadrat rozdeleni s 5 stupni volnosti
qchisq(0.1,df=5)

# ukazka vypoctu 40% kvantilu t rozdeleni s 10 stupni volnosti
qt(0.4,df=10)

# ukazka vypoctu 95% kvantilu F rozdeleni s 12 a 17 stupni volnosti
qf(0.95,df1=12,df2=17)