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!