$C$$O_2$ Mauna Loa Gaussian Process Regression¶

import warnings
from IPython.display import display_html

%matplotlib inline
import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sns, gpflow

x = np.array(data1.time)
y = np.array(data1.value)
X = x.reshape(-1, 1)
Y = y.reshape(-1, 1)

k1_ = gpflow.kernels.RBF(1, variance=(66.0 ** 2), lengthscales=67.0)
k2_exp_sine_squred_gamma= 2.0 / 1.3 ** 2.0
k2_exp_sine_squred_period = 1.0
k2_ = gpflow.kernels.RBF(1, variance=(2.4 ** 2.0), lengthscales=90.0) * gpflow.kernels.Periodic(1, period=k2_exp_sine_squred_period, variance=1.0,  lengthscales=1.0/k2_exp_sine_squred_gamma)
# k3 = how to do a rational quadratic term in GPflow?
k4_ = gpflow.kernels.RBF(1, variance=(0.18 ** 2), lengthscales=1.6) + gpflow.kernels.White(1, variance=0.19)
k1 = k1_
k2 = k2_
k4 = k4_
kernel = k1 + k2 + k4
k_ = kernel

xx = np.linspace(1990, 2000, 1000).reshape(-1,1)

def plot(m, with_data=True):
    if not isinstance(m, list):
        m = [m]
    predictions = [model.predict_y(xx) for model in m]
    
    p = bokeh.plotting.figure(title="Fit to the Mauna Loa Data", #x_axis_type='datetime', 
                          plot_width=900, plot_height=600)
    p.yaxis.axis_label = 'CO2 [ppm]'
    p.xaxis.axis_label = 'Date'
   
    if with_data:
        # true value
        p.circle(x, y, color="black", legend="Observed data")
    for mean, var in predictions:
        xpts = xx[:,0]
        p.line(xpts, mean[:,0],  line_width=1)
        lowerband = mean[:,0] - 2*np.sqrt(var[:,0])
        upperband = mean[:,0] + 2*np.sqrt(var[:,0])
        band_x = np.append(xpts, xpts[::-1])
        band_y = np.append(lowerband, upperband[::-1])
        p.patch(band_x, band_y, color='#7570B3', fill_alpha=0.2)
    p.legend.location = "top_left"
    bokeh.plotting.show(p)

import functools

# http://www.cs.toronto.edu/~duvenaud/cookbook/index.html

class GPCovarianceFunctionWhiteNoise:
    def __init__(self, sigma):
        self.sigma = float(sigma)

    def covariance(self, x1, x2):
        lr = np.zeros(shape=(len(x1), len(x2)))
        np.fill_diagonal(lr, self.sigma ** 2)
        return lr

class GPCovarianceFunctionSquaredExponential:
    def __init__(self, l, sigma):
        self.l = float(l)
        self.sigma = float(sigma)

    def covariance(self, x1, x2):
        return (self.sigma**2) * np.exp( -0.5 * (np.subtract.outer(x1, x2)**2) / (self.l ** 2))

class GPCovarianceFunctionExpSine2Kernel:
    def __init__(self, l, period, sigma):
        self.l = float(l)
        self.period = float(period)
        self.sigma = float(sigma)

    def covariance(self, x1, x2):
        return (self.sigma**2) * np.exp( -2.0 * np.power(np.sin((np.pi / self.period) * np.abs(np.subtract.outer(x1, x2))),
                                                         2) / (self.l ** 2))

class GPCovarianceFunctionSum:
    def __init__(self, k1, k2):
        self.k1 = k1
        self.k2 = k2

    def covariance(self, x1, x2):
        s1 = self.k1.covariance(x1, x2)
        s2 = self.k2.covariance(x1, x2)
        s = s1 + s2
        return  s

class GPCovarianceFunctionProduct:
    def __init__(self, k1, k2):
        self.k1 = k1
        self.k2 = k2

    def covariance(self, x1, x2):
        p1 = self.k1.covariance(x1, x2)
        p2 = self.k2.covariance(x1, x2)
        p = p1 * p2
        return  p

def conditional1(x_new, x, y, cov, sigma_n=0):
    if not isinstance(cov, list):
        cov = [cov]

    if len(cov) < 2:
        total_covariance_function = cov[0]
    else:
        total_covariance_function = functools.reduce(lambda a, x: GPCovarianceFunctionSum(a, x), cov)

    A = total_covariance_function.covariance(x_new, x_new)
    C = total_covariance_function.covariance(x_new, x)
    B = total_covariance_function.covariance(x, x) + np.power(sigma_n,2)*np.diag(np.ones(len(x)))

    mu = [np.linalg.inv(B).dot(C.T).T.dot(y).squeeze()]
    sigma = [(A - C.dot(np.linalg.inv(B).dot(C.T))).squeeze()]

    for i in range(0, len(cov)):
        partial_covariance_function = cov[i]
        C_ = partial_covariance_function.covariance(x_new, x)
        mu_ = np.linalg.inv(B).dot(C_.T).T.dot(y).squeeze()
        mu.append(mu_)

    return (mu, sigma)

def predict1(x_new, x, y, cov, sigma_n=0):
    l_y_pred, l_sigmas = conditional1(x_new, x, y, cov=cov, sigma_n=sigma_n)
    if len(l_sigmas[0].shape) > 1:
        return l_y_pred, [np.diagonal(ls) for ls in l_sigmas]
    else:
        return l_y_pred, l_sigmas

xx = np.linspace(1990, 2000, 1000)
k1 = GPCovarianceFunctionSquaredExponential(l=m.kern.rbf_1.lengthscales.value, 
                                            sigma=np.sqrt(m.kern.rbf_1.variance.value))
k2 = GPCovarianceFunctionProduct(
    GPCovarianceFunctionSquaredExponential(l=m.kern.product.rbf.lengthscales.value,
                                           sigma=np.sqrt(m.kern.product.rbf.variance.value)),
    GPCovarianceFunctionExpSine2Kernel(l=m.kern.product.periodic.lengthscales.value, 
                                       period=m.kern.product.periodic.period.value, 
                                       sigma=np.sqrt(m.kern.product.periodic.variance.value)))

k4 = GPCovarianceFunctionSum(GPCovarianceFunctionSquaredExponential(l=m.kern.rbf_2.lengthscales.value, 
                                                                    sigma=np.sqrt(m.kern.rbf_2.variance.value)),
                             GPCovarianceFunctionWhiteNoise(sigma=np.sqrt(m.kern.white.variance.value)))

k = GPCovarianceFunctionSum(k1, GPCovarianceFunctionSum(k2, k4))

y_pred, sigmas = predict1(xx, x, y, cov=[k1, k2, k4], sigma_n=np.sqrt(m.likelihood.variance.value))

# plot several years

def dates_to_idx(timelist):
    reference_time = pd.to_datetime('1959-01-01')
    t = (timelist - reference_time) / pd.Timedelta(1, "Y")
    return np.asarray(t)+1959

p = bokeh.plotting.figure(title="Several years of the seasonal component",
           plot_width=900, plot_height=600)
p.yaxis.axis_label = 'Δ CO2 [ppm]'
p.xaxis.axis_label = 'Month'

colors = bokeh.palettes.brewer['Paired'][5]
years = ["1960", "1970", "1980", "1990", "2000"]

XX = xx.reshape(-1, 1)

def conditional2(x_new, x, y, cov, sigma_n=0):
    if not isinstance(cov, list):
        cov = [cov]

    if len(cov) < 2:
        total_covariance_function = cov[0]
    else:
        # total_covariance_function = None
        total_covariance_function = functools.reduce(lambda a, x: a + x, cov)

    A = total_covariance_function.compute_K_symm(x_new)
    C = total_covariance_function.compute_K(x_new, x)
    B = total_covariance_function.compute_K_symm(x) + np.power(sigma_n,2)*np.diag(np.ones(len(x)))

    mu = [np.linalg.inv(B).dot(C.T).T.dot(y).squeeze()]
    sigma = [(A - C.dot(np.linalg.inv(B).dot(C.T))).squeeze()]

    for i in range(0, len(cov)):
        partial_covariance_function = cov[i]
        C_ = partial_covariance_function.compute_K(x_new, x)
        mu_ = np.linalg.inv(B).dot(C_.T).T.dot(y).squeeze()
        mu.append(mu_)

    return (mu, sigma)

def predict2(x_new, x, y, cov, sigma_n=0):
    l_y_pred, l_sigmas = conditional2(x_new, x, y, cov=cov, sigma_n=sigma_n)
    if len(l_sigmas[0].shape) > 1:
        return l_y_pred, [np.diagonal(ls) for ls in l_sigmas]
    else:
        return l_y_pred, l_sigmas

y_pred2, sigmas2 = predict2(XX, X, Y, cov=[k1_, k2_, k4_], sigma_n=np.sqrt(m.likelihood.variance.value))

import GPy

def conditional3(x_new, x, y, cov, sigma_n=0):
    if not isinstance(cov, list):
        cov = [cov]

    if len(cov) < 2:
        total_covariance_function = cov[0]
    else:
        # total_covariance_function = None
        total_covariance_function = functools.reduce(lambda a, x: a + x, cov)

    A = total_covariance_function.K(x_new) # total_covariance_function.K(x_new,x_new)
    C = total_covariance_function.K(x_new, x)
    B = total_covariance_function.K(x) + np.power(sigma_n,2)*np.diag(np.ones(len(x))) # total_covariance_function.K(x,x)

    mu = [np.linalg.inv(B).dot(C.T).T.dot(y).squeeze()]
    sigma = [(A - C.dot(np.linalg.inv(B).dot(C.T))).squeeze()]

    for i in range(0, len(cov)):
        partial_covariance_function = cov[i]
        C_ = partial_covariance_function.K(x_new, x)
        mu_ = np.linalg.inv(B).dot(C_.T).T.dot(y).squeeze()
        mu.append(mu_)
        
        #         kxX = C_
        #         K_1 = np.linalg.inv(B)
        #         mean = C_.dot(np.linalg.inv(B)).dot(Y)


    return (mu, sigma)

def predict3(x_new, x, y, cov, sigma_n=0):
    l_y_pred, l_sigmas = conditional3(x_new, x, y, cov=cov, sigma_n=sigma_n)
    if len(l_sigmas[0].shape) > 1:
        return l_y_pred, [np.diagonal(ls) for ls in l_sigmas]
    else:
        return l_y_pred, l_sigmas

def conditional(X_new, X, Y, cov, sigma_n=0, framework='gpy'):
    if not isinstance(cov, list):
        cov = [cov]

    if len(cov) < 2:
        total_covariance_function = cov[0]
    else:
        # total_covariance_function = None
        total_covariance_function = functools.reduce(lambda a, x: a + x, cov)

    if framework == 'gpflow':
        A = total_covariance_function.compute_K_symm(X_new)
        C = total_covariance_function.compute_K(X_new, X)
        B = total_covariance_function.compute_K_symm(X) + np.power(sigma_n,2)*np.diag(np.ones(len(X)))
    elif framework == 'gpy':
        A = total_covariance_function.K(X_new) # total_covariance_function.K(x_new,x_new)
        C = total_covariance_function.K(X_new, X)
        B = total_covariance_function.K(X) + np.power(sigma_n,2)*np.diag(np.ones(len(X))) # total_covariance_function.K(x,x)
    else:
        raise ValueError('Unknown framework parameter: {}'.format(framework))

    # mu = [np.linalg.inv(B).dot(C.T).T.dot(y).squeeze()]
    BY = np.linalg.solve(B, Y)
    mu = [C.dot(BY).squeeze()]
    sigma = [(A - C.dot(np.linalg.inv(B).dot(C.T))).squeeze()]

    for i in range(0, len(cov)):
        partial_covariance_function = cov[i]
        C_ = None
        if framework == 'gpflow':
            C_ = partial_covariance_function.compute_K(X_new, X)
        elif framework == 'gpy':
            C_ = partial_covariance_function.K(X_new, X)
        else:
            raise ValueError('Unknown framework parameter: {}'.format(framework))
        # mu_ = np.linalg.inv(B).dot(C_.T).T.dot(y).squeeze()
        mu_ = C_.dot(BY).squeeze()
        mu.append(mu_)

    return (mu, sigma)

def predict(X_new, X, Y, cov, sigma_n=0, framework='gpy'):
    l_y_pred, l_sigmas = conditional(X_new, X, Y, cov=cov, sigma_n=sigma_n, framework=framework)
    if len(l_sigmas[0].shape) > 1:
        return l_y_pred, [np.diagonal(ls) for ls in l_sigmas]
    else:
        return l_y_pred, l_sigmas

with warnings.catch_warnings():
    warnings.filterwarnings("ignore")
    y_pred2, sigmas2 = predict(XX, X, Y, cov=[k1_ , k2_ , k4_] , sigma_n=np.sqrt(m.likelihood.variance.value)   , framework='gpflow')
    y_pred3, sigmas3 = predict(XX, X, Y, cov=[k1__, k2__, k4__], sigma_n=np.sqrt(m__.Gaussian_noise.variance[0]), framework='gpy')