# optimExample.R
# Author: Markus Gesmann
# Demonstrate the use of the optim function.
# We create data and try to fit a linear curve to it, via minimsing 
# the residual sum of squares

# Create data
df=data.frame(y=c(1,3,5,6,8,12), x=c(1,2,3,4,5,6))

# Create a function which will calculate the residual some of squares
# for a given data and parameter set.

min.RSS <- function(data, par){
 y <- data$y
 x <- data$x
 sum( ( par[1]+par[2]*x - y)^2 )
}

# Use the optim function:
result <- optim(par=c(0,1), min.RSS, data=df)

# Look at the result, you'll find the optimised parameters in result$par
# the minimised RSS is in result$value

result

plot(y ~ x, data=df)
abline(a=result$par[1], b=result$par[2], col="red")


# Let's compare this against linear regression:

lm(y ~ x, data=df)


## Second example:
## Count data, so fit a Poisson distribution

x=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 17, 42, 43)
y=c(1392, 1711,  914,  468, 306, 192, 96, 56, 35, 17, 15, 6, 2, 2, 1, 1)
plot(y ~ x, t='h')
y.rel =y/sum(y)

## Let's try the curve fitting approach to estimate the parameter lambda
pois.opt <- function(data, par){ sum( ( dpois(data$x, par) - data$y )^2 )}
result = optim(par=2, pois.opt, data=data.frame(x,y.rel) )
result

plot(y.rel ~ x, t='h', xlim=c(0,13))
lines(dpois(x, result$par), lwd=2, col="red")

## Compare the result with fitdistr which uses maximum liklihood
library(MASS)
fitdistr(rep(x,y), "Poisson")

# No surprise it gives the mean back
mean( rep(x,y) )

# Let's do the maximum liklihood approach with optim.
# The log-liklihood function is
poisson <- function(x, lambda) lambda^x/factorial(x)*exp(-lambda)

# and volia, we get the fitdistr result (Beware of the minus sign in front of sum,
# as we want the MAXIMUM liklihood):
optim( par=2, function(x, par){-sum(log(poisson(x,par)))}, x=rep(x,y) )