# R code for examples accompanying lecture 10, 36-350, Fall 2013
# See lecture slides for context

# Load data, create missing column
gmp <- read.table("http://www.stat.cmu.edu/~cshalizi/statcomp/13/lectures/04/gmp.dat")
gmp$pop <- gmp$gmp/gmp$pcgmp


############### Example of using optim() ############################

# MSE function for the West et al. model
  # Input: a length-two vector of parameters (theta)
  # Output: the MSE at those parameters
mse <- function(theta) { mean((gmp$pcgmp - theta[1]*gmp$pop^theta[2])^2) }
require(numDeriv) # Numerical derivative package for calculating gradient
# Gradient of the West et al. MSE
  # Input: a length-two vector of parameters (theta)
  # Output: a length-two vector of partial derivatives
grad.mse <- function(theta) { grad(func=mse,x=theta) }
# An initial parameter guess
theta0=c(5000,0.15)
# Minimize the MSE, using the Newton-ish BFGS method
system.time(fit1 <- optim(theta0,mse,grad.mse,method="BFGS",hessian=TRUE))
  # Running time 0.027 seconds
fit1  # Display the output of the minimization: argmin ("par"), min ("value"),
      # information about the optimization, Hessian matrix

# Exercise: run this again, without giving a gradient argument to optim() ---
# how different is it?
# Exercise: change the starting value theta0 and run it again --- how different
# is it?


################# Example of using nls() ###############################
# Fit the same model with nls()
  # order of arguments: formula, data frame, starting values
  # note: different starting values than before
system.time(fit2 <- nls(pcgmp~y0*pop^a,data=gmp,start=list(y0=5000,a=0.1)))
  # Running time 0.049 seconds
summary(fit2)  # Display information about the fitted model
  # Relative difference between the two sets of parameter estiates is < 10^-4



################ Visual comparison of the two fits ######################
# Plot the data, with the horizontal axis (population) on a log scale
plot(pcgmp~pop,data=gmp,log="x")
pop.order <- order(gmp$pop) # Sort the cities in order of population
# Add a solid black line based on fit2
lines(gmp$pop[pop.order],fitted(fit2)[pop.order])
  # What happens if you run
  #### lines(gmp$pop,fitted(fit2))
  # and why?
# Add a dashed blue line based on fit1
curve(fit1$par[1]*x^fit1$par[2],add=TRUE,lty="dashed",col="blue")