How to calculate Jensen-Shannon divergence?

Question

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]

Answer 1

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]