# Foundations Lab, Day 2: Using R for Maximum Likelihood

###  Writing functions  ###

plus7 <- function(x) return(x + 7)

plus7alt <- function(y) {
  y + 7
}

plus7(5)
plus7(x=4)


histdensity <- function(x) {
  hist(x, prob=T)
  lines(density(x))
} 
# Question: What is returned?

my.data <- rgamma(1000, 1)      # Some fake data
histdensity(my.data)            # Seems to work

hist(my.data, breaks=50)        # Wouldn't more bars be nice?
histdensity(my.data, breaks=50) # Can't just add this to my function

histdensity <- function(x, ...) {
  hist(x, ..., prob=T)
  lines(density(x))
}
histdensity(my.data, breaks=50)                               # Sweet
histdensity(my.data, breaks=50, main="My Graph", ylim=c(0,1)) # Even sweeter


###  Looping with "for"  ###

for(i in 1:3) {
  fake.data <- rnorm(1000, mean=i, sd=i)
  histdensity(fake.data, breaks=50, main=paste("Mean and SD =", i), xlim=c(-10,15))
}
# Where did they go?

par(mfrow=c(1,3))       # Creates a 1x3 array of plots
for(i in 1:3) {
  fake.data <- rnorm(1000, mean=i, sd=i)
  histdensity(fake.data, breaks=50, main=paste("Mean and SD =", i), xlim=c(-10,20) )
}

###  Conditional evaluation with 'if'  ###

match.penny <- function(guess) {
  # Given guess in ("H", "T") returns T for win, F for loss
  
  my.coin <- sample(c("H", "T"), 1)
  
  cat("I guessed ", my.coin, ".\n",
    "You guessed ", guess, ".\n", sep="")
  
  if( guess == my.coin ) {
    cat("They matched, so you win.\n")
    return(T)
  } else {
    cat("They did not match, so you lose.\n")
    return(F)
  }
}

match.penny("H")


###  Finding a max with 'optim'  ###

f <- function(x) return( x * (4-x) )
plot(f, xlim=c(0,4))

max.uncon <- optim(0, f,
  control=list(fnscale=-1),
  hessian=T
)
max.uncon

abline(v=1, lty=2)

max.con <- optim(0, f,
  control=list(fnscale=-1),
  method="L-BFGS-B", lower=-Inf, upper=1,
  hessian=T
)
max.con


g <- function(v) return( -sqrt(sum((v-c(2,3))^2)) )

max.g <- optim(
  rep(0,2),                     # starting value for parameter vector
  g,                            # function to optimize
  control=list(fnscale=-1),     # maximize, don't minimize
  hessian=T                     # get Information matrix
)
max.g

###  Estimating an MLE: bivariate normal  ###

like.one <- function(lambda, x) {
  return( lambda * exp(-lambda * x) )
}

like <- function(lambda, X) {
  # lambda = 1 / mean
  # X = vector of exponential observations
  n <- length(X)
  out <- 1
  for(i in 1:n) {
    out <- out * like.one(lambda, X[i])
  }
  return( out ) 
}

bogus <- rnorm(100, mean=10, sd=2)

like(.1,bogus)

mle.exp <- optim(par=2, 
  fn=like,
  control=list(fnscale=-1),     # maximize, don't minimize
  hessian=T,                    # get Information matrix
  X=bogus
)
mle.exp

# Hmm, try again with log-likelihood
llike.one <- function(lambda, x) {
  return( log(lambda) - lambda * x )
}

llike <- function(lambda, X) {
  # lambda = 1 / mean
  # X = vector of exponential observations
  n <- length(X)
  out <- 0
  for(i in 1:n) {
    out <- out + llike.one(lambda, X[i])
  }
  return( out ) 
}

llike(.1,bogus)

mle.exp.log <- optim(par=2, 
  fn=llike,
  control=list(fnscale=-1),     # maximize, don't minimize
  hessian=T,                    # get Information matrix
  X=bogus
)
mle.exp.log

SE <- sqrt(diag(solve(-1 * mle.exp.log$hessian)))