tests/testthat/kde_functionality.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

test_that("univariate default case works as density()", {
  xx = rnorm(1000, sd=1.0)
  kdx = kdeDensity(xx)
  dx = density(xx, kernel="gaussian", from=-5, to=5, n = 101)
  ys = kdx(dx$x)

  # plot(dx)
  # curve(kdx, from=-4, to=4, col=2, add=TRUE)

  expect_equal(sum(ys)*(dx$x[2]-dx$x[1]), 1, tolerance=1e-2)
  expect_equal(ys, dx$y, tolerance = 1e-2)
})



test_that("multivariate default case works as expected", {
  # perhaps change to a Gauss-Quadrature?
  set.seed(4445)
  xx = cbind(rnorm(100, sd=2.5), rnorm(100, sd=0.5), rnorm(100, sd=0.25))
  kdxy = kdeDensity(xx)
  # check that the integral is approx 1:
  x0 = seq(from=-2.5, to=2.5, length.out=25)
  xys = expand.grid(x0*2.5, x0*0.5, x0*0.25)
  zs =  kdxy(xys)
  dx0 = x0[2]-x0[1]
  A = sum(zs)*dx0^3*2.5*0.5*0.25 # bad approx but within expectations given low precision
  expect_equal(A, 1, tolerance=1e-1)

  # compare with kde from stats::density for one marginal:
  dim(zs) = c(25,25,25)
  dx1 = density(xx[,1], kernel="gaussian", from=-2.5*2.5, to=2.5*2.5, n=25)$y
  dx2 = apply(zs,1,sum)*dx0^2*0.5*0.25/A
  expect_equal(dx1, dx2, tolerance=1e-1)

  dx1 = density(xx[,2], kernel="gaussian", from=-2.5*0.5, to=2.5*0.5, n=25)$y
  dx2 = apply(zs,2,sum)*dx0^2*2.5*0.25/A
  expect_equal(dx1, dx2, tolerance=1e-1)

})