![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp-optimize-eigen.png\")
\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/two_point_sample012.svg\")
\n",
"\n",
"Figure: The joint Gaussian over $\\mappingFunction_1$ and\n",
"$\\mappingFunction_2$ along with the conditional distribution of\n",
"$\\mappingFunction_2$ given $\\mappingFunction_1$\n",
"\n",
"- The single contour of the Gaussian density represents the\n",
" joint distribution,\n",
" $p(\\mappingFunction_1, \\mappingFunction_2)$\n",
"\n",
". . .\n",
"\n",
"- We observe that $\\mappingFunction_1=?$\n",
"\n",
". . .\n",
"\n",
"- Conditional density:\n",
" $p(\\mappingFunction_2|\\mappingFunction_1=?)$\n",
"\n",
"- Prediction of $\\mappingFunction_2$ from $\\mappingFunction_1$\n",
" requires *conditional density*.\n",
"\n",
"- Conditional density is *also* Gaussian. $$\n",
" p(\\mappingFunction_2|\\mappingFunction_1) = {\\mathcal{N}\\left(\\mappingFunction_2|\\frac{\\kernelScalar_{1, 2}}{\\kernelScalar_{1, 1}}\\mappingFunction_1,\\kernelScalar_{2, 2} - \\frac{\\kernelScalar_{1,2}^2}{\\kernelScalar_{1,1}}\\right)}\n",
" $$ where covariance of joint density is given by $$\n",
" \\kernelMatrix= \\begin{bmatrix} \\kernelScalar_{1, 1} & \\kernelScalar_{1, 2}\\\\ \\kernelScalar_{2, 1} & \\kernelScalar_{2, 2}\\end{bmatrix}\n",
" $$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pods\n",
"from ipywidgets import IntSlider"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"pods.notebook.display_plots('two_point_sample{sample:0>3}.svg', \n",
" '../slides/diagrams/gp', \n",
" sample=IntSlider(13, 13, 17, 1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/two_point_sample013.svg\")
\n",
"\n",
"Figure: Sample from the joint Gaussian model, points indexed by 1 and\n",
"8 highlighted.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/two_point_sample017.svg\")
\n",
"\n",
"Figure: The joint Gaussian over $\\mappingFunction_1$ and\n",
"$\\mappingFunction_8$ along with the conditional distribution of\n",
"$\\mappingFunction_8$ given $\\mappingFunction_1$\n",
"\n",
"- The single contour of the Gaussian density represents the\n",
" joint distribution,\n",
" $p(\\mappingFunction_1, \\mappingFunction_8)$\n",
"\n",
". . .\n",
"\n",
"- We observe a value for\n",
" $\\mappingFunction_1=-?$\n",
"\n",
". . .\n",
"\n",
"- Conditional density:\n",
" $p(\\mappingFunction_5|\\mappingFunction_1=?)$.\n",
"\n",
"- Prediction of $\\mappingFunctionVector_*$ from\n",
" $\\mappingFunctionVector$ requires multivariate *conditional\n",
" density*.\n",
"\n",
"- Multivariate conditional density is *also* Gaussian. ![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/neal-infinite-priors.png\")
\n",
"\n",
"Figure: Page 37 of [Radford Neal's 1994\n",
"thesis](http://www.cs.toronto.edu/~radford/ftp/thesis.pdf)\n",
"\n",
"In loose terms, what Radford considers is what happens to the elements\n",
"of the covariance function, $$\n",
" \\begin{align*}\n",
" \\kernel_\\mappingFunction\\left(\\inputVector_i, \\inputVector_j\\right) & = \\alpha \\activationVector\\left(\\mappingMatrix_1, \\inputVector_i\\right)^\\top \\activationVector\\left(\\mappingMatrix_1, \\inputVector_j\\right)\\\\\n",
" & = \\alpha \\sum_k \\activationScalar\\left(\\mappingVector^{(1)}_k, \\inputVector_i\\right) \\activationScalar\\left(\\mappingVector^{(1)}_k, \\inputVector_j\\right)\n",
" \\end{align*}\n",
" $$ if instead of considering a finite number you sample infinitely\n",
"many of these activation functions, sampling parameters from a prior\n",
"density, $p(\\mappingVectorTwo)$, for each one, $$\n",
"\\kernel_\\mappingFunction\\left(\\inputVector_i, \\inputVector_j\\right) = \\alpha \\int \\activationScalar\\left(\\mappingVector^{(1)}, \\inputVector_i\\right) \\activationScalar\\left(\\mappingVector^{(1)}, \\inputVector_j\\right) p(\\mappingVector^{(1)}) \\text{d}\\mappingVector^{(1)}\n",
"$$ And that's not *only* for Gaussian $p(\\mappingVectorTwo)$. In fact\n",
"this result holds for a range of activations, and a range of prior\n",
"densities because of the *central limit theorem*.\n",
"\n",
"To write it in the form of a probabilistic program, as long as the\n",
"distribution for $\\phi_i$ implied by this short probabilistic program,\n",
"$$\n",
" \\begin{align*}\n",
" \\mappingVectorTwo & \\sim p(\\cdot)\\\\\n",
" \\phi_i & = \\activationScalar\\left(\\mappingVectorTwo, \\inputVector_i\\right), \n",
" \\end{align*}\n",
" $$ has finite variance, then the result of taking the number of hidden\n",
"units to infinity, with appropriate scaling, is also a Gaussian process.\n",
"\n",
"## Further Reading\n",
"\n",
"To understand this argument in more detail, I highly recommend reading\n",
"chapter 2 of Neal's thesis [@Neal:bayesian94], which remains easy to\n",
"read and clear today. Indeed, for readers interested in Bayesian neural\n",
"networks, both Raford Neal's and David MacKay's PhD thesis\n",
"[@MacKay:bayesian92] remain essential reading. Both theses embody a\n",
"clarity of thought, and an ability to weave together threads from\n",
"different fields that was the business of machine learning in the 1990s.\n",
"Radford and David were also pioneers in making their software widely\n",
"available and publishing material on the web.\n",
"\n",
"## Gaussian Process\n",
"\n",
"In our \\refnotes{session on Bayesian regression}{bayesian-regression} we\n",
"sampled from the prior over paraemters. Through the properties of\n",
"multivariate Gaussian densities this prior over parameters implies a\n",
"particular density for our data observations, $\\dataVector$. In this\n",
"session we sampled directly from this distribution for our data,\n",
"avoiding the intermediate weight-space representation. This is the\n",
"approach taken by *Gaussian processes*. In a Gaussian process you\n",
"specify the *covariance function* directly, rather than *implicitly*\n",
"through a basis matrix and a prior over parameters. Gaussian processes\n",
"have the advantage that they can be *nonparametric*, which in simple\n",
"terms means that they can have *infinite* basis functions. In the\n",
"lectures we introduced the *exponentiated quadratic* covariance, also\n",
"known as the RBF or the Gaussian or the squared exponential covariance\n",
"function. This covariance function is specified by $$\n",
"\\kernelScalar(\\inputVector, \\inputVector^\\prime) = \\alpha \\exp\\left( -\\frac{\\left\\Vert \\inputVector-\\inputVector^\\prime\\right\\Vert^2}{2\\ell^2}\\right),\n",
"$$ where $\\left\\Vert\\inputVector - \\inputVector^\\prime\\right\\Vert^2$ is\n",
"the squared distance between the two input vectors $$\n",
"\\left\\Vert\\inputVector - \\inputVector^\\prime\\right\\Vert^2 = (\\inputVector - \\inputVector^\\prime)^\\top (\\inputVector - \\inputVector^\\prime) \n",
"$$ Let's build a covariance matrix based on this function. First we\n",
"define the form of the covariance function,"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%load -s eq_cov mlai.py"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use this to compute *directly* the covariance for\n",
"$\\mappingFunctionVector$ at the points given by `x_pred`. Let's define a\n",
"new function `K()` which does this,"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%load -s Kernel mlai.py"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can image the resulting covariance,"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"kernel = Kernel(function=eq_cov, variance=1., lengthscale=10.)\n",
"K = kernel.K(x_pred, x_pred)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To visualise the covariance between the points we can use the `imshow`\n",
"function in matplotlib."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(8,8))\n",
"im = ax.imshow(K, interpolation='none')\n",
"fig.colorbar(im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we can sample functions from the marginal likelihood."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize(8, 5))\n",
"for i in range(10):\n",
" y_sample = np.random.multivariate_normal(mean=np.zeros(x_pred.size), cov=K)\n",
" ax.plot(x_pred.flatten(), y_sample.flatten())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"**Moving Parameters** Have a play with the parameters for this\n",
"covariance function (the lengthscale and the variance) and see what\n",
"effects the parameters have on the types of functions you observe.\n",
"\n",
"### Write your answer to Exercise here"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Use this box for any code you need for the exercise\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bayesian Inference by Rejection Sampling \\[edit\\]\n",
"\n",
"One view of Bayesian inference is to assume we are given a mechanism for\n",
"generating samples, where we assume that mechanism is representing on\n",
"accurate view on the way we believe the world works.\n",
"\n",
"This mechanism is known as our *prior* belief.\n",
"\n",
"We combine our prior belief with our observations of the real world by\n",
"discarding all those samples that are inconsistent with our prior. The\n",
"*likelihood* defines mathematically what we mean by inconsistent with\n",
"the prior. The higher the noise level in the likelihood, the looser the\n",
"notion of consistent.\n",
"\n",
"The samples that remain are considered to be samples from the\n",
"*posterior*.\n",
"\n",
"This approach to Bayesian inference is closely related to two sampling\n",
"techniques known as *rejection sampling* and *importance sampling*. It\n",
"is realized in practice in an approach known as *approximate Bayesian\n",
"computation* (ABC) or likelihood-free inference.\n",
"\n",
"In practice, the algorithm is often too slow to be practical, because\n",
"most samples will be inconsistent with the data and as a result the\n",
"mechanism has to be operated many times to obtain a few posterior\n",
"samples.\n",
"\n",
"However, in the Gaussian process case, when the likelihood also assumes\n",
"Gaussian noise, we can operate this mechanism mathematically, and obtain\n",
"the posterior density *analytically*. This is the benefit of Gaussian\n",
"processes."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pods\n",
"from ipywidgets import IntSlider"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"pods.notebook.display_plots('gp_rejection_sample{sample:0>3}.png', \n",
" directory='../slides/diagrams/gp', \n",
" sample=IntSlider(1,1,5,1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp_rejection_sample003.png\")
\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp_rejection_sample004.png\")
\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp_rejection_sample005.png\")
\n",
"\n",
"Figure: One view of Bayesian inference is we have a machine for\n",
"generating samples (the *prior*), and we discard all samples\n",
"inconsistent with our data, leaving the samples of interest (the\n",
"*posterior*). The Gaussian process allows us to do this\n",
"analytically.\n",
"\n",
"## Gaussian Process\n",
"\n",
"The Gaussian process perspective takes the marginal likelihood of the\n",
"data to be a joint Gaussian density with a covariance given by\n",
"$\\kernelMatrix$. So the model likelihood is of the form, $$\n",
"p(\\dataVector|\\inputMatrix) =\n",
"\\frac{1}{(2\\pi)^{\\frac{\\numData}{2}}|\\kernelMatrix|^{\\frac{1}{2}}}\n",
"\\exp\\left(-\\frac{1}{2}\\dataVector^\\top \\left(\\kernelMatrix+\\dataStd^2\n",
"\\eye\\right)^{-1}\\dataVector\\right)\n",
"$$ where the input data, $\\inputMatrix$, influences the density through\n",
"the covariance matrix, $\\kernelMatrix$ whose elements are computed\n",
"through the covariance function,\n",
"$\\kernelScalar(\\inputVector, \\inputVector^\\prime)$.\n",
"\n",
"This means that the negative log likelihood (the objective function) is\n",
"given by, $$\n",
"\\errorFunction(\\boldsymbol{\\theta}) = \\frac{1}{2} \\log |\\kernelMatrix|\n",
"+ \\frac{1}{2} \\dataVector^\\top \\left(\\kernelMatrix +\n",
"\\dataStd^2\\eye\\right)^{-1}\\dataVector\n",
"$$ where the *parameters* of the model are also embedded in the\n",
"covariance function, they include the parameters of the kernel (such as\n",
"lengthscale and variance), and the noise variance, $\\dataStd^2$. Let's\n",
"create a class in python for storing these variables."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%load -s GP mlai.py"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Making Predictions\n",
"\n",
"We now have a probability density that represents functions. How do we\n",
"make predictions with this density? The density is known as a process\n",
"because it is *consistent*. By consistency, here, we mean that the model\n",
"makes predictions for $\\mappingFunctionVector$ that are unaffected by\n",
"future values of $\\mappingFunctionVector^*$ that are currently\n",
"unobserved (such as test points). If we think of\n",
"$\\mappingFunctionVector^*$ as test points, we can still write down a\n",
"joint probability density over the training observations,\n",
"$\\mappingFunctionVector$ and the test observations,\n",
"$\\mappingFunctionVector^*$. This joint probability density will be\n",
"Gaussian, with a covariance matrix given by our covariance function,\n",
"$\\kernelScalar(\\inputVector_i, \\inputVector_j)$. $$\n",
"\\begin{bmatrix}\\mappingFunctionVector \\\\ \\mappingFunctionVector^*\\end{bmatrix} \\sim \\gaussianSamp{\\zerosVector}{\\begin{bmatrix} \\kernelMatrix & \\kernelMatrix_\\ast \\\\\n",
"\\kernelMatrix_\\ast^\\top & \\kernelMatrix_{\\ast,\\ast}\\end{bmatrix}}\n",
"$$ where here $\\kernelMatrix$ is the covariance computed between all the\n",
"training points, $\\kernelMatrix_\\ast$ is the covariance matrix computed\n",
"between the training points and the test points and\n",
"$\\kernelMatrix_{\\ast,\\ast}$ is the covariance matrix computed betwen all\n",
"the tests points and themselves. To be clear, let's compute these now\n",
"for our example, using `x` and `y` for the training data (although `y`\n",
"doesn't enter the covariance) and `x_pred` as the test locations."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# set covariance function parameters\n",
"variance = 16.0\n",
"lengthscale = 8\n",
"# set noise variance\n",
"sigma2 = 0.05\n",
"\n",
"kernel = Kernel(eq_cov, variance=variance, lengthscale=lengthscale)\n",
"K = kernel.K(x, x)\n",
"K_star = kernel.K(x, x_pred)\n",
"K_starstar = kernel.K(x_pred, x_pred)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we use this structure to visualise the covariance between test data\n",
"and training data. This structure is how information is passed between\n",
"test and training data. Unlike the maximum likelihood formalisms we've\n",
"been considering so far, the structure expresses *correlation* between\n",
"our different data points. However, just like the\n",
"\\refnotes{naive Bayes approach}{naive-bayes} we now have a *joint\n",
"density* between some variables of interest. In particular we have the\n",
"joint density over\n",
"$p(\\mappingFunctionVector, \\mappingFunctionVector^*)$. The joint density\n",
"is *Gaussian* and *zero mean*. It is specified entirely by the\n",
"*covariance matrix*, $\\kernelMatrix$. That covariance matrix is, in\n",
"turn, defined by a covariance function. Now we will visualise the form\n",
"of that covariance in the form of the matrix, $$\n",
"\\begin{bmatrix} \\kernelMatrix & \\kernelMatrix_\\ast \\\\ \\kernelMatrix_\\ast^\\top\n",
"& \\kernelMatrix_{\\ast,\\ast}\\end{bmatrix}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(8,8))\n",
"im = ax.imshow(np.vstack([np.hstack([K, K_star]), np.hstack([K_star.T, K_starstar])]), interpolation='none')\n",
"# Add lines for separating training and test data\n",
"ax.axvline(x.shape[0]-1, color='w')\n",
"ax.axhline(x.shape[0]-1, color='w')\n",
"fig.colorbar(im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are four blocks to this color plot. The upper left block is the\n",
"covariance of the training data with itself, $\\kernelMatrix$. We see\n",
"some structure here due to the missing data from the first and second\n",
"world wars. Alongside this covariance (to the right and below) we see\n",
"the cross covariance between the training and the test data\n",
"($\\kernelMatrix_*$ and $\\kernelMatrix_*^\\top$). This is giving us the\n",
"covariation between our training and our test data. Finally the lower\n",
"right block The banded structure we now observe is because some of the\n",
"training points are near to some of the test points. This is how we\n",
"obtain 'communication' between our training data and our test data. If\n",
"there is no structure in $\\kernelMatrix_*$ then our belief about the\n",
"test data simply matches our prior.\n",
"\n",
"## Conditional Density\n",
"\n",
"Just as in naive Bayes, we first defined the joint density (although\n",
"there it was over both the labels and the inputs,\n",
"$p(\\dataVector, \\inputMatrix)$ and now we need to define *conditional*\n",
"distributions that answer particular questions of interest. In\n",
"particular we might be interested in finding out the values of the\n",
"function for the prediction function at the test data given those at the\n",
"training data, $p(\\mappingFunctionVector_*|\\mappingFunctionVector)$. Or\n",
"if we include noise in the training observations then we are interested\n",
"in the conditional density for the prediction function at the test\n",
"locations given the training observations,\n",
"$p(\\mappingFunctionVector^*|\\dataVector)$.\n",
"\n",
"As ever all the various questions we could ask about this density can be\n",
"answered using the *sum rule* and the *product rule*. For the\n",
"multivariate normal density the mathematics involved is that of *linear\n",
"algebra*, with a particular emphasis on the *partitioned inverse* or\n",
"[*block matrix\n",
"inverse*](http://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion),\n",
"but they are beyond the scope of this course, so you don't need to worry\n",
"about remembering them or rederiving them. We are simply writing them\n",
"here because it is this *conditional* density that is necessary for\n",
"making predictions.\n",
"\n",
"The conditional density is also a multivariate normal, $$\n",
"\\mappingFunctionVector^* | \\mappingFunctionVector \\sim \\gaussianSamp{\\meanVector_\\mappingFunction}{\\mathbf{C}_\\mappingFunction}\n",
"$$ with a mean given by $$\n",
"\\meanVector_\\mappingFunction = \\kernelMatrix_*^\\top \\left[\\kernelMatrix + \\dataStd^2\n",
"\\eye\\right]^{-1} \\dataVector\n",
"$$ and a covariance given by $$\n",
"\\mathbf{C}_\\mappingFunction\n",
"= \\kernelMatrix_{*,*} - \\kernelMatrix_*^\\top \\left[\\kernelMatrix + \\dataStd^2\n",
"\\eye\\right]^{-1} \\kernelMatrix_\\ast.\n",
"$$ Let's compute what those posterior predictions are for the olympic\n",
"marathon data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%load -s posterior_f mlai.py"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# attach the new method to class GP():\n",
"GP.posterior_f = posterior_f"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"model = GP(x, y, sigma2, exponentiated_quadratic, variance=variance, lengthscale=lengthscale)\n",
"mu_f, C_f = model.posterior_f(x_pred)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where for convenience we've defined\n",
"\n",
"$$\n",
"\\mathbf{A} = \\left[\\kernelMatrix + \\dataStd^2\\eye\\right]^{-1}\\kernelMatrix_*.\n",
"$$\n",
"\n",
"We can visualize the covariance of the *conditional*,"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(8,8))\n",
"im = ax.imshow(C_f, interpolation='none')\n",
"fig.colorbar(im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and we can plot the mean of the conditional"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot(x, y, 'rx')\n",
"plt.plot(x_pred, mu_f, 'b-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"as well as the associated error bars. These are given (similarly to the\n",
"Bayesian parametric model from the last lab) by the standard deviations\n",
"of the marginal posterior densities. The marginal posterior variances\n",
"are given by the diagonal elements of the posterior covariance,"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"var_f = np.diag(C_f)[:, None]\n",
"std_f = np.sqrt(var_f)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"They can be added to the underlying mean function to give the error\n",
"bars,"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot(x, y, 'rx')\n",
"plt.plot(x_pred, mu_f, 'b-')\n",
"plt.plot(x_pred, mu_f+2*std_f, 'b--')\n",
"plt.plot(x_pred, mu_f-2*std_f, 'b--')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This gives us a prediction from the Gaussian process. Remember machine\n",
"learning is $$\n",
"\\text{data} + \\text{model} \\rightarrow \\text{prediction}.\n",
"$$ Here our data is from the olympics, and our model is a Gaussian\n",
"process with two parameters. The assumptions about the world are encoded\n",
"entirely into our Gaussian process covariance. The GP covariance assumes\n",
"that the function is highly smooth, and that correlation falls off with\n",
"distance (scaled according to the length scale, $\\ell$). The model\n",
"sustains the uncertainty about the function, this means we see an\n",
"increase in the size of the error bars during periods like the 1st and\n",
"2nd World Wars when no olympic marathon was held.\n",
"\n",
"### Exercise\n",
"\n",
"Now try changing the parameters of the covariance function (and the\n",
"noise) to see how the predictions change.\n",
"\n",
"Then try sampling from this conditional density to see what your\n",
"predictions look like. What happens if you sample from the conditional\n",
"density in regions a long way into the future or the past? How does this\n",
"compare with the results from the polynomial model?\n",
"\n",
"### Write your answer to Exercise here"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Use this box for any code you need for the exercise\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Importance of the Covariance Function\n",
"\n",
"The covariance function encapsulates our assumptions about the data. The\n",
"equations for the distribution of the prediction function, given the\n",
"training observations, are highly sensitive to the covariation between\n",
"the test locations and the training locations as expressed by the matrix\n",
"$\\kernelMatrix_*$. We defined a matrix $\\mathbf{A}$ which allowed us to\n",
"express our conditional mean in the form, $$\n",
"\\meanVector_\\mappingFunction = \\mathbf{A}^\\top \\dataVector,\n",
"$$ where $\\dataVector$ were our *training observations*. In other words\n",
"our mean predictions are always a linear weighted combination of our\n",
"*training data*. The weights are given by computing the covariation\n",
"between the training and the test data ($\\kernelMatrix_*$) and scaling\n",
"it by the inverse covariance of the training data observations,\n",
"$\\left[\\kernelMatrix + \\dataStd^2 \\eye\\right]^{-1}$. This inverse is the\n",
"main computational object that needs to be resolved for a Gaussian\n",
"process. It has a computational burden which is $O(\\numData^3)$ and a\n",
"storage burden which is $O(\\numData^2)$. This makes working with\n",
"Gaussian processes computationally intensive for the situation where\n",
"$\\numData>10,000$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from IPython.lib.display import YouTubeVideo\n",
"YouTubeVideo('ewJ3AxKclOg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Improving the Numerics\n",
"\n",
"In practice we shouldn't be using matrix inverse directly to solve the\n",
"GP system. One more stable way is to compute the *Cholesky\n",
"decomposition* of the kernel matrix. The log determinant of the\n",
"covariance can also be derived from the Cholesky decomposition."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%load -s update_inverse mlai.py"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"GP.update_inverse = update_inverse"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Capacity Control\n",
"\n",
"Gaussian processes are sometimes seen as part of a wider family of\n",
"methods known as kernel methods. Kernel methods are also based around\n",
"covariance functions, but in the field they are known as Mercer kernels.\n",
"Mercer kernels have interpretations as inner products in potentially\n",
"infinite dimensional Hilbert spaces. This interpretation arises because,\n",
"if we take $\\alpha=1$, then the kernel can be expressed as $$\n",
"\\kernelMatrix = \\basisMatrix\\basisMatrix^\\top \n",
"$$ which imples the elements of the kernel are given by, $$\n",
"\\kernelScalar(\\inputVector, \\inputVector^\\prime) = \\basisVector(\\inputVector)^\\top \\basisVector(\\inputVector^\\prime).\n",
"$$ So we see that the kernel function is developed from an inner product\n",
"between the basis functions. Mercer's theorem tells us that any valid\n",
"*positive definite function* can be expressed as this inner product but\n",
"with the caveat that the inner product could be *infinite length*. This\n",
"idea has been used quite widely to *kernelize* algorithms that depend on\n",
"inner products. The kernel functions are equivalent to covariance\n",
"functions and they are parameterized accordingly. In the kernel modeling\n",
"community it is generally accepted that kernel parameter estimation is a\n",
"difficult problem and the normal solution is to cross validate to obtain\n",
"parameters. This can cause difficulties when a large number of kernel\n",
"parameters need to be estimated. In Gaussian process modelling kernel\n",
"parameter estimation (in the simplest case proceeds) by maximum\n",
"likelihood. This involves taking gradients of the likelihood with\n",
"respect to the parameters of the covariance function.\n",
"\n",
"## Gradients of the Likelihood\n",
"\n",
"The easiest conceptual way to obtain the gradients is a two step\n",
"process. The first step involves taking the gradient of the likelihood\n",
"with respect to the covariance function, the second step involves\n",
"considering the gradient of the covariance function with respect to its\n",
"parameters.\n",
"\n",
"## Overall Process Scale\n",
"\n",
"In general we won't be able to find parameters of the covariance\n",
"function through fixed point equations, we will need to do gradient\n",
"based optimization.\n",
"\n",
"## Capacity Control and Data Fit\n",
"\n",
"The objective function can be decomposed into two terms, a capacity\n",
"control term, and a data fit term. The capacity control term is the log\n",
"determinant of the covariance. The data fit term is the matrix inner\n",
"product between the data and the inverse covariance.\n",
"\n",
"Can we determine covariance parameters from the data?\n",
"\n",
"$$\n",
"\\gaussianDist{\\dataVector}{\\mathbf{0}}{\\kernelMatrix}=\\frac{1}{(2\\pi)^\\frac{\\numData}{2}{\\det{\\kernelMatrix}^{\\frac{1}{2}}}}{\\exp\\left(-\\frac{\\dataVector^{\\top}\\kernelMatrix^{-1}\\dataVector}{2}\\right)}\n",
"$$\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" \\gaussianDist{\\dataVector}{\\mathbf{0}}{\\kernelMatrix}=\\frac{1}{(2\\pi)^\\frac{\\numData}{2}{\\color{black} \\det{\\kernelMatrix}^{\\frac{1}{2}}}}{\\color{black}\\exp\\left(-\\frac{\\dataVector^{\\top}\\kernelMatrix^{-1}\\dataVector}{2}\\right)}\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" \\log \\gaussianDist{\\dataVector}{\\mathbf{0}}{\\kernelMatrix}=&{\\color{black}-\\frac{1}{2}\\log\\det{\\kernelMatrix}}{\\color{black}-\\frac{\\dataVector^{\\top}\\kernelMatrix^{-1}\\dataVector}{2}} \\\\ &-\\frac{\\numData}{2}\\log2\\pi\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"$$\n",
"\\errorFunction(\\parameterVector) = {\\color{black} \\frac{1}{2}\\log\\det{\\kernelMatrix}} + {\\color{black} \\frac{\\dataVector^{\\top}\\kernelMatrix^{-1}\\dataVector}{2}}\n",
"$$\n",
"\n",
"The parameters are *inside* the covariance function (matrix).\n",
"\\normalsize\n",
"$$\\kernelScalar_{i, j} = \\kernelScalar(\\inputVals_i, \\inputVals_j; \\parameterVector)$$\n",
"\n",
"\\Large\n",
"$$\\kernelMatrix = \\rotationMatrix \\eigenvalueMatrix^2 \\rotationMatrix^\\top$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"gpoptimizePlot1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"| \n",
"\n",
" | \n",
"\n", "$\\eigenvalueMatrix$ represents distance on axes. $\\rotationMatrix$ gives\n", "rotation.\n", " | \n", "
\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp-optimise006.svg\") \n",
" | \n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp-optimise010.svg\") \n",
" | \n",
"
\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp-optimise016.svg\") \n",
" | \n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp-optimise021.svg\") \n",
" | \n",
"
\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/kern/eq_covariance.svg\") \n",
" | \n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/kern/eq_covariance.gif\") \n",
" | \n",
"
| \n", "- Gold medal times for Olympic Marathon since 1896.\n", "- Marathons before 1924 didn't have a standardised distance.\n", "- Present results using pace per km.\n", "- In 1904 Marathon was badly organised leading to very slow times.\n", "\n", " | \n", "\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/Stephen_Kiprotich.jpg\") \n",
"Image from Wikimedia Commons | \n",
"
![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/datasets/olympic-marathon.svg\")
\n",
"\n",
"Figure: Olympic marathon pace times since 1892.\n",
"\n",
"Things to notice about the data include the outlier in 1904, in this\n",
"year, the olympics was in St Louis, USA. Organizational problems and\n",
"challenges with dust kicked up by the cars following the race meant that\n",
"participants got lost, and only very few participants completed.\n",
"\n",
"More recent years see more consistently quick marathons.\n",
"\n",
"## Alan Turing \\[edit\\]\n",
"\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/turing-times.gif\") \n",
" | \n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/turing-run.jpg\") \n",
" | \n",
"
![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/olympic-marathon-gp.svg\")
\n",
"\n",
"Figure: Gaussian process fit to the Olympic Marathon data. The error\n",
"bars are too large, perhaps due to the outlier from 1904.\n",
"\n",
"## Fit Quality\n",
"\n",
"In the fit we see that the error bars (coming mainly from the noise\n",
"variance) are quite large. This is likely due to the outlier point in\n",
"1904, ignoring that point we can see that a tighter fit is obtained. To\n",
"see this making a version of the model, `m_clean`, where that point is\n",
"removed."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x_clean=np.vstack((x[0:2, :], x[3:, :]))\n",
"y_clean=np.vstack((y[0:2, :], y[3:, :]))\n",
"\n",
"m_clean = GPy.models.GPRegression(x_clean,y_clean)\n",
"_ = m_clean.optimize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gene Expression Example \\[edit\\]\n",
"\n",
"We now consider an example in gene expression. Gene expression is the\n",
"measurement of mRNA levels expressed in cells. These mRNA levels show\n",
"which genes are 'switched on' and producing data. In the example we will\n",
"use a Gaussian process to determine whether a given gene is active, or\n",
"we are merely observing a noise response.\n",
"\n",
"## Della Gatta Gene Data \\[edit\\]\n",
"\n",
"- Given given expression levels in the form of a time series from\n",
" @DellaGatta:direct08."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pods"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=937)\n",
"\n",
"x = data['X']\n",
"y = data['Y']\n",
"\n",
"offset = y.mean()\n",
"scale = np.sqrt(y.var())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import teaching_plots as plot\n",
"import mlai"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"xlim = (-20,260)\n",
"ylim = (5, 7.5)\n",
"yhat = (y-offset)/scale\n",
"\n",
"fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n",
"_ = ax.plot(x, y, 'r.',markersize=10)\n",
"ax.set_xlabel('time/min', fontsize=20)\n",
"ax.set_ylabel('expression', fontsize=20)\n",
"ax.set_xlim(xlim)\n",
"ax.set_ylim(ylim)\n",
"\n",
"mlai.write_figure(figure=fig, \n",
" filename='../slides/diagrams/datasets/della-gatta-gene.svg', \n",
" transparent=True, \n",
" frameon=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/datasets/della-gatta-gene.svg\")
\n",
"\n",
"Figure: Gene expression levels over time for a gene from data\n",
"provided by @DellaGatta:direct08. We would like to understand whethere\n",
"there is signal in the data, or we are only observing noise.\n",
"\n",
"- Want to detect if a gene is expressed or not, fit a GP to each gene\n",
" @Kalaitzis:simple11.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/health/1471-2105-12-180_1.png\")
\n",
"\n",
"Figure: The example is taken from the paper \"A Simple Approach to\n",
"Ranking Differentially Expressed Gene Expression Time Courses through\n",
"Gaussian Process Regression.\" @Kalaitzis:simple11.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/della-gatta-gene-gp.svg\")
\n",
"\n",
"Figure: Result of the fit of the Gaussian process model with the time\n",
"scale parameter initialized to 50 minutes.\n",
"\n",
"Now we try a model initialized with a longer length scale."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m_full2 = GPy.models.GPRegression(x,yhat)\n",
"m_full2.kern.lengthscale=2000\n",
"_ = m_full2.optimize() # Optimize parameters of covariance function"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import teaching_plots as plot"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n",
"plot.model_output(m_full2, scale=scale, offset=offset, ax=ax, xlabel='time/min', ylabel='expression', fontsize=20, portion=0.2)\n",
"ax.set_xlim(xlim)\n",
"ax.set_ylim(ylim)\n",
"ax.set_title('log likelihood: {ll:.3}'.format(ll=m_full2.log_likelihood()), fontsize=20)\n",
"mlai.write_figure(figure=fig,\n",
" filename='../slides/diagrams/gp/della-gatta-gene-gp2.svg', \n",
" transparent=True, frameon=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/della-gatta-gene-gp2.svg\")
\n",
"\n",
"Figure: Result of the fit of the Gaussian process model with the time\n",
"scale parameter initialized to 2000 minutes.\n",
"\n",
"Now we try a model initialized with a lower noise."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m_full3 = GPy.models.GPRegression(x,yhat)\n",
"m_full3.kern.lengthscale=20\n",
"m_full3.likelihood.variance=0.001\n",
"_ = m_full3.optimize() # Optimize parameters of covariance function"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import teaching_plots as plot"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=plot.big_wide_figsize)\n",
"plot.model_output(m_full3, scale=scale, offset=offset, ax=ax, xlabel='time/min', ylabel='expression', fontsize=20, portion=0.2)\n",
"ax.set_xlim(xlim)\n",
"ax.set_ylim(ylim)\n",
"ax.set_title('log likelihood: {ll:.3}'.format(ll=m_full3.log_likelihood()), fontsize=20)\n",
"mlai.write_figure(figure=fig,\n",
" filename='../slides/diagrams/gp/della-gatta-gene-gp3.svg', \n",
" transparent=True, frameon=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/della-gatta-gene-gp3.svg\")
\n",
"\n",
"Figure: Result of the fit of the Gaussian process model with the\n",
"noise initialized low (standard deviation 0.1) and the time scale\n",
"parameter initialized to 20 minutes.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/multiple-optima000.svg\")
\n",
"\n",
"Figure: \n",
"\n",
"\n",
"## Example: Prediction of Malaria Incidence in Uganda \\[edit\\]\n",
"\n",
"[![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/people/2013_03_28_180606.JPG\")
]{style=\"text-align:right\"}\n",
"\n",
"As an example of using Gaussian process models within the full pipeline\n",
"from data to decsion, we'll consider the prediction of Malaria incidence\n",
"in Uganda. For the purposes of this study malaria reports come in two\n",
"forms, HMIS reports from health centres and Sentinel data, which is\n",
"curated by the WHO. There are limited sentinel sites and many HMIS\n",
"sites.\n",
"\n",
"The work is from Ricardo Andrade Pacheco's PhD thesis, completed in\n",
"collaboration with John Quinn and Martin Mubangizi\n",
"[@Andrade:consistent14; @Mubangizi:malaria14]. John and Martin were\n",
"initally from the AI-DEV group from the University of Makerere in\n",
"Kampala and more latterly they were based at UN Global Pulse in Kampala.\n",
"\n",
"Malaria data is spatial data. Uganda is split into districts, and health\n",
"reports can be found for each district. This suggests that models such\n",
"as conditional random fields could be used for spatial modelling, but\n",
"there are two complexities with this. First of all, occasionally\n",
"districts split into two. Secondly, sentinel sites are a specific\n",
"location within a district, such as Nagongera which is a sentinel site\n",
"based in the Tororo district.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/health/uganda-districts-2006.png\")
\n",
"\n",
"Figure: Ugandan districs. Data SRTM/NASA from\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/health/Kapchorwa_District_in_Uganda.svg\")
\n",
"\n",
"Figure: The Kapchorwa District, home district of Stephen\n",
"Kiprotich.\n",
"\n",
"Stephen Kiprotich, the 2012 gold medal winner from the London Olympics,\n",
"comes from Kapchorwa district, in eastern Uganda, near the border with\n",
"Kenya.\n",
"\n",
"The common standard for collecting health data on the African continent\n",
"is from the Health management information systems (HMIS). However, this\n",
"data suffers from missing values [@Gething:hmis06] and diagnosis of\n",
"diseases like typhoid and malaria may be confounded.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/health/Tororo_District_in_Uganda.svg\")
\n",
"\n",
"Figure: The Tororo district, where the sentinel site, Nagongera, is\n",
"located.\n",
"\n",
"[World Health Organization Sentinel Surveillance\n",
"systems](https://www.who.int/immunization/monitoring_surveillance/burden/vpd/surveillance_type/sentinel/en/)\n",
"are set up \"when high-quality data are needed about a particular disease\n",
"that cannot be obtained through a passive system\". Several sentinel\n",
"sites give accurate assessment of malaria disease levels in Uganda,\n",
"including a site in Nagongera.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/health/sentinel_nagongera.png\")
\n",
"\n",
"Figure: Sentinel and HMIS data along with rainfall and temperature\n",
"for the Nagongera sentinel station in the Tororo district.\n",
"\n",
"In collaboration with the AI Research Group at Makerere we chose to\n",
"investigate whether Gaussian process models could be used to assimilate\n",
"information from these two different sources of disease informaton.\n",
"Further, we were interested in whether local information on rainfall and\n",
"temperature could be used to improve malaria estimates.\n",
"\n",
"The aim of the project was to use WHO Sentinel sites, alongside rainfall\n",
"and temperature, to improve predictions from HMIS data of levels of\n",
"malaria.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/health/Mubende_District_in_Uganda.svg\")
\n",
"\n",
"Figure: The Mubende District.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/health/mubende.png\")
\n",
"\n",
"Figure: Prediction of malaria incidence in Mubende.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gpss/1157497_513423392066576_1845599035_n.jpg\")
\n",
"\n",
"Figure: The project arose out of the Gaussian process summer school\n",
"held at Makerere in Kampala in 2013. The school led, in turn, to the\n",
"Data Science Africa initiative.\n",
"\n",
"## Early Warning Systems\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/health/Kabarole_District_in_Uganda.svg\")
\n",
"\n",
"Figure: The Kabarole district in Uganda.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/health/kabarole.gif\")
\n",
"\n",
"Figure: Estimate of the current disease situation in the Kabarole\n",
"district over time. Estimate is constructed with a Gaussian process with\n",
"an additive covariance funciton.\n",
"\n",
"Health monitoring system for the Kabarole district. Here we have fitted\n",
"the reports with a Gaussian process with an additive covariance\n",
"function. It has two components, one is a long time scale component (in\n",
"red above) the other is a short time scale component (in blue).\n",
"\n",
"Monitoring proceeds by considering two aspects of the curve. Is the blue\n",
"line (the short term report signal) above the red (which represents the\n",
"long term trend? If so we have higher than expected reports. If this is\n",
"the case *and* the gradient is still positive (i.e. reports are going\n",
"up) we encode this with a *red* color. If it is the case and the\n",
"gradient of the blue line is negative (i.e. reports are going down) we\n",
"encode this with an *amber* color. Conversely, if the blue line is below\n",
"the red *and* decreasing, we color *green*. On the other hand if it is\n",
"below red but increasing, we color *yellow*.\n",
"\n",
"This gives us an early warning system for disease. Red is a bad\n",
"situation getting worse, amber is bad, but improving. Green is good and\n",
"getting better and yellow good but degrading.\n",
"\n",
"Finally, there is a gray region which represents when the scale of the\n",
"effect is small.\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/health/monitor.gif\")
\n",
"\n",
"Figure: The map of Ugandan districts with an overview of the Malaria\n",
"situation in each district.\n",
"\n",
"These colors can now be observed directly on a spatial map of the\n",
"districts to give an immediate impression of the current status of the\n",
"disease across the country.\n",
"\n",
"## Additive Covariance \\[edit\\]\n",
"\n",
"An additive covariance function is derived from considering the result\n",
"of summing two Gaussian processes together. If the first Gaussian\n",
"process is $g(\\cdot)$, governed by covariance\n",
"$\\kernelScalar_g(\\cdot, \\cdot)$ and the second process is $h(\\cdot)$,\n",
"governed by covariance $\\kernelScalar_h(\\cdot, \\cdot)$ then the combined\n",
"process $f(\\cdot) = g(\\cdot) + h(\\cdot)$ is govererned by a covariance\n",
"function, $$\n",
"\\kernelScalar_f(\\inputVector, \\inputVector^\\prime) = \\kernelScalar_g(\\inputVector, \\inputVector^\\prime) + \\kernelScalar_h(\\inputVector, \\inputVector^\\prime)\n",
"$$\n",
"\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/kern/add_covariance.svg\") \n",
" | \n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/kern/add_covariance.gif\") \n",
" | \n",
"
![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/ml/bialik-fridaythe13th-1.png\")
\n",
"\n",
"Figure: This is a retrospective analysis of US births by Aki Vehtari.\n",
"The challenges of forecasting. Even with seasonal and weekly effects\n",
"removed there are significant effects on holidays, weekends, etc.\n",
"\n",
"There's a nice analysis of US birth rates by Gaussian processes with\n",
"additive covariances in @Gelman:bayesian13. A combination of covariance\n",
"functions are used to take account of weekly and yearly trends. The\n",
"analysis is summarized on the cover of the book.\n",
"\n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/ml/bda_cover_1.png\") \n",
" | \n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/ml/bda_cover.png\") \n",
" | \n",
"
\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/kern/basis_covariance.svg\") \n",
" | \n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/kern/basis_covariance.gif\") \n",
" | \n",
"
\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/kern/brownian_covariance.svg\") \n",
" | \n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/kern/brownian_covariance.gif\") \n",
" | \n",
"
\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/kern/mlp_covariance.svg\") \n",
" | \n",
"\n",
"![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/kern/mlp_covariance.gif\") \n",
" | \n",
"
![](\"http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gpy.png\")
\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