Also, instead of doing the Bernoulli layer, have you thought about just doing a simple mean value imputation before doing factor analysis, similar to what was done here:
https://pymc-devs.github.io/pymc3/notebooks/pmf-pymc.html
# Perform mean value imputation
nan_mask = np.isnan(self.data)
self.data[nan_mask] = self.data[~nan_mask].mean()