# Calculating Jensen-Shannon and Kullback-Leibler for density splines

I need to find a way to calculate the JS distance (and, by necessity, the KL divergence) between continuous, non-normal (generating distribution unknown) distributions.

I have many classes, and many variables for each class, and what I want to do is measure the distance between each class' distribution of each variable.

My data is sparse for many of my classes, and is truncated. The best method I've found for getting densities for the distributions, is to calculate log splines using the oldlogspline function from the R logspline package, since it supports truncated distributions. The output is, for each distribution, a vector of the locations of spline knots and a set of coefficients for the spline.

What I need to do is (a) calculate the KL divergence between splines representing the distributions, and (b) estimate a spline density as the average of two other spline densities.

Both (a) and (b) exceed my meager calculus skills.

I would greatly appreciate if anyone can tell me how to do this (preferrably with R code! ;) ).

Here is some sample code to generate an example:

library(logspline)

ab1 <- c(0.5, 3.6) + 5
ab2 <- c(0.3, 2.8) + 5
set1 <- rnorm(n = 5000, mean = sample(ab1, size=5000, replace = T))
hist(set1, breaks=50)

set2 <- rnorm(n = 5000, mean = sample(ab2, size=5000, replace = T))
hist(set2, breaks=50)

breakpoint <- 3

set1 <- set1[set1 > breakpoint]
set2 <- set2[set2 > breakpoint]

spline1 <- oldlogspline(set1, lbound=breakpoint)
plot(spline1)

spline2 <- oldlogspline(set2, lbound=breakpoint)
plot(spline2)

Edit:


I could solve this if I could find the antiderivative of:

$\int e^{b_0 + b_1x + \sum_{i \in k_i < x} b_i (x - k_i)^3} \log\frac{e^{b_0 + b_1x + \sum_{i \in k_i < x} b_i (x - k_i)^3}}{e^{c_0 + c_1x + \sum_{i \in j_i < x} c_i (x - j_i)^3}}dx$

Which simplifies to:

$\int e^{b_0 + b_1x + \sum_{i \in k_i < x} b_i (x - k_i)^3} \cdot \left((b_0 - c_0) + (b_1 - c_1)x + \sum_{i \in k_i < x} b_i (x - k_i)^3 - \sum_{i \in j_i < x} c_i (x - j_i)^3\right)dx$

• Your data is numerical and finite, so why the interest in solving the integral in closed form. Once you've create via splines a function, you can do numerical integration. In addition, are you familiar with the Kolmogorov-Smirnov test/theorem ? It gives a way to compare two distributions. It is intrinsically symmetric and may be closer in spirit to what you want. – aginensky May 21 '18 at 18:31
• You're correct that numerical integration is an option. That's what I'm doing at the moment. But the densities can be somewhat unusual. Multi-modal, etc. There are some 14,000 cases across the whole dataset, and I'm concerned that some of them (perhaps many) will be sufficiently pathological that numerical integration with default hyperparameters will be inaccurate. So if I can find an analytic solution, that would be preferable. I'm somewhat familiar with KS, but hadn't considered it for this. If I'm unable to solve it, I may consider alternative distance metrics such as Wasserstein or KS. – Bob May 21 '18 at 20:14
• To elaborate on that a bit, one approach I took was to try to estimate using a gaussian kernel and intervals. The problem that arose was that in many of the cases, there are large blocks of the density range where the interval density is 0. Using splines helps address that issue, but I'm concerned that numerical integration may be hiding the problem rather than solving it. – Bob May 21 '18 at 20:18