2
from sympy import *
x, y, mu, sigma, density1, density2 = symbols('x y mu sigma density1 density2')
eq1 = Eq(density1, 1/(sqrt(2*pi)*sigma)
                         *exp(-(x-mu)**2/(2*sigma**2)))     # normal 
eq2 = Eq(y, exp(x))                                         # substitution
eq3 = Eq(density2, 1/(y*sqrt(2*pi)*sigma)
                         *exp(-(ln(y)-mu)**2/(2*sigma**2))) # lognormal
[eq1, eq2, eq3]

Output: LaTeX output

How can I make SymPy take the normal density (eq1), apply the x to y substitution (eq2) and output the lognormal density (eq3)?

(I received no answer to this question at https://stats.stackexchange.com/q/55353/14202 .)

winerd
  • 1,813
  • 1
  • 14
  • 12

1 Answers1

1

When we change a variable in a probability density function, it is also necessary to multiply the density by the derivative of the function that performs the substitution. This is how substitution works in integrals, and probability density must have integral 1, so we have to respect that.

Let's call the substitution-performing function f (it's log(y) in your example). Here is how the process works, starting from your setup: 

f = solve(eq2, x)[0]
new_density = eq1.rhs.subs(x, f) * f.diff()
# test it now
Eq(new_density, eq3.rhs)  #  True