Why won't SymPy integrate a standard Log-Normal PDF to 1?
I'm running the following code in Python 3.x and SymPy 1.0.1:
from sympy.stats import density, LogNormal
from sympy import Symbol, integrate, oo
mu, sigma = 0, 1
z = Symbol('z')
X = LogNormal('x', mu, sigma)
f = density(X)(z)
integrate(f, (z, 0, oo))
which should(?) return 1 but outputs:
sqrt(2)*Integral(exp(-log(z)**2/2)/z, (z, 0, oo))/(2*sqrt(pi))
Does anyone know what's going on here?
Apparently, Sympy fails to find the closed form solution of this integral.
You can, however, help Sympy perform the integration. One approach is to perform a transformation of the integration variable with the hope that it will result in a simpler integrand expression that Sympy can handle. Sympy offers a convenient transform() method for this purpose.
import sympy as sp
import sympy.stats
mu, sigma = 0, 1
z = sp.Symbol('z', nonnegative=True)
X = sympy.stats.LogNormal('x', mu, sigma)
f = sympy.stats.density(X)(z)
I = sp.Integral(f, (z, 0, sp.oo))
print(I)
This is the original integral form, which Sympy fails to evaluate. (Note the use of sympy.Integral which returns an unevaluated integral.) One (obvious?) transformation of the integration variable is z -> exp(z), which results in a new integral as follows
I2 = I.transform(z,sp.exp(z))
print(I2)
Now, we may call the doit() method to evaluate the transformed integral:
I2.doit()
1
