# How to calculate Jensen-Shannon divergence?

How to write the function (or any other methods) to calculate Jensen-Shannon divergence (JSD) for two (p and q) discrete probability distributions? I need to calculate JSD for too many different (over 1000) probability distributions so it should be run able by for loop. Thank you for your help. I wrote the below code but it is not working. You can see the part of the probability distribution (PD) of p and q. There are a lot 0 in the both p and q PD.

p={0.13253, 0.0361446, 0.0120482, 0.0240964, 0., 0., 0., 0., 0., 0.,0., 0., 0.,0,0}

q={0.0282448, 0.0163522, 0.010406, 0.0163522, 0.0208119, 0.0118925,
0.00891941, 0.00297314, 0.0044597, 0.00148657, 0., 0., 0., 0.,0.00148657}

Function[{p, q},

N[1/2 (p*Log2 p + q*Log2 q) - ((p + q)/2*Log2 (p + q)/2)]][

p@DjointProbaility, q@DmarjinalProbaility]

• I think you'll get more help with your questions if you would show some programming effort and/or definitions of the statistics of interest.
– JimB
Feb 9, 2022 at 0:58
• I tried many different functions and method before asking the question. Feb 9, 2022 at 1:04
• That's much better but without examples of the p and q, how can we know what "it is not working" means.
– JimB
Feb 9, 2022 at 1:06
• For the user's convenience, see Wiki for info. Feb 9, 2022 at 6:53
• Jensen-Shannon divergence is only well-defined for pairs of probability distributions that are absolutely continuous with respect to each other. Meaning: both p/q and q/p must exists (without division by zero). Your example vectors just don't satisfy this, so it takes no wonder that this does not "work". This metric is just not applicable to your data. Feb 9, 2022 at 8:58

The function cJSdiv below is Listable, so when called on two vectors p and q goes through the entries and returns a vector of same length. It tests whether both are zero or both are nonzero. That's the two cases in which the summand for the Jensen-Shannon divergence is defined; then cJSdiv returns this summand at the corresponding position of the output vector. Otherwise the maximum double precision number 1.7976931348623157*^308 is returned. The function JSdiv serves as a wrapper and performs the summation.

cJSdiv = Compile[{{p, _Real}, {q, _Real}},
Block[{minv},
If[p > 0. && q > 0.,
minv = 2./(p + q);
0.5 (p (Log[p minv]) + q (Log[q minv])),
If[p == 0. && q == 0.,
0.,
1.7976931348623157*^308
]
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable}
];

JSdiv = {p, q} |-> Total[cJSdiv[p, q]];


Just call it like this:

JSdiv[p,q]

• Executing your code, I obtain "CCompilerDriverCreateLibrary::nocomp: A C compiler cannot be found on your system. Please consult the documentation to learn how to set up suitable compilers." and "Compile::nogen: A library could not be generated from the compiled function.". Feb 9, 2022 at 10:38
• Needs["SymbolicC"] and Needs["CCompilerDriver"] do not help. Feb 9, 2022 at 11:06
• Well, then you have to install a C/C++ compiler. Natural choices would be the following: Apple clang for macos; gcc or clang for Linux; Microsoft Virtual Studio for Windows. See here for instructions: reference.wolfram.com/language/CCompilerDriver/tutorial/… Feb 9, 2022 at 11:09
• Alternatively, you can change CompilationTarget -> "C", to CompilationTarget -> "WVM",. Then no external compiler is required, but the performance is not as nearly good. Feb 9, 2022 at 11:12
• Thank you. I used CompilationTarget -> "WVM" and obtained 1.258385194403621*10^309.BTW, this differs from 1.7976931348623157*^308`. Feb 9, 2022 at 13:01