Gaussian Processes

• http://www.gaussianprocess.org

import sys

import matplotlib.pyplot as pl
%matplotlib inline

import numpy as np
np.warnings.simplefilter('ignore')
import scipy.optimize as op

from scipy.linalg import cho_factor, cho_solve

def lnlike_gp(r, K):
    """
    The multivariate Gaussian ln-likelihood (up to a constant) for the
    vector ``r`` given a covariance matrix ``K``.
    
    :param r: ``(N,)``   The residual vector with ``N`` points.
    :param K: ``(N, N)`` The square (``N x N``) covariance matrix.
    
    :returns lnlike: ``float`` The Gaussian ln-likelihood. 
    
    """
    alpha = np.linalg.solve(K, r)
    s, d = np.linalg.slogdet(K)
    return -0.5 * (np.dot(r, alpha) + d)


def lnlike_gp_faster(r, K):
    """
    The multivariate Gaussian ln-likelihood (up to a constant) for the
    vector ``r`` given a covariance matrix ``K``.
    
    :param r: ``(N,)``   The residual vector with ``N`` points.
    :param K: ``(N, N)`` The square (``N x N``) covariance matrix.
    
    :returns lnlike: ``float`` The Gaussian ln-likelihood. 
    
    """
    L, lower = cho_factor(K)
    Lt_alpha = cho_solve((L,lower), r)
    ln_det = 2 * sum(np.log(np.diag(L)))
    return -0.5 * (np.dot(r, Lt_alpha) + ln_det)

def expsq_kernel(theta, dx):
    """
    The exponential-squared kernel function. The difference matrix
    can be an arbitrarily shaped numpy array so make sure that you
    use functions like ``numpy.exp`` for exponentiation.
    
    :param theta: ``(2,)`` The parameter vector ``(amp, ell)``.
    :param dx: ``numpy.array`` The difference matrix. This can be
        a numpy array with arbitrary shape.
    
    :returns K: The kernel matrix (should be the same shape as the
        input ``dx``).
        
    """
    amp, tau = theta
    return amp ** 2 * np.exp( -1 * dx**2 / (2 * tau ** 2) )

def dx(x):
    return x[None, :] - x[:, None]

fig, axs = pl.subplots(2, 3, figsize=(11,7), sharex=True, sharey=True)
x = np.linspace(-10, 10, 500)
thetas = [1,1], [1,-1], [1,2], [2,1], [2,-1], [2,2]
for i, theta in enumerate(thetas):
    K = expsq_kernel(np.exp(theta), dx(x))
    y = np.random.multivariate_normal(np.zeros_like(x), K)
    axs.flat[i].plot(x, y.T, lw=1.5)
    axs.flat[i].set_title(r"$\alpha = {},\ \tau = {}$".format(*theta), fontsize=20)

fig.tight_layout()

def expsinesq_kernel(theta, dx):
    """
    The exponential-sine-squared kernel function. The difference matrix
    can be an arbitrarily shaped numpy array so make sure that you
    use functions like ``numpy.exp`` for exponentiation.
    
    :param theta: ``(2,)`` The parameter vector ``(amp, ell)``.
    :param dx: ``numpy.array`` The difference matrix. This can be
        a numpy array with arbitrary shape.
    
    :returns K: The kernel matrix (should be the same shape as the
        input ``dx``).
        
    """    
    gamma, period = theta
    return np.exp(-gamma * np.sin(np.pi * dx / period) ** 2)


fig, axs = pl.subplots(2, 3, figsize=(11,7), sharex=True, sharey=True)
x = np.linspace(-10, 10, 500)
thetas = [1,1], [1,-1], [1,2], [2,1], [2,-1], [2,2]
for i, theta in enumerate(thetas):
    K = expsinesq_kernel(np.exp(theta), dx(x))
    y = np.random.multivariate_normal(np.zeros_like(x), K)
    axs.flat[i].plot(x, y.T, lw=1.5)
    axs.flat[i].set_title(r"$\Gamma = {},\ P = {}$".format(*theta), fontsize=20)

fig.tight_layout()

def least_squares(x, y, yerr):

    A = np.vander(x, 2)
    AT = A.T

    cov = np.linalg.inv(np.dot(AT, A / yerr[:, None] ** 2))
    mu = np.dot(cov, np.dot(AT, y / yerr ** 2))

    return mu, cov


def compute(theta, kernel, x, y, yerr, xs):

    rxx = x[:, None] - x[None, :]
    rss = xs[:, None] - xs[None, :]
    rxs = x[None, :] - xs[:, None]
    ye2 = yerr ** 2
    
    # Compute the covariance matrices.
    Kxx = kernel(theta, rxx)
    Kxx[np.diag_indices_from(Kxx)] += ye2
    Kss = kernel(theta, rss)
    Kxs = kernel(theta, rxs)

    # Compute the residuals.
    mu, _ = least_squares(x, y, yerr)
    m, b = mu
    resid = y - (m*x+b)
    model = m*xs+b
    a = np.linalg.solve(Kxx, resid)

    # Compute the predictive mean.
    mu = np.dot(Kxs, a) + model

    # Compute the predictive covariance.
    cov = Kss - np.dot(Kxs, np.linalg.solve(Kxx, Kxs.T))

    return mu, cov


def sample_conditional(theta, kernel, x, y, yerr, xs, n):
    
    mu, cov = compute(theta, kernel, x, y, yerr, xs)
    
    ys = np.random.multivariate_normal(mu, cov, n)

    return ys

n = 10
t = np.sort(20*np.random.rand(n)-10)
sig = 1
y = np.random.randn(n) * sig
yerr = np.ones_like(y) * sig

t_pred = np.linspace(-10, 10, 500)

theta = 1, 1

fig, axs = pl.subplots(1, 3, figsize=(15,3), sharex=True, sharey=True)

amp, tau = np.exp(theta)
mu = np.zeros_like(t_pred)
K = expsq_kernel(np.exp(theta), dx(t_pred))    
for ys in np.random.multivariate_normal(mu, K, 10):
    axs[0].plot(t_pred, ys, color='C9', lw=1.5, alpha=0.5)
axs[0].set_title('GP samples from prior')

theta = amp, tau
kernel = expsq_kernel
mu, cov = compute(theta, kernel, t, y, yerr, t_pred)
for ys in np.random.multivariate_normal(mu, cov, 20):
    axs[1].plot(t_pred, ys, color='C9', lw=1.5, alpha=0.5, zorder=0)
axs[1].errorbar(t, y, yerr, marker='o', color='C0', linestyle='none', zorder=1)
axs[1].set_title('GP conditional samples')

pred_mean = mu
pred_std = np.sqrt(cov.diagonal())
axs[2].plot(t_pred, pred_mean, color='C9', zorder=0)
axs[2].fill_between(t_pred, pred_mean+pred_std, pred_mean-pred_std, color='C9', alpha=0.5, edgecolor="none", zorder=0)
axs[2].errorbar(t, y, yerr, marker='o', color='C0', linestyle='none', zorder=1)
axs[2].set_title('GP conditional mean and sigma')

fig.tight_layout()

n = 10
t = np.sort(20*np.random.rand(n)-10)
sig = 1
y = np.random.randn(n) * sig
yerr = np.ones_like(y) * sig

t_pred = np.linspace(-10, 10, 100)

fig, axs = pl.subplots(2, 3, figsize=(15,8), sharex=False, sharey=False)

thetas = list(zip([1,5,10],[-1,2,5]))
for i, theta in enumerate(thetas):

    mu = np.zeros_like(t_pred)
    K = expsq_kernel(np.exp(theta), dx(t_pred))    

    axs[0,i].imshow(K, interpolation='none', cmap=pl.cm.gray)
    axs[0,i].xaxis.set_visible(False)
    axs[0,i].yaxis.set_visible(False)
    
    kernel = expsq_kernel
    mu, cov = compute(theta, kernel, t, y, yerr, t_pred)

    pred_mean = mu
    pred_std = np.sqrt(cov.diagonal())

    axs[1,i].plot(t_pred, pred_mean, color='C9', zorder=0)
    axs[1,i].fill_between(t_pred, pred_mean+pred_std, pred_mean-pred_std, color='C9', alpha=0.5, edgecolor="none", zorder=0)
    axs[1,i].errorbar(t, y, yerr, marker='o', color='C0', linestyle='none', zorder=1)
    axs[1,i].set_title(r"$\alpha = {},\ \tau = {}$".format(*theta), fontsize=20)

fig.tight_layout()