How to interpolate values within pymc Model?

Questions v5
1

I’m developing a model using PyMC to estimate the relative efficiency of a detector. This efficiency is dependent on the angular position of the radiation source. To do this, I need to compute some intermediate values using sample values from Uniform distribution about the source’s location. Could you please provide guidance on how to accomplish this within the PyMC framework? Unfortunately, I’m not familiar with all features of pymc. My current code is shown below

import pytensor.tensor as tt
import numpy as np
import pymc as pm


def calculate_angles_tensor(
        x_position: tt.TensorVariable,
        y_position: tt.TensorVariable,
        src_x: tt.TensorVariable,
        src_y: tt.TensorVariable
) -> tt.TensorVariable:
    """
    Calculates relative angles of incidence using  tensors.

    :param x_position: Pytensor tensor of x coordinates
    :param y_position: Pytensor tensor of y coordinates
    :param src_x: Pytensor tensor x position of the orphan source
    :param src_y: Pytensor tensor y position of the orphan source
    :return: Pytensor tensor of relative angles expressed in degrees
    """
    # This function has been tested separately with tt 

def mean_angular_count_rate(
        src_x,
        src_y,
        activity,
        bkg,
        rel_eff
):
    dist = tt.sqrt((X_POS - src_x) ** 2 + (Y_POS - src_y) ** 2)
    angular_count_rate = (activity * SCALE * BRANCH_RATIO * DET_EFF * rel_eff * tt.exp(-mu_air * dist)) / (
            4 * np.pi * dist ** 2)
    return angular_count_rate + bkg

data_x = at.as_tensor_variable(X_POS)
data_y = at.as_tensor_variable(Y_POS)

#Example of interpolator 
x = np.linspace(-180,180, num = 100)
y = np.sin(np.deg2rad(x))
interpolator = interp1d(x,y)


with pm.Model() as cps_angl_model:
    x = pm.Uniform("src_x", lower=-1000, upper=1000)
    y = pm.Uniform("src_y", lower=-1000, upper=1000)
    angles = calculate_angles_tensor(data_x, data_y, x, y)  
    eff_rel = interpolator(angles) # <---- HERE I GET AN ERROR 
    
    act = pm.Gamma("act", alpha=1.2, beta=0.001)
    bkg_ = pm.Normal("bkg", mu=MEAN_BKG_CPS, sigma=3)

    count_rate = mean_angular_count_rate(x, y, act, bkg_, eff_rel)
    cps_obs = pm.Poisson("cps_obs", mu=count_rate, observed=cps)

Are there any ways to implement such intermediate interpolation? How can TensorVariable be interpolated to a vector of efficiency values?

2

If you want to use interpolation, your best best bet might be bspline_interpolation from pymc_extras. There will be some subtle differences between how that works and how interp1d works so be sure to carefully read the documentation.

If you want something more Bayesian, you can use Bayesian BSplines. See this example for details: Splines — PyMC example gallery.

2 Likes
3

Hi

Thanks for the tip. I’m trying to import bspline_interpolation from pymc_extras.utils.spline and getting this error “cannot import name ‘MeasurableOp’ from ‘pymc.logprob.abstract’” Couldn’t find anything useful about this issue. Advice?

4

The only advice I could give is to make sure PyMC and pymc_extras are up to date. If that’s not the issue, then I’d suspect there’s a bug. I can’t be much help there, sorry!

1 Like
5

Simple interpolation a la interp1d is sorely missing in pytensor. I opened a PR today to try to address this. If you’re antsy to try it, you can copy/paste some of the code – it’s quite straight forward to implement.

1 Like
6

Hi

Thank you for the link! This worked for me.

7

Glad this helped. I just merged the PR, so the next release of pytensor will have this functionality.

@dreycenfoiles The Bsplines from pymc-extras are also a great solution, but be aware that you can’t backprop through them. So they are good for generating data that doesn’t depend on parameters, but not if you want something that depends on data + parameters (which I thought was the case here).

Based on the posted error, they might be out of date. If one of you could confirm that they don’t work and open an issue about it on pymc-extras I’d be very grateful.