Convolution in pymc3 with numpy and/or theano

Questions
1

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)
2

For those curious, i have been trying to define a theano function that does 1D convolution like in the last example here Link but i am running into problems with theano. Firstly conv2d = T.signal.conv.conv2d

no longer seems to be supported only

 from theano.tensor.nnet import conv
conv2d = conv.conv2d

seems to work and also the x.shape used to define the vector length gives an error related to theano test values

ValueError: Cannot compute test value: input 0 (<TensorType(float64, matrix)>) of Op Shape(<TensorType(float64, matrix)>) missing default value.
3

You should use theano convolution for sure, otherwise the performance would be pretty bad.
How are you calling the conv2d in this case?

4

I first define conv2d as

conv2d = theano.tensor.nnet.conv2d

and then i define the two input matrices together with the length of the 1-d vector for the inputshape

x = tt.matrix()
y = tt.matrix()
veclen = x.shape[1]

I can then define the expression for the 1-d convolution

conv1d_expr = conv2d(x, y, input_shape=(1, veclen), border_mode='full')
conv1d = theano.function([x, y], outputs=conv1d_expr)

I fixed the ValueError related to the lack of test values, but now i seem to have a problem with how veclen is defined since i get the error

ValueError: imshp should be None or a tuple of constant int values

I seem to be misunderstanding how to define the image shape for the 1-d case.

I also found out that conv2d does not have the “boder_mode=‘same’” like np.convolve has from here and here, which is what i am looking for.

5

You dont need to define a theano.function, just pass the conv1d_expr around (this is the symbolic tensor represents the result after convolution).

6

ahh i see what you mean i can use the expression inside the model instead of defining a function and then using it inside the model. So i can write the following inside the model

conv2d = theano.tensor.nnet.conv2d

conv1d_expr = conv2d(transfer,muJ, input_shape=(len(muJ),1), border_mode='full')

likelihood = pm.Normal('y', mu=conv1d_expr, sigma=sdK, observed=muK)

tracetransfer = pm.sample(2000)

the expression seems to work now i just need to pass the function and data in correctly. I am getting an error right now because cov2d wants a 4D tensor as input

TypeError                                 Traceback (most recent call last)
<ipython-input-14-7f47da93c2db> in <module>
 43     conv2d = theano.tensor.nnet.conv2d
 44 
---> 45     conv1d_expr = conv2d(transfer,muJ, input_shape=(len(muJ),1), border_mode='full')
 46 
 47     #conv= np.convolve(transfer,muJ,'same')

~\Anaconda3\lib\site-packages\theano\tensor\nnet\__init__.py in conv2d(input, filters, input_shape, filter_shape, border_mode, subsample, filter_flip, image_shape, filter_dilation, num_groups, unshared, **kwargs)
167     return abstract_conv2d(input, filters, input_shape, filter_shape,
168                            border_mode, subsample, filter_flip,
--> 169                            filter_dilation, num_groups, unshared)
170 
171 

~\Anaconda3\lib\site-packages\theano\tensor\nnet\abstract_conv.py in conv2d(input, filters, input_shape, filter_shape, border_mode, subsample, filter_flip, filter_dilation, num_groups, unshared)
642                              num_groups=num_groups,
643                              unshared=unshared)
--> 644     return conv_op(input, filters)
645 
646 

~\Anaconda3\lib\site-packages\theano\gof\op.py in __call__(self, *inputs, **kwargs)
613         """
614         return_list = kwargs.pop('return_list', False)
--> 615         node = self.make_node(*inputs, **kwargs)
616 
617         if config.compute_test_value != 'off':

~\Anaconda3\lib\site-packages\theano\tensor\nnet\abstract_conv.py in make_node(self, img, kern)
2333 
2334         if img.type.ndim != 2 + self.convdim:
-> 2335             raise TypeError('img must be %dD tensor' % (2 + self.convdim))
2336 
2337         if self.unshared:

TypeError: img must be 4D tensor

My guess is i will have to format my input in such a way that the input is a 4D tensor with two dimensions of length 1 and then match that for the input_shape.

7

Exactly :wink:

8

Okay so after struggling with dimensions for a bit i created a simple example to test whether the convolutions works

import theano.tensor.signal.conv
from scipy import signal
xtest=np.linspace(0,10,100)
filt=2.0*xtest
conv=signal.convolve(xtest, filt, mode='full')

with pm.Model() as model:

A =pm.Uniform('A',lower=0.0,upper=4.0)#pm.Normal('A', mu=2.0, sigma=1.0)
filtt = A*xtest
convol=theano.tensor.signal.conv.conv2d(xtest[None,:],filtt[None,:],(1,100),(1,100),border_mode='full')

likelihood = pm.Normal('conv', mu=convol[0], observed=conv)

trace = pm.sample(1000)

I tried both using a uniform and normal prior to see how robust the solution was. While the convolution seems to work i am suprised that the sampler seems to have a high acceptance probability and a small number of effective samples

The acceptance probability does not match the target. It is 0.9979773108121525, 
but should be close to 0.8. Try to increase the number of tuning steps.
The estimated number of effective samples is smaller than 200 for some parameters.

Though the statictics seem to be fine for the constant

|mean|sd|mc_error|hpd_2.5|hpd_97.5|n_eff|Rhat|
|A|2.000464|0.004782|0.000403|1.999902|2.000096|95.187888|1.012574|

convtest
convtest872×152 8.45 KB

Now i just need a way to use the ‘same’ border mode that the scipy and numpy convolve functions have. My idea is to obtain the middle of the convolution output and obtain an array of the same length as the input. The trouble there will be to write it in a way theano understands.

1 Like
9

Okay i finally managed to do a 1d convolution with the right border mode. For those having problems in the future here is a simple example to show it works

import theano
import theano.tensor as tt
from scipy import signal

xtest=np.linspace(0.1,10,101)
filt=2.0*xtest
conv=signal.convolve(xtest, filt, mode='same')

with pm.Model() as model:

 A = pm.Normal('A', mu=2.0, sigma=1.0)#also works for 
 pm.Uniform('A',lower=0.0,upper=5.0)

 xtestT=np.reshape(xtest,(1,1,101,1))#np.reshape(xtest,(1,1,1,101)) should also work

 filtt = A*xtestT

 convol=theano.tensor.nnet.conv2d(xtestT,filtt,border_mode='half')
 #'half' work the same way as 'same' for odd sized filters. 

 likelihood = pm.Normal('conv', mu=convol[0,0,:,0], observed=conv)

 trace = pm.sample(1000)

The important thing to remember is how the data is formatted nnet.conv2d wants a 4d tensor, so we have to reshape the input and filter (since the filter is calculated from the input data in this case i only have to do it for the input). The input has the shape (batch size, input channels, input rows, input columns) and the important ones in our case is the last two. The sampling result are similar to the ones above for signal.conv.conv2d but we can now use the ‘half’ border mode! Just remember your filter has to be oddly sized for it to work like the ‘same’ border mode from numpy and scipy.

1 Like