# How to implement Kullback-Leibler divergence using Mathematica's probability and distribution functions?

Last year Jens answered a question about how to implement Kullback-Leibler divergence in the discrete case. He gave this answer:

pmfA = {1/6, 1/6, 1/6, 1/6, 1/6, 1/6};
pmfB = {1/4, 1/4, 1/2, 0, 0, 0};

Function[
{p, q},
Limit[p*Log[(p + ϵ)/(q + ϵ)], ϵ -> 0]][pmfB, pmfA]
Total[%]


However, I was curious as to if there was a way to do this with continuous distributions using Mathematica's built-in probability related functions.

P = NormalDistribution[];
Q = GumbelDistribution[];

Integrate[
Probability[x, x \[Distributed] P] Log[
Probability[x, x \[Distributed] P]/
Probability[x, x \[Distributed] Q]], {x, -Infinity, Infinity}]


But this doesn't yield anything usable...

Similarly

\[Integral]P Log[P/Q]


doesn't do much either...

I have always had a hard time using Mathematica's Probability so I would appreciate your clarification.

In the continuous case the Kullback-Leibler-Divergence from distribution $Q$ to distribution $P$ is defined as

$$D_{KL}( P ||Q) = \int \limits_{-\infty}^{+\infty} p(x)\cdot \log \frac{p(x)}{q(x)} dx$$

So what you need is the probability density function which is PDF in Mathematica:

distP = NormalDistribution[]; (* avoid capital letters *)
distQ = GumbelDistribution[];

klDivergenceContinuous = Function[ { dist1, dist2 },
NIntegrate[
PDF[ dist1, x ] × ( Log @ PDF[ dist1, x ] - Log @ PDF[ dist2, x ] ),
{ x, -∞, +\[Infinity] }
]
];

klDivergenceContinuous[ distP, distQ ]


0.229783

## Update: More general solution

The definition of the expected value in the continuous case is given as

$$\mathbb{E}[ f(X) ] = \int \limits_{-\infty}^{+\infty} f(x)\cdot p(x) dx$$

We can use this definition to find the KL-Divergence. Since Mathematica's implementation of Expectation will handle discrete distributions also (effectivly summing over all values instead of integrating over them), the following routine might turn out to be more general:

Options[ klDivergence ] = Options @ Expectation;

klDivergence[ distP_ , distQ_, opts:OptionsPattern[ klDivergence ] ] := Expectation[
Log @ PDF[ distP, \[FormalX] ] - Log @ PDF[ distQ, \[FormalX] ],
\[FormalX] \[Distributed] distP,
opts
]

klDivergence[ 1, distP_, distQ_, opts:OptionsPattern[ klDivergence ] ] := klDivergence[
distP,
distQ,
opts
]

klDivergence[ n_Integer, distP_, distQ_, opts:OptionsPattern[ klDivergence ] ]
/; n > 0 := Module[
{
vars = Table[ Unique[ \[FormalX] ], n ]
},
Expectation[
Log @ PDF[ distP, vars ] - Log @ PDF[ distQ, vars ],
vars \[Distributed] distP,
opts
]
]


This might (in principle) work for discrete as well as for continuous distributions:

Continuous Distributions (OP)

klDivergence[ distP, distQ ]


$-\frac{1}{2}+ \sqrt{\mathrm{e}} - \frac{1}{2} \log 2\pi$

N[%]


0.229783

Discrete Distributions

We may use the example given in (104506) to check how this works out in the discrete case:

pmfA = EmpiricalDistribution[ {1/6, 1/6, 1/6, 1/6, 1/6, 1/6} -> Range[6] ];
pmfB = EmpiricalDistribution[ {1/4, 1/4, 1/8, 1/8, 1/8, 1/8} -> Range[6] ];

klDivergence[ pmfA, pmfB ]


1/3 (Log[4] - 3 Log[6] + 2 Log[8])

The problems arising from a discrete distribution containing zero probabilities (cf. the question linked above) may be solved by dropping the zero-entries from the distribution:

(* pmfB = {1/4, 1/4, 1/2, 0, 0, 0} *)
pmfBzero = EmpiricalDistribution[ {1/4, 1/4, 1/2} -> Range[3] ];

klDivergence[ pmfBzero, pmfA ]


Log[ 3/Sqrt[2] ]

• thank you for this. I do have a few questions though. Stylistically, why do you use Function[{args},code] verse function[args_]:=? Also what is the internal difference between NormalDistribution[] and PDF[NormalDistribution[]]? Why wouldn't functions just work on the NormalDistribution[]? What is the difference from PDF[dist] and Probability[x, x ~ dist]? Feb 17, 2017 at 13:21
• It is clear why the discrete case works (as in my original post), but does PDF work for discrete distributions / or should one use ProbabilityDistribution? Feb 17, 2017 at 13:32
• @SumNeuron You do need to understand the difference between a (continuous) probability distribution, its density function, a probability (which makes use of the density) and an expectation. Start here maybe. All of this is concisely implemented in Mathematica.
– gwr
Feb 17, 2017 at 15:27
• Just a small suggestion: I'd localize the variable x with Module[{x},... (+1).
– Jens
Feb 17, 2017 at 17:05
• I have now expanded the solution to include the multivariate case as well. It is quite simple to extend this further by copying the definitions to have another function nklDivergence which will use NExpectation instead of Expectation.
– gwr
Mar 8, 2017 at 12:16

I know that this question already has a very good answer, but I still wanted to present the Wolfram Function Repository function I wrote specifically for this purpose (before I even found this question).

Like the solution by gwr, it uses (N)Expectation to compute the result, though I use LogLikelihood instead of Log @ PDF[...]. The function also checks the domains of the distributions for you.

edit The current version will complain if you use an EmpiricalDistribution, but I submitted a new version that fixes that problem.