R defining a function using optim

I wish to write a function in R that fits a quantile function to two given tertiles, t1, t2. If I know what the quantile function is, qgamma, say, then I can write:

 gapfn <- function(par, t1, t2) {
  a <- t1 - qgamma(1 / 3, shape = par[1], rate = par[2])
  b <- t2 - qgamma(2 / 3, shape = par[1], rate = par[2])
  return(a ^ 2 + b ^ 2)
}

result <- optim(par = c(10,1), gapfn, t1, t2)

and then result gives me what I want.

But what I really want to do is to write a function that looks like this:

tertile_fit <- function(t1, t2, qfun, par, par_start)

where qfun can be any quantile function defined by parameters par, par_start are starting values for the optimisation, and where the output of tertile_fit is a list as produced by optim.

The reason I want such a function is that I wish to try fitting a number of different quantile functions and want to avoid writing multiple slightly varying pieces of code.

I have tried to write such a function using

gapfn <- function(par, t1, t2, qfun) {
  a <- t1 - qfun(1 / 3, par)
  b <- t2 - qfun(2 / 3, par)
  return(a ^ 2 + b ^ 2)
}

to measure the closeness of fit, and then use

result <- optim(par = par_start, 
                  gapfn, 
                  T1 = T1, 
                  T2 = T2,
                  qfun = qfun)

But that does not work:

tertile_fit <- function(t1, t2, qfun, par_start) {
  gapfn <- function(par, t1, t2, qfun) {
    a <- t1 - qfun(1 / 3, par)
    b <- t2 - qfun(2 / 3, par)
    return(a ^ 2 + b ^ 2)
  }
  result <- optim(
    par = par_start,
    gapfn,
    t1,
    t2,
    qfun = qfun
  )
  return(result)
}

tertile_fit(t1 = 0.5, t2 = 2, qfun = qgamma, par_start = c(1, 10))

delivers:

Error in fn(par, ...) : argument "t2" is missing, with no default

What should I do?

You may be looking for do.call:

tertile_fit <- function(t1, t2, qfun, par_start) {
  gapfn <- function(par, t1, t2, qfun) {
    a <- t1 - do.call(qfun, as.list(c(1/3, par)))
    b <- t2 - do.call(qfun, as.list(c(2/3, par)))
    a ^ 2 + b ^ 2
  }
  optim(par = par_start, gapfn, t1 = t1, t2 = t2, qfun = qfun)
}

par <- tertile_fit(t1 = 0.5, t2 = 2, qfun = qgamma, par_start = c(1, 10))$par
#> Warning in (function (p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE, :
#> NaNs produced
#> Warning in (function (p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE, :
#> NaNs produced
#> Warning in (function (p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE, :
#> NaNs produced
#> Warning in (function (p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE, :
#> NaNs produced
qgamma((1:2)/3, par[1], par[2])
#> [1] 0.4998295 2.0000464

The warnings are from searching negative numbers. This can be fixed by optimizing on the log of the parameters:

tertile_fit <- function(t1, t2, qfun, par_start) {
  gapfn <- function(par, t1, t2, qfun) {
    par <- exp(par)
    a <- t1 - do.call(qfun, as.list(c(1/3, par)))
    b <- t2 - do.call(qfun, as.list(c(2/3, par)))
    a ^ 2 + b ^ 2
  }
  optim(par = log(par_start), gapfn, t1 = t1, t2 = t2, qfun = qfun)
}

par <- exp(tertile_fit(t1 = 0.5, t2 = 2, qfun = qgamma, par_start = c(1, 10))$par)
qgamma((1:2)/3, par[1], par[2])
#> [1] 0.5000493 1.9999386

Or if qfun is vectorized:

tertile_fit <- function(t1, t2, qfun, par_start) {
  gapfn <- function(par, t1, t2, qfun) {
    sum((c(t1, t2) - do.call(qfun, c(list((1:2)/3), as.list(exp(par)))))^2)
  }
  optim(par = log(par_start), gapfn, t1 = t1, t2 = t2, qfun = qfun)
}

par <- exp(tertile_fit(t1 = 0.5, t2 = 2, qfun = qgamma, par_start = c(1, 10))$par)
qgamma((1:2)/3, par[1], par[2])
#> [1] 0.5000493 1.9999386