setwd(' ')

xx<-read.table("kostky.csv",header=TRUE,sep=";")
x<-xx$sestky # data ulozime do vektoru x

#1a) prozkoumame data
x<-xx$sestky # data ulozime do vektoru x
table(x)
barplot(table(x),xlab = 'pocet sestek', ylab = 'absolutni cetnost',col=3,density=50)
plot(table(x),xlab = 'pocet sestek', ylab = 'absolutni cetnost')

#1bc)
# zkusime binomicky model; odhadneme jeho parametr p odpovidajici pravdepodobnosti padnuti sestky
n<-length(x)
p<-sum(x)/(10*n) # pocet padnutych sestek/pocet vsech hodu
p # odhadnuta pravdepodobnost padnuti sestky
1/6 # porovnejme s teoretickou hodnotou

#1d)
# muzeme porovnat relativni cetnosti a odhadnute pravdepodobnosti binomickym modelem
prop.table(table(x))
round(dbinom(0:10,10,p),4)

# toto neni prilis prehledne, proto se spise porovnavaji ocekavane o pozorovane cetnosti
table(x) # pozorovane cetnosti
round(n*dbinom(0:10,10,p)) # ocekavane cetnosti

# ulozme tyto hodnoty do vektoru a vykresleme prislusne grafy
poz_cetnosti<-as.numeric(table(x))
oce_cetnosti<-round(n*dbinom(0:10,10,p))

plot(0:10, oce_cetnosti, type = 'h', lty = 2, col = 'black',xlab = 'Pocet sestek', ylab = 'Cetnosti')
points(0:10, oce_cetnosti, pch = 21, col = 'black', bg = 'black')
lines(0:4, poz_cetnosti, type = 'h', col = 'red')
points(0:4, poz_cetnosti, pch = 21, col = 'darkred', bg = 'red')

legend('topright', pch = c(21, 21), col = c('darkred', 'black'), pt.bg = c('red', 'black'),
       legend = c('pozorovane', 'ocekavane'), bty = 'n')

#1e)
# pro popis modelu vykreslime jeste graf odhadnute pravdepodobnostni funkce
plot(0:10,dbinom(0:10,10,p),type='h',xlab = 'Pocet sestek', ylab = 'Pravdepodobnostni funkce',main='Binomicke rozdeleni')
points(0:10,dbinom(0:10,10,p), col = 'red', pch = 19)

# a jeste graf odhadnute distribucni funkce
plot(0:10, pbinom(0:10,10,p), type = 's', main = 'Binomicke rozdeleni',
     xlab = "Pocet sestek", ylab = 'Distribucni funkce', col = 'blue')

# komu se nelibi, muzeme jej malinko vylepsit
plot(0:10, pbinom(0:10,10,p), type = 'n', main = 'Binomicke rozdeleni',
     xlab = "Pocet sestek", ylab = 'Distribucni funkce', col = 'blue')
segments (-1:10, c(0,pbinom(0:10,10,p)) , 0:11, c(0,pbinom(0:10,10,p)))
points (0:10, c(0,pbinom(0:9,10,p)), pch = 21, col=1)
points (0:10, pbinom(0:10,10,p), pch = 19, col=2)

#1f)
# pocitejme (odhadujme) pravdepodobnosti
# pripomenme graf pravdepodobnostni funkce
plot(0:10,dbinom(0:10,10,p),type='h',xlab = 'Pocet sestek', ylab = 'Pravdepodobnostni funkce',main='Binomicke rozdeleni')
points(0:10,dbinom(0:10,10,p), col = 'red', pch = 19)

dbinom(1,10,p) # presne 1 sestka

#1g)
sum(dbinom(0:2,10,p))
pbinom(2,10,p)

#1h)
sum(dbinom(6:10,10,p))
1-pbinom(5,10,p)

########################################################################
########################################################################
#2

x<-xx$neuspechy # data ulozime do vektoru x
table(x)
barplot(table(x),xlab = 'pocet neuspechu', ylab = 'absolutni cetnost',col=3,density=50)
plot(table(x),xlab = 'pocet neuspechu', ylab = 'absolutni cetnost')

#2bc)
# zkusime geometricky model; odhadneme jeho parametr p odpovidajici pravdepodobnosti padnuti sestky
n<-length(x)
p<-1/mean(x+1) # prevracena hodnota prumerneho poctu vsech hodu
p # odhadnuta pravdepodobnost padnuti sestky
1/6 # porovnejme s teoretickou hodnotou

#2d)
# muzeme porovnat relativni cetnosti a odhadnute pravdepodobnosti geometrickym modelem
prop.table(table(x))
round(dgeom(0:15,p),4)

# toto neni prilis prehledne, proto se spise porovnavaji ocekavane o pozorovane cetnosti
table(x) # pozorovane cetnosti
round(n*dgeom(0:15,p)) # ocekavane cetnosti

# ulozme tyto hodnoty do vektoru a vykresleme prislusne grafy
poz_cetnosti<-as.numeric(table(x))
oce_cetnosti<-round(n*dgeom(0:15,p))

plot(0:15, oce_cetnosti, type = 'h', lty = 2, col = 'black',xlab = 'Pocet neuspechu', ylab = 'Cetnosti')
points(0:15, oce_cetnosti, pch = 21, col = 'black', bg = 'black')
lines(c(0,1,2,4,5,6,7,8,9,10,11,14), poz_cetnosti, type = 'h', col = 'red')
points(c(0,1,2,4,5,6,7,8,9,10,11,14), poz_cetnosti, pch = 21, col = 'darkred', bg = 'red')

legend('topright', pch = c(21, 21), col = c('darkred', 'black'), pt.bg = c('red', 'black'),
       legend = c('pozorovane', 'ocekavane'), bty = 'n')

#2e)
# pro popis modelu vykreslime jeste graf odhadnute pravdepodobnostni funkce
plot(0:15,dgeom(0:15,p),type='h',xlab = 'Pocet neuspechu', ylab = 'Pravdepodobnostni funkce',main='
Geometricke rozdeleni',ylim=c(0,0.2))
points(0:15,dgeom(0:15,p), col = 'red', pch = 19)

# a jeste graf odhadnute distribucni funkce
plot(0:15, pgeom(0:15,p), type = 's', main = 'Geometricke rozdeleni',
xlab = "Pocet neuspechu", ylab = 'Distribucni funkce',col = 'blue', ylim=c(0,1))

# komu se nelibi, muzeme jej malinko vylepsit
plot(0:15, pgeom(0:15,p), type = 'n', main = 'Geometricke rozdeleni',
     xlab = "Pocet neuspechu", ylab = 'Distribucni funkce', col = 'blue',ylim=c(0,1))
segments (-1:15, c(0,pgeom(0:15,p)) , 0:16, c(0,pgeom(0:15,p)))
points (0:15, c(0,pgeom(0:14,p)), pch = 21, col=1)
points (0:15, pgeom(0:15,p), pch = 19, col=2)

#2f)
# pocitejme (odhadujme) pravdepodobnosti
# pripomenme graf pravdepodobnostni funkce
plot(0:15,dgeom(0:15,p),type='h',xlab = 'Pocet neuspechu', ylab = 'Pravdepodobnostni funkce',main='Geometricke rozdeleni',ylim=c(0,0.2))
points(0:15,dgeom(0:15,p), col = 'red', pch = 19)

dgeom(4,p) # presne 4 neuspechy

#2g)
dgeom(0,p)

#2h)
1-pgeom(6,p)
sum(dgeom(7:1000,p)) # priblizny vypocet

########################################################################
########################################################################
#3

x<-xx$sestky20s # data ulozime do vektoru x
table(x)
barplot(table(x),xlab = 'pocet sestek', ylab = 'absolutni cetnost',col=3,density=50)
plot(table(x),xlab = 'pocet sestek', ylab = 'absolutni cetnost')

#3bc)
# zkusime Poissonuv model; odhadneme jeho parametr lambda odpovidajici ocekavanemu poctu sestek
n<-length(x)
lambda<-mean(x)
lambda # odhadnuty parametr lambda

#3d)
# muzeme porovnat relativni cetnosti a odhadnute pravdepodobnosti Poissonovym modelem
prop.table(table(x))
round(dpois(0:10,lambda),4)

# toto neni prilis prehledne, proto se spise porovnavaji ocekavane o pozorovane cetnosti
table(x) # pozorovane cetnosti
round(n*dpois(0:10,lambda)) # ocekavane cetnosti

# ulozme tyto hodnoty do vektoru a vykresleme prislusne grafy
poz_cetnosti<-as.numeric(table(x))
oce_cetnosti<-round(n*dpois(0:6,lambda))

plot(0:6, oce_cetnosti, type = 'h', lty = 2, col = 'black',xlab = 'Pocet sestek', ylab = 'Cetnosti')
points(0:6, oce_cetnosti, pch = 21, col = 'black', bg = 'black')
lines(0:5, poz_cetnosti, type = 'h', col = 'red')
points(0:5, poz_cetnosti, pch = 21, col = 'darkred', bg = 'red')

legend('topright', pch = c(21, 21), col = c('darkred', 'black'), pt.bg = c('red', 'black'),
       legend = c('pozorovane', 'ocekavane'), bty = 'n')

#3e)
# pro popis modelu vykreslime jeste graf odhadnute pravdepodobnostni funkce
plot(0:10,dpois(0:10,lambda),type='h',xlab = 'Pocet sestek', ylab = 'Pravdepodobnostni funkce',main='
Poissonovo rozdeleni')
points(0:10,dpois(0:10,lambda), col = 'red', pch = 19)

# a jeste graf odhadnute distribucni funkce
plot(0:10, ppois(0:10,lambda), type = 's', main = 'Poissonovo rozdeleni',
     xlab = "Pocet sestek", ylab = 'Distribucni funkce', col = 'blue',ylim=c(0,1))

# komu se nelibi, muzeme jej malinko vylepsit
plot(0:10, ppois(0:10,lambda), type = 'n', main = 'Poissonovo rozdeleni',
     xlab = "Pocet sestek", ylab = 'Distribucni funkce', col = 'blue',ylim=c(0,1))
segments (-1:10, c(0,ppois(0:10,lambda)) , 0:11, c(0,ppois(0:10,lambda)))
points (0:10, c(0,ppois(0:9,lambda)), pch = 21, col=1)
points (0:10, ppois(0:10,lambda), pch = 19, col=2)

#3f)
# pocitejme (odhadujme) pravdepodobnosti
# pripomenme graf pravdepodobnostni funkce
plot(0:10,dpois(0:10,lambda),type='h',xlab = 'Pocet sestek', ylab = 'Pravdepodobnostni funkce',main='Poissonovo rozdeleni')
points(0:10,dpois(0:10,lambda), col = 'red', pch = 19)

dpois(1,lambda) # presne 1 sestka

#3g)
sum(dpois(0:2,lambda))
ppois(2,lambda)

#3h)
1-ppois(5,lambda)
sum(dpois(6:1000,lambda)) #priblizny vypocet


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

data<-read.table("newborns.txt",header=TRUE)
x<-data$prch.N
x<-na.omit(x)

barplot(table(x),xlab = 'pocet starsich sourozencu', ylab = 'absolutni cetnost',col=3,density=50)
plot(table(x),xlab = 'pocet starsich sourozencu', ylab = 'absolutni cetnost')

#4b)

n<-length(x)
lambda<-mean(x)
lambda # odhadnuty parametr lambda Poissonova rozdeleni

p<-1/mean(x+1) # 
p # odhadnuty parametr p geometrickeho rozdeleni

#4c)

table(x) # pozorovane cetnosti
round(n*dpois(0:9,lambda)) # ocekavane cetnosti pro Poissonovo rozdeleni
round(n*dgeom(0:9,p)) # ocekavane cetnosti pro geometricke rozdeleni

# ulozme tyto hodnoty do vektoru a vykresleme prislusne grafy
poz_cetnosti<-as.numeric(table(x))
oce_cetnostiP<-round(n*dpois(0:9,lambda))
oce_cetnostiG<-round(n*dgeom(0:9,p))

plot(0:9, oce_cetnostiP, type = 'h', lty = 2, col = 'black',xlab = 'Pocet starsich sourozencu', ylab = 'Cetnosti',ylim=c(0,700))
points(0:9, oce_cetnostiP, pch = 21, col = 'black', bg = 'black')

lines(0:9, oce_cetnostiG, type = 'h', lty = 2, col = 'blue',xlab = 'Pocet starsich sourozencu', ylab = 'Cetnosti')
points(0:9, oce_cetnostiG, pch = 21, col = 'blue', bg = 'blue')

lines(0:9, poz_cetnosti, type = 'h', col = 'red')
points(0:9, poz_cetnosti, pch = 21, col = 'darkred', bg = 'red')

legend('topright', pch = c(21, 21,21), col = c('darkred', 'black','blue'), pt.bg = c('red', 'black','blue'),
       legend = c('pozorovane', 'ocekavane Pois','ocekavane Geom'), bty = 'n')


#4d)
dpois(0,lambda) 
dgeom(0,p)

#4e)
1-dpois(0,lambda) 
1-dgeom(0,p)

#4f)
ppois(2,lambda) 
pgeom(2,p)