Problem when using shared variable and linear mean function with GP conditional

Questions
1

Hi,
Thanks to all devs for creating this wonderful tool!

I run into a problem when trying to do prediction using GP. First, I have no issue when not using shared variable. I can run my model and sample posterior predictive using gp conditionnal for new X. Here’s my model…

with pm.Model() as gp_model_1:

#Mean Function
mean_func = pm.gp.mean.Linear(coeffs=-0.2, intercept=0.5)
#mean_func = pm.gp.mean.Constant(0.3)
#custom Mean function
#mean_func = InvLogit(a = -5, x0 = 0.5, c = 0.27, d = 0.3)

#Cov Function
cov_base = pm.gp.cov.WhiteNoise(0.15)
cov_func = pm.gp.cov.ScaledCov(1, scaling_func=logistic, args=(a,x0,c,d), cov_func=cov_base)

# Specify the GP
gp = pm.gp.Latent(mean_func = mean_func, cov_func=cov_func)
# make gp prior
f = gp.prior("f", X=df.x.values.reshape(-1,1), shape= df.x.values.shape)

sigma = pm.HalfNormal('sigma', sd =0.1)
p = pm.Deterministic("p", pm.math.invlogit(f))
y = pm.Beta("y", mu = p, sigma = sigma, observed = df.y_DH.values)
mp = pm.find_MAP()
#gp_trace1 = pm.sample(1000, tune=1000, chains=2)

#set up log-normal parameter on new X for the shared variable
zeta = np.sqrt(np.log(1 + 0.4**2))
lam = np.log(0.4) - (zeta**2)/2
n_pixel = 625

with gp_model_1:
    x_depth = stats.lognorm.rvs(scale = np.exp(lam),s = zeta, loc = 0, size = n_pixel).reshape(-1,1)
    x_tensor = shared(value = x_depth)
    print(x_tensor.get_value().shape)
    f_pred = gp.conditional(f'y_pred', x_tensor, shape = n_pixel)
    
"""
 end GP Model definition
 begin Sampling with shared variable
"""

But when I tried with shared variable so I can sample through a for loop and update new Xs, I get this error message in the GP Conditionnal about the mean function. If I use a constant mean function, no problem, but I get the same error if I use the linear mean function or my custom mean function. I was able to use my custom or the Linear mean function just fine before using shared variable… The shape of the Teano tensor is (625,1).

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-20-dd75cbc07b09> in <module>
     68     x_tensor = shared(value = x_depth)
     69     print(x_tensor.get_value().shape)
---> 70     f_pred = gp.conditional(f'y_pred', x_tensor, shape = n_pixel)
     71 
     72 """

~/anaconda3/envs/bayes/lib/python3.7/site-packages/pymc3/gp/gp.py in conditional(self, name, Xnew, given, **kwargs)
    218         """
    219         givens = self._get_given_vals(given)
--> 220         mu, cov = self._build_conditional(Xnew, *givens)
    221         shape = infer_shape(Xnew, kwargs.pop("shape", None))
    222         return pm.MvNormal(name, mu=mu, cov=cov, shape=shape, **kwargs)

~/anaconda3/envs/bayes/lib/python3.7/site-packages/pymc3/gp/gp.py in _build_conditional(self, Xnew, X, f, cov_total, mean_total)
    183         A = solve_lower(L, Kxs)
    184         v = solve_lower(L, f - mean_total(X))
--> 185         mu = self.mean_func(Xnew) + tt.dot(tt.transpose(A), v)
    186         Kss = self.cov_func(Xnew)
    187         cov = Kss - tt.dot(tt.transpose(A), A)

~/anaconda3/envs/bayes/lib/python3.7/site-packages/theano/tensor/var.py in __add__(self, other)
    129     def __add__(self, other):
    130         try:
--> 131             return theano.tensor.basic.add(self, other)
    132         # We should catch the minimum number of exception here.
    133         # Otherwise this will convert error when Theano flags

~/anaconda3/envs/bayes/lib/python3.7/site-packages/theano/gof/op.py in __call__(self, *inputs, **kwargs)
    672                 thunk.outputs = [storage_map[v] for v in node.outputs]
    673 
--> 674                 required = thunk()
    675                 assert not required  # We provided all inputs
    676 

~/anaconda3/envs/bayes/lib/python3.7/site-packages/theano/gof/op.py in rval()
    860 
    861         def rval():
--> 862             thunk()
    863             for o in node.outputs:
    864                 compute_map[o][0] = True

~/anaconda3/envs/bayes/lib/python3.7/site-packages/theano/gof/cc.py in __call__(self)
   1737                 print(self.error_storage, file=sys.stderr)
   1738                 raise
-> 1739             reraise(exc_type, exc_value, exc_trace)
   1740 
   1741 

~/anaconda3/envs/bayes/lib/python3.7/site-packages/six.py in reraise(tp, value, tb)
    701             if value.__traceback__ is not tb:
    702                 raise value.with_traceback(tb)
--> 703             raise value
    704         finally:
    705             value = None

ValueError: Input dimension mis-match. (input[0].shape[1] = 1, input[1].shape[1] = 625)

Is this a bug or I’m simply feeding wrong shape or params into the model?
Any help would be appreciated.

Versions
pymc3 - 3.9.3
theano 1.0.5
python - 3.7.3
Ubuntu - 18.04.3

2

I am not entirely sure why is it working for the mean functions other that the Linear mean function, but I think the problem is with passing the shape parameter in the gp.prior and gp.conditional method. Does it raise an exception if you don’t pass the shape parameter? If not, then can you please remove it and try again?

Other reason could be that your custom mean function doesn’t output a tensor of correct shape. Make sure that the mean function you create or choose outputs a tensor of shape (X.shape[0], ) ((625, ) in your case). Shapes other that this may lead to problems…

If none of the above comment resolve your issue, can you please share the code for custom mean functions that you are using?

3

Thanks for the reply,

  1. I’ve tried not specifying a shape parameter like you said but the gp.conditional gives an error (Cannot infer shape). The prior doesn’t seem to care but I think I’ve seen in another post you need to specify the shape if you’re using shared variable… So, I do have to pass a shape to the conditional…

  2. I’m not sure how I can specify the shape of the mean function like you ask (any hint)… Here’s the code for my custom mean function and I added the code for the linear function. I’ve tried to replicate the linear function output to do mine.

     class Linear(Mean):
 R"""
 Linear mean function for Gaussian process. (pymc3)
 Parameters
 ----------
 coeffs: variables
     Linear coefficients
 intercept: variable, array or integer
     Intercept for linear function (Defaults to zero)
 """

 def __init__(self, coeffs, intercept=0):
     Mean.__init__(self)
     self.b = intercept
     self.A = coeffs

 def __call__(self, X):
     return tt.squeeze(tt.dot(X, self.A) + self.b)

 class InvLogit(pm.gp.mean.Mean):
 """
 InvLogit function for Gaussian process. Customn Mean function (JulienM)
 """

 def __init__(self, a, x0, c, d):
     pm.gp.mean.Mean.__init__(self)
     self.a = a
     self.x0 = x0
     self.c = c
     self.d = d

 def __call__(self, X):
     return tt.squeeze(self.d * pm.math.invlogit(self.a*(X - self.x0)) + self.c)

I still don’t understand why it work without shared variable… but I need to sample ppc for a big number of parameters with the x_depth tensor … So I kind of need shared variable.

4

Could this be link to issue 3345 and the init of the mean_func shape?

5

Input data to test code
input_GP_error.csv (10.5 KB)

6

I didn’t check in great detail your code or request, please tell me if I missed something. From what I can see, the problem is on the tt.dot inside of the mean function because you are sending it a scalar and asking for it to be broadcasted. I just did the following and it runs for me.

mean_func = Linear(coeffs=np.array([-0.2]), intercept=np.array([0.5]))
7

thanks! It does work for the linear mean function! It doesn’t work for my custom invLogit function but at least it work for the linear. I didn’t want to use the constant function for my application…

I will have a deeper look into my custom function, it’s certainly my bad not the code…

thanks again!

1 Like
8

Glad it helped. If you get stuck, feel free to post the new code. Handling the shapes in the custom mean funcs of GPs is tricky.

2 Likes
9

I’m trying to do a mean function that follows an inverse logic based on this scaling function from the GP covariance documentation.

d * pm.math.invlogit(a*(x - x0)) + c

I’ve define my mean function like this:

def __init__(self, a, x0, c, d):
pm.gp.mean.Mean.__init__(self)
self.a = a
self.x0 = x0
self.c = c
self.d = d

def __call__(self, X):
return tt.squeeze(self.d * pm.math.invlogit(self.a*(X - self.x0)) + self.c)

and use the input like you said

mean_func = InvLogit(a = np.array([-5]), 
                     x0 = np.array([0.5]), 
                     c = np.array([0.27]), 
                     d = np.array([0.3]))

I still get the same error

ValueError: Input dimension mis-match. (input[0].shape[1] = 1, input[1].shape[1] = 625)

I’ve replaced the multiplication of both constant (self.a and self.d) with the pm.math.dot function but it doesn’t work…

10

I can’t understand why, but tt.squeeze seems to not be doing what it is supposed to. You can make it work by using another way to remove the second dimension.

(self.d * pm.math.invlogit(self.a*(X - self.x0)) + self.c).reshape((-1,))
2 Likes
11

Thanks so much!

1 Like