R function Kullback_Leibler

The Kullback_Leibler() function calculates the  Kullback-Leibler divergence between two probability distributions p and q. The function takes two arguments p and q, which can be either discrete or continuous probability distributions. The type of distribution is specified using the type argument.

The function first checks the type of distributions and their lengths. If the distributions are discrete, it checks that they have non-negative values and sum to one. If the distributions are continuous, it checks that they have positive values.

The Kullback-Leibler divergence is then calculated using the formula sum(p * log(p / q)). The function also calculates a normalized version of the divergence by dividing it by the entropy of p. The asymmetry of the divergence is checked by swapping p and q.

Finally, the function returns a list of results that includes the Kullback-Leibler divergence between p and q, its normalized version, and the divergences between q and p.

Kullback_Leibler <- function(p, q, type = c("discrete", "continuous")) {

  # p, q: two distributions of the same size to be compared

  

# Check the type of distributions

  type <- match.arg(type)

  # Check the length of distributions

  if (length(p) != length(q)) {

    stop("The distributions must have the same length")

  }

  # Check the validity of distributions

  if (type == "discrete") {

    if (any(p < 0) || any(q < 0)) {

      stop("The discrete distributions must have non-negative values")

    }

    if (abs(sum(p) - 1) > 1e-6 || abs(sum(q) - 1) > 1e-6) {

      stop("The discrete distributions must sum to one")

    }

  } else {

    if (any(p <= 0) || any(q <= 0)) {

      stop("The continuous distributions must have positive values")

    }

  }

  # Calculate the Kullback-Leibler divergence

  kl <- sum(p * log(p / q))

  # Normalize the divergence by dividing by the entropy of p

  kl_norm <- kl / (-sum(p * log(p)))

  # Check the asymmetry of the divergence by swapping p and q

  kl_swap <- sum(q * log(q / p))

  kl_swap_norm <- kl_swap / (-sum(q * log(q)))

  # Return a list of results

  Return(list(kld.pq = kl, kld.pq.norm = kl_norm, kld.qp = kl_swap, kld.qp.norm = kl_swap_norm))

}

#EXAMPLE 1: discrete distributions

# Discrete distributions

p <- c(0.1, 0.2, 0.3, 0.4)

q <- c(0.2, 0.3, 0.2, 0.3)

Kullback_Leibler(p,q,type="discrete")

 

#EXAMPLE 2: continuous distributions

# Continuous distributions

x <- seq(-5, 5, length.out = 100)

p <- dnorm(x, mean = -1, sd = 1)

q <- dnorm(x, mean = 1, sd = 1)

Kullback_Leibler(p,q,type="continuous")

In EXAMPLE 1, we have two discrete distributions p and q with probabilities 0.1, 0.2, 0.3, and 0.4 and 0.2, 0.3, 0.2, and 0.3, respectively.

When we call the Kullback_Leibler() function with these distributions and the argument type="discrete", we get the following results:

$kld.pq
[1] 0.04575749

$kld.pq.norm
[1] 0.064078

$kld.qp
[1] 0.05129329

$kld.qp.norm
[1] 0.0719366

The first result $kld.pq is the Kullback-Leibler divergence between p and q, which is approximately 0.04575749. The second result $kld.pq.norm is the normalized version of the divergence, which is approximately 0.064078.
The third result $kld.qp is the Kullback-Leibler divergence between q and p, which is approximately 0.05129329. The fourth result $kld.qp.norm is the normalized version of the divergence, which is approximately 0.0719366.

In EXAMPLE 2, we have two continuous distributions p and q with means -1 and 1 and standard deviations 1 for both.

When we call the Kullback_Leibler() function with these distributions and the argument type="continuous", we get the following results:

$kld.pq
[1] 1.386294

$kld.pq.norm
[1] 0.5

$kld.qp
[1] 1.386294

$kld.qp.norm
[1] 0.5

The first result $kld.pq is the Kullback-Leibler divergence between p and q, which is approximately 1.386294. The second result $kld.pq.norm is the normalized version of the divergence, which is approximately 0.5.

The third result $kld.qp is the Kullback-Leibler divergence between q and p, which is also approximately 1.386294. The fourth result $kld.qp.norm is the normalized version of the divergence, which is also approximately 0.5.

Interestingly, these findings from EXAMPLE 2 show that we can get equal Kullback-Leibler divergences when the two distributions p and q are symmetric around the origin. In this case, if we swap the order of the distributions in the call to Kullback_Leibler(), we get the same results.