Posterior for custom likelihoods

Questions
1

Hello, I am trying to implement custom likelihoods using GP for the posterior inference for 13 parameters with normal priors. And I have been following this github(Using a "black box" likelihood function in PyMC3 | Samplers Demo) using theano.tensor.Op.

Instantiate a Gaussian Process model (using sklearn)

kernel = C(2, (1e-3, 1e3)) * RBF((10,10,10,10,10,10,10,10,10,10,10,10,10), (1e-2, 1e3))+ WhiteKernel()
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)

class LogLike2(tt.Op):

"""
Specify what type of object will be passed and returned to the Op when it is
called. In our case we will be passing it a vector of values (the parameters
that define our model) and returning a single "scalar" value (the
log-likelihood)
"""
itypes = [tt.dvector] # expects a vector of parameter values when called
otypes = [tt.dscalar] # outputs a single scalar value (the log likelihood)

def __init__(self,x,y):
    """
    Initialise the Op with various things that our log-likelihood function
    requires. Below are the things that are needed in this particular
    example.

    Parameters
    ----------
    x:
        The dependent variable (aka 'x') that our model requires
    y:
        The dependent variable (aka 'y') that our model requires
    """

    # add inputs as class attributes
    self.x = x
    self.y = y
    
def perform(self, node, inputs, outputs):
    # the method that is used when calling the Op
    #theta, = inputs  # this will contain my variables
    theta, = inputs
    # call the log-likelihood function
    gp.fit(self.x,self.y)
    logl = (-1/2)*gp.predict(theta)

    outputs[0][0] = np.array(logl) # output the log-likelihood
    #outputs = np.array(logl) # output the log-likelihood

create our Op

logl2 = LogLike2(x_train, y_train)

priors info

priors_dictionary ={“theta”: [0.5, 0.2]}

with pm.Model() as full_bayes:

## Priors definition
# Theta prior
theta_1 = pm.Normal("theta_1", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_2 = pm.Normal("theta_2", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_3 = pm.Normal("theta_3", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_4 = pm.Normal("theta_4", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_5 = pm.Normal("theta_5", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_6 = pm.Normal("theta_6", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_7 = pm.Normal("theta_7", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_8 = pm.Normal("theta_8", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_9 = pm.Normal("theta_9", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_10 = pm.Normal("theta_10", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_11 = pm.Normal("theta_11", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_12 = pm.Normal("theta_12", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_13 = pm.Normal("theta_13", mu = priors_dictionary["theta"][0],sd= priors_dictionary["theta"][1])
theta_random= tt.as_tensor_variable([theta_1,theta_2,theta_3,theta_4,theta_5,theta_6,theta_7,theta_8,theta_9,theta_10,theta_11,theta_12,theta_13])

# Likelihood of the data

pm.DensityDist('obs',  lambda v:logl2(v), observed = {'v':theta_random})
trace = pm.sample(n_sample, tune= n_tune, chains = 1, discard_tuned_samples = False)

But I got following errors, which I am not sure how to fix. I believe it is something to do with gp.predict. and I tried ( logl = (-1/2)*gp.predict(theta.reshape(1,-1)) ) but I keep getting the same error. Could anyone help?

pm.DensityDist(‘obs’, lambda v:logl2(v), observed = {‘v’:theta_random})
File “/mnt/home/sonmooky/anaconda3/lib/python3.8/site-packages/pymc3/distributions/distribution.py”, line 47, in new
return model.Var(name, dist, data, total_size)
File “/mnt/home/sonmooky/anaconda3/lib/python3.8/site-packages/pymc3/model.py”, line 940, in Var
var = MultiObservedRV(name=name, data=data, distribution=dist,
File “/mnt/home/sonmooky/anaconda3/lib/python3.8/site-packages/pymc3/model.py”, line 1543, in init
self.logp_elemwiset = distribution.logp(**self.data)
File “MCMC_posterior_Normal_Fayan.py”, line 288, in
pm.DensityDist(‘obs’, lambda v:logl2(v), observed = {‘v’:theta_random})
File “/mnt/home/sonmooky/anaconda3/lib/python3.8/site-packages/theano/gof/op.py”, line 674, in call
required = thunk()
File “/mnt/home/sonmooky/anaconda3/lib/python3.8/site-packages/theano/gof/op.py”, line 892, in rval
r = p(n, [x[0] for x in i], o)
File “MCMC_posterior_Normal_Fayan.py”, line 231, in perform
logl = (-1/2)*gp.predict(theta)
File “/mnt/home/sonmooky/anaconda3/lib/python3.8/site-packages/sklearn/gaussian_process/_gpr.py”, line 323, in predict
X = check_array(X, ensure_2d=True, dtype=“numeric”)
File “/mnt/home/sonmooky/anaconda3/lib/python3.8/site-packages/sklearn/utils/validation.py”, line 63, in inner_f
return f(*args, **kwargs)
File “/mnt/home/sonmooky/anaconda3/lib/python3.8/site-packages/sklearn/utils/validation.py”, line 637, in check_array
raise ValueError(
ValueError: Expected 2D array, got 1D array instead:
array=[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5].
Reshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample.