Gaussian process simple example

Gaussian_process_example.py

# Import required libraries
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')

# Define GP Class

class GaussianProcess:

    def __init__(self, param, data, num_of_line):
        self.parameter = param
        self.data = data
        self.num_of_line = num_of_line
        self.gram_matrix = self.build_gram_matrix(self.data)

    def exponential_of_quad_kernel(self,X_n, X_m):

        theta = self.parameter
        K_Xn_Xm = theta[0] * np.exp(-0.5 * theta[1] * (X_n - X_m) ** 2) + theta[2] + theta[3] * np.dot(X_n, X_m)

        return K_Xn_Xm

    def build_gram_matrix(self,data):

        gram_matrix = np.zeros((len(data),len(data)))

        for i, X_n in enumerate(data):
            for j, X_m in enumerate(data):
                gram_matrix[i,j] = self.exponential_of_quad_kernel(X_n,X_m)

        return gram_matrix

    def distribution_of_y(self):

        mean = np.zeros(len(self.data))
        y = np.random.multivariate_normal(mean,self.gram_matrix)

        return y

    def plot(self):
        plt.figure(figsize=(6,4))
        plt.title("parameter = {}".format(str(self.parameter)))
        for i in range(self.num_of_line):
            y = self.distribution_of_y()
            plt.plot(self.data,y)

        plt.show()

    def plot_gaussian(self,scale):

        x, y = np.mgrid[-scale:scale:30j, -scale:scale:30j]
        xy = np.column_stack([x.flat, y.flat])
        mu = np.array([0.0, 0.0])
        covariance = self.build_gram_matrix(np.array([1,-1]))
        z = multivariate_normal.pdf(xy, mean=mu, cov=covariance)

        z = z.reshape(x.shape)

        fig = plt.figure()
        plt.title("parameter = {}".format(str(self.parameter)))
        ax = fig.add_subplot(111, projection='3d')

        ax.plot_surface(x, y, z)
#         ax.plot_wireframe(x,y,z)

        plt.show()

def main():
    data = np.linspace(-1, 1, num=1000)

    GP_1 = GaussianProcess([1, 4, 0, 0], data, 6)
    GP_1.plot()
    GP_1.plot_gaussian(scale=10)

    GP_2 = GaussianProcess([9, 4, 0, 0], data, 6)
    GP_2.plot()
    GP_2.plot_gaussian(scale=10)

    GP_3 = GaussianProcess([1, 64, 0, 0], data, 6)
    GP_3.plot()
    GP_3.plot_gaussian(scale=10)

    GP_4 = GaussianProcess([1, 0.25, 0, 0], data, 6)
    GP_4.plot()
    GP_4.plot_gaussian(scale=10)

    GP_5 = GaussianProcess([1, 4, 10, 0], data, 6)
    GP_5.plot()
    GP_5.plot_gaussian(scale=10)

    ## Extreme case

    GP_5_E = GaussianProcess([1, 4, 300, 0], data, 6)
    GP_5_E.plot()
    GP_5_E.plot_gaussian(scale=10)

    GP_6 = GaussianProcess([1, 4, 0, 5], data, 6)
    GP_6.plot()
    GP_6.plot_gaussian(scale=10)

    # Extreme case

    GP_6_E = GaussianProcess([1, 4, 0, 300], data, 6)
    GP_6_E.plot()
    GP_6_E.plot_gaussian(scale=10)

if __name__ == "__main__":
    main()