p <- seq(from=0.00001, to=0.99999, length=2048)
n <- 10
N <- 20
  likely <- function(p, x, N){
    x*log(p)+(N-x)*log(1-p)
  }
  poptim <- optimize(likely, c(0.001, 0.999), x=n, N=N, maximum=T)
  pmax   <- poptim$maximum
  
  plot(p, likely(p, n, N), type='l', ylim=c(-100, -5),
       ylab='l(p|x)', xlab=paste('pravdepodobnost p \n  maximum v bode p=', n/N, sep=''))
  abline(v=pmax, lty=2, col='red')




l <- function(p, x, N){x*log(p) + (N-x)*log(1-p)}
S <- function(p, x, N){x/p-(N-x)/(1-p)}
I <- function(p, x, N){N/(p*(1-p))}

NewtonRaphson<-function(p0, treshold, x, N){
  k         <- 0
  kriterium <- 1
  p1        <- NULL
  for(i in 1:100){
    p1[i] <- p0 + S(p0,x,N)/I(p0,x,N)
    kriterium <- abs(l(p1[i],x,N)-l(p0,x,N))
    p0 <- p1[i]
    k  <- k+1
    if(kriterium < treshold){
      break}
    pmax <- p1[i]
    return(pmax)
  }  
}
pmax <- NewtonRaphson(0.001, 0.00005, n, N)

plot(p, l(p, n, N), type='l', ylim=c(-100, -5),
     ylab='l(p|x)', xlab=paste('pravdepodobnost p \n  maximum v bode p=', round(pmax,2), sep=''))
abline(v=pmax, lty=2, col='red')