Translate matlab fminsearchbnd in R?

Viewed 135

I am trying to translate a Matlab code to R, finding the minmum error of two functions.

Here is the Matlab code:

erro=@(x) sum(((f14.*x(2)./(((x(1))/0.9924)+x(5))-(f14.*x(2)./(((x(1))/0.9924)+x(5))-x(3)).*exp(-(((x(1))/0.9924)+x(5))*(t-t1)))-mo).^2)./std(mo)+ sum (((f15.*x(2)./((x(1)+x(5)))-(f15.*x(2)./((x(1)+x(5)))-x(4)).*exp(-((x(1)+x(5))).*(t-t1)))-no).^2)/std(no); 

x=fminsearchbnd(erro,xo,[0 0 200 0.9 0],[])

% x = 0.0123128476186327    13.5174782795704    368.926386686419    76.4422946745777    8.28109088192835e-09

So i tried the pracma-package first, unfortunately I could find fminbnd() and fminsearch() function, but no fminsearchbnd(). Is there a way to combine these two? Or can anyone think of different solutions? I tried the optim() function with "method = c("L-BFGS-B") and it gave me the best results so far. However I wonder if there is a way to get exact the same results in Matlab and R.

My translation to R:

t <- c(0.5, 1, 2, 3)
mo <- c(368.54, 374.04, 381.75, 390.92)
no <- c(76.46, 75.96, 75.14, 74.23)
f15 <- c(0.003431)  
f14 <- c(1-f15)    
xo <-  c(0.01211952, 13.30265, 368.54, 76.46, -0.000504064)


erro <- function(x){sum(((f14*x[2]/(((x[1])/0.9924)+x[5])-(f14*x[2]/(((x[1])/0.9924)+x[5])-x[3])*exp(-(((x[1])/0.9924)+x[5])*(t-t[1])))-mo)^2)/std(mo) 
  + sum (((f15*x[2]/((x[1]+x[5]))-(f15*x[2]/((x[1]+x[5]))-x[4])*exp(-((x[1]+x[5]))*(t-t[1])))-no)^2)/std(no)} 

x <- optim(xo, erro, lower = c(0, 0, 200, 0.9, 0), method = c("L-BFGS-B"))
# $par
#[1]   0.01230393  13.30265160 368.54000000  76.44231900   0.00000000

Thanks for your help and please forgive any mistakes, this is my first post and I have no idea about matlab! ;)

Update: I wanted to compare the optim-results with those from MATLAB for different values, so I tried with:

mo <- c(322.79, 327.06, 323.11, 322.60)
no <- c(70.01023605, 64.18709552, 59.03694208, 55.83444661)
lo <- c(0.2573, 0.2525, 0.2398, 0.2388)
xo <- c(0.07856273, 24.87679, 322.79, 70.01024, 0.03182682)

I got the error 'L-BFGS-B needs finite values of 'fn'' so I set my first bound up from 0 to 0.00001, resulting in:

[1,] 1e-05 24.91 322.79 68.59785 0.09070217

However, my first value always sticks to the lower bound, whereas Matlab results are:

0.05460 28.88570 324.74040 68.58710 0.03620
1 Answers

The region where the minimum lives is very flat. Therefore, to get good results in R, one could add an explicit numerical gradient, as the gradient used by optim() may not be accurate enough. Here we employ the gradient function from the 'numDeriv' package.

f15 <- c(0.003431); f14 <- c(1-f15)
t <- c(0.5, 1, 2, 3)

mo <- c(322.79, 327.06, 323.11, 322.60)
no <- c(70.01023605, 64.18709552, 59.03694208, 55.83444661)
lo <- c(0.2573, 0.2525, 0.2398, 0.2388)

# starting parameters
xo <- c(0.07856273, 24.87679, 322.79, 70.01024, 0.03182682)

erro <- function(x){
    sum(((f14*x[2]/(((x[1])/0.9924)+x[5])-(f14*x[2]/(((x[1])/0.9924)+x[5])-x[3])*exp(-(((x[1])/0.9924)+x[5])*(t-t[1])))-mo)^2)/sd(mo) + sum (((f15*x[2]/((x[1]+x[5]))-(f15*x[2]/((x[1]+x[5]))-x[4])*exp(-((x[1]+x[5]))*(t-t[1])))-no)^2)/sd(no)
} 

# minimum value at MATLAB solution xm
xm <- c(0.05460, 28.88570, 324.74040, 68.58710, 0.03620)
erro(xm)
## [1] 6.367856

# numerical gradient
erro_gr <- function(x) numDeriv::grad(erro, x)

s <- optim(xo, erro, gr = erro_gr,
           lower = c(1.e-07, 0, 200, 0.9, 0),
           method = c("L-BFGS-B"))
s$par;s$value
## [1] 0.08281183 28.97220 324.7395 68.58882 0.008001717
## [1] 6.367579

This is even smaller than the MATLAB solution 6.367856, see above.

Related