Using custom python function within PyMC3

I want to do bayesian regression in Python using PyMC3.
The regression is in the following form:

y = intercept + beta*function(data, parameters) + error

Where the function embeds parameters (L, P, D) that are random variables I want to estimate and for which I initialize a prior distribution. Let say:

L~Uniform(0, 10)
P~Uniform(0, 10)
D~Beta(3, 3)

The function in a mathematical form is the following:

enter image description here

That can be translate in python as:

def apply_adstock(x, L, P, D):
    x: original media variable, array
    L: length
    P: peak, delay in effect
    D: decay, retain rate
    array, adstocked media variable
    x = np.append(np.zeros(L-1), x)

    weights = np.zeros(L)
    for l in range(L):
        weight = D**((l-P)**2)
        weights[L-1-l] = weight

    adstocked_x = []
    for i in range(L-1, len(x)):
        x_array = x[i-L+1:i+1]
        xi = sum(x_array * weights)/sum(weights)
    adstocked_x = np.array(adstocked_x)
    return adstocked_x

Alternatively, this function can also be rewritten in a much more synthetic form:

def apply_adstock(x, L, P, D):
    return np.convolve(x, D**((np.arange(0, L, 1) - P)**2))[:-(L-1)] / sum(D**((np.arange(0, L, 1) - P)**2))

The problem that I am encountering, is that I want to estimate L, P, D that are random variables and that enter in a Python function using bayesian Inference. Is there a way to do so?

I have written the following code:

with Model() as model:  
    # Define priors
    sigma = HalfCauchy("sigma", beta=10, testval=1.0)
    intercept = Normal("Intercept", 0, sigma=20)
    beta = Normal("x", 0, sigma=20)
    L = pm.Uniform('L', lower=0, upper=10)
    P = pm.Uniform('P', lower=0, upper=10)
    D = pm.Beta('D', 3, 3)

    # Define likelihood
    likelihood = Normal("y", mu=intercept + beta * apply_adstock(x, L, P, D), sigma=sigma, observed=y)

    # Inference!
    # draw 3000 posterior samples using NUTS sampling
    trace = sample(300, return_inferencedata=True)

But I get the following error:

ValueError: setting an array element with a sequence.

I have searched for a solution online but I have no clue on how to estimate the posterior ditribution of the parameters within the PyMC3 model.

Thank you in advance

Hey! It should be relatively easy to port your adstock function from numpy to aesara. See for example Media Effect Estimation with PyMC: Adstock, Saturation & Diminishing Returns - Dr. Juan Camilo Orduz :slightly_smiling_face:


Thank you! This is a great source for addressing my problem :slight_smile: