source("newton-raphson.R") make.functions <- function(successes,n.trials) { # Log of the likelihood (without the constant) log.likelihood <- function(pi) { successes * log(pi) + (n.trials - successes) * log(1-pi) } # First derivative of log likelihood d1.log.likelihood <- function(pi) { successes/pi - (n.trials - successes)/(1-pi) } # Second derivative of log likelihood d2.log.likelihood <- function(pi) { -successes/pi^2 - (n.trials - successes)/(1-pi)^2 } list(log.likelihood,d1.log.likelihood,d2.log.likelihood) } successes <- 21 n.trials <- 50 fns <- make.functions(successes,n.trials) mle.via.newton.raphson <- find.mle(fns,0.05,0.01,TRUE,c(0,1),TRUE) mle.via.newton.raphson mle.via.analyticals <- successes/n.trials mle.via.analyticals