What is the correct way to implement Jensen-Shannon Distance? I'm trying to use this code to compute the Jensen-Shannon distance:
def js_dist(P, Q):
    """Compute the Jensen-Shannon distance between two probability distributions.

    Input
    -----
    P, Q : array-like
        Probability distributions of equal length that sum to 1
    """

    def _kldiv(A, B):
        # Calculate Kullback-Leibler divergence
        
        return np.sum([v for v in A * np.log2(A/B) if not np.isnan(v)])

    P = np.array(P)
    Q = np.array(Q)

    M = 0.5 * (P + Q)

    # Get the JS DIVERGENCE
    result = 0.5 * (_kldiv(P, M) +_kldiv(Q, M))
    # Take sqrt to get the JS DISTANCE
    return np.sqrt(result)

but when comparing it to Scipy's implementation:
from scipy.spatial.distance import jensenshannon

#See here for the Scipy source code: 
https://github.com/scipy/scipy/blob/v1.5.2/scipy/spatial/distance.py#L1230-L1287 

I get different results:
a = [0.2, 0.4, 0.4]
b = [0.3, 0.2, 0.5]

[In]: js_dist(a,b)
[Out]: 0.18918

[In]: jensenshannon(a,b)
[Out]: 0.15750

Where am I going wrong??
 A: You are using logarithm base 2 by default.
scipy.spatial.distance.jensenshannon uses the default base of scipy.stats.entropy.
https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.jensenshannon.html
If you set the optional base parameter in jensenshannon(a,b, base=2.0), function will use log2 too, and you will obtain the same result as in your function 0.18918.
Your functions is well implemented.
A: If you want calculate "jensen shannon divergence", you could use following code:
from scipy.stats import entropy
from numpy.linalg import norm
import numpy as np

def JSD(P, Q):
    _P = P / norm(P, ord=1)
    _Q = Q / norm(Q, ord=1)
    _M = 0.5 * (_P + _Q)
    return 0.5 * (entropy(_P, _M) + entropy(_Q, _M))

but if you want " jensen-shanon distance", you can take square root of JSD or use scipy.spatial.distance.jensenshannon
A: Here is a minimal example based on two normal distributions (built based on the answers already exist in this thread):

import scipy.stats
import scipy.spatial
import numpy as np


# Start with two normally distributed samples
# with identical standard deviations;
# The two means are 2 standard deviations away from each other
x1 = scipy.stats.norm.rvs(loc=0, scale=1, size=1000)
x2 = scipy.stats.norm.rvs(loc=2, scale=1, size=1000)

# Construct empirical PDF with these two samples
hist1 = np.histogram(x1, bins=10)
hist1_dist = scipy.stats.rv_histogram(hist1)
hist2 = np.histogram(x2, bins=10)
hist2_dist = scipy.stats.rv_histogram(hist2)

X = np.linspace(-4, 6, 10)
Y1 = hist1_dist.pdf(X)
Y2 = hist2_dist.pdf(X)

# Obtain point-wise mean of the two PDFs Y1 and Y2, denote it as M
M = (Y1 + Y2) / 2

# Compute Kullback-Leibler divergence between Y1 and M
d1 = scipy.stats.entropy(Y1, M, base=2)
# d1 = 0.406
# Compute Kullback-Leibler divergence between Y2 and M
d2 = scipy.stats.entropy(Y2, M, base=2)
# d2 = 0.300

# Take the average of d1 and d2 
# we get the symmetric Jensen-Shanon divergence
js_dv = (d1 + d2) / 2
# js_dv = 0.353
# Jensen-Shanon distance is the square root of the JS divergence
js_distance = np.sqrt(js_dv)
# js_distance = 0.594


# Check it against scipy's calculation
js_distance_scipy = scipy.spatial.distance.jensenshannon(Y1, Y2)
# js_distance_scipy = 0.493


The difference between the KL-divergence-derived JS distance and scipy's JS distance may have been caused by the very coarse binning used to construct the empirical PDF.
