diff --git a/GPy/GPyCrashCourse.ipynb b/GPy/GPyCrashCourse.ipynb
index 40829d4..16844d1 100644
--- a/GPy/GPyCrashCourse.ipynb
+++ b/GPy/GPyCrashCourse.ipynb
@@ -16,6 +16,34 @@
     "In this tutorial I will run through a few basic features of GPy."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n",
+    "try:\n",
+    "    import pods\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", 'pods']\n",
+    "    )\n",
+    "    import pods\n",
+    "\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -26,9 +54,8 @@
     "import pods\n",
     "from matplotlib import pyplot as plt\n",
     "np.random.seed(1)\n",
-    "%matplotlib inline\n",
     "#%config InlineBackend.figure_format = 'svg'\n",
-    "#import matplotlib;matplotlib.rcParams['figure.figsize'] = (8,5)"
+    "#import matplotlib;matplotlib.rcParams['figure.figsize'] = (8,5)\n"
    ]
   },
   {
diff --git a/GPy/Integral_kernel.ipynb b/GPy/Integral_kernel.ipynb
index 6335457..3a36ab8 100644
--- a/GPy/Integral_kernel.ipynb
+++ b/GPy/Integral_kernel.ipynb
@@ -8,6 +8,24 @@
     "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/SheffieldML/notebook/blob/main/GPy/Integral_kernel.ipynb)\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -19,7 +37,6 @@
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "\n",
-    "%matplotlib inline\n",
     "\n",
     "def plot_hist(bins,bin_ranges,errs,alpha):\n",
     "    \"\"\"Plots a histogram\"\"\"\n",
@@ -46,7 +63,7 @@
     "    simple_inputs = np.array(simple_inputs)\n",
     "\n",
     "    original_outputs = np.array(outputs)\n",
-    "    return inputs,simple_inputs,original_outputs"
+    "    return inputs,simple_inputs,original_outputs\n"
    ]
   },
   {
@@ -240,13 +257,6 @@
    "source": [
     "The above curve is based on the histogram 'binned' data."
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/GPy/MagnificationFactor.ipynb b/GPy/MagnificationFactor.ipynb
index 3cac298..767e273 100644
--- a/GPy/MagnificationFactor.ipynb
+++ b/GPy/MagnificationFactor.ipynb
@@ -20,6 +20,24 @@
     "Note: Requires GPy Version >= 0.8"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -29,11 +47,10 @@
     "import numpy as np\n",
     "if not hasattr(np, \"NaN\"):\n",
     "    np.NaN = np.nan  # pods mocap still uses np.NaN\n",
-    "%matplotlib inline\n",
     "from matplotlib import pyplot as plt\n",
     "import GPy # import GPy package\n",
     "np.random.seed(12345)\n",
-    "GPy.plotting.change_plotting_library('matplotlib')"
+    "GPy.plotting.change_plotting_library('matplotlib')\n"
    ]
   },
   {
diff --git a/GPy/ParametricNonParametricInference.ipynb b/GPy/ParametricNonParametricInference.ipynb
index d0e978e..905467d 100644
--- a/GPy/ParametricNonParametricInference.ipynb
+++ b/GPy/ParametricNonParametricInference.ipynb
@@ -74,6 +74,24 @@
     "$$\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -82,8 +100,7 @@
    "source": [
     "import GPy, numpy as np, pandas as pd\n",
     "from GPy.kern import LinearSlopeBasisFuncKernel, DomainKernel, ChangePointBasisFuncKernel\n",
-    "%matplotlib inline\n",
-    "from matplotlib import pyplot as plt"
+    "from matplotlib import pyplot as plt\n"
    ]
   },
   {
diff --git a/GPy/Poisson regression tutorial.ipynb b/GPy/Poisson regression tutorial.ipynb
index 7ed664d..5c113dd 100644
--- a/GPy/Poisson regression tutorial.ipynb	
+++ b/GPy/Poisson regression tutorial.ipynb	
@@ -1,390 +1,407 @@
 {
-  "cells": [
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "# Count Data Tutorial\n",
-        "\n",
-        "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/SheffieldML/notebook/blob/main/GPy/Poisson%20regression%20tutorial.ipynb)\n"
-      ]
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "## 28th May 2014"
-      ]
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "### by Alan Saul, with edits by Neil Lawrence"
-      ]
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "In this tutorial we will look at using a Gaussian Process with a Poisson likelihood to model a response variable that is count data, i.e $ Y \\in \\mathbb{N}^{0} = \\{0, 1, 2, 3, ..\\}$.\n",
-        "\n",
-        "This allows us to incorporate an approximate Poisson Process regression model into our Gaussian Process framework, indeed as the Poisson rate parameter, $\\lambda$, becomes large the Poisson Process asymptotically approaches a Gaussian Process. This is important as if your counts are large, the rate, $\\lambda$ at which events is happening is likely to be large, and it may not be required to use the approximation routine used here. Instead it may be more appropriate simply use a Gaussian Process, this will drastically reduce the computation required without a significant loss in accuracy.\n",
-        "\n",
-        "The approximation we will use here is known as the Laplace approximation. The Laplace approximation was applied to Gaussian processes for classification by [Williams and Barber, 1999](http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=735807&url=http%3A%2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D735807). More details are given in [Rasmussen and Williams, 2006](http://www.gaussianprocess.org/gpml/),  Section 3.4 and Section 5.5.1. It has also been suggested for Gaussian random fields by [Rue, Martino and Chopin, 2009](http://hachamam-inla.googlecode.com/hg/r-inla.org/papers/inla-rss.pdf). The Laplace approximation is a local approximation. It proceeds by finding the mode of the posterior distribution and computing the curvature to find the fit.\n",
-        "\n",
-        "First we perform some setup."
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "%pylab inline\n",
-        "import numpy as np\n",
-        "import scipy as sp\n",
-        "from scipy import stats\n",
-        "import pylab as pb"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "Next we show the Poisson distribution for several different rate parameters."
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "K = np.arange(0, 30, 1)[:, None]\n",
-        "rates = np.arange(0, 40, 1)[:, None]\n",
-        "\n",
-        "#Make a matrix with PMF for (rate x count) combinations\n",
-        "rates_count = np.array([sp.stats.poisson(r).pmf(K) for r in rates]).reshape(rates.shape[0], K.shape[0])\n",
-        "\n",
-        "rates_plot = [1,3,8,20]\n",
-        "\n",
-        "#Plot each rate as a function of counts K\n",
-        "for r in rates_plot:\n",
-        "    pb.plot(K, rates_count[r, :], label='rate = {}'.format(float(rates[r].item())))\n",
-        "\n",
-        "pb.title('Poisson probability mass distribution for different rates')\n",
-        "pb.ylabel('PMF')\n",
-        "pb.xlabel('Counts')\n",
-        "pb.legend()"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "As the rate increases, the Poisson distribution over discrete values begins to take on the form of a Gaussian density, over continuous values, as can be seen below."
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "Ks = [1,5,10,25]\n",
-        "pb.figure()\n",
-        "for k in Ks:\n",
-        "    pb.plot(rates, rates_count[:, k], label='K = {}'.format(int(K[k].item())))\n",
-        "\n",
-        "pb.xlabel('rate')\n",
-        "pb.ylabel('PMF')\n",
-        "pb.legend()"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "The variance of the Poisson distribution is equal to its mean. The Poisson only has support for non-negative integers, whereas a Gaussian has support over all real numbers (including negatives). This can complicate matters because if we place a Gaussian process prior directly over the rate of the Poisson, then the implication is that the rate can go negative. In practice we normally consider a Gaussian process over the logarithm of the rate. This is known as a log Gaussian process. It is formed by combining a logarithmic link function with a Gaussian model. "
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "small_rate = 2\n",
-        "K = np.arange(0, 20, 1)[:, None]\n",
-        "Kcont = np.linspace(0, 20, 100)[:, None]\n",
-        "gauss20 = sp.stats.norm(loc=small_rate, scale=np.sqrt(small_rate)).pdf(Kcont)\n",
-        "poisson20 = sp.stats.poisson(small_rate).pmf(K)\n",
-        "pb.plot(K, poisson20, label='Poisson')\n",
-        "pb.plot(Kcont, gauss20, label='Gaussian')\n",
-        "pb.title('Gaussian and Poisson small rate')\n",
-        "pb.ylabel('PDF/PMF')\n",
-        "pb.xlabel('Counts')\n",
-        "pb.legend()\n",
-        "\n",
-        "pb.figure()\n",
-        "large_rate = 40\n",
-        "K = np.arange(0, 80, 1)[:, None]\n",
-        "Kcont = np.linspace(0, 80, 200)[:, None]\n",
-        "gauss40 = sp.stats.norm(loc=large_rate, scale=np.sqrt(large_rate)).pdf(Kcont)\n",
-        "poisson40 = sp.stats.poisson(large_rate).pmf(K)\n",
-        "pb.plot(K, poisson40, label='Poisson')\n",
-        "pb.plot(Kcont, gauss40, label='Gaussian')\n",
-        "pb.title('Gaussian and Poisson large rate')\n",
-        "pb.ylabel('PDF/PMF')\n",
-        "pb.xlabel('Counts')\n",
-        "pb.legend()"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "rates = np.arange(0, 20, 0.01)[:, None]\n",
-        "small_counts = 1\n",
-        "gauss3 = sp.stats.norm(loc=rates, scale=np.sqrt(rates)).pdf(small_counts)\n",
-        "#gauss3 = np.array([sp.stats.norm(loc=r, scale=np.sqrt(r)).pdf(small_counts) for r in rates])\n",
-        "poisson3 = sp.stats.poisson(rates).pmf(small_counts)\n",
-        "pb.plot(rates, poisson3, label='Poisson')\n",
-        "pb.plot(rates, gauss3, label='Gaussian')\n",
-        "pb.title('Gaussian and Poisson small count (K = {})'.format(small_counts))\n",
-        "pb.ylabel('PDF/PMF')\n",
-        "pb.xlabel('Rate')\n",
-        "pb.legend()\n",
-        "\n",
-        "pb.figure()\n",
-        "rates = np.arange(0, 100, 0.01)[:, None]\n",
-        "large_counts = 40\n",
-        "gauss40 = sp.stats.norm(loc=rates, scale=np.sqrt(rates)).pdf(large_counts)\n",
-        "poisson40 = sp.stats.poisson(rates).pmf(large_counts)\n",
-        "pb.plot(rates, poisson40, label='Poisson')\n",
-        "pb.plot(rates, gauss40, label='Gaussian')\n",
-        "pb.title('Gaussian and Poisson large count (K = {})'.format(large_counts))\n",
-        "pb.ylabel('PDF/PMF')\n",
-        "pb.xlabel('Rate')\n",
-        "pb.legend()\n"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "If the counts are believed to follow a Poisson distribution and are very large, then the underlying rate is also likely to be large and we might reasonably approximate the Poisson likelihood with a Gaussian likelihood by matching the mean and variance to the observed counts. \n",
-        "\n",
-        "If, however, counts are relatively small, and we wish to fit Poisson process regression into a Gaussian Process framework, we may need a more involved approximation to the posterior.\n",
-        "\n",
-        "The approximation used in this tutorial is called the Laplace approximation, it uses the location of the mode of the marginal distribution and the curvature of this mode (the hessian) to find a more appropriate Gaussian posterior approximation. As we shall see, in practice this is quite effective, and future tutorials will show us how it can be used to make Gaussian Process approximations to a range of different likelihoods, rather than just the Poisson, with relative ease.\n",
-        "\n"
-      ]
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "Approximate Poisson Process Regression\n",
-        "======================================\n",
-        "\n",
-        "Poisson process regression, contrasting to the above example, allows the rate to change, for example over time. That is the intensity of the counts can change, initially perhaps the data suggests there were few counts due to a low number of counts, but later this intensity is assumed to rise as we start observing larger more frequent counts. \n",
-        "\n",
-        "An example of data that might follow this change in 'intensity' is the times at which buses come throughout the day, we will choose this as our toy modelling example.\n",
-        "\n",
-        "First we will generate some toy data, we will say that there are very few buses overnight, thus low counts, then during the morning rush-hour the number of observed counts rises, it then reduces and plateus until the evening rush-hour, and then finally as the evening draws on we observe again a reduced number "
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "X = np.linspace(0,23.50,48)[:, None]\n",
-        "intensities = 4 + 3*np.sin(10 + 0.6*X) + np.sin(0.1*X)\n",
-        "pb.plot(X, intensities)\n",
-        "pb.title('Real underlying intensities (rate at which buses are issued)')\n",
-        "\n",
-        "Y = np.random.poisson(intensities.ravel())[:, None]\n",
-        "pb.figure()\n",
-        "pb.bar(X.ravel(), Y.ravel(), width=0.5)\n",
-        "pb.xlabel('Time')\n",
-        "pb.ylabel('Buses in half an hour')\n",
-        "pb.title('Observed (noisey) counts of buses throughout day')"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "We will now fit a Gaussian Process to this data using a Poisson likelihood and a Laplace approximation using GPy, as we shall see this is relatively simple, however to use Gaussian Processes we must make some assumptions about the intensity, we do this by defining a kernel."
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "import GPy\n",
-        "kernel = GPy.kern.RBF(1, variance=1.0, lengthscale=1.0)"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "By defining a radial basis function (RBF) kernel (also known as a exponentiated quadratic, squared exponential or Gaussian kernel), we are making *a priori* a smoothness assumption about the intensity parameter throughout time in this case. This is a very general assumption that is required in many learning algorithms, that data 'close' to one another shares similar properties. In this case we are saying that times that are close to one another have a similar number of buses arriving. \n",
-        "\n",
-        "The measure of closeness is defined by the lengthscale (or timescale) parameter, $\\ell$, initially we set this to be 2.0, roughly speaking this says that as we have moved two time-steps (1 hour in this case) the intensity shares almost nothing with the intensity a hour previous.\n",
-        "\n",
-        "The variance parameter describes our prior assumptions about the scale of the function being learnt. Both of these parameters although initialised here are usually optimized to find the most appropriate parameters, as such we set them to something that seems sensible but do not need to worry about having them absolutely correct.\n",
-        "\n",
-        "There are a variety of likelihoods implemented in GPy, type 'GPy.likelihoods.' then press tab to list them.\n",
-        "\n",
-        "There is also a number of alternative inference methods implemented, a subject for another tutorial, however it is worth noting that not all likelihoods are yet compatible with all inference methods. \n",
-        "\n",
-        "We now choose the type of likelihood and inference method we would like to use, in this case the Poisson likelihood and Laplace inference method. "
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "poisson_likelihood = GPy.likelihoods.Poisson()\n",
-        "laplace_inf = GPy.inference.latent_function_inference.Laplace()"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "Finally we are ready to create our model Gaussian Process model"
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "m = GPy.core.GP(X=X, Y=Y, likelihood=poisson_likelihood, inference_method=laplace_inf, kernel=kernel)\n",
-        "print(m)\n",
-        "\n",
-        "#Predictions of counts of buses, including noise arising from Poisson noise\n",
-        "m.plot()\n",
-        "pb.xlabel('Time')\n",
-        "pb.ylabel('Buses counts per 30 minutes')"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "As we can see, the model as it stands is not all that great, the lengthscale parameter is soo small. We noted before that the kernels hyperparameters we chose were not the correct ones, only an initial guess. Next we try maximizing the likelihood of the model with respect to the parameters. Optimization within GPy is straightforward."
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "m.optimize()\n",
-        "print(m)\n",
-        "\n",
-        "#Predictions of counts of buses, including noise arising from Poisson noise\n",
-        "m.plot()\n",
-        "pb.xlabel('Time')\n",
-        "pb.ylabel('Buses counts per 30 minutes')"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "Now we should have a clearly better fit to the data, we have made predictions of the number of counts we would see throughout the day with our model and it seems to fit the training data we have provided quite well.\n",
-        "\n",
-        "To get a better indication of the intensity we have infered we can plot the latent function $f$ of the model, this is the models prediction, before we have added Poisson noise to it. In fact the model is modelling the log intensity, this is what is displayed below"
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "#Posterior distribution of the intensity of buses, and true intensity in red\n",
-        "m.plot_f()\n",
-        "pb.plot(X, np.log(intensities), '--r', linewidth=2)\n",
-        "pb.title('Real log intensity vs posterior belief')\n",
-        "pb.xlabel('Time')\n",
-        "pb.ylabel('Log intensity')\n",
-        "plt.ylim(-2,3)"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "To get the predictions from the model is straightforward with the plot functionality."
-      ]
-    },
-    {
-      "cell_type": "code",
-      "metadata": {},
-      "source": [
-        "pred_points = np.linspace(0,24,100)[:, None]\n",
-        "#Predictive GP for log intensity mean and variance\n",
-        "f_mean, f_var = m._raw_predict(pred_points)\n",
-        "f_upper, f_lower = f_mean + 2*np.sqrt(f_var), f_mean - 2.*np.sqrt(f_var)\n",
-        "pb.plot(X, intensities, '--r', linewidth=2, label='true intensity')\n",
-        "#Plotting Y on an exponential scale as we are now looking at intensity rather than log intensity\n",
-        "from GPy.plotting import Tango\n",
-        "Tango.reset()\n",
-        "b = Tango.nextMedium()\n",
-        "pb.plot(pred_points, np.exp(f_mean), color=b, lw=2)\n",
-        "pb.fill_between(pred_points[:,0], np.exp(f_lower[:,0]), np.exp(f_upper[:,0]), color=b, alpha=.1)\n",
-        "pb.title('Real intensity vs posterior belief')\n",
-        "pb.xlabel('Time')\n",
-        "pb.ylabel('intensity')\n",
-        "pb.legend()"
-      ],
-      "execution_count": null,
-      "outputs": []
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": []
-    }
-  ],
-  "metadata": {
-    "kernelspec": {
-      "display_name": "Python 3",
-      "language": "python",
-      "name": "python3"
-    },
-    "language_info": {
-      "codemirror_mode": {
-        "name": "ipython",
-        "version": 3
-      },
-      "file_extension": ".py",
-      "mimetype": "text/x-python",
-      "name": "python",
-      "nbconvert_exporter": "python",
-      "pygments_lexer": "ipython3",
-      "version": "3.12.0"
-    }
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Count Data Tutorial\n",
+    "\n",
+    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/SheffieldML/notebook/blob/main/GPy/Poisson%20regression%20tutorial.ipynb)\n"
+   ]
   },
-  "nbformat": 4,
-  "nbformat_minor": 0
-}
\ No newline at end of file
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 28th May 2014"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### by Alan Saul, with edits by Neil Lawrence"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this tutorial we will look at using a Gaussian Process with a Poisson likelihood to model a response variable that is count data, i.e $ Y \\in \\mathbb{N}^{0} = \\{0, 1, 2, 3, ..\\}$.\n",
+    "\n",
+    "This allows us to incorporate an approximate Poisson Process regression model into our Gaussian Process framework, indeed as the Poisson rate parameter, $\\lambda$, becomes large the Poisson Process asymptotically approaches a Gaussian Process. This is important as if your counts are large, the rate, $\\lambda$ at which events is happening is likely to be large, and it may not be required to use the approximation routine used here. Instead it may be more appropriate simply use a Gaussian Process, this will drastically reduce the computation required without a significant loss in accuracy.\n",
+    "\n",
+    "The approximation we will use here is known as the Laplace approximation. The Laplace approximation was applied to Gaussian processes for classification by [Williams and Barber, 1999](http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=735807&url=http%3A%2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D735807). More details are given in [Rasmussen and Williams, 2006](http://www.gaussianprocess.org/gpml/),  Section 3.4 and Section 5.5.1. It has also been suggested for Gaussian random fields by [Rue, Martino and Chopin, 2009](http://hachamam-inla.googlecode.com/hg/r-inla.org/papers/inla-rss.pdf). The Laplace approximation is a local approximation. It proceeds by finding the mode of the posterior distribution and computing the curvature to find the fit.\n",
+    "\n",
+    "First we perform some setup."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matplotlib import pyplot as plt\n",
+    "import numpy as np\n",
+    "import scipy as sp\n",
+    "from scipy import stats\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we show the Poisson distribution for several different rate parameters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "K = np.arange(0, 30, 1)[:, None]\n",
+    "rates = np.arange(0, 40, 1)[:, None]\n",
+    "\n",
+    "#Make a matrix with PMF for (rate x count) combinations\n",
+    "rates_count = np.array([sp.stats.poisson(r).pmf(K) for r in rates]).reshape(rates.shape[0], K.shape[0])\n",
+    "\n",
+    "rates_plot = [1,3,8,20]\n",
+    "\n",
+    "#Plot each rate as a function of counts K\n",
+    "for r in rates_plot:\n",
+    "    plt.plot(K, rates_count[r, :], label='rate = {}'.format(float(rates[r].item())))\n",
+    "\n",
+    "plt.title('Poisson probability mass distribution for different rates')\n",
+    "plt.ylabel('PMF')\n",
+    "plt.xlabel('Counts')\n",
+    "plt.legend()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As the rate increases, the Poisson distribution over discrete values begins to take on the form of a Gaussian density, over continuous values, as can be seen below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Ks = [1,5,10,25]\n",
+    "plt.figure()\n",
+    "for k in Ks:\n",
+    "    plt.plot(rates, rates_count[:, k], label='K = {}'.format(int(K[k].item())))\n",
+    "\n",
+    "plt.xlabel('rate')\n",
+    "plt.ylabel('PMF')\n",
+    "plt.legend()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The variance of the Poisson distribution is equal to its mean. The Poisson only has support for non-negative integers, whereas a Gaussian has support over all real numbers (including negatives). This can complicate matters because if we place a Gaussian process prior directly over the rate of the Poisson, then the implication is that the rate can go negative. In practice we normally consider a Gaussian process over the logarithm of the rate. This is known as a log Gaussian process. It is formed by combining a logarithmic link function with a Gaussian model. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "small_rate = 2\n",
+    "K = np.arange(0, 20, 1)[:, None]\n",
+    "Kcont = np.linspace(0, 20, 100)[:, None]\n",
+    "gauss20 = sp.stats.norm(loc=small_rate, scale=np.sqrt(small_rate)).pdf(Kcont)\n",
+    "poisson20 = sp.stats.poisson(small_rate).pmf(K)\n",
+    "plt.plot(K, poisson20, label='Poisson')\n",
+    "plt.plot(Kcont, gauss20, label='Gaussian')\n",
+    "plt.title('Gaussian and Poisson small rate')\n",
+    "plt.ylabel('PDF/PMF')\n",
+    "plt.xlabel('Counts')\n",
+    "plt.legend()\n",
+    "\n",
+    "plt.figure()\n",
+    "large_rate = 40\n",
+    "K = np.arange(0, 80, 1)[:, None]\n",
+    "Kcont = np.linspace(0, 80, 200)[:, None]\n",
+    "gauss40 = sp.stats.norm(loc=large_rate, scale=np.sqrt(large_rate)).pdf(Kcont)\n",
+    "poisson40 = sp.stats.poisson(large_rate).pmf(K)\n",
+    "plt.plot(K, poisson40, label='Poisson')\n",
+    "plt.plot(Kcont, gauss40, label='Gaussian')\n",
+    "plt.title('Gaussian and Poisson large rate')\n",
+    "plt.ylabel('PDF/PMF')\n",
+    "plt.xlabel('Counts')\n",
+    "plt.legend()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rates = np.arange(0, 20, 0.01)[:, None]\n",
+    "small_counts = 1\n",
+    "gauss3 = sp.stats.norm(loc=rates, scale=np.sqrt(rates)).pdf(small_counts)\n",
+    "#gauss3 = np.array([sp.stats.norm(loc=r, scale=np.sqrt(r)).pdf(small_counts) for r in rates])\n",
+    "poisson3 = sp.stats.poisson(rates).pmf(small_counts)\n",
+    "plt.plot(rates, poisson3, label='Poisson')\n",
+    "plt.plot(rates, gauss3, label='Gaussian')\n",
+    "plt.title('Gaussian and Poisson small count (K = {})'.format(small_counts))\n",
+    "plt.ylabel('PDF/PMF')\n",
+    "plt.xlabel('Rate')\n",
+    "plt.legend()\n",
+    "\n",
+    "plt.figure()\n",
+    "rates = np.arange(0, 100, 0.01)[:, None]\n",
+    "large_counts = 40\n",
+    "gauss40 = sp.stats.norm(loc=rates, scale=np.sqrt(rates)).pdf(large_counts)\n",
+    "poisson40 = sp.stats.poisson(rates).pmf(large_counts)\n",
+    "plt.plot(rates, poisson40, label='Poisson')\n",
+    "plt.plot(rates, gauss40, label='Gaussian')\n",
+    "plt.title('Gaussian and Poisson large count (K = {})'.format(large_counts))\n",
+    "plt.ylabel('PDF/PMF')\n",
+    "plt.xlabel('Rate')\n",
+    "plt.legend()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If the counts are believed to follow a Poisson distribution and are very large, then the underlying rate is also likely to be large and we might reasonably approximate the Poisson likelihood with a Gaussian likelihood by matching the mean and variance to the observed counts. \n",
+    "\n",
+    "If, however, counts are relatively small, and we wish to fit Poisson process regression into a Gaussian Process framework, we may need a more involved approximation to the posterior.\n",
+    "\n",
+    "The approximation used in this tutorial is called the Laplace approximation, it uses the location of the mode of the marginal distribution and the curvature of this mode (the hessian) to find a more appropriate Gaussian posterior approximation. As we shall see, in practice this is quite effective, and future tutorials will show us how it can be used to make Gaussian Process approximations to a range of different likelihoods, rather than just the Poisson, with relative ease.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Approximate Poisson Process Regression\n",
+    "======================================\n",
+    "\n",
+    "Poisson process regression, contrasting to the above example, allows the rate to change, for example over time. That is the intensity of the counts can change, initially perhaps the data suggests there were few counts due to a low number of counts, but later this intensity is assumed to rise as we start observing larger more frequent counts. \n",
+    "\n",
+    "An example of data that might follow this change in 'intensity' is the times at which buses come throughout the day, we will choose this as our toy modelling example.\n",
+    "\n",
+    "First we will generate some toy data, we will say that there are very few buses overnight, thus low counts, then during the morning rush-hour the number of observed counts rises, it then reduces and plateus until the evening rush-hour, and then finally as the evening draws on we observe again a reduced number "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "X = np.linspace(0,23.50,48)[:, None]\n",
+    "intensities = 4 + 3*np.sin(10 + 0.6*X) + np.sin(0.1*X)\n",
+    "plt.plot(X, intensities)\n",
+    "plt.title('Real underlying intensities (rate at which buses are issued)')\n",
+    "\n",
+    "Y = np.random.poisson(intensities.ravel())[:, None]\n",
+    "plt.figure()\n",
+    "plt.bar(X.ravel(), Y.ravel(), width=0.5)\n",
+    "plt.xlabel('Time')\n",
+    "plt.ylabel('Buses in half an hour')\n",
+    "plt.title('Observed (noisey) counts of buses throughout day')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We will now fit a Gaussian Process to this data using a Poisson likelihood and a Laplace approximation using GPy, as we shall see this is relatively simple, however to use Gaussian Processes we must make some assumptions about the intensity, we do this by defining a kernel."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import GPy\n",
+    "kernel = GPy.kern.RBF(1, variance=1.0, lengthscale=1.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "By defining a radial basis function (RBF) kernel (also known as a exponentiated quadratic, squared exponential or Gaussian kernel), we are making *a priori* a smoothness assumption about the intensity parameter throughout time in this case. This is a very general assumption that is required in many learning algorithms, that data 'close' to one another shares similar properties. In this case we are saying that times that are close to one another have a similar number of buses arriving. \n",
+    "\n",
+    "The measure of closeness is defined by the lengthscale (or timescale) parameter, $\\ell$, initially we set this to be 2.0, roughly speaking this says that as we have moved two time-steps (1 hour in this case) the intensity shares almost nothing with the intensity a hour previous.\n",
+    "\n",
+    "The variance parameter describes our prior assumptions about the scale of the function being learnt. Both of these parameters although initialised here are usually optimized to find the most appropriate parameters, as such we set them to something that seems sensible but do not need to worry about having them absolutely correct.\n",
+    "\n",
+    "There are a variety of likelihoods implemented in GPy, type 'GPy.likelihoods.' then press tab to list them.\n",
+    "\n",
+    "There is also a number of alternative inference methods implemented, a subject for another tutorial, however it is worth noting that not all likelihoods are yet compatible with all inference methods. \n",
+    "\n",
+    "We now choose the type of likelihood and inference method we would like to use, in this case the Poisson likelihood and Laplace inference method. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "poisson_likelihood = GPy.likelihoods.Poisson()\n",
+    "laplace_inf = GPy.inference.latent_function_inference.Laplace()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Finally we are ready to create our model Gaussian Process model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m = GPy.core.GP(X=X, Y=Y, likelihood=poisson_likelihood, inference_method=laplace_inf, kernel=kernel)\n",
+    "print(m)\n",
+    "\n",
+    "#Predictions of counts of buses, including noise arising from Poisson noise\n",
+    "m.plot()\n",
+    "plt.xlabel('Time')\n",
+    "plt.ylabel('Buses counts per 30 minutes')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As we can see, the model as it stands is not all that great, the lengthscale parameter is soo small. We noted before that the kernels hyperparameters we chose were not the correct ones, only an initial guess. Next we try maximizing the likelihood of the model with respect to the parameters. Optimization within GPy is straightforward."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m.optimize()\n",
+    "print(m)\n",
+    "\n",
+    "#Predictions of counts of buses, including noise arising from Poisson noise\n",
+    "m.plot()\n",
+    "plt.xlabel('Time')\n",
+    "plt.ylabel('Buses counts per 30 minutes')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we should have a clearly better fit to the data, we have made predictions of the number of counts we would see throughout the day with our model and it seems to fit the training data we have provided quite well.\n",
+    "\n",
+    "To get a better indication of the intensity we have infered we can plot the latent function $f$ of the model, this is the models prediction, before we have added Poisson noise to it. In fact the model is modelling the log intensity, this is what is displayed below"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Posterior distribution of the intensity of buses, and true intensity in red\n",
+    "m.plot_f()\n",
+    "plt.plot(X, np.log(intensities), '--r', linewidth=2)\n",
+    "plt.title('Real log intensity vs posterior belief')\n",
+    "plt.xlabel('Time')\n",
+    "plt.ylabel('Log intensity')\n",
+    "plt.ylim(-2,3)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To get the predictions from the model is straightforward with the plot functionality."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pred_points = np.linspace(0,24,100)[:, None]\n",
+    "#Predictive GP for log intensity mean and variance\n",
+    "f_mean, f_var = m._raw_predict(pred_points)\n",
+    "f_upper, f_lower = f_mean + 2*np.sqrt(f_var), f_mean - 2.*np.sqrt(f_var)\n",
+    "plt.plot(X, intensities, '--r', linewidth=2, label='true intensity')\n",
+    "#Plotting Y on an exponential scale as we are now looking at intensity rather than log intensity\n",
+    "from GPy.plotting import Tango\n",
+    "Tango.reset()\n",
+    "b = Tango.nextMedium()\n",
+    "plt.plot(pred_points, np.exp(f_mean), color=b, lw=2)\n",
+    "plt.fill_between(pred_points[:,0], np.exp(f_lower[:,0]), np.exp(f_upper[:,0]), color=b, alpha=.1)\n",
+    "plt.title('Real intensity vs posterior belief')\n",
+    "plt.xlabel('Time')\n",
+    "plt.ylabel('intensity')\n",
+    "plt.legend()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/GPy/SVI.ipynb b/GPy/SVI.ipynb
index c441ab9..b51dff2 100644
--- a/GPy/SVI.ipynb
+++ b/GPy/SVI.ipynb
@@ -34,6 +34,34 @@
     "This notebook needs [climin](https://github.com/BRML/climin) for stochastic optimisation. Climin is an **optional extra**, not a GPy dependency. Install it with `pip install climin` before running the cells below."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n",
+    "try:\n",
+    "    import ipywidgets\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", 'ipywidgets']\n",
+    "    )\n",
+    "    import ipywidgets\n",
+    "\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -68,9 +96,8 @@
     "            for i in idxs:\n",
     "                yield slices[i]\n",
     "climin.util.draw_mini_slices = _draw_mini_slices\n",
-    "%matplotlib inline\n",
     "%config InlineBackend.figure_format = 'svg'\n",
-    "#plt.rcParams['figure.figsize'] = 8,8"
+    "#plt.rcParams['figure.figsize'] = 8,8\n"
    ]
   },
   {
@@ -208,13 +235,6 @@
     "\n",
     "fig1.tight_layout()"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/GPy/Space_time_effects.ipynb b/GPy/Space_time_effects.ipynb
index 737a090..aff5543 100644
--- a/GPy/Space_time_effects.ipynb
+++ b/GPy/Space_time_effects.ipynb
@@ -17,6 +17,24 @@
     "$y_{s,t} | λ(s, t) ∼ Poisson(\\lambda)$\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -27,8 +45,7 @@
     "import numpy as np\n",
     "import GPy\n",
     "from scipy.stats import poisson\n",
-    "np.random.seed(100)\n",
-    "%matplotlib inline"
+    "np.random.seed(100)\n"
    ]
   },
   {
diff --git a/GPy/basic_classification.ipynb b/GPy/basic_classification.ipynb
index 1c015e9..858350c 100644
--- a/GPy/basic_classification.ipynb
+++ b/GPy/basic_classification.ipynb
@@ -8,6 +8,24 @@
     "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/SheffieldML/notebook/blob/main/GPy/basic_classification.ipynb)\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -15,7 +33,6 @@
    "outputs": [],
    "source": [
     "#configure plotting\n",
-    "%matplotlib inline\n",
     "%config InlineBackend.figure_format = 'svg'\n",
     "import matplotlib;matplotlib.rcParams['figure.figsize'] = (8,5)\n",
     "import matplotlib;matplotlib.rcParams['text.usetex'] = True\n",
@@ -24,7 +41,7 @@
     "import GPy\n",
     "import numpy as np\n",
     "from matplotlib import pyplot as plt\n",
-    "np.random.seed(1)"
+    "np.random.seed(1)\n"
    ]
   },
   {
@@ -248,13 +265,6 @@
     "plt.ylabel('$f(x)$');plt.xlabel('$x$')\n",
     "plt.legend()"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/GPy/basic_gp.ipynb b/GPy/basic_gp.ipynb
index 71d58be..880f30d 100644
--- a/GPy/basic_gp.ipynb
+++ b/GPy/basic_gp.ipynb
@@ -8,6 +8,24 @@
     "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/SheffieldML/notebook/blob/main/GPy/basic_gp.ipynb)\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
diff --git a/GPy/basic_kernels.ipynb b/GPy/basic_kernels.ipynb
index e2a4a8d..e3147af 100644
--- a/GPy/basic_kernels.ipynb
+++ b/GPy/basic_kernels.ipynb
@@ -8,6 +8,24 @@
     "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/SheffieldML/notebook/blob/main/GPy/basic_kernels.ipynb)\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -15,10 +33,9 @@
    "outputs": [],
    "source": [
     "#configure plotting\n",
-    "%matplotlib inline\n",
     "%config InlineBackend.figure_format = 'svg'\n",
     "import matplotlib;matplotlib.rcParams['figure.figsize'] = (8,5)\n",
-    "from matplotlib import pyplot as plt"
+    "from matplotlib import pyplot as plt\n"
    ]
   },
   {
@@ -217,13 +234,6 @@
     "k = k1 * k2\n",
     "plot_sample(k)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/GPy/coregionalized_regression_tutorial.ipynb b/GPy/coregionalized_regression_tutorial.ipynb
index 45a15aa..829a005 100644
--- a/GPy/coregionalized_regression_tutorial.ipynb
+++ b/GPy/coregionalized_regression_tutorial.ipynb
@@ -44,6 +44,24 @@
     "The first thing to do is to set the plots to be interactive and to import GPy."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -54,10 +72,9 @@
    },
    "outputs": [],
    "source": [
-    "%pylab inline\n",
-    "import pylab as pb\n",
-    "pylab.ion()\n",
-    "import GPy"
+    "import numpy as np\n",
+    "from matplotlib import pyplot as plt\n",
+    "import GPy\n"
    ]
   },
   {
@@ -112,7 +129,7 @@
    "outputs": [],
    "source": [
     "xlim = (0,100); ylim = (0,50)\n",
-    "fig = pb.figure(figsize=(12,8))\n",
+    "fig = plt.figure(figsize=(12,8))\n",
     "ax1 = fig.add_subplot(211)\n",
     "ax1.set_xlim(xlim)\n",
     "ax1.set_title('Output 1')\n",
@@ -124,7 +141,7 @@
     "ax2.set_title('Output 2')\n",
     "ax2.plot(X2[:,:1],Y2,'kx',mew=1.5,label='Train set')\n",
     "ax2.plot(Xt2[:,:1],Yt2,'rx',mew=1.5,label='Test set')\n",
-    "ax2.legend()"
+    "ax2.legend()\n"
    ]
   },
   {
@@ -145,7 +162,7 @@
    "outputs": [],
    "source": [
     "def plot_2outputs(m,xlim,ylim):\n",
-    "    fig = pb.figure(figsize=(12,8))\n",
+    "    fig = plt.figure(figsize=(12,8))\n",
     "    #Output 1\n",
     "    ax1 = fig.add_subplot(211)\n",
     "    ax1.set_xlim(xlim)\n",
@@ -157,7 +174,7 @@
     "    ax2.set_xlim(xlim)\n",
     "    ax2.set_title('Output 2')\n",
     "    m.plot(plot_limits=xlim,fixed_inputs=[(1,1)],which_data_rows=slice(100,200),ax=ax2)\n",
-    "    ax2.plot(Xt2[:,:1],Yt2,'rx',mew=1.5)"
+    "    ax2.plot(Xt2[:,:1],Yt2,'rx',mew=1.5)\n"
    ]
   },
   {
@@ -343,7 +360,7 @@
     "m1.optimize()\n",
     "m2 = GPy.models.GPRegression(X2,Y2,kernel=K.copy())\n",
     "m2.optimize()\n",
-    "fig = pb.figure(figsize=(12,8))\n",
+    "fig = plt.figure(figsize=(12,8))\n",
     "#Output 1\n",
     "ax1 = fig.add_subplot(211)\n",
     "m1.plot(plot_limits=xlim,ax=ax1)\n",
@@ -353,7 +370,7 @@
     "ax2 = fig.add_subplot(212)\n",
     "m2.plot(plot_limits=xlim,ax=ax2)\n",
     "ax2.plot(Xt2[:,:1],Yt2,'rx',mew=1.5)\n",
-    "ax2.set_title('Output 2')"
+    "ax2.set_title('Output 2')\n"
    ]
   },
   {
diff --git a/GPy/heteroscedastic_regression.ipynb b/GPy/heteroscedastic_regression.ipynb
index 8d5172c..219705b 100644
--- a/GPy/heteroscedastic_regression.ipynb
+++ b/GPy/heteroscedastic_regression.ipynb
@@ -51,13 +51,40 @@
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "6",
    "metadata": {},
    "outputs": [],
    "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n",
+    "try:\n",
+    "    import ipywidgets\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", 'ipywidgets']\n",
+    "    )\n",
+    "    import ipywidgets\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matplotlib import pyplot as plt\n",
     "import numpy as np\n",
-    "import pylab as pb\n",
-    "import GPy \n",
-    "%pylab inline"
+    "import GPy\n"
    ]
   },
   {
@@ -76,9 +103,9 @@
     "def f(X):\n",
     "    return 10. + .1*X + 2*np.sin(X)/X\n",
     "\n",
-    "fig,ax = pb.subplots()\n",
+    "fig,ax = plt.subplots()\n",
     "ax.plot(np.linspace(-15,25),f(np.linspace(-10,20)),'r-')\n",
-    "ax.grid()"
+    "ax.grid()\n"
    ]
   },
   {
@@ -99,10 +126,10 @@
     "X = np.hstack([np.random.uniform(-1,1,1),X]) #Prepend a point between -1 and 1  (just for illustration)\n",
     "error = np.random.normal(0,.2,X.size)\n",
     "Y = f(X) + error\n",
-    "fig,ax = pb.subplots()\n",
+    "fig,ax = plt.subplots()\n",
     "ax.plot(np.linspace(-15,25),f(np.linspace(-10,20)),'r-')\n",
     "ax.plot(X,Y,'kx',mew=1.5)\n",
-    "ax.grid()"
+    "ax.grid()\n"
    ]
   },
   {
@@ -141,9 +168,9 @@
     "m.optimize()\n",
     "\n",
     "m.plot_f() #Show the predictive values of the GP.\n",
-    "pb.errorbar(X,Y,yerr=np.array(m.likelihood.flattened_parameters).flatten(),fmt=\"none\",ecolor='r',zorder=1)\n",
-    "pb.grid()\n",
-    "pb.plot(X,Y,'kx',mew=1.5)"
+    "plt.errorbar(X,Y,yerr=np.array(m.likelihood.flattened_parameters).flatten(),fmt=\"none\",ecolor='r',zorder=1)\n",
+    "plt.grid()\n",
+    "plt.plot(X,Y,'kx',mew=1.5)\n"
    ]
   },
   {
@@ -165,10 +192,10 @@
     "    m.optimize()\n",
     " \n",
     "    m.plot_f() \n",
-    "    pb.errorbar(X.flatten(),Y.flatten(),yerr=np.array(m.likelihood.flattened_parameters).flatten(),fmt=\"none\",ecolor='r',zorder=1)        \n",
-    "    pb.plot(X[1:],Y[1:],'kx',mew=1.5)\n",
-    "    pb.plot(X[:1],Y[:1],'ko',mew=.5)\n",
-    "    pb.grid()"
+    "    plt.errorbar(X.flatten(),Y.flatten(),yerr=np.array(m.likelihood.flattened_parameters).flatten(),fmt=\"none\",ecolor='r',zorder=1)        \n",
+    "    plt.plot(X[1:],Y[1:],'kx',mew=1.5)\n",
+    "    plt.plot(X[:1],Y[:1],'ko',mew=.5)\n",
+    "    plt.grid()\n"
    ]
   },
   {
@@ -214,14 +241,14 @@
     "m2.optimize()\n",
     "\n",
     "m1.plot_f()\n",
-    "pb.title('Homoscedastic model')\n",
+    "plt.title('Homoscedastic model')\n",
     "m2.plot_f()\n",
-    "pb.title('Heteroscedastic model')\n",
+    "plt.title('Heteroscedastic model')\n",
     "\n",
     "print(\"Kernel parameters (optimized) in the heteroscedastic model\")\n",
     "print(m1.kern)\n",
     "print(\"\\nKernel parameters (optimized) in the homoscedastic model\")\n",
-    "print(m2.kern)"
+    "print(m2.kern)\n"
    ]
   },
   {
@@ -263,7 +290,7 @@
     "m.plot_data(ax=ax)\n",
     "m.plot_errorbars_trainset(ax=ax, alpha=1)\n",
     "fig.tight_layout()\n",
-    "pb.grid()"
+    "plt.grid()\n"
    ]
   },
   {
@@ -301,15 +328,8 @@
     "m.plot_data(ax=ax)\n",
     "m.plot_errorbars_trainset(ax=ax, alpha=1)\n",
     "fig.tight_layout()\n",
-    "pb.grid()"
+    "plt.grid()\n"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/GPy/models_basic.ipynb b/GPy/models_basic.ipynb
index fe36b7a..e36e4e3 100644
--- a/GPy/models_basic.ipynb
+++ b/GPy/models_basic.ipynb
@@ -8,6 +8,24 @@
     "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/SheffieldML/notebook/blob/main/GPy/models_basic.ipynb)\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -16,11 +34,10 @@
    "source": [
     "#import necessary modules, set up the plotting\n",
     "import numpy as np\n",
-    "%matplotlib inline\n",
     "%config InlineBackend.figure_format = 'svg'\n",
     "import matplotlib;matplotlib.rcParams['figure.figsize'] = (8,6)\n",
     "from matplotlib import pyplot as plt\n",
-    "import GPy"
+    "import GPy\n"
    ]
   },
   {
@@ -413,13 +430,6 @@
     "fig = m.plot(plot_density=True)\n",
     "GPy.plotting.show(fig, filename='gpy_sparse_gp_example')"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/GPy/multiple outputs.ipynb b/GPy/multiple outputs.ipynb
index 0e11068..a7a3b52 100644
--- a/GPy/multiple outputs.ipynb	
+++ b/GPy/multiple outputs.ipynb	
@@ -18,18 +18,45 @@
     "Just as in the first lab, we firstly specify to include plots in the notebook and to import relevant libraries.\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n",
+    "try:\n",
+    "    import pods\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", 'pods']\n",
+    "    )\n",
+    "    import pods\n",
+    "\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "%matplotlib inline\n",
     "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
     "import GPy\n",
     "import pods\n",
-    "from IPython.display import display"
+    "from IPython.display import display\n"
    ]
   },
   {
diff --git a/GPy/optimizer-implementation.ipynb b/GPy/optimizer-implementation.ipynb
index 19bf82d..f1d108d 100644
--- a/GPy/optimizer-implementation.ipynb
+++ b/GPy/optimizer-implementation.ipynb
@@ -8,6 +8,24 @@
     "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/SheffieldML/notebook/blob/main/GPy/optimizer-implementation.ipynb)\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -15,8 +33,7 @@
    "outputs": [],
    "source": [
     "import GPy\n",
-    "from GPy.inference.optimization import Optimizer\n",
-    "%matplotlib inline"
+    "from GPy.inference.optimization import Optimizer\n"
    ]
   },
   {
diff --git a/GPy/sampling_hmc.ipynb b/GPy/sampling_hmc.ipynb
index ac726fd..f638e99 100644
--- a/GPy/sampling_hmc.ipynb
+++ b/GPy/sampling_hmc.ipynb
@@ -13,14 +13,41 @@
   {
    "cell_type": "code",
    "execution_count": null,
+   "id": "1",
    "metadata": {},
    "outputs": [],
    "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n",
+    "try:\n",
+    "    import seaborn\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", 'seaborn']\n",
+    "    )\n",
+    "    import seaborn\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matplotlib import pyplot as plt\n",
     "import numpy as np\n",
     "import GPy\n",
-    "%matplotlib inline\n",
-    "#%config InlineBackend.figure_format = 'svg'\n",
-    "from pylab import *"
+    "#%config InlineBackend.figure_format = 'svg'\n"
    ]
   },
   {
@@ -45,8 +72,8 @@
    "source": [
     "# Let's make some synthetic data\n",
     "x = np.linspace(0.,2*np.pi,100)[:,None]\n",
-    "y = -cos(x)+np.random.randn(*x.shape)*0.3+1\n",
-    "_ = plot(x,y,'.')"
+    "y = -np.cos(x)+np.random.randn(*x.shape)*0.3+1\n",
+    "_ = plt.plot(x,y,'.')\n"
    ]
   },
   {
@@ -109,7 +136,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot(s)"
+    "plt.plot(s)\n"
    ]
   },
   {
@@ -133,8 +160,8 @@
     "xs = np.linspace(xmin,xmax,100)\n",
     "for i in range(samples.shape[1]):\n",
     "    kernel = stats.gaussian_kde(samples[:,i])\n",
-    "    plot(xs,kernel(xs),label=labels[i])\n",
-    "_ = legend()"
+    "    plt.plot(xs,kernel(xs),label=labels[i])\n",
+    "_ = plt.legend()\n"
    ]
   },
   {
@@ -150,7 +177,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig = figure(figsize=(14,4))\n",
+    "fig = plt.figure(figsize=(14,4))\n",
     "ax = fig.add_subplot(131)\n",
     "_=ax.plot(samples[:,0],samples[:,1],'.')\n",
     "ax.set_xlabel(labels[0]); ax.set_ylabel(labels[1])\n",
@@ -159,7 +186,7 @@
     "ax.set_xlabel(labels[1]); ax.set_ylabel(labels[2])\n",
     "ax = fig.add_subplot(133)\n",
     "_=ax.plot(samples[:,0],samples[:,2],'.')\n",
-    "ax.set_xlabel(labels[0]); ax.set_ylabel(labels[2])"
+    "ax.set_xlabel(labels[0]); ax.set_ylabel(labels[2])\n"
    ]
   },
   {
@@ -229,7 +256,7 @@
     "x_new,mi = m.infer_newX(y_new)\n",
     "print(mi)\n",
     "m.plot()\n",
-    "plot(x_new,y_new,'or')"
+    "plt.plot(x_new,y_new,'or')\n"
    ]
   },
   {
@@ -262,7 +289,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "_ = plot(s_new[:,:])"
+    "_ = plt.plot(s_new[:,:])\n"
    ]
   },
   {
@@ -285,7 +312,7 @@
     "xs = np.linspace(xmin,xmax,100)\n",
     "for i in range(samples_new.shape[1]):\n",
     "    kernel = stats.gaussian_kde(samples_new[:,i])\n",
-    "    plot(xs,kernel(xs))"
+    "    plt.plot(xs,kernel(xs))\n"
    ]
   },
   {
@@ -345,7 +372,7 @@
     "#hmc = GPy.inference.optimization.HMC(m, stepsize=5e-1)\n",
     "#t = hmc.sample(m_iters=50000,hmc_iters=20)\n",
     "\n",
-    "_=plot(t)"
+    "_=plt.plot(t)\n"
    ]
   },
   {
@@ -356,7 +383,7 @@
    "source": [
     "print(t.mean(axis=0))\n",
     "print(t.std(axis=0))\n",
-    "_=hist(t[:,:2],50)"
+    "_=plt.hist(t[:,:2],50)\n"
    ]
   },
   {
@@ -408,13 +435,6 @@
    "source": [
     "df"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/GPy/sparse_gp_regression.ipynb b/GPy/sparse_gp_regression.ipynb
index 3da2c31..d458298 100644
--- a/GPy/sparse_gp_regression.ipynb
+++ b/GPy/sparse_gp_regression.ipynb
@@ -18,6 +18,24 @@
     "This example shows the variational compression effect of so-called 'sparse' Gaussian processes. In particular we show how using the variational free energy framework of [Titsias, 2009](http://jmlr.csail.mit.edu/proceedings/papers/v5/titsias09a/titsias09a.pdf) we can compress a Gaussian process fit. First we set up the notebook with a fixed random seed, and import GPy."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    import GPy\n",
+    "except ImportError:\n",
+    "    import subprocess\n",
+    "    import sys\n",
+    "    subprocess.check_call(\n",
+    "        [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"gpy[plotting]>=1.15.0\"]\n",
+    "    )\n",
+    "    import GPy\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -28,11 +46,10 @@
    },
    "outputs": [],
    "source": [
-    "%matplotlib inline\n",
     "%config InlineBackend.figure_format = 'svg'\n",
     "import GPy\n",
     "import numpy as np\n",
-    "np.random.seed(101)"
+    "np.random.seed(101)\n"
    ]
   },
   {
@@ -180,13 +197,6 @@
    "source": [
     "This time, we have enough inducing points and the fit resembles that of the GP. This is verified by the fact that the bound on the marginal likelihood is tight, which means that our variational approximation must be good (the difference between the bound and the true likelihood is the Kullback Leibler divergence between the approximation and the truth). "
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/backlog/features/2026-10-24_gpy-tutorials-colab-gpy-install.md b/backlog/features/2026-10-24_gpy-tutorials-colab-gpy-install.md
new file mode 100644
index 0000000..30f3109
--- /dev/null
+++ b/backlog/features/2026-10-24_gpy-tutorials-colab-gpy-install.md
@@ -0,0 +1,49 @@
+---
+id: "2026-10-24_gpy-tutorials-colab-gpy-install"
+title: "Install GPy from tutorial notebooks when it is missing"
+status: "Completed"
+priority: "Medium"
+created: "2026-10-24"
+last_updated: "2026-10-24"
+category: "features"
+related_cips: ["0006"]
+owner: "Neil Lawrence"
+dependencies: []
+tags:
+- backlog
+- gpy
+- colab
+---
+
+# Task: Colab GPy install check
+
+> **Note**: Backlog tasks are DOING the work defined in CIPs (HOW).
+> Use `related_cips` to link to CIPs. Don't link directly to requirements (bottom-up pattern).
+
+## Description
+
+Insert a first code cell in each executable `GPy/*.ipynb` that imports GPy if present, otherwise `pip install gpy[plotting]>=1.15.0`. Add the same pattern for `pods`, `seaborn`, and `ipywidgets` only where those notebooks already import them. Skip `index.ipynb` and `config.ipynb`. Do not auto-install climin.
+
+## Acceptance Criteria
+
+- [x] Canonical ensure-GPy cell is the first code cell in each executable GPy tutorial
+- [x] Extra ensure lines for pods / seaborn / ipywidgets where needed
+- [x] CIP-0003 execute stays green (import-hit path; no pip in CI)
+
+## Implementation Notes
+
+Use `subprocess` + `sys.executable`, not `%pip`. Do not upgrade an already-importable GPy.
+
+## Related
+
+- CIP: 0006
+- PRs:
+- Documentation: `cip/cip0006.md`
+
+## Progress Updates
+
+### 2026-10-24
+
+Task created as In Progress after CIP-0006 was accepted.
+
+Inserted the ensure-GPy cell (plus pods/seaborn/ipywidgets where needed). Local execute: 19 passed (import-hit path).
diff --git a/backlog/features/2026-10-24_gpy-tutorials-remove-pylab.md b/backlog/features/2026-10-24_gpy-tutorials-remove-pylab.md
new file mode 100644
index 0000000..75b2f51
--- /dev/null
+++ b/backlog/features/2026-10-24_gpy-tutorials-remove-pylab.md
@@ -0,0 +1,49 @@
+---
+id: "2026-10-24_gpy-tutorials-remove-pylab"
+title: "Replace pylab in GPy tutorials with matplotlib.pyplot"
+status: "Completed"
+priority: "Medium"
+created: "2026-10-24"
+last_updated: "2026-10-24"
+category: "features"
+related_cips: ["0005"]
+owner: "Neil Lawrence"
+dependencies: []
+tags:
+- backlog
+- gpy
+- plotting
+---
+
+# Task: Replace pylab in GPy tutorials
+
+> **Note**: Backlog tasks are DOING the work defined in CIPs (HOW).
+> Use `related_cips` to link to CIPs. Don't link directly to requirements (bottom-up pattern).
+
+## Description
+
+Remove `pylab`, `%pylab inline`, and `%matplotlib inline` from `GPy/*.ipynb`. Use `matplotlib.pyplot as plt` and qualify star-imported names in `sampling_hmc.ipynb`.
+
+## Acceptance Criteria
+
+- [x] No `pylab` or `%matplotlib` in `GPy/*.ipynb`
+- [x] `sampling_hmc.ipynb` uses `np` / `plt` instead of `from pylab import *`
+- [x] CIP-0003 execute still passes
+
+## Implementation Notes
+
+Do not add `%matplotlib inline`. GPy `plotting_library().figure` / `show` calls stay.
+
+## Related
+
+- CIP: 0005
+- PRs:
+- Documentation: `cip/cip0005.md`
+
+## Progress Updates
+
+### 2026-10-24
+
+Task created as In Progress after CIP-0005 was accepted.
+
+Replaced pylab with pyplot, removed matplotlib inline magics, and added `import numpy as np` in the coregionalized tutorial (it previously came from `%pylab`). Local execute: 19 passed.
diff --git a/cip/cip0005.md b/cip/cip0005.md
new file mode 100644
index 0000000..41878dd
--- /dev/null
+++ b/cip/cip0005.md
@@ -0,0 +1,122 @@
+---
+author: Neil Lawrence
+created: "2026-10-24"
+id: "0005"
+last_updated: "2026-10-24"
+status: "Closed"
+compressed: false
+related_requirements: []
+related_cips:
+- "0001"
+- "0002"
+- "0003"
+tags:
+- cip
+- gpy
+- plotting
+- matplotlib
+title: "Replace pylab in GPy tutorials with matplotlib.pyplot"
+---
+
+# CIP-0005: Replace pylab in GPy tutorials with matplotlib.pyplot
+
+> **Note**: CIPs describe HOW to achieve requirements (WHAT).
+> Use `related_requirements` to link to the requirements this CIP implements.
+
+## Status
+
+- [x] Proposed - Initial idea documented
+- [x] Accepted - Approved, ready to start work
+- [x] In Progress - Actively being implemented
+- [x] Implemented - Work complete, awaiting verification
+- [x] Closed - Verified and complete
+- [ ] Rejected - Will not be implemented
+- [ ] Deferred - Postponed
+
+## Summary
+
+Stop using Matplotlib's deprecated `pylab` interface in `GPy/*.ipynb`. Import
+`matplotlib.pyplot as plt` (and `numpy as np` where needed), delete
+`%pylab inline` (do **not** replace it with `%matplotlib inline`), and
+qualify names that today come from `from pylab import *`. Also delete
+existing `%matplotlib inline` magics: Jupyter and Colab show figures
+without that line.
+
+## Motivation
+
+GPy 1.15.0 already moved the **library** plotting path from `pylab` to
+`pyplot`. Four tutorials still teach the old notebook style:
+
+| Notebook | Current pylab use |
+|----------|-------------------|
+| `Poisson regression tutorial.ipynb` | `%pylab inline`, `import pylab as pb`, mixed `pb.*` / `plt.ylim` |
+| `heteroscedastic_regression.ipynb` | `import pylab as pb`, `%pylab inline`, mixed `pb.*` / `plt.subplots` |
+| `coregionalized_regression_tutorial.ipynb` | `%pylab inline`, `import pylab as pb`, `pylab.ion()` |
+| `sampling_hmc.ipynb` | `from pylab import *`, then bare `cos`, `plot`, … |
+
+`pylab` combines NumPy and pyplot into one namespace. Matplotlib has
+discouraged that for years; `%pylab` is an IPython compatibility shim. The
+other GPy tutorials already use `matplotlib.pyplot as plt`.
+
+## Detailed Description
+
+**Scope:** `GPy/*.ipynb` in this repository. Not `lab_classes/`, `compbio/`,
+`GPy-phil/`, and not the GPy library (that change shipped in 1.15.0).
+
+**House style:**
+
+```python
+import numpy as np
+from matplotlib import pyplot as plt
+```
+
+No IPython matplotlib magics. Current Jupyter and Colab display figures
+without `%matplotlib inline`; `%pylab inline` is the old combined import.
+
+Replacements:
+
+1. Delete `%pylab inline`. Do not write `%matplotlib inline` in its place.
+2. Delete existing `%matplotlib inline` in `GPy/*.ipynb` (most tutorials still have it).
+3. `import pylab as pb` / `import pylab` → `from matplotlib import pyplot as plt` (or `import matplotlib.pyplot as plt`). Rename `pb.` to `plt.`.
+4. `pylab.ion()` → drop.
+5. `from pylab import *` in `sampling_hmc.ipynb` → explicit `np` / `plt` names (`np.cos`, `plt.plot`, …). Do not leave a `pb` alias; the rest of the tree uses `plt`.
+6. Keep GPy's own plotting API (`GPy.plotting.change_plotting_library('matplotlib')`, `GPy.plotting.show`, `plotting_library().figure`) as it is. That is not `pylab`. Rewriting those calls to raw `plt.subplots` is out of scope.
+
+After the edits, `GPy/` must not contain `pylab`, `%pylab`, `from pylab`, or
+`%matplotlib`.
+
+## Implementation Plan
+
+1. Edit the four pylab notebooks; grep `GPy/` for `pylab`.
+2. Strip `%matplotlib inline` from every `GPy/*.ipynb` that still has it.
+3. In `sampling_hmc.ipynb`, convert star-imported names to `np` / `plt`.
+4. Run CIP-0003 (`scripts/run_gpy_tutorials.sh`) on the pylab four, then the full `GPy/` set.
+
+## Backward Compatibility
+
+Readers still see the same figures. Cells that relied on `%pylab` injecting
+`array` / `linspace` into the notebook namespace will need `np.` prefixes
+where they do not already import NumPy (they do, except via the pylab star
+import in `sampling_hmc.ipynb`).
+
+## Testing Strategy
+
+- Grep: no `pylab` or `%matplotlib` in `GPy/*.ipynb`.
+- CIP-0003 `execute` green against `gpy[plotting]==1.15.0`.
+
+## Related Requirements
+
+None yet.
+
+## Implementation Status
+
+- [x] Replace pylab in the four notebooks
+- [x] Remove `%matplotlib inline` from `GPy/*.ipynb`
+- [x] Confirm `sampling_hmc.ipynb` has no star-imported pylab names
+- [x] CIP-0003 execute green
+
+## References
+
+- CIP-0002 (Chart Studio → matplotlib library choice)
+- GPy 1.15.0 changelog: plotting `pylab` → `pyplot`
+- Matplotlib: [pylab](https://matplotlib.org/stable/api/index.html#module-pylab) is a historical namespace
diff --git a/cip/cip0006.md b/cip/cip0006.md
new file mode 100644
index 0000000..7b42993
--- /dev/null
+++ b/cip/cip0006.md
@@ -0,0 +1,152 @@
+---
+author: Neil Lawrence
+created: "2026-10-24"
+id: "0006"
+last_updated: "2026-10-24"
+status: "Closed"
+compressed: false
+related_requirements: []
+related_cips:
+- "0003"
+- "0004"
+- "0005"
+tags:
+- cip
+- gpy
+- colab
+- notebooks
+title: "Install GPy from the tutorial notebooks when it is missing"
+---
+
+# CIP-0006: Install GPy from the tutorial notebooks when it is missing
+
+> **Note**: CIPs describe HOW to achieve requirements (WHAT).
+> Use `related_requirements` to link to the requirements this CIP implements.
+
+## Status
+
+- [x] Proposed - Initial idea documented
+- [x] Accepted - Approved, ready to start work
+- [x] In Progress - Actively being implemented
+- [x] Implemented - Work complete, awaiting verification
+- [x] Closed - Verified and complete
+- [ ] Rejected - Will not be implemented
+- [ ] Deferred - Postponed
+
+## Summary
+
+Each executable GPy tutorial gets a first code cell that imports GPy if it is
+already present, and otherwise installs `gpy[plotting]>=1.15.0` (the floor
+these tutorials target) then imports it. Colab is the main case: CIP-0004
+badges open a notebook in an empty runtime.
+
+## Motivation
+
+Open-in-Colab badges land on a machine with NumPy and Matplotlib, not GPy.
+The first `import GPy` fails unless the reader knows to `%pip install gpy`.
+Local clones and CIP-0003 already have GPy, so the same cell must be a
+no-op when the import works.
+
+## Detailed Description
+
+**Scope:** executable `GPy/*.ipynb` (every tutorial with a code cell that
+imports GPy). Skip `index.ipynb` and `config.ipynb` (markdown only). Not
+`lab_classes/`, `compbio/`, `GPy-phil/`.
+
+**Canonical cell** (identical text in every in-scope notebook, first code
+cell, after the Colab badge markdown):
+
+```python
+try:
+    import GPy
+except ImportError:
+    import subprocess
+    import sys
+    subprocess.check_call(
+        [sys.executable, "-m", "pip", "install", "-q", "gpy[plotting]>=1.15.0"]
+    )
+    import GPy
+```
+
+Rationale:
+
+- **Check then install**, not an unconditional `%pip`. CIP-0003 already
+  installed GPy; a always-on pip would hit the network on every notebook.
+- **`subprocess` + `sys.executable`**, not `%pip`, so nbclient and Colab
+  both run it. `%pip` is Colab-friendly but is a Jupyter magic.
+- **`gpy[plotting]>=1.15.0`**: tutorials need the 1.15 matplotlib path;
+  Colab should get current PyPI, not be stuck on 1.14.2. Do not pin `==`
+  in the notebook (that fights CIP-0003's `execute-latest` canary).
+- **Do not upgrade** an already-importable GPy. A local 1.14.2 editable
+  install must keep working for library development; only a missing
+  module triggers pip.
+- **Do not detect Colab only.** Missing GPy in any Jupyter session gets
+  the same install. If import fails for a reason other than a missing
+  package, the `ImportError` still surfaces after pip (second `import GPy`).
+
+Place this cell **before** other imports that assume GPy (`GPy.plotting`,
+climin guards, `import pylab` / `pyplot`). Existing `import GPy` lines in
+later cells can stay (re-import is cheap) or be dropped if they become
+redundant; prefer leaving them so a reader who runs a mid-notebook cell
+still has `GPy` in that cell's story.
+
+**Per-notebook extras** (same try/install pattern, only where that notebook
+already depends on them):
+
+| Extra | Notebooks |
+|-------|-----------|
+| `pods` | `GPyCrashCourse.ipynb`, `multiple outputs.ipynb` (and any other `import pods`) |
+| `seaborn` | `sampling_hmc.ipynb` |
+| `ipywidgets` | `heteroscedastic_regression.ipynb` |
+| `climin` | keep the existing ImportError skip in `SVI.ipynb` and `optimizer-implementation.ipynb`; do **not** auto-install climin (optional extra, Python 3.11 `imp` constraint). Colab on 3.12 would install a package that cannot import. |
+
+East Timor in `GPyCrashCourse` stays on `datasets_mirror` (CIP-0003). Do not
+wait for a pods release for that file.
+
+**GPy folder vs `import GPy`:** the CIP-0003 runner uses the notebook
+directory as cwd so the `GPy/` tutorial folder does not shadow the library.
+This bootstrap cell does not change that. Colab downloads a single file, so
+there is no shadowing there.
+
+## Implementation Plan
+
+1. Insert the canonical cell as the first code cell in each in-scope
+   notebook (or prepend it to the existing first code cell if that cell is
+   only imports).
+2. Add the extra try/install lines only in notebooks that import those
+   packages, in the same cell or the cell that first needs them.
+3. Leave `index.ipynb` and `config.ipynb` without a pip cell.
+4. Run CIP-0003: every notebook must still pass with GPy already present
+   (the `except` branch must not run in CI). Optionally execute one
+   notebook in a venv without GPy to confirm the install branch.
+
+## Backward Compatibility
+
+Readers with GPy installed see one extra successful import. Readers without
+GPy need network access on first run (Colab has it). The cell writes into
+the current Python environment; that is the point of the Colab path.
+
+## Testing Strategy
+
+- CIP-0003 `execute` and `execute-latest` stay green (import-hit path).
+- A throwaway venv without GPy: run one short notebook (`config` has no
+  code; use `basic_kernels.ipynb` or `index` skip) and confirm pip then
+  `import GPy` succeeds.
+- Colab: open one badge after merge, run the first cell, confirm GPy
+  installs and the next cell imports.
+
+## Related Requirements
+
+None yet.
+
+## Implementation Status
+
+- [x] Canonical ensure-GPy cell in each executable `GPy/*.ipynb`
+- [x] Extra ensure lines for pods / seaborn / ipywidgets where needed
+- [x] CIP-0003 green on the import-hit path
+
+## References
+
+- CIP-0004 Colab badges
+- CIP-0003 runner and GPy 1.15 pin
+- GPy 1.15.0: `pip install 'gpy[plotting]>=1.15.0'`