Marginal Likelihood optimization with python

Question

I'm trying to optimize the marginal likelihood to estimate parameters for a Gaussian process regression. So i defined the marginal log likelihood this way:

def marglike(par,X,Y):
l,sigma_n = par
n = len(X)
dist_X = (X.T - X)**2
k = np.exp(-(1/(2*(l**2)))*dist_X)
inverse = np.linalg.inv(k + (sigma_n**2)*np.eye(len(k))) 
ml = (1/2)*np.dot(np.dot(Y.T,inverse),Y) + (1/2)*np.log(np.linalg.det(k + (sigma_n**2)*np.eye(len(k)))) + (n/2)*np.log(2*np.pi)
return ml

Where the parameter to be optimized are " l " and " sigma_n ". With some initial values and data, the function give some value back:

X = np.linspace(1,10,20)
F = np.sin(X)
start = np.array([1,0.05]) #initial parameters values

marglike(start,X,F)

marglike(start,X,F)
Out[75]: array([[1872.6511786]])

But when i try to optimize the parameters with " minimize ", i get this:

re = minimize(marglike,start,args=(X,F),method="BFGS",options = {'disp':True})

re = minimize(marglike,start,args=(X,F),method="BFGS",options = {'disp':True})
Optimization terminated successfully.
     Current function value: 22.863446
     Iterations: 8
     Function evaluations: 60
     Gradient evaluations: 15

re.x
Out[89]: array([1.        , 0.70845989])

I dont know why but the parameter "l" don't seem to be optimized, but it matches the starting value that i fixed.

Any suggest ?

Answer 1

You need to reshape X to 2d first for X.T-X to work. Also, you need to add one more parameter called variance (var in the code below) in optimization. Let me know if the below code solves your problem.

from scipy.optimize import minimize

def marglike(par,X,Y):
  # print(par)
  l,var,sigma_n = par
  n = len(X)
  dist_X = (X - X.T)**2
  # print(dist_X)
  k = var*np.exp(-(1/(2*(l**2)))*dist_X)
  inverse = np.linalg.inv(k + (sigma_n**2)*np.eye(len(k))) 
  ml = (1/2)*np.dot(np.dot(Y.T,inverse),Y) + (1/2)*np.log(np.linalg.det(k + (sigma_n**2)*np.eye(len(k)))) + (n/2)*np.log(2*np.pi)
  return ml

X = np.linspace(1,10,20).reshape(-1,1) # Reshaping
F = np.sin(X)
start = np.array([1.1,1.6,0.05]) #initial parameters values
print(marglike(start,X,F))

re = minimize(marglike,start,args=(X,F),method="L-BFGS-B",options = {'disp':True})
re.x