Unknown joint distribution, known conditional marginals

Questions
2

This looks like a job for DensityDist.

def myfun_f(y, x, theta):
    ...  # p(y1 | x, theta)

def myfun_g(y, x):
    ... # p(y2 | x)

with pm.Model() as model:
    X = pm.Normal('X', 0., 1.)  # or whatever prior
    theta = pm.HalfNormal('theta', 2.)  # or whatever prior
    def logp_f(y_):
        return log(myfun_f(y_, X, theta))
    def logp_g(y_):
        return log(myfun_g(y_, X))
    Y1 = pm.DensityDist('Y1', logp_f, observed={'y_': y1})
    Y2 = pm.DensityDist('Y2', logp_g, observed={'y_': y2})
2 Likes