Using custom python function within PyMC3

Viewed 137

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):
    '''
    params:
    x: original media variable, array
    L: length
    P: peak, delay in effect
    D: decay, retain rate
    returns:
    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.append(xi)
    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

1 Answers
Related