I am wondering about the most effective way of incorporating the above improper prior for eta in a hierarchical model as a statial random effect. The Q matrix is singular and has rank (n-1) where n is the number of data points. The rho has a prior following a Gamma(0.5, 0.005). The reason why the constraint is added onto the improper prior joint distribution is so it can be identifiable. This is what i’ve tried to do.
`with pm.Model() as m: .... # other parameter definitions before this line tau = pm.Gamma('tau', alpha=0.5, beta=0.005) #eta = pm.MvNormal('eta', mu=0, tau=tau*Q, shape=(1, n)) eta = pm.DensityDist('eta', lambda value: 0.5*(n-1)*tt.log(tau)-0.5*tau* tt.dot(tt.transpose(value),tt.dot(Q, value)) ) #eta_c = pm.Deterministic('eta_c', eta-tt.mean(eta)) eta = eta - tt.mean(eta) .... #line using the generated eta vector as an input to expression for the distribution of another parameter`
AT first it seems to run and shortly after while give this output to the console:
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [eta, tau__log__, alpha, beta] Could not pickle model, sampling singlethreaded. Sequential sampling (2 chains in 1 job) NUTS: [eta, tau__log__, alpha, beta]
A then shortly after will throw this error:
ValueError: Mass matrix contains zeros on the diagonal. Some derivatives might always be zero.
From the output it seems to be generating log(tau) instead of tau, plus the eta posterior samples are not needed.
My question now is, what am I doing wrong in the implementation? I cant use the available pyMC3 distributions since the prior is improper.
BTW this formulation could be used as an extension to a model like the one on this link (Problem with hierachical occupancy model) where the eta is added as a spatial random effect on the x*beta expression modelling psi.