| 
  • If you are citizen of an European Union member nation, you may not use this service unless you are at least 16 years old.

  • Want to get organized in 2022? Let Dokkio put your cloud files (Drive, Dropbox, and Slack and Gmail attachments) and documents (Google Docs, Sheets, and Notion) in order. Try Dokkio (from the makers of PBworks) for free. Available on the web, Mac, and Windows.

View
 

RSouls20October2006

Page history last edited by PBworks 15 years, 3 months ago

 
# 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) )

Comments (0)

You don't have permission to comment on this page.