#SCRIPT K BAKALARSKE PRACI
#kod je naimplementovan v programu RStudio
#Nazev prace: Dvourozmerne negativne binomicke rozdeleni
#Jmeno autora: David Sir

###########################################################################
#JSOU ZDE NAIMPLEMENTOVANY ODHADY PRO NEZNAME PARAMETRY DVOUROZMERNEHO NEGATIVNE BINOMICKEHO ROZDELENI
#Odvozeni odhadu viz kapitola 3.5 bakalarske prace.
###########################################################################
#Predpokladame, ze mame data pochazejici z dvourozmerneho negativne binomickeho rozdeleni.



      #MOMENTOVA METODA#
#Funkce mom.odhad nalezne momentovy odhad pro data z dvourozmerneho negativne binomickeho rozdeleni.
#Vstupni parametr: data, pro nez chceme nalezt odhady (matice o dvou sloupcich s nahodnymi vybery X a Y).
mom.odhad = function(data){
  X <- data[,1]
  Y <- data[,2]
  n <- length(X)    #rozsah vyberu
  
  Xn <- sum(X)/n    #vyberovy prumer X
  Yn <- sum(Y)/n    #vyberovy prumer Y
  Sxy <- (1/(n-1))*sum((X-Xn)*(Y-Yn))    #vyberova kovariance X a Y
  
  MOMaa <- (Xn*Yn)/Sxy             #odhad parametru a
  MOMpp1 <- Xn/(MOMaa + Xn + Yn)   #odhad parametru p1
  MOMpp2 <- Yn/(MOMaa + Xn + Yn)   #odhad parametru p2
  MOMpp0 <- 1 - MOMpp1 - MOMpp2    #odhad parametru p0
  
  odhad <- c(MOMaa,MOMpp1,MOMpp2)
  names(odhad) <- c("a","p1","p2")
  return(odhad)
}



      #METODA MAXIMALNI VEROHODNOSTI#
#Funkce verohodnost nalezne logaritmickou verohodnost l_n(a,p1,p2).
#Funkce je nezbytna pro funkci mle.odhad, ktera hleda maximalne verohodny odhad.
verohodnost = function(param,data){
  a <- param[1]
  p1 <- param[2]
  p2 <- param[3]
  X <- data[,1]
  Y <- data[,2]
  n <- nrow(data)
  
  if((p1<0)|(p2<0)|(p1+p2>1)|(a<0)) verohodnost = -Inf  else
    verohodnost <- a*n*log(1-p1-p2)+sum(X)*log(p1)+sum(Y)*log(p2)+sum(log(gamma(a+X+Y)))-n*log(gamma(a))  
  return(-verohodnost) # vraci minus verohodnost, protoze funkce optim hleda minimum
}



#Funkce mle.odhad nalezne maximalne verohodny odhad pro data z dvourozmerneho negativne binomickeho rozdeleni.
#Vstupni parametr: data, pro nez chceme nalezt odhady (matice o dvou sloupcich s nahodnymi vybery X a Y).
mle.odhad = function(data){
  init <- mom.odhad(data)    #pocatecni hodnotu do optimalizace vezmeme momentove odhady
  odhad <- optim(par=init,fn=verohodnost,data=data)$par   #nalezneme maximum logaritmicke verohodnosti
  names(odhad) <- c("a","p1","p2")
  return(odhad)
}