library(VGAM)
library(VGAMdata)
library(Rfast)

#Pripomenme si nase data
data=belcap[,1]
plot(table(data),main='Graf cetnosti',ylab='cetnost')

# 1.)

#a) maximalizujme funkci verohodnosti
verohodnost=function(x){
a=1
for (i in 1:length(data)){
a=a*dgeom(data[i],x)
}
return(a)
}

#maximalizace nefunguje
optimize(verohodnost,c(0,1),maximum=TRUE)
curve(verohodnost,0,1,xlab="p",ylab=expression(l(p)))

#b) maximalizujme logaritmickou verohodnost
logverohodnost=function(x){
a=0
for (i in 1:length(data)){
a=a+(dgeom(data[i],x,log=TRUE))
}
return(a)
}

p_hat=optimize(logverohodnost,c(0,1),maximum=TRUE)$max
p_hat

curve(logverohodnost,0,1,xlab="p",ylab=expression(l(p)))
points(p_hat,logverohodnost(p_hat),col=2,pch=19,cex=2)

#c)
1/(1+mean(data)) #MLE a odhad metodou momentu geometrickeho rozdeleni

#d)
x=(0:8)
par(mfrow=c(1,2))
plot(table(data)/length(data), type='p',ylab='pravdepodobnost',ylim=c(0,0.25),main='Metoda momentu')
points(dgeom(x,p_hat)~x,col=2,type='h')

plot(table(data)/length(data), type='p',ylab='pravdepodobnost',ylim=c(0,0.25),main='Metoda maximalni verohodnosti')
points(dgeom(x,p_hat)~x,col=3,type='h')

#############################################
# 2.)

#a) maximalizujme funkci verohodnosti

verohodnost=function(x){
a=1
for (i in 1:length(data)){
a=a*dzipois(data[i],x[2],x[1])
}
return(a)
}

#opet nefunguje
optim(c(0.2,4),verohodnost,control=list(fnscale=-1))

#b) maximalizujme logaritmickou verohodnost

logverohodnost=function(x){
a=0
for (i in 1:length(data)){
a=a+(dzipois(data[i],x[2],x[1],log=TRUE))
}
return(a)
}

opt=optim(c(0.2,4),logverohodnost,control=list(fnscale=-1))
pihat=opt$par[1]
lambdahat=opt$par[2]

pihat
lambdahat

#nakreslime graf logaritmicke verohodnosti
x=seq(0.1,0.3,by=0.01)
y=seq(3,6,by=0.1)
z=matrix(rep(0,length(x)*length(y)),nrow=length(x))

for (i in 1:length(x)){
for (j in 1:length(y)){
z[i,j]=logverohodnost(c(x[i],y[j]))
}}


persp(x,y,z,xlim = c(0.1, 0.3), ylim=c(3, 6),xlab=expression(pi),ylab=expression(lambda))
library(rgl) 
persp3d(x,y,z,xlim = c(0.1, 0.3), ylim=c(3, 6),xlab=expression(pi),
ylab=expression(lambda),main='Log-likelihood function')

contour(x,y,z,nlevels=50,xlim = c(0.1, 0.3), ylim=c(3, 6),xlab=expression(pi),
ylab=expression(lambda),main='Contours of log-likelihood function')


#c)
zip.mle(data) #funkce z knihovy Rfast

M1=mean(data)
M2=mean(data^2)
lambdatilde=M2/M1-1      
pitilde=1-M1^2/(M2-M1)

print(c("Odhad metodou momentu:",pitilde,lambdatilde))
print(c("Odhad metodou max. verohodnosti:",pihat,lambdahat))


#d)
x=(0:8)
par(mfrow=c(1,2))
plot(table(data)/length(data), type='p',ylab='pravdepodobnost',ylim=c(0,0.25),main='Metoda momentu')
points(dzipois(x,lambdatilde,pitilde)~x,col=2,type='h')

plot(table(data)/length(data), type='p',ylab='pravdepodobnost',ylim=c(0,0.25),main='Metoda maximalni verohodnosti')
points(dzipois(x,lambdahat,pihat)~x,col=3,type='h')