# 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)    # generate draws with i as mean and standard deviation
  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) {
  # Play "matching pennies" against a fair coin
  # Given guess in ("H", "T") returns TRUE for win, FALSE 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(TRUE)
  } else {
    cat("They did not match, so you lose.\n")
    return(FALSE)
  }
}

match.penny("H")


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

f <- function(x) return( x * (4-x) )
plot(f, xlim=c(0,4))      # take a look at the graph of the function

max.uncon <- optim(0, f,
  control=list(fnscale=-1),     # scales the function by -1 so we're maximizing instead of minimizing
  hessian=T                     # handy later for finding standard errors
)
max.uncon                       # LOOK AT THE OUTPUT.  I warned you.

# What about a constrained optimum?  Limit f to (-Inf,1].
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

# Optimize a function of two parameters.
# NOTE: both are passed in a *single* vector
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: exponential  ###
########################################

like.exp.one <- function(lambda, x) {
  # Density of exponential distribution at x given mean lambda
  return( 1/lambda * exp(-x/lambda) )
}

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

bogus <- rnorm(100, mean=10, sd=2)      # Yes, I know this is not exponential
histdensity(bogus, breaks=20, main="Normal Draws with Mean=10, SD=2")

like.exp(lambda=5, X=bogus)

mle.exp <- optim(par=5, 
  fn=like.exp,
  control=list(fnscale=-1),     # maximize, don't minimize
  method="L-BFGS-B", lower=0, upper=Inf,
  hessian=T,                    # get Information matrix
  X=bogus
)
mle.exp
# Hmm, that seems odd. Notice that the output $par is the starting value.
# It didn't go anywhere!

# Let's try again with the log-likelihood
llike.exp.one <- function(lambda, x) {
  # Log-density of exponential distribution at x given mean lambda
  return( -log(lambda) - x/lambda )
}

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

llike.exp(lambda=.1, X=bogus)

mle.exp.log <- optim(par=5, 
  fn=llike.exp,
  control=list(fnscale=-1),     # maximize, don't minimize
  method="L-BFGS-B", lower=0.1, upper=Inf,
  hessian=T,                    # get Information matrix
  X=bogus
)
mle.exp.log

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


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

install.packages("mvtnorm")
library(mvtnorm)                # gives us multivariate normal and student-t draws

true.params <- list(mean=c(10, 20), sd=c(2,5), cor=.75)
bogus.2d <- rmvnorm(
  n=200, 
  mean=c(10, 20), 
  sigma=matrix(c(true.params$sd[1]^2, 
    true.params$sd[1]*true.params$sd[2] * true.params$cor,
    true.params$sd[1]*true.params$sd[2] * true.params$cor,
    true.params$sd[2]^2), nrow=2)
)
plot(bogus.2d)

llike.mvn.one <- function(par, x) {
  # Log-density of multivariate normal at x given par, where
  # par[1:2] is the mean vector
  # par[3] = log(standard deviation of first variable)
  # par[4] = log(standard deviation of second variable)
  # par[5] = correlation between first and second variables
  this.mean <- par[1:2]
  this.sigma <- matrix(c(par[3]^2, 
    par[3]*par[4]*par[5],
    par[3]*par[4]*par[5],
    par[4]^2), nrow=2)
  return( dmvnorm(x=x, mean=this.mean, sigma=this.sigma) )
}
llike.mvn.one(par=c(1,2,3,4,.5), c(6,7))

llike.mvn <- function(par, X) {
  # par[1:2] is the mean vector
  # par[3] = standard deviation of first variable
  # par[4] = standard deviation of second variable
  # par[5] = correlation between first and second variables
  # X = nx2 matrix of observations
  n <- nrow(X)
  out <- 0
  for(i in 1:n) {
    out <- out + llike.mvn.one(par, X[i,])
  }
  return( out ) 
}
llike.mvn(par=c(15,15,1,1,0), X=bogus.2d)

mle.mvn <- optim(
  par=c(15,15, 1,1, 0),         # starting values
  fn=llike.mvn,
  control=list(fnscale=-1),     # maximize, don't minimize
  method="L-BFGS-B", 
    lower=c(-Inf, -Inf, 0, 0, -1), 
    upper=c(Inf, Inf, Inf, Inf, 1),
  hessian=T,                    # get Information matrix
  X=bogus.2d
)
mle.mvn

SE.mean <- sqrt(diag(solve(-1 * mle.mvn$hessian)))[1:2]
SE.mean