Thanks for the sharing the detail! I understand now how it works when Y is 1-dimensional. The problem I am working on having Y as multi-dimensional (think of Y as multi-variate normal) and for each row of Y, it could be missing for just some columns or the entire row. So would pymc3 use all the non-missing-at-all rows to infer \mu then sample missing Y’s or incorporate partially missing rows as well? If so, could you comment on how it works briefly under this framework? Would there be any conditional distributions estimated as well in case which I would image it’d be quite computationally pricy?