# Jensen-Shannon Divergence for multiple probability distributions?

What is the correct mathematical expression for computing Jensen-Shannon divergence ($$JSD$$) between multiple probability distributions?

I found the following expression on Wikipedia, but I did not find any official reference:

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

where $$H$$ is the Shannon-entropy.

• What is the intuition of divergence in multiple distribution case?
• In the case of two univariate distributions $$P$$ & $$Q$$, it was interpreted as a distance between them $$JSD(P||Q)$$. So in the multi-distribution case, is it measuring the distance between two multivariate distributions?
• Is there an efficient Python code available for it?

I found one paper that uses the JS-Divergence of multiple distributions to estimate the hardness of a query (in the area of Information Retrieval). The paper can be found here and they themselves refer a paper called "Divergence measures based on the shannon entropy (1991)". They also give a little bit different mathematical expression for it.

As for the interpretation, they explain it as follows:

Given a set of distributions thus obtained, we employ the well known Jensen-Shannon divergence [8] to measure the diversity of the distributions corresponding [...] So the Jensen-Shannon divergence can be seen to measure the overall diversity between all the probability distributions.

As for the Python code, I couldn't find any package that implements the JSD for more than two distributions. But there is already one quite straightforward code example on crossvalidated (see here)

.

• @ Bobipuegi Python code you referred is for 2 distributions only while I was looking for a general case. Anyway, I can write it myself as I have done for 2 distributions. Also, the paper you referred does not give exactly the same expression as given in my question, so I am not sure if it is correct (though I have seen on Stack Exchange someone using above expression). I want to make sure which expression is the correct one. More importantly, my question is about understanding intuition behind multi-distribution JS divergence. Between which distributions does it compute distance? – Hello World Sep 20 '17 at 17:15
• 1. I cannot do the math right now, but it seems to me that the formula you give and the one in the paper are equivalent after a bit of reshaping. 2. the Python code in the link I attached is able to handle multiple distributions 3. as cited in the paper, the js-divergence of multiple distributions isn't a distance measure but an overall diversity measure – Bobipuegi Sep 20 '17 at 17:41
• @ Bobipuegi So is it JS(P1 vs. P2, P3,...,Pn) or JS(P1 vs. P2 vs. P3 vs....vs., Pn) or something like this? – Hello World Sep 20 '17 at 18:00
• @ Bobipuegi For 2 distributions, it is JS(P vs. Q) – Hello World Sep 20 '17 at 18:07
• I understand it as first calculating a "centroid" distribution and then calculating something like JS(P1vs. Centroid, P2 vs. Centroid .... Pn vs. Centroid) – Bobipuegi Sep 20 '17 at 18:13

The Wikipedia has the formula on this, and the references about Information Radius (I think it is the first 2) should give you the formula. Another reference is work by Bob Williamson at ANU (Australian National University). You should be able to find him on Scholar.google.com. Good luck and thanks for asking the question.

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