Use arbitrary (numpy) functions on pytensor for posterior estimation with MCMC

Questions v5
1

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:

  1. 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.

  2. 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.