library(VGAM) library(VGAMdata) library(Rfast) #Pripomenme si nase data data<-belcap[,1] plot(table(data),main='Graf cetnosti',ylab='cetnost') # 1.) #a) maximalizujme verohodnostni funkci verohodnost<-function(x){ verohodnost<-prod(dgeom(data,x)) } #maximalizace nefunguje optimize(verohodnost,c(0,1),maximum=TRUE) #vykreslime graf maximalizovane verohodnosti na mrizce xx<-seq(0.01,0.99, by=0.01) plot(Vectorize(verohodnost)(xx)~xx,xlab="p",ylab=expression(L(p)), type="l") #b) maximalizujme logaritmickou verohodnost logverohodnost<-function(x){ logverohodnost<-sum(dgeom(data,x,log=TRUE)) } p_hat<-optimize(logverohodnost,c(0,1),maximum=TRUE)$max p_hat #vykreslime graf maximalizovane logaritmicke verohodnosti na mrizce plot(Vectorize(logverohodnost)(xx)~xx,xlab="p",ylab=expression(l(p)), type="l") 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 verohodnostni funkci verohodnost<-function(x){ verohodnost<-prod(dzipois(data,x[2],x[1])) } #opet nefunguje optim(c(0.2,4),verohodnost,control=list(fnscale=-1)) #b) maximalizujme logaritmickou verohodnost logverohodnost<-function(x){ logverohodnost<-sum(dzipois(data,x[2],x[1],log=TRUE)) } opt<-optim(c(0.2,4),logverohodnost, method= "BFGS",control=list(fnscale=-1)) opt pihat<-opt$par[1] lambdahat<-opt$par[2] #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')