0

I'm trying to get the parameters of a bi-exponential model by minimizing the Kullback-Leibler divergence using optim. The function I'm using have 3 parameters, but when I pass them to optim using par = par it throws the error "list object cannot be coerced to type 'double'", but I'm not even using lists.

Here is the code I'm using:

library(logKDE) # for kernel density of positive distributions
# Simulate rv of the bi-exponential
p <- 0.7
n <- 50
w <- 2
b <- 0.2
delt <- 0.01

biexp_data <- (p * rexp(n, 1/w) + (1 - p) * rexp(n, 1/b)) -  delt

# define kld to optimize
kld_optim <- function(x, par, from_a, to_b) {
  par <- unlist(par)
  w <- par[1]
  b <- par[2]
  p <- par[3]
  d <- 0.002
  
  integrand <- function(x, w, b, p, d, t) {
    
    denx <- logdensity(x, bw = 'logG', from = from_a, to = to_b)
    f.y <- approx(unlist(denx$x), unlist(denx$y), t)$y
    f.x <- p * dexp(t - d, rate = 1/w) + (1 - p) * dexp(t - d, rate = 1/b)
    tmpRatio <- f.x * (log2(f.x) - log2(f.y))
    # Return
    ifelse(is.infinite(tmpRatio), 0, ifelse(is.na(tmpRatio), 0, tmpRatio))

  }
  
  integrate(integrand, 
            from_a, to_b, 
            x = x,
            w = w, b = b, p = p, d = d)
  
}

optim(par = c(2, 0.1, 0.6),
      fn = kld_optim, 
      from_a = 0.01, 
      to_b = 20,
      x = biexp_data)

Why is this happening?

Thanks!

1 Answer 1

0

The function you pass to optim needs to return a scalar value, and your kld_optim function returns the result of a call to integrate() which according to the ?integrate help page, returns a list, not a numeric value. The "value" in contained in that list under the name "value". So change your integrate() call to

  integrate(integrand, 
            from_a, to_b, 
            x = x,
            w = w, b = b, p = p, d = d)$value
2
  • Yes, you're right. Now I had another problem. It say "extremely bad integrating behavior". Now, I guess thats a problem with my procedure. Do you know of other methods for intergation than could work?
    – jealcalat
    Aug 11, 2020 at 5:20
  • @jealcalat Try massively increasing the number of subdivisions in integrate.
    – Roland
    Aug 11, 2020 at 6:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.