Hello everyone

I have been trying to find a way to do convolution in pymc3 as part of a regression problem. My situation is this, i have a set of data and a function and i want to fit the convolution of those two to an other set of data to obtain the parameters for the function. I have been trying to use numpy.convolve, but it does not seem to work with the functions defined in pymc3 (theano), as it gives a ValueError because i am setting an array element with a sequence. Is there a way to do this type of convolution in theano?

Here is my attempt so far using numpy.convolve

```
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import theano
import theano.tensor as tt
from theano.tensor import fft
muJ=np.asarray([-0.03425363, -0.03351397, -0.02664699, -0.02273714, -0.03174768,
-0.03957353, -0.04738784, -0.05586893, -0.06363427, -0.07391645,
-0.08178534, -0.08926872, -0.09593623, -0.09036217, -0.07807825,
-0.06753805, -0.05767188, -0.06243789, -0.06809364, -0.06405182,
-0.06114183, -0.04508844, -0.03290343, -0.02404407, -0.02079736,
-0.03644591, -0.03598996, -0.03445805, -0.03331908, -0.01566788,
-0.01656279, -0.01754197, -0.0130746 , 0.00193958, 0.01056539,
0.01698442, 0.0382919 , 0.05958532, 0.06605687, 0.07316159,
0.07929714, 0.08111757, 0.08236758, 0.08409051, 0.07132862,
0.06410184, 0.06383741, 0.0638774 , 0.06268115, 0.0642296 ,
0.07699274, 0.09465466, 0.10719381, 0.07496593, 0.04823132,
0.02464346, 0.00900288, -0.00262479, -0.01581547, -0.02713472,
-0.04211868, -0.06695968, -0.0843794 , -0.0838455 , -0.09070336,
-0.09610623, -0.09881969, -0.09027169, -0.08033855, -0.06809048,
-0.04661056, -0.01563495, 0.0156263 , 0.03325961, 0.04721029,
0.0665903 , 0.09242785, 0.10892494, 0.12218285, 0.1282392 ,
0.13511139, 0.14373963, 0.15130526, 0.15932193, 0.1431794 ,
0.08651746, 0.06932233, 0.05781805, 0.05321184, 0.04833902,
0.04359331, 0.03972242, 0.03617235, 0.0324834 , 0.04958008,
0.05654622, 0.05459009, 0.07174552, 0.09021518, 0.10966458])
XK_new=np.linspace(57690.160596,58541.349392,100,dtype=np.float64)
muK = np.asarray([-0.0188423 , -0.01101513, -0.02182555, -0.04066994, -0.0372273 ,
-0.04142979, -0.04643643, -0.05618733, -0.06305947, -0.07055236,
-0.07793845, -0.08585502, -0.06063842, -0.04144011, -0.04482767,
-0.05074738, -0.05828078, -0.06202257, -0.0469031 , -0.03229114,
-0.01682147, -0.00109132, 0.01100408, 0.02217135, 0.00313528,
-0.01391616, -0.02103778, -0.00784842, -0.00958049, -0.00958786,
0.01109297, -0.00763489, -0.02091672, -0.01876265, 0.00469278,
0.0188317 , 0.03078942, 0.04751188, 0.06810743, 0.07498746,
0.08292102, 0.08732753, 0.10046523, 0.11947633, 0.13786948,
0.14400703, 0.13987333, 0.13111406, 0.12287179, 0.11626273,
0.11051379, 0.09606653, 0.03146849, -0.03145527, -0.03731231,
0.05864158, 0.10667121, 0.12383762, 0.11461967, 0.08906796,
0.06559065, 0.04109558, 0.02030753, -0.00256419, -0.02360279,
-0.04425853, -0.07565115, -0.11281323, -0.158842 , -0.15194773,
-0.13758974, -0.12698046, -0.11453399, -0.09043922, -0.06368814,
-0.03801246, -0.01361853, 0.01826796, 0.04772783, 0.0783482 ,
0.10210635, 0.11343715, 0.11997221, 0.12888862, 0.14047052,
0.1510065 , 0.16296829, 0.15034888, 0.12048199, 0.10888384,
0.10030352, 0.09039116, 0.08457157, 0.07565612, 0.06965578,
0.06232899, 0.05684617, 0.07264807, 0.09185271, 0.10704618])
with pm.Model() as convmodel:
# Define priors
sigma_DT=pm.Uniform('sigma_DT', lower=-7.0, upper=2.0)
sigma_AD=pm.Uniform('sigma_AD', lower=-7.0, upper=2.0)
mu_DT=pm.Uniform('mu_DT', lower=2.7, upper=4.0)
mu_AD=pm.Uniform('mu_AD', lower=0.0, upper=2.3)
A_T=pm.Uniform('A_T', lower=0.0, upper=1.0)
T=pm.Uniform('T', lower=1300.0, upper=1700.0)
wav=pm.Uniform('wav', lower=0.0, upper=10000.0)
wav_0=pm.Uniform('wav_0', lower=0.0, upper=10000.0)
index=pm.Uniform('index', lower=0.0, upper=10.0)
sigma=pm.HalfNormal("sigma", sigma=np.sqrt(np.pi/(np.pi-2)))
h = 6.626e-34#Plancks constant
c = 299792458#speed of light
k = 1.38e-23#Boltzmanns constant
a = 2.0*h*c**2
b = h*c/(wav*k*T)
BB = a/ ( (wav**5) * (np.exp(b) - 1.0) )
exp_DT = -((tt.log(XK_new)-mu_DT)**2/(2*sigma_DT**2))
front_DT = A_T/(XK_new*sigma_DT*np.sqrt(2*np.pi))
Psi_DT = BB*front_DT*np.exp(exp_DT)
powr = (wav/wav_0)**(index)
exp_AD = -((tt.log(XK_new)-mu_AD)**2/(2*sigma_AD**2))
front_AD = (1.0-A_T)/(XK_new*sigma_AD*np.sqrt(2*np.pi))
Psi_AD = powr*front_AD*tt.exp(exp_AD)
transfer = Psi_DT + Psi_AD
conv= np.convolve(transfer,muJ,'same')
likelihood = pm.Normal('y', mu=conv, sigma=sigma, observed=muK)
tracetransfer = pm.sample(2000)
```