Hi all,
For one of my projects, I wanted to estimate the posterior distribution of some kinetic parameters using MCMC. (Find my code below)
In a nutshell, I have 3 kinetic parameters DVR, R1 and k2, for which I know the prior distribution. These 3 parameters are given to a function create_activity_curve
that performs multiplications/divisions, interpolations, and convolutions to produce a time activity curve (TAC), which is a 54x1 vector. I happen to have (several) TAC data sample for which I know exactly what DVR, R1 and k2 were used, and I use one sample as observed variable to estimate the posterior distribution of the DVR, R1 and k2. I guess I have two problems there:

My main problem is that the function I use to produce the TAC vector is kind of a black box and is coded with numpy, and therefore won’t take pytensor symbolic variables. I tried using @as_op but the results do not seem to make sense to me. Is this the correct way? I also tried using pm.draw() and feeding that to my function that produce the TACS, but realized it breaks the symbolic links.

The 2nd problem is more conceptual: In a first step I would like to consider that the produced TACs is noiseless. But I feel like I need the Normal distribution to pass the “observed variable” as argument. Should I create a potential instead that calculate the logp between the TAC and the observed variable?
Here is my code:
import pymc as pm
import pytensor as ptt
from pytensor.compile.ops import as_op
create_activity_curve = ... # numpy function that creates TACs from kinetic parameters
y_obs = ... # 54x1 vector data sample for which DVR, R1, and k2 is known
@as_op(itypes=[ptt.tensor.dscalar, ptt.tensor.dscalar, ptt.tensor.dscalar],
otypes=[ptt.tensor.dvector])
def ptt_create_activity_curve(k2,R1,k2a):
kin_params = {"k2": k2, "R_1": R1, "k2a": k2a}
sn = create_activity_curve(kin_params) # Input 3 scalars
return sn # Output 54x1 vector
iter_mcmc = 10000
burn_mcmc = 5000
MCMC_model = pm.Model()
# define the model
with MCMC_model:
dist_DVR = pm.Normal.dist(mu=1.25, sigma=0.05)
dist_k2 = pm.Normal.dist(mu=0.03, sigma=0.05)
dist_R1 = pm.Normal.dist(mu=0.51, sigma=0.05)
varDVR = pm.Truncated("var_DVR", dist_DVR, lower=0, upper=None)
vark2 = pm.Truncated("var_k2", dist_k2, lower=0, upper=None)
varR1 = pm.Truncated("var_R1", dist_R1, lower=0, upper=None)
sn = ptt_create_activity_curve(vark2, varR1, vark2 / varDVR)
dist_sn = pm.Normal.dist(mu=sn, sigma=1)
varsn = pm.Truncated("var_sn", dist_sn, lower=0, upper=None, observed=y_obs)
step = pm.Metropolis()
idata = pm.sample(draws=iter_mcmc, tune=burn_mcmc, step=step)
Thanks in advance.