# Example of fitting a non-linear model with several parameters by # maximum likelihood, using "deriv" and "nlm". # # The model is for the data on sunspot counts supplied with R. We # hypothesize that the counts (actually integers, but modelled as # non-negative continuous values) come from taking the absolute value # of a normal variate whose mean varies with time according to a sine # wave. The model parameters are a, b, f, and M, which define the # sine wave, as M*(a+sin(b+f*t), and the log of the standard deviation # of the observation, lsigma. # # We use "nlm" to find the maximum likelihood estimates, using the # gradient values automatically produced using "deriv". Note that # "deriv" only knows about "dnorm" with one argument (the standard normal # density), not with threee arguments. # # We'll fit only the first 500 points. # # Question: What happens if we try to fit all the points? How about # using the estimate from the first 500 as an initial value? What # does the result trying to fit all the points say about the adequacy # of this model? # Function computing the log probability of each observations, with # gradient. (Using the hessian too seems to produce worse results.) logp <- deriv (quote ( -lsigma + log (dnorm((x-M*(a+sin(b+f*t)))/exp(lsigma)) + dnorm((-x-M*(a+sin(b+f*t)))/exp(lsigma)))), c("a", "b", "f", "M", "lsigma"), fun=TRUE) # Function for nlm to minimize - minus the log likelihood summing # log probabilities (and gradient and hessian) over data points. logl <- function (p) { lp <- logp(p[1],p[2],p[3],p[4],p[5]) ll <- -sum(lp) attr(ll,"gradient") <- -colSums(attr(lp,"gradient")) ll } # The data. x <- sunspots[1:500] t <- 1:length(x) # Estimation, from a carefully chosen starting point. est <- nlm (logl, c(0, -1, 2*pi/200, 80, 4)) print(est) mle <- est$estimate # Plot the mean from the model with the MLE parameter estimates together # with the data points. a <- mle[1] b <- mle[2] f <- mle[3] M <- mle[4] plot(t,x,pch=20) lines(t,abs(M*(a+sin(b+f*t))),col="red")