For simplicity, you can use betas = pm.beta(1, 2, size=100) instead.
With that, as @ricardoV94 mentioned, you can define the observation distribution as mnorm = pm.MvNormal('mnorm', mu=betas, cov=cov, observed=obs) given an array of observations, obs. From there you can use pm.sample.