##################################################
## R code for section 4
## 2-19-2014
## Stephen Pettigrew


# Writing a Poisson log likelihood function into R:

loglike <- function(parameters, outcome, covariates){
  # Add a column of 1's to your covariate matrix 
  # to account for the intercept
  cov <- as.matrix(cbind(1, covariates))
  
  # To make our code easier to read let's multiply 
  # the covariate matrix, X, by the parameters, beta
  xb <- cov %*% parameters
  
  sum(outcome * xb - exp(xb))
}


# Let's look test whether our function works with the 
# curling data. The x covariate is whether the US
# had the hammer in that end

x <- c(0,1,0,1,0,1,0,1,0,1)
y <- c(1,0,1,2,0,1,0,2,0,1)



# Let's look at whether our function works by plugging
# the two parameters values into the function.

loglike(pars = c(1.2, .57), # these were arbitrarily picked
        outcome = y,
        covariates = x) # It doesn't give us an error so it works!




## Next week in section we'll talk about how to optimize functions,
## but just as a preview, let's find the MLE for the betas based on
##  the data we generated

results <- optim(par = c(0,0),
                 outcome = y, 
                 covariates = x, 
                 fn = loglike,
                 method = "BFGS",
                 control = list(fnscale = -1, maxit=100000), 
                 hessian = TRUE)

results$par ## There are our beta0 and beta1 estimates!