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:
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
