library(nleqslv)
library(VGAM)
library(VGAMdata)

belcap

data<-belcap[,1]
plot(table(data),main='Graf cetnosti',ylab='cetnost')


#a) Geometricke rozdeleni
#bude stacit prvni moment

M1<-mean(data)

#Definujeme rovnici, jejiz koren budeme hledat.
#Nutno prevest do tvaru, kdy na prave strane je 0.

f<-function(x){
f<-(1-x)/x-M1
}

curve(f,0,1)
abline(h=0,col=2)

uniroot(f,c(0,1))
ptilde<-uniroot(f,c(0,1))$root

#vzorecek
ptilde<-1/(1+M1)
ptilde


#b) ZIP rozdeleni
#nyni budeme potrebovat prvni dva momenty

M1<-mean(data)
M2<-mean(data^2)

#Definujeme soustavu, kterou budeme resit.
#Nutno prevest do tvaru, kdy na prave strane je (0,0).

f<-function(x){
a<-(1-x[1])*x[2]-M1
b<-(1-x[1])*(x[2]+x[2]^2)-M2
return(c(a,b))
}

f(c(0.2,5)) #Ukazka vystupu funkce f

nleqslv(c(0.2,5),f) #Resime soustavu rovnic. Zvolime vhodne pocatecni podminky.

pihat<-nleqslv(c(0.2,5),f)$x[1]      #Odhad parametru pi
lambdahat<-nleqslv(c(0.2,5),f)$x[2]  #Odhad parametru lambda

M2/M1-1         #presny vzorec
1-M1^2/(M2-M1)  #presny vzorec

#c) graf
x<-(0:8)
par(mfrow=c(1,2))

plot(table(data)/length(data), type='p',ylim=c(0,0.25),main='Geometricky model',ylab='pravdepodobnost')
points(dgeom(x,ptilde)~x,col=2,type='h')

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