Obliczanie dywergencji Jensena-Shannona dla 3 rozkładów probów: Czy to w porządku?

Chciałbym obliczyć dywergencję Jensen-Shannon dla 3 następujących rozkładów. Czy poniższe obliczenia są prawidłowe? (Postępowałem zgodnie ze wzorem JSD z wikipedii ):

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

Z góry dziękuję...

EDYCJA Oto prosty, brudny kod Pythona, który również to oblicza:

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

$(5/18, 28/90, 37/90)$$(1/6, 1/5, 9/30)$

Podam szczegóły jednego obliczenia:

W podobny sposób pozostałe terminy to 0,325083 i 1,098612. Zatem końcowy wynik to 1,084503 - (0,6931472 + 0,325083 + 1,098612) / 3 = 0,37889

Pyton:

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

Jawa:

/**
 * 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
}

Podałeś odniesienie do Wikipedii. Podaję tutaj pełne wyrażenie rozbieżności Jensen-Shannon z wieloma rozkładami prawdopodobieństwa:

Pierwotne pytanie zostało zadane bez matematycznego wyrażenia rozbieżności JS w wielu dystrybucjach, co prowadzi do nieporozumień w zakresie rozumienia dostarczonych obliczeń. weightUżyto również terminu, który ponownie powoduje zamieszanie, że sposób wyboru odpowiednich wag do mnożenia. Powyższe wyrażenie wyjaśnia te zamieszania. Jak wynika z powyższego wyrażenia, wagi są wybierane automatycznie w zależności od liczby dystrybucji.

Wersja Scala rozbieżności JS dwóch dowolnych sekwencji długości:

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

Sprawdź tę odpowiedź z kodem w sekcji edycji pytania:

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

Wersja ogólna, dla n rozkładów prawdopodobieństwa, w pythonie opartym na formule Wikipedii i komentarzach w tym poście z wektorem wag ( pi ) jako parametrem i niestandardową bazą danych :

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