Solved – How to use optim function in R on the custom Residual Sum of Squares function

by

I'm trying to use the optim function in R to find the optimal values for $alpha$ and $l_0$ for my custom residual sum of squares (RSS) function, but the values I'm getting back aren't even close to those picked by R's ses function, which were $alpha$ = 0.2971 and $l_0$ = 77260, for the pigs dataset in Rob Hyndman's fpp2 library.

Here's my custom RSS function:

custom_rss <- function(ts, alpha, l0) {   N <- length(ts)   res <- vector()   for (j in 1:N-1) {     y_hat <- custom_ses(ts[1:j], alpha, l0)     res[j] <- (ts[j+1] - y_hat)^2   }   return(sum(res)) }

My custom function for simple exponential smoothing is:

custom_ses <- function(ts, alpha, l0) {   N <- length(ts)   y_hat = 0   for (j in 1:N-1) {     y_hat <- y_hat + alpha * ((1 - alpha)^j) * ts[N-j] + ((1 - alpha)^N)*l0   }   return(y_hat) }

And finally, here's how I used the optim function:

optim(par=c(0,1), fn=custom_rss, alpha=0.2971, l0=77260.0561)

So, I'm not sure where the problem lies: my SES function, my RSS function, or the way I'm using the optim function.

Contents hide
1 Best Answer
2 Similar Posts:

Best Answer

I think your objective and smoothing functions are not correct.

More precisely, the problems I found are:

  • The way you pass the parameters to the optimizer
  • The absence of constraints on the parameters ($alpha in [0, 1])$
  • The formula for y_hat. I used the Weighted Average Form in https://otexts.org/fpp2/ses.html#eq:7-ses

The parameters to optim() should be passed through the vector par. Then, from this vector, you can extract alpha and l0 in your objective function definition. Second, I think it is appropriate to estimate those parameter in such a way you are sure they will respect the necessary constraints. This can be achieved by transforming their values within the optimization routine.

Below you will find a modified example inspired by your R codes above

rm(list = ls()); cat("14"); graphics.off() library(fpp2)  # Objective function custom_rss = function(par, ts)  {   # Be seure to have alpha in [0,1]   alpha = plogis(par[1])   l0 = par[2]    N = length(ts)   fit = custom_ses(ts, alpha, l0)   res = ts - fit    out = sum(res^2)/(N-1)   return(out) }  # Smoother function       custom_ses = function(ts, alpha, l0)  {   N = length(ts)   y_hat = c(l0, 0 * (2:N))   for (j in 1:(N-1))   {     new = alpha * ts[j] + (1-alpha) * y_hat[j]     y_hat[j+1] = new     # print(y_hat)   }    return(y_hat) }  # The data data("pigs") ts = pigs  # optim-based estimates  est = optim(par = c(0.1, 5000), fn = custom_rss, ts = ts) alpha = plogis(est$par[1]) l0 = est$par[2] fit = custom_ses(ts, alpha = alpha, l0 = l0)  # Ses estimates sesMod = ses(ts)  # Check cbind(ses = sesMod$model$par, custom = c(alpha, l0))  # Plots plot(1:length(ts), ts, type = "l", lwd = 1, col = 8) points(1:length(ts), ts, pch = 16, col = 8, cex = 0.5) lines(1:length(ts), sesMod$fitted, col = 2, type = "l", lwd = 2) lines(1:length(ts), fit, col = 3, type = "l", lty = 2,  lwd = 2) legend("bottomright", legend = c("Data", "ses", "custom"), col = c(8, 2, 3),         lty = c(1,1,2), lwd = 2)

enter image description here

Similar Posts:

Rate this post

Leave a Comment