Help with Multivariate Normal with nested prior information

Currently, in your model there is no free parameter, in another word, there is nothing to be estimated. You are basically sampling from a MvNormal(mu0, V0).

To represent prior information in the model, you need to represent them as distribution, and let them interact with each other, for example, you can do:

with pm.Model() as Maher_model:
    theta = pm.MvNormal('theta', mu=mu0, cov=V0, shape=6)
    H0 = pm.Uniform('H0', 0, 1, shape=(4, 6))
    g = tt.dot(H0, theta)
    likelihood = pm.MvNormal('likelihood', mu=g, cov=sigma, shape=4, observed=g)

where the prior information about g is represented using theta and H0.