Kullback Leibler divergence between two Gaussians

kullback_leibler.py

import sympy as sp 
x = sp.symbols('x', real=True)
p, q = sp.symbols('p q', positive=True)
KL = sp.Integral(p*sp.log(p) - p*sp.log(q), (x, -sp.oo, sp.oo))
mu, mu1, mu2 = sp.symbols('mu mu1 mu2', real=True)
sig, sig1, sig2 = sp.symbols('sig sig1 sig2', positive=True)
n = 1/sp.sqrt(2*sp.pi*sig**2)*sp.exp(-(x - mu)**2/(2*sig**2))
KL_n = sp.simplify(
    KL.subs({p: n.subs({mu: mu1, sig: sig1}), q: n.subs({mu: mu2, sig: sig2})}).doit()    
)
KL_n = sp.collect(sp.logcombine(sp.expand(KL_n)), sig2)
KL_n

# Out[]: log(sig2/sig1) - 1/2 + (mu1**2/2 - mu1*mu2 + mu2**2/2 + sig1**2/2)/sig2**2