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

12

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]])
kanzen_master
źródło
2
Przy okazji ładny kod Pythona!
gui11aume

Odpowiedzi:

13

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

Podam szczegóły jednego obliczenia:

H(1/2,1/2,0)=1/2log(1/2)1/2log(1/2)+0=0.6931472

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

gui11aume
źródło
3
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)}pqp <- 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)0.378889334/1551/9213/45714/453737/90
1
Nie takie brudne ... ;-)
gui11aume
4
(1) Ponów matematykę. (2) Entropię można mierzyć przy użyciu dowolnej podstawy logarytmu, o ile ci się podoba, o ile jesteś spójny. Naturalne, pospolite i logi base-2 są konwencjonalne. (3) To naprawdę średnia rozbieżność między rozkładami a ich średnią. Jeśli myślisz o każdym rozkładzie jako punkcie, tworzą one chmurę. Patrzysz na średnią „odległość” między środkiem chmury a jej punktami, trochę jak średni promień. Intuicyjnie mierzy rozmiar chmury.
whuber
1
@ Legend Myślę, że masz rację. Nie przetestowałem wystarczająco po stwierdzeniu, że jeden wynik zgadza się z odpowiedzią uzyskaną w inny sposób (z Mathematica ).
whuber
1
@dmck W moim komentarzu rzeczywiście są literówki: (1) fraza h <- function(x) {została wklejona dwukrotnie. Wystarczy usunąć: wszystko inne działa i daje wyniki, które zacytowałem. Następnie zmodyfikuj apply(p, 2, h)do, apply(p, 1, h)jak wskazano w komentarzu Legend .
whuber
6

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
}
Renaud
źródło
0

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

JSmetric(p1,...,pm)=H(p1+...+pmm)j=1mH(pj)m

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.

Witaj świecie
źródło
Jest to automatycznie oznaczane jako niskiej jakości, prawdopodobnie dlatego, że jest tak krótkie. Obecnie jest to raczej komentarz niż odpowiedź według naszych standardów. Czy możesz to rozwinąć? Możemy również zamienić to w komentarz.
gung - Przywróć Monikę
To brzmi jak komentarz wyjaśniający, a nie odpowiedź. Czy powinna to być edycja pytania?
gung - Przywróć Monikę
@ gung, zmodyfikowałem moją odpowiedź. Mam nadzieję, że to pomoże.
Hello World,
0

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
Mageswaran
źródło
0

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

alemol
źródło