Minimal fitting example

Read in data from some text file

d <- read.table("peaks.txt", head = TRUE)
head(d)

##   wavenumber intensity
## 1        400     98.32
## 2        401     98.49
## 3        402     98.86
## 4        403     98.51
## 5        404     99.28
## 6        405     99.32

str(d)

## 'data.frame':    401 obs. of  2 variables:
##  $ wavenumber: int  400 401 402 403 404 405 406 407 408 409 ...
##  $ intensity : num  98.3 98.5 98.9 98.5 99.3 ...

plot(d)

Define the model and objective (cost) functions

model  <- function(p, x) { # lorentzian
  p[3] / pi * p[2] / ((x - p[1])^2 + p[2]^2) + p[4]
}

objective  <- function(p, data=NULL, # pass a data.frame, or alternatively x and y
                       x=data[,1], y=data[,2], ...){
    predicted <- model(p, x, ...)
    sum((predicted - y)^2)
}

Fitting using \Sexpr{optim()}

guess  <-  c(500, 5, 500, 100) # initial search point

fit  <-  optim(guess, objective, 
  data = subset(d, wavenumber < 550 & # fit only over a limited range
                wavenumber > 450))

# combine fit and guess into a data.frame
res <- transform(d, guess = model(guess, wavenumber),
                 fit = model(fit$par, wavenumber))

Plotting the results against the data and initial guess

with(res, matplot(wavenumber, cbind(intensity, guess, fit), t = "n", lty = c(1, 
    2, 1), ylab = "Intensity", xlab = expression(Wavelength/nm)))
rect(450, -1000, 550, 1000, col = "grey95", border = NA)  # fit region
with(res, matlines(wavenumber, cbind(intensity, guess, fit), t = "l", lty = c(1, 
    2, 1)))