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

belcap

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


#a)Poissonovo rozdeleni

lambdatilde=mean(data) #odhad parametru lambda metodou momentu
lambdatilde

#b)ZIP rozdeleni

M1=mean(data)
M2=mean(data^2)

M1
M2

#Definujeme soustavu, kterou budeme resit. Nutno prevest do tvaru, kdy na prave strane je 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))
}

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='Poissonuv model',ylab='pravdepodobnost')
points(dpois(x,lambdatilde)~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')