Hello!
I am trying to create an R optimization routine for a task that's
currently being done using Excel (lots of tables, formulas, and
Solver).
However, otpim seems to be finding a local minimum.
Example data, functions, and comparison with the solution found in
Excel are below.
I am not experienced in optimizations so thanks a lot for your advice!
Dimitri
### 2 Inputs:
IV<-data.frame(IV=c(0,6672895.687,13345791.37,20018687.06,26691582.75,33364478.44,40037374.12,46710269.81,53383165.5,60056061.18,66728956.87,73401852.56,80074748.24,86747643.93,93420539.62,100093435.3,106766331,113439226.7,120112122.4,126785018.1,133457913.7,140130809.4,146803705.1,153476600.8,160149496.5,166822392.2,173495287.9,180168183.5,186841079.2,193513974.9,200186870.6))
DV<-data.frame(DV=c(0,439.8839775,829.7360945,1176.968757,1487.732038,1767.147276,2019.49499,2248.366401,2456.78592,2647.310413,2822.109854,2983.033036,3131.661246,3269.352233,3397.276321,3516.446162,3627.741311,3731.928591,3829.679009,3921.581866,4008.156537,4089.862363,4167.106955,4240.253215,4309.625263,4375.513474,4438.178766,4497.856259,4554.75841,4609.077705,4660.988983))
## Function "transformIV" transforms a data frame column
"IV" using
parameters .alpha & .beta:
## It returns a data frame column IV_transf:
transformIV = function(.alpha,.beta) {
IV_transf <- as.data.frame(1 - (1/exp((IV/.beta)^.alpha)))
return(IV_transf)
}
### Function "mysum" calculates the sum of absolute residuals after a
regression with a single predictor:
mysum<- function(myIV,myDV){
regr<-lm(myDV[[1]] ~ 0 + myIV[[1]])
mysum<-sum(abs(regr$resid))
return(mysum)
}
### Function to be optimized;
### param is a vector of 2 values (.alpha and .beta)
myfunc <- function(param){
myalpha<-param[1]
mybeta<-param[2]
IVtransf<-transformIV(myalpha, mybeta)
sumofdevs<-mysum(myIV=IVtransf,myDV=DV)
return(sumofdevs)
}
# Optimizing using optim:
myopt <- optim(fn=myfunc, par=c(0.1,max(IV)), method="L-BFGS-B",
lower=0)
(myopt)
myfunc(myopt$par)
## Comparing this solution to Excel Solver solution:
myfunc(c(0.888452533990788,94812732.0897449))
--
Dimitri Liakhovitski
marketfusionanalytics.com
Just to add: I also experimented with the starting parameters (par) under optim, especially with the second one. I tried 1, 10, 100, 1000, etc. When I tried 100,000,000 then I got a somewhat better solution (but still not as good as in Excel). However, under message it said: "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH" Dimitri On Thu, Nov 10, 2011 at 1:50 PM, Dimitri Liakhovitski <dimitri.liakhovitski at gmail.com> wrote:> Hello! > > I am trying to create an R optimization routine for a task that's > currently being done using Excel (lots of tables, formulas, and > Solver). > However, otpim seems to be finding a local minimum. > Example data, functions, and comparison with the solution found in > Excel are below. > I am not experienced in optimizations so thanks a lot for your advice! > Dimitri > > ### 2 Inputs: > IV<-data.frame(IV=c(0,6672895.687,13345791.37,20018687.06,26691582.75,33364478.44,40037374.12,46710269.81,53383165.5,60056061.18,66728956.87,73401852.56,80074748.24,86747643.93,93420539.62,100093435.3,106766331,113439226.7,120112122.4,126785018.1,133457913.7,140130809.4,146803705.1,153476600.8,160149496.5,166822392.2,173495287.9,180168183.5,186841079.2,193513974.9,200186870.6)) > DV<-data.frame(DV=c(0,439.8839775,829.7360945,1176.968757,1487.732038,1767.147276,2019.49499,2248.366401,2456.78592,2647.310413,2822.109854,2983.033036,3131.661246,3269.352233,3397.276321,3516.446162,3627.741311,3731.928591,3829.679009,3921.581866,4008.156537,4089.862363,4167.106955,4240.253215,4309.625263,4375.513474,4438.178766,4497.856259,4554.75841,4609.077705,4660.988983)) > > ## Function "transformIV" transforms a data frame column "IV" using > parameters .alpha & .beta: > ## It returns a data frame column IV_transf: > transformIV = function(.alpha,.beta) { > ?IV_transf <- as.data.frame(1 - (1/exp((IV/.beta)^.alpha))) > ?return(IV_transf) > } > > ### Function "mysum" calculates the sum of absolute residuals after a > regression with a single predictor: > mysum<- function(myIV,myDV){ > ?regr<-lm(myDV[[1]] ~ 0 + myIV[[1]]) > ?mysum<-sum(abs(regr$resid)) > ?return(mysum) > } > > ### Function to be optimized; > ### param is a vector of 2 values (.alpha and .beta) > myfunc <- function(param){ > ?myalpha<-param[1] > ?mybeta<-param[2] > ?IVtransf<-transformIV(myalpha, mybeta) > ?sumofdevs<-mysum(myIV=IVtransf,myDV=DV) > ?return(sumofdevs) > } > > # Optimizing using optim: > myopt <- optim(fn=myfunc, par=c(0.1,max(IV)), method="L-BFGS-B", lower=0) > (myopt) > myfunc(myopt$par) > ## Comparing this solution to Excel Solver solution: > myfunc(c(0.888452533990788,94812732.0897449)) > > -- > Dimitri Liakhovitski > marketfusionanalytics.com >-- Dimitri Liakhovitski marketfusionanalytics.com
I won't requote all the other msgs, but the latest (and possibly a bit glitchy) version of optimx on R-forge 1) finds that some methods wander into domains where the user function fails try() (new optimx runs try() around all function calls). This includes L-BFGS-B 2) reports that the scaling is such that you really might not expect to get a good solution then 3) Actually gets a better result than the> xlf<-myfunc(c(0.888452533990788,94812732.0897449)) > xlf[1] 334.607>with Kelley's variant of Nelder Mead (from dfoptim package), with> myoptxmethod par fvalues fns grs itns conv KKT1 4 LBFGSB NA, NA 8.988466e+307 NA NULL NULL 9999 NA 2 Rvmmin 0.1, 200186870.6 25593.83 20 1 NULL 0 FALSE 3 bobyqa 6.987875e-01, 2.001869e+08 1933.229 44 NA NULL 0 FALSE 1 nmkb 8.897590e-01, 9.470163e+07 334.1901 204 NA NULL 0 FALSE KKT2 xtimes meths 4 NA 0.01 LBFGSB 2 FALSE 0.11 Rvmmin 3 FALSE 0.24 bobyqa 1 FALSE 1.08 nmkb But do note the terrible scaling. Hardly surprising that this function does not work. I'll have to delve deeper to see what the scaling setup should be because of the nature of the function setup involving some of the data. (optimx includes parscale on all methods). However, original poster DID include code, so it was easy to do a quick check. Good for him. JN> ## Comparing this solution to Excel Solver solution: > myfunc(c(0.888452533990788,94812732.0897449)) > > -- Dimitri Liakhovitski marketfusionanalytics.com
Hi Dimitri, Your problem has little to do with local versus global optimum. You can convince yourself that the solution you got is not even a local optimum by checking the gradient at the solution. The main issue is that your objective function is not differentiable everywhere. So, you have 2 options: either you use a smooth objective function (e.g. squared residuals) or you use an optimization algorithm than can handle non-smooth objective function. Here I show that your problem is well solved by the `nmkb' function (a bound-constraints version of Nelder-Mead simplex method) from the "dfoptim" package. library(dfoptim)> myopt2 <- nmkb(fn=myfunc, par=c(0.1,max(IV)), lower=0) > myopt2$par [1] 8.897590e-01 9.470163e+07 $value [1] 334.1901 $feval [1] 204 $restarts [1] 0 $convergence [1] 0 $message [1] "Successful convergence" Then, there is also the issue of properly scaling your function, because it is poorly scaled. Look how different the 2 parameters are - they are 7 orders of magnitude apart. You are really asking for trouble here. Hope this is helpful, Ravi. ------------------------------------------------------- Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvaradhan@jhmi.edu<mailto:rvaradhan@jhmi.edu> [[alternative HTML version deleted]]