Jensen-Shannon divergence calculation for 3 prob distributions: Is this ok?

I would like to calculate the jensen-shannon divergence for he following 3 distributions. Is the calculation below correct? (I followed the JSD formula from wikipedia):

P1  a:1/2  b:1/2    c:0
P2  a:0    b:1/10   c:9/10
P3  a:1/3  b:1/3    c:1/3
All distributions have equal weights, ie 1/3.

JSD(P1, P2, P3) = H[(1/6, 1/6, 0) + (0, 1/30, 9/30) + (1/9,1/9,1/9)] -
[1/3*H[(1/2,1/2,0)] + 1/3*H[(0,1/10,9/10)] + 1/3*H[(1/3,1/3,1/3)]]

JSD(P1, P2, P3) = H[(1/6, 1/5, 9/30)] - [0 + 1/3*0.693 + 0] = 1.098-0.693 = 0.867


Thanks in advance...

EDIT Here's some simple dirty Python code that calculates this as well:

    def entropy(prob_dist, base=math.e):
return -sum([p * math.log(p,base) for p in prob_dist if p != 0])

def jsd(prob_dists, base=math.e):
weight = 1/len(prob_dists) #all same weight
js_left = [0,0,0]
js_right = 0
for pd in prob_dists:
js_left[0] += pd[0]*weight
js_left[1] += pd[1]*weight
js_left[2] += pd[2]*weight
js_right += weight*entropy(pd,base)
return entropy(js_left)-js_right

usage: jsd([[1/2,1/2,0],[0,1/10,9/10],[1/3,1/3,1/3]])

• Nice Python code by the way! – gui11aume Jun 1 '12 at 23:55

6 Answers

There is mistake in the mixture distribution. It should be $(5/18, 28/90, 37/90)$ instead of $(1/6, 1/5, 9/30)$ which does not sum up to 1. The entropy (with natural log) of that is 1.084503. Your other entropy terms are wrong.

I will give the detail of one computation:

$$H(1/2,1/2,0) = -1/2*\log(1/2) - 1/2*\log(1/2) + 0 = 0.6931472$$

In a similar way, the other terms are 0.325083 and 1.098612. So the final result is 1.084503 - (0.6931472 + 0.325083 + 1.098612)/3 = 0.378889

• +1. Quick and dirty R calculation: h <- function(x) {h <- function(x) {y <- x[x > 0]; -sum(y * log(y))}; jsd <- function(p,q) {h(q %*% p) - q %*% apply(p, 2, h)}. Argument p is a matrix whose rows are the distributions and argument q is the vector of weights. E.g., p <- matrix(c(1/2,1/2,0, 0,1/10,9/10, 1/3,1/3,1/3), ncol=3, byrow=TRUE); q <- c(1/3,1/3,1/3); jsd(p,q) returns $0.378889$ (which approximates the log of $3^{34/15} 5^{1/9} 2^{-13/45} 7^{-14/45} 37^{-37/90}$). – whuber May 31 '12 at 19:05
• Not so dirty... ;-) – gui11aume May 31 '12 at 22:19
• (1) Redo the math. (2) Entropy can be measured using any base of logarithm you like, as long as you are consistent. Natural, common, and base-2 logs are all conventional. (3) It's really a mean discrepancy between the distributions and their average. If you think of each distribution as a point, they form a cloud. You're looking at the average "distance" between the center of the cloud and its points, kind of like a mean radius. Intuitively, it measures the cloud's size. – whuber May 31 '12 at 23:51
• @Legend I think you're right. I did not test sufficiently after finding that one result agreed with the answer I obtained in another way (with Mathematica). – whuber Jan 17 '13 at 3:26
• @dmck There are indeed typos in my comment: (1) the phrase h <- function(x) { was pasted twice. Just delete it: everything else works and produces the results I quote. Then modify the apply(p, 2, h) to apply(p, 1, h) as pointed out in the comment by Legend. – whuber Dec 15 '14 at 22:52

Python:

import numpy as np
# @author: jonathanfriedman

def jsd(x,y): #Jensen-shannon divergence
import warnings
warnings.filterwarnings("ignore", category = RuntimeWarning)
x = np.array(x)
y = np.array(y)
d1 = x*np.log2(2*x/(x+y))
d2 = y*np.log2(2*y/(x+y))
d1[np.isnan(d1)] = 0
d2[np.isnan(d2)] = 0
d = 0.5*np.sum(d1+d2)
return d

jsd(np.array([0.5,0.5,0]),np.array([0,0.1,0.9]))


Java:

/**
* Returns the Jensen-Shannon divergence.
*/
public static double jensenShannonDivergence(final double[] p1,
final double[] p2) {
assert (p1.length == p2.length);
double[] average = new double[p1.length];
for (int i = 0; i < p1.length; ++i) {
average[i] += (p1[i] + p2[i]) / 2;
}
return (klDivergence(p1, average) + klDivergence(p2, average)) / 2;
}

public static final double log2 = Math.log(2);

/**
* Returns the KL divergence, K(p1 || p2).
*
* The log is w.r.t. base 2.
* <p>
* *Note*: If any value in <tt>p2</tt> is <tt>0.0</tt> then the
* KL-divergence is <tt>infinite</tt>. Limin changes it to zero instead of
* infinite.
*/
public static double klDivergence(final double[] p1, final double[] p2) {
double klDiv = 0.0;
for (int i = 0; i < p1.length; ++i) {
if (p1[i] == 0) {
continue;
}
if (p2[i] == 0.0) {
continue;
} // Limin

klDiv += p1[i] * Math.log(p1[i] / p2[i]);
}
return klDiv / log2; // moved this division out of the loop -DM
}


You gave a Wikipedia reference. Here I give the complete expression for Jensen-Shannon divergence with multiple probability distributions:

$$JSmetric(p^1,...,p^m)=H(\frac{p^1+...+p^m}{m})-\frac{\sum_{j=1}^{m} H(p^j)}{m}$$

Original question was posted without mathematical expression of multi-distribution JS divergence which lead to a confusion in understanding the computation provided. Also, term weight was used which again causes confusion that how you select appropriate weights for multiplication. Above expression clarifies these confusions. As clear from above expression, weights are automatically chosen depending on the number of distribution.

• This is being automatically flagged as low quality, probably because it is so short. At present it is more of a comment than an answer by our standards. Can you expand on it? We can also turn it into a comment. – gung - Reinstate Monica Sep 18 '17 at 18:20
• That sounds like a clarifying comment, rather than an answer. Should this be an edit into the question? – gung - Reinstate Monica Sep 18 '17 at 18:27
• @ gung, modified my answer. Hope it helps. – Hello World Sep 18 '17 at 18:29

Scala version of JS divergence of two arbitrary length sequences:

def entropy(dist: WrappedArray[Double]) = -(dist.filter(_ != 0.0).map(i => i * Math.log(i)).sum)

val jsDivergence = (dist1: WrappedArray[Double], dist2: WrappedArray[Double]) => {
val weights = 0.5 //since we are considering inly two sequences
val left = dist1.zip(dist2).map(x => x._1 * weights + x._2 * weights)
// println(left)
// println(entropy(left))
val right = (entropy(dist1) * weights) + (entropy(dist2) * weights)
// println(right)
entropy(left) - right

}

jsDivergence(Array(0.5,0.5,0), Array(0,0.1,0.9))

res0: Double = 0.557978817900054


Cross check this answer with code in the question edit section:

jsd([np.array([0.5,0.5,0]), np.array([0,0.1,0.9])])
0.55797881790005399


A general version, for n probability distributions, in python based in Wikipedia formula and comments in this post with vector of weights (pi) as parameter and custom logbase:

import numpy as np
from scipy.stats import entropy as H

def JSD(prob_distributions, weights, logbase=2):
# left term: entropy of mixture
wprobs = weights * prob_distributions
mixture = wprobs.sum(axis=0)
entropy_of_mixture = H(mixture, base=logbase)

# right term: sum of entropies
entropies = np.array([H(P_i, base=logbase) for P_i in prob_distributions])
wentropies = weights * entropies
# wentropies = np.dot(weights, entropies)
sum_of_entropies = wentropies.sum()

divergence = entropy_of_mixture - sum_of_entropies
return(divergence)

# From the original example with three distributions:
P_1 = np.array([1/2, 1/2, 0])
P_2 = np.array([0, 1/10, 9/10])
P_3 = np.array([1/3, 1/3, 1/3])

prob_distributions = np.array([P_1, P_2, P_3])
n = len(prob_distributions)
weights = np.empty(n)
weights.fill(1/n)

print(JSD(prob_distributions, weights))


0.546621319446

I could not find any library so I created following Python code. I cross verified with scipy implementation for 2 distributions and get the same result source.

# Calculates the JSD for multiple probability distributions
def jsd(self, prob_dists):
weight = 1/len(prob_dists)              # Set weights to be uniform
js_left = np.zeros(len(prob_dists[0]))
js_right = 0
for pd in prob_dists:
js_left += pd * weight
js_right += weight * self.entropy(pd, normalize=False)

jsd = self.entropy(js_left, normalize=False)-js_right
return jsd

# Entropy function
def entropy(self, prob_dist, normalize=True):
entropy = -sum([p * math.log2(p) for p in prob_dist if p != 0])
if normalize:
max_entropy = math.log2(prob_dist.shape[0])
return entropy/max_entropy
return entropy