{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Inducing Variable Demo\n", "\n", "### 2021-05-19" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "plt.rcParams.update({'font.size': 22})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## pods\n", "\n", "\\[edit\\]\n", "\n", "In Sheffield we created a suite of software tools for ‘Open Data\n", "Science.’ Open data science is an approach to sharing code, models and\n", "data that should make it easier for companies, health professionals and\n", "scientists to gain access to data science techniques.\n", "\n", "You can also check this blog post on [Open Data\n", "Science](http://inverseprobability.com/2014/07/01/open-data-science).\n", "\n", "The software can be installed using" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install --upgrade git+https://github.com/lawrennd/ods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "from the command prompt where you can access your python installation.\n", "\n", "The code is also available on github: \n", "\n", "Once `pods` is installed, it can be imported in the usual manner." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import pods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## mlai\n", "\n", "\\[edit\\]\n", "\n", "The `mlai` software is a suite of helper functions for teaching and\n", "demonstrating machine learning algorithms. It was first used in the\n", "Machine Learning and Adaptive Intelligence course in Sheffield in 2013.\n", "\n", "The software can be installed using" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install --upgrade git+https://github.com/lawrennd/mlai.git" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "from the command prompt where you can access your python installation.\n", "\n", "The code is also available on github: \n", "\n", "Once `mlai` is installed, it can be imported in the usual manner." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mlai" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install gpy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## GPy: A Gaussian Process Framework in Python\n", "\n", "\\[edit\\]\n", "\n", "Gaussian processes are a flexible tool for non-parametric analysis with\n", "uncertainty. The GPy software was started in Sheffield to provide a easy\n", "to use interface to GPs. One which allowed the user to focus on the\n", "modelling rather than the mathematics.\n", "\n", "\n", "\n", "Figure: GPy is a BSD licensed software code base for implementing\n", "Gaussian process models in Python. It is designed for teaching and\n", "modelling. We welcome contributions which can be made through the Github\n", "repository \n", "\n", "GPy is a BSD licensed software code base for implementing Gaussian\n", "process models in python. This allows GPs to be combined with a wide\n", "variety of software libraries.\n", "\n", "The software itself is available on\n", "[GitHub](https://github.com/SheffieldML/GPy) and the team welcomes\n", "contributions.\n", "\n", "The aim for GPy is to be a probabilistic-style programming language,\n", "i.e. you specify the model rather than the algorithm. As well as a large\n", "range of covariance functions the software allows for non-Gaussian\n", "likelihoods, multivariate outputs, dimensionality reduction and\n", "approximations for larger data sets.\n", "\n", "The documentation for GPy can be found\n", "[here](https://gpy.readthedocs.io/en/latest/)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A Simple Regression Problem\n", "\n", "\\[edit\\]\n", "\n", "Here we set up a simple one dimensional regression problem. The input\n", "locations, $\\mathbf{X}$, are in two separate clusters. The response\n", "variable, $\\mathbf{ y}$, is sampled from a Gaussian process with an\n", "exponentiated quadratic covariance." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import GPy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.random.seed(101)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 50\n", "noise_var = 0.01\n", "X = np.zeros((50, 1))\n", "X[:25, :] = np.linspace(0,3,25)[:,None] # First cluster of inputs/covariates\n", "X[25:, :] = np.linspace(7,10,25)[:,None] # Second cluster of inputs/covariates\n", "\n", "# Sample response variables from a Gaussian process with exponentiated quadratic covariance.\n", "k = GPy.kern.RBF(1)\n", "y = np.random.multivariate_normal(np.zeros(N),k.K(X)+np.eye(N)*np.sqrt(noise_var)).reshape(-1,1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we perform a full Gaussian process regression on the data. We\n", "create a GP model, `m_full`, and fit it to the data, plotting the\n", "resulting fit." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m_full = GPy.models.GPRegression(X,y)\n", "_ = m_full.optimize(messages=True) # Optimize parameters of covariance function" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import mlai\n", "import mlai.plot as plot \n", "from mlai.gp_tutorial import gpplot" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "plot.model_output(m_full, ax=ax, xlabel='$x$', ylabel='$y$', fontsize=20, portion=0.2)\n", "xlim = ax.get_xlim()\n", "ylim = ax.get_ylim()\n", "mlai.write_figure(figure=fig,\n", " filename='sparse-demo-full-gp.svg',\n", " directory='./gp/')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: Full Gaussian process fitted to the data set.\n", "\n", "Now we set up the inducing variables, $\\mathbf{u}$. Each inducing\n", "variable has its own associated input index, $\\mathbf{Z}$, which lives\n", "in the same space as $\\mathbf{X}$. Here we are using the true covariance\n", "function parameters to generate the fit." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kern = GPy.kern.RBF(1)\n", "Z = np.hstack(\n", " (np.linspace(2.5,4.,3),\n", " np.linspace(7,8.5,3)))[:,None]\n", "m = GPy.models.SparseGPRegression(X,y,kernel=kern,Z=Z)\n", "m.noise_var = noise_var\n", "m.inducing_inputs.constrain_fixed()\n", "display(m)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "plot.model_output(m, ax=ax, xlabel='$x$', ylabel='$y$', fontsize=20, portion=0.2, xlim=xlim, ylim=ylim)\n", "mlai.write_figure(figure=fig,\n", " filename='sparse-demo-constrained-inducing-6-unlearned-gp.svg', \n", " directory='./gp/')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: Sparse Gaussian process fitted with six inducing variables,\n", "no optimization of parameters or inducing variables." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "_ = m.optimize(messages=True)\n", "display(m)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "plot.model_output(m, ax=ax, xlabel='$x$', ylabel='$y$', fontsize=20, portion=0.2, xlim=xlim, ylim=ylim)\n", "mlai.write_figure(figure=fig,\n", " filename='sparse-demo-constrained-inducing-6-learned-gp.svg',\n", " directory='./gp/')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: Gaussian process fitted with inducing variables fixed and\n", "parameters optimized" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m.randomize()\n", "m.inducing_inputs.unconstrain()\n", "_ = m.optimize(messages=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "plot.model_output(m, ax=ax, xlabel='$x$', ylabel='$y$', fontsize=20, portion=0.2,xlim=xlim, ylim=ylim)\n", "mlai.write_figure(figure=fig,\n", " filename='sparse-demo-unconstrained-inducing-6-gp.svg', \n", " directory='./gp/')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: Gaussian process fitted with location of inducing variables\n", "and parameters both optimized\n", "\n", "Now we will vary the number of inducing points used to form the\n", "approximation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "m.num_inducing=8\n", "m.randomize()\n", "M = 8\n", "m.set_Z(np.random.rand(M,1)*12)\n", "\n", "_ = m.optimize(messages=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n", "plot.model_output(m, ax=ax, xlabel='$x$', ylabel='$y$', fontsize=20, portion=0.2, xlim=xlim, ylim=ylim)\n", "mlai.write_figure(figure=fig,\n", " filename='sparse-demo-sparse-inducing-8-gp.svg', \n", " directory='./gp/')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "Figure: Comparison of the full Gaussian process fit with a sparse\n", "Gaussian process using eight inducing varibles. Both inducing variables\n", "and parameters are optimized.\n", "\n", "And we can compare the probability of the result to the full model." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(m.log_likelihood(), m_full.log_likelihood())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References" ] } ], "nbformat": 4, "nbformat_minor": 5, "metadata": {} }