3

I would like to minimize a function with two variables.

First I have made a function (rba), which is needed inside the function (kvasum) that I need to minimize. The values to minimized over are a part of rba.

# Data
vpk = data.frame(V1 =c(3650000000, 19233, 2211.2, 479.47, 168.46, 83.447, 52.349, 38.738,
                     32.34, 29.588), V2 = 1:10)
n = nrow(vpk)

# functions to minimize 

# This function returns a vector with 10 values

rba = function(par){
  v <- matrix(ncol = 1, nrow = 10)
  for (p in 1:10){
    k<- ifelse (par[1] < 1-1/p && par[1]>0 && p > par[2] && 
par[2]>0 && par[2]<2, par[2]*p,
                    ifelse(par[1] < 1-1/par[2] && par[1] > 0 && 
p < par[2] && par[2]>0 && par[2]<2, -1+(par[1]+1/par[2]),
                           ifelse(par[1] > (1 - 1 / max(p,par[2])) && 
par[2]>0 && par[2]<2, -1+p, "error")))
    v[p] <- k 
  }
  return(v)
}

# This function uses the function rba, and returns a value

kvasum = function(par){
  sum( (log(vpk$V1)/log(1/n) - rba(par) )^2)
}

# what I would I to do is to find par[1] and par[2] such that kvasum is minimized

m1 = optim(par=c(0.1,0.4),kvasum, lower=0)

I've tried to use the optim function, but I can't get it to work. I get a non-numeric argument, and have tried everything I can think of. Any help is appreciated.

Mathy
  • 167
  • 1
  • 8
  • replace `&&` with `&`. See [here](http://stackoverflow.com/questions/6558921/r-boolean-operators-and) for the difference. – user227710 Jun 18 '15 at 19:45
  • Does "tried everything I can think of" include: using `control=list(trace=...))` arg to optim, or printing out the parameters in the function, or checking what happens if the parameters are negative... – Spacedman Jun 18 '15 at 19:49

1 Answers1

5

There are a few problems with your overall process that cause trouble.

First of all as @user227710 mentions in the comments you should replace && with &. These have different meanings.

Now for the optimiser

It looks like you want to set limits to your parameters (i.e. what is known as box constraints). In order to do this and thus use the lower argument, you need to use the L-BFGS-B method. When you use this you also need to specify the upper argument as well.

The error that you receive, you receive it because your ifelse statement only works when values are between 0 and 1 roughly. Otherwise, the k variable gets the value error (which is the value returned if all conditions in the ifelse statements are FALSE) which is why you get the

Error in log(vpk$V1)/log(1/n) - rba(par) : 
  non-numeric argument to binary operator

error.

Therefore, if you specify your box constraints accordingly (or maybe have a look at your ifelse statement because you might have coded it wrong) this seems to work perfect:

# Data
vpk = data.frame(V1 =c(3650000000, 19233, 2211.2, 479.47, 168.46, 83.447, 52.349, 38.738,
                       32.34, 29.588), V2 = 1:10)
n = nrow(vpk)

# functions to minimize 

# This function returns a vector with 10 values

rba = function(par){
  v <- matrix(ncol = 1, nrow = 10)
  for (p in 1:10){
    k<- ifelse (par[1] < 1-1/p & par[1]>0 & p > par[2] & 
                  par[2]>0 & par[2]<2, par[2]*p,
                ifelse(par[1] < 1-1/par[2] & par[1] > 0 & 
                         p < par[2] & par[2]>0 & par[2]<2, -1+(par[1]+1/par[2]),
                       ifelse(par[1] > (1 - 1 / max(p,par[2])) & 
                                par[2]>0 & par[2]<2, -1+p, "error")))
    #I am adding a line here so that you know why the optim failed
    if(k=='error') stop('your ifelse function returned an error')
    v[p] <- k 
  }
  return(v)
}

# This function uses the function rba, and returns a value

kvasum = function(par){
  sum( (log(vpk$V1)/log(1/n) - rba(par) )^2)
}

# what I would I to do is to find par[1] and par[2] such that kvasum is minimized

m1 = optim(par=c(0.1,0.4),kvasum, method='L-BFGS-B', lower= c(0.1,0.1), upper=c(0.9,0.9))

Output:

> m1
$par
[1] 0.1 0.1

$value
[1] 171.5774

$counts
function gradient 
       2        2 

$convergence
[1] 0

$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"
LyzandeR
  • 34,139
  • 12
  • 63
  • 78