1D parameter estimation using maximum likelihood estimation

This example will cover:

• Use maximum likelihood estimation to optimise kernel parameters

from gptide import cov
from gptide import GPtideScipy
import numpy as np
import matplotlib.pyplot as plt

Generate some data

Start off with the same kernel as Example 1 and generate some data.

####
# These are our kernel input parameters
np.random.seed(1)
noise = 0.5
η = 1.5
ℓ = 100
covfunc = cov.expquad_1d

###
# Domain size parameters
dx = 25.
N = 100
covparams = (η, ℓ)

# Input data points
xd = np.arange(0,dx*N,dx)[:,None]

GP = GPtideScipy(xd, xd, noise, covfunc, covparams)

# Use the .prior() method to obtain some samples
yd = GP.prior(samples=1)

Inference

We now use the gptide.mle function do the parameter estimation. This calls the scipy.optimize.minimize routine.

# mle function import
from gptide import mle

# Initial guess of the noise and covariance parameters (these can matter)
noise_ic = 0.01
covparams_ic = [1., 10.]

# There is no mean function in this case
# meanfunc = None
# meanparams_ic = ()


soln = mle(
    xd, yd,
    covfunc,
    covparams_ic,
    noise_ic,
    verbose=False)

print('Noise (true): {:3.2f}, |Noise| (mle): {:3.2f}'.format(noise, abs(soln['x'][0])))
print('ℓ (true): {:3.2f}, ℓ (mle): {:3.2f}'.format(covparams[0], soln['x'][1]))
print('η (true): {:3.2f}, η (mle): {:3.2f}'.format(covparams[1], soln['x'][2]))

Noise (true): 0.50, |Noise| (mle): 0.45
ℓ (true): 1.50, ℓ (mle): 1.21
η (true): 100.00, η (mle): 95.02