My mentors were not very convinced about creating a JointDistributionHandmade object out of each compound distributions, and marginalising over the latent distribution. Finally, after some discussion the Gitter channel, and on #14989, the final API that was decided gives handles compound distribution to give results as in the follows:

N1 = Normal('N1', 0, 1)
N2 = Normal('N2', N1, 2)
assert density(N2)(0).doit() == sqrt(10)/(10*sqrt(pi))
assert simplify(density(N2, Eq(N1, 1))(x)) == \
    sqrt(2)*exp(-(x - 1)**2/8)/(4*sqrt(pi))

This result doesn’t break the previous API, and gives the correct results in terms of density. This is done by creating a CompoundDistribution object instead of SingleDistribution objects, and calculating the PDFs using MarginalDistribution. There are checks in place that see to it that the CompoundDistribution objects indeed have a RV as an argument, and return single distribution objects otherwise. Though compound distributions took a longer time than expected, the PR was merged into the master branch, and compound distributions are implemented in the development version.