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]])

 A: 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
}

A: 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
A: 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

A: 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.
A: 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: 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

