source("newton-raphson.R")

make.functions <- function(data) {

  n <- length(data)

  sumy  <- function(x) sum(data-x)
  sumy2 <- function(x) sum((data-x)^2)

  # Log of the likelihood (without the constant)
  log.likelihood <- function(x) {
    -n/2.0*log(2*pi) - n/2.0*log(x) - 1.0/2.0*x^(-1)*sumy2(x)
  }

  # First derivative of log likelihood
  d1.log.likelihood <- function(x) {
    -n/2.0*x^(-1) + 1.0/2.0*x^(-2)*sumy2(x) + x^(-1)*sumy(x)
  }

  # Second derivative of log likelihood
  d2.log.likelihood <- function(x) {
    -n*x^(-1) + (-2.0*sumy(x)+n/2.0)*x^(-2) - x^(-3)*sumy2(x)
  }

  list(log.likelihood,d1.log.likelihood,d2.log.likelihood)
}

data <-  c(3.08, 0.68, 2.09, 0.87, -0.02, 0.25, 1.98, 1.47, 1.95, 0.99)
mean(data)
var(data)

fns <- make.functions(data)

mle.via.newton.raphson <- find.mle(fns[[1]],fns[[2]],fns[[3]],0.95,0.01,FALSE)
mle.via.newton.raphson

mle.via.newton.raphson <- find.mle(fns[[1]],fns[[2]],fns[[3]],1.334,0.01,FALSE)
mle.via.newton.raphson

mle.via.newton.raphson <- find.mle(fns[[1]],fns[[2]],fns[[3]],0.911,0.01,FALSE)
mle.via.newton.raphson

mle.via.newton.raphson <- find.mle(fns[[1]],fns[[2]],fns[[3]],0.2,0.01,FALSE)
mle.via.newton.raphson