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