```julia julia> using Symbolics julia> @variables μp σp μq σq x; julia> logdensity_rel(Normal(μp,σp), Normal(μq,σq), x) |> simplify (1//2)*(((x - μq) / σq)^2) - (1//2)*(((x - μp) / σp)^2) ``` Note that it simplifies more than this (Factor out the 1/2, then it's a difference of squares). We can make this efficient by adding `logdensity_def` methods