Hello pymc3 community!
I research analytical methods for protein structure modeling and have recently used pymc3 to estimate uncertainty in experimentally data obtained modeling parameters. Now, I want to make the model that we are using more general and have a question (actually two) concerning a custom sampling step and was hoping maybe some of you could provide some input.
Okay, here is the problem: The whole model is pretty complex and has quite a lot of parameters, so I don’t really have a minimal example to post right now. But basically what I want to do is to have a custom sampling step for parameter (vector) P. P is very similar to a multivariate Gaussian but must be constrained to nonnegative values. We implemented this step (but not the full model) in MATLAB in previous works (for those curious, here is the link) and want to port this to python and pymc3 for obvious reasons.
Here are my two questions

I am struggling a little with telling pymc3 to do a custom step for one of the parameters while using NUTS for all the others. I found a few examples online, but those are either old and don’t work with the current pymc3 version or are not what I want to do.

What is the best approach for the custom draw for P. It is quite a bit of code that involves matrix inverses and a fast nonnegative least squares solution (fnnls), and also includes a while condition. Converting the sampler into python code was pretty straight forward, but I am having issues with making it compatible with theano, in particular the fnnls part with the while loop. Is there an easier way of doing this? This draw also depends on a hyperparameter δ that is sampled via NUTS.
On a side note: What I am doing currently is to sample P with
P = pm.MvNormal("P", mu=P0, chol = CL, shape = nr)
which runs but doesn’t give acceptable results. I also played around with restraining this
P = pm.Bound(pm.MvNormal, lower=0.0, upper=5.0)("P", mu=P0, chol = CL, shape = nr)
but this leads to issues with the gradient.
I’d appreciate any help and I’ll try to come up with a MWE (or minimal non working example ) in the meanwhile.
Thanks,
Stephan