Laplace kernel density ratio estimation

laplace_densityratio

# density ratio estimation using ulsif with laplace kernel
library(kernlab)
library(rmutil)
set.seed(45)

# two samples
x <- matrix(rlaplace(10000))
y <- matrix(rnorm(10000, sd = 0.707))

# kernel stuff
kern <- laplacedot(sigma = 1)
nquantile <- function(x, n) quantile(x, cumsum(rep(1/(n+1), n)))
n_centers <- 100
centers <- matrix(nquantile(x, n_centers))

# create kernel gram matrices
phi_x <- kernelMatrix(kern, x, centers)
phi_y <- kernelMatrix(kern, y, centers)

# compute using ulsif
H <- crossprod(phi_x) / length(x)
h <- crossprod(phi_y, rep(1, length(y))) / length(y)
lambda <- 1e-2
theta_hat <- solve(H + lambda * diag(n_centers), h)

# plot
plot(
  x = seq(-3, 3, 0.01),
  y = pmax(0, kernelMult(kern, matrix(seq(-3, 3, 0.01)), centers, theta_hat)),
  type = "l",
  lty = 2,
  main = "Density ratio estimation",
  xlab = "x",
  ylab = "Density",
  bty = "L"
)
curve(dnorm(x, sd = 0.707) / dlaplace(x), add = TRUE)
mtext("Lambda = 1e-2, sigma = 1, Laplace kernel", line = 0.5)
legend("topright",
       lty = c(1, 2),
       legend = c("True ratio", "estimated ratio"))