{ "metadata": { "name": "", "signature": "sha256:0a44ee41e4ba70b93c51a1ba3c59e49ed69bb33ce5be371bb7175befa45c5bbf" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "DRAFT\n", "==\n", "# Differential Equations and Gaussian Processes\n", "\n", "# Gaussian Process Summer School, Melbourne, Australia\n", "### 25th-27th February 2015\n", "\n", "### Cristian Guarnizo Lemus, Mu Niu and Neil D. Lawrence" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gaussian process models can deal with linear operations and retain their nature as a Gaussian process. One example of this type of linear operation is a convolution. Such convolutions are often the result of driving a differential equaition with a function or setting that function as the initial condition. In this notebook we demonstrate the combination of Gaussian processes with differential equations, an approach that we refer to as *latent force models*. We consider a second order dynamical system and a spatial differential equation\n", "\n", "## Simulation Study\n", "\n", "First, we generate some data points from a second order differential equations. We start by importing some libraries" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "from scipy import integrate\n", "from IPython.display import display\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we define a function that allow us to simulate a dynamical system" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def calc_deri(yvec, time, damper, spring, s):\n", " return (yvec[1], 1.*s - damper*yvec[1] - spring*yvec[0])" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the previous function, We simulate the dynamical system regarding the following equation:\n", "$$\\frac{\\text{d}^2 y_d(t)}{\\text{d} t^2} + C_d \\frac{\\text{d} y_d(t)}{\\text{d} t} + B_d y_d(t) = \\sum_{q=1}^{Q} S_{d,q}u_q(t)$$\n", "In the function coded before, u(t) is the unit step." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# System parameters:\n", "# Spring constants\n", "B=[1.,3.]\n", "# Damper Constants\n", "C=[4.,1.]\n", "# Sensitivities\n", "S=[[1.],[1.]]\n", "# Numer of inputs\n", "q=1\n", "# Number of outputs\n", "d=2\n", "# Initialization of output values\n", "Nd = 100\n", "ti = 0.\n", "tf = 10.\n", "t = np.linspace(ti, tf, Nd)\n", "start = 0\n", "Y = np.empty([d*Nd, 1])\n", "T = np.empty([d*Nd, 1])\n", "index = np.empty(Y.shape, dtype = np.int16)\n", "for i in range(d):\n", " T[start:start+Nd, 0] = t\n", " yarr = integrate.odeint(calc_deri, (0, 0), t, args=(C[i], B[i], S[i][0]))\n", " Y[start:start+Nd, 0] = yarr[:, 0] + 0.01*np.random.randn(yarr.shape[0])\n", " index[start:start+Nd, 0] = i*np.ones((Nd,), dtype = np.int16)\n", " start = start + Nd\n", "\n", "X = np.hstack((T, index))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the minimum and maximum value of time inputs we generate the inducing inputs (here we make use of a trick to differentiate between ouput and latent time inputs)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "Mq = 50\n", "indexz = np.empty([Mq*q, 1], dtype = np.int16)\n", "z = np.empty([Mq*q, 1])\n", "start = 0\n", "zi = np.linspace(ti, tf, Mq)\n", "for i in range(q):\n", " z[start:start+Mq] = zi.reshape((zi.size, 1))\n", " indexz[start:start+Mq] = i*np.ones((zi.size, 1), dtype = np.int16)\n", " start = start + Mq\n", "#Index number greater than the number of outputs are considered as latent time inputs\n", "indexz += d\n", "Z = np.hstack((z, indexz))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model identification using GPy's sparse inference\n", "\n", "Kernel is created using the second order differential equation covariance function. Then, the model is trained using GPy's SparseGPRegression model." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import GPy\n", "kern = GPy.kern.EQ_ODE2(2, output_dim=d, rank=q)\n", "model = GPy.models.SparseGPRegression(X=X, Y=Y, Z=Z, kernel=kern)\n", "model.inducing_inputs.fix()\n", "model.optimize('scg', max_iters = 100)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we review the outcome from the sparse inference." ] }, { "cell_type": "code", "collapsed": false, "input": [ "display(model)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\n", "Name : sparse gp mpi\n", "Log-likelihood : 586.791497527\n", "Number of Parameters : 108\n", "Parameters:\n", " sparse_gp_mpi. | Value | Constraint | Prior | Tied to\n", " \u001b[1minducing inputs \u001b[0;0m | (50, 2) | fixed | | \n", " \u001b[1meq_ode2.lengthscale \u001b[0;0m | 1.61296828862 | +ve | | \n", " \u001b[1meq_ode2.C \u001b[0;0m | (2,) | +ve | | \n", " \u001b[1meq_ode2.B \u001b[0;0m | (2,) | +ve | | \n", " \u001b[1meq_ode2.W \u001b[0;0m | (2, 1) | | | \n", " \u001b[1mGaussian_noise.variance\u001b[0;0m | 0.000112833006706 | +ve | | \n" ] } ], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additionally we print the parameters." ] }, { "cell_type": "code", "collapsed": false, "input": [ "display(model.kern.B)\n", "display(model.kern.C)\n", "display(model.kern.W)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " Index | \u001b[1msparse_gp_mpi.eq_ode2.B\u001b[0;0m | Constraint | Prior | Tied to\n", " [0] | 0.53967786 | +ve | | N/A \n", " [1] | 2.3962038 | +ve | | N/A Index | \u001b[1msparse_gp_mpi.eq_ode2.C\u001b[0;0m | Constraint | Prior | Tied to\n", " [0] | 2.0954422 | +ve | | N/A \n", " [1] | 1.0402882 | +ve | | N/A Index | \u001b[1msparse_gp_mpi.eq_ode2.W\u001b[0;0m | Constraint | Prior | Tied to\n", " [0 0] | 0.35545259 | | | N/A \n", " [1 0] | 0.53738722 | | | N/A \n" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we perform predictions from the training data, and compare them with the real output values." ] }, { "cell_type": "code", "collapsed": false, "input": [ "Yp_mean = m.predict(X)\n", "Yp = Yp_mean[0][:, 0]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "plt.plot(X[:100, 0], Yp[:100], X[:100, 0], Y[:100, 0])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 9, "text": [ "[,\n", " ]" ] }, { "metadata": {}, "output_type": "display_data", "svg": [ "\n", "\n", "\n", "\n" ], "text": [ "" ] } ], "prompt_number": 9 }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(X[100:, 0], Yp[100:], X[100:, 0], Y[100:, 0])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 10, "text": [ "[,\n", " ]" ] }, { "metadata": {}, "output_type": "display_data", "svg": [ "\n", "\n", "\n", "\n" ], "text": [ "" ] } ], "prompt_number": 10 } ], "metadata": {} } ] }