import pods
pods.notebook.display_google_book(id='1YQPAAAAQAAJ', page='PR17-IA2')

import pods
pods.notebook.display_google_book(id='1YQPAAAAQAAJ', page='PR17-IA4')

import pods
from ipywidgets import IntSlider

pods.notebook.display_plots('gp_rejection_sample{sample:0>3}.png', 
                            directory='../slides/diagrams/gp', 
                            sample=IntSlider(1,1,5,1))

%load -s Kernel mlai.py

%load -s polynomial_cov mlai.py

%load -s exponentiated_quadratic mlai.py

import pods
from ipywidgets import IntSlider

pods.notebook.display_plots('two_point_sample{sample:0>3}.svg', '../slides/diagrams/gp', sample=IntSlider(0, 0, 8, 1))

import pods
from ipywidgets import IntSlider

pods.notebook.display_plots('two_point_sample{sample:0>3}.svg', 
                            '../slides/diagrams/gp', 
                            sample=IntSlider(9, 9, 12, 1))

import pods
from ipywidgets import IntSlider

pods.notebook.display_plots('two_point_sample{sample:0>3}.svg', 
                            '../slides/diagrams/gp', 
                            sample=IntSlider(13, 13, 17, 1))

$$\mappingFunction_D(\inputVector_*) = \kernelVector(\inputVector_*, \inputMatrix) \boldsymbol{\alpha}$$

import numpy as np
import pods

data = pods.datasets.olympic_marathon_men()
x = data['X']
y = data['Y']

offset = y.mean()
scale = np.sqrt(y.var())

import matplotlib.pyplot as plt
import teaching_plots as plot
import mlai


xlim = (1875,2030)
ylim = (2.5, 6.5)
yhat = (y-offset)/scale

fig, ax = plt.subplots(figsize=plot.big_wide_figsize)
_ = ax.plot(x, y, 'r.',markersize=10)
ax.set_xlabel('year', fontsize=20)
ax.set_ylabel('pace min/km', fontsize=20)
ax.set_xlim(xlim)
ax.set_ylim(ylim)

mlai.write_figure(figure=fig, 
                  filename='../slides/diagrams/datasets/olympic-marathon.svg', 
                  transparent=True, 
                  frameon=True)

import GPy

m_full = GPy.models.GPRegression(x,yhat)
_ = m_full.optimize() # Optimize parameters of covariance function

xt = np.linspace(1870,2030,200)[:,np.newaxis]
yt_mean, yt_var = m_full.predict(xt)
yt_sd=np.sqrt(yt_var)

import teaching_plots as plot

fig, ax = plt.subplots(figsize=plot.big_wide_figsize)
plot.model_output(m_full, scale=scale, offset=offset, ax=ax, xlabel='year', ylabel='pace min/km', fontsize=20, portion=0.2)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
mlai.write_figure(figure=fig,
                  filename='../slides/diagrams/gp/olympic-marathon-gp.svg', 
                  transparent=True, frameon=True)

x_clean=np.vstack((x[0:2, :], x[3:, :]))
y_clean=np.vstack((y[0:2, :], y[3:, :]))

m_clean = GPy.models.GPRegression(x_clean,y_clean)
_ = m_clean.optimize()

gpoptimizePlot1


\only<1>{\input{../../../gp/tex/diagrams/gpOptimiseQuadratic1.svg}}\only<2>{\input{../../../gp/tex/diagrams/gpOptimiseQuadratic2.svg}}\only<3>{\input{../../../gp/tex/diagrams/gpOptimiseQuadratic3.svg}}-->
<table>
<tr>
<td width="50%">
<img src="http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp-optimise006.svg" class="" align="100%" style="vertical-align:middle;">
</td>
<td width="50%">
<img src="http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp-optimise010.svg" class="" align="100%" style="vertical-align:middle;">
</td>
</tr>
</table>
<table>
<tr>
<td width="50%">
<img src="http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp-optimise016.svg" class="" align="100%" style="vertical-align:middle;">
</td>
<td width="50%">
<img src="http://inverseprobability.com/talks/slides/../slides/diagrams/gp/gp-optimise021.svg" class="" align="100%" style="vertical-align:middle;">
</td>
</tr>
</table>
Figure: <i>Variation in the data fit term, the capacity term and the
negative log likelihood for different lengthscales.</i>

## Gene Expression Example \[<a href="https://github.com/lawrennd/snippets/edit/main/_gp/includes/della-gatta-gene-gp.md" target="_blank" >edit</a>\]

We now consider an example in gene expression. Gene expression is the
measurement of mRNA levels expressed in cells. These mRNA levels show
which genes are 'switched on' and producing data. In the example we will
use a Gaussian process to determine whether a given gene is active, or
we are merely observing a noise response.

## Della Gatta Gene Data \[<a href="https://github.com/lawrennd/snippets/edit/main/_ml/includes/della-gatta-gene-data.md" target="_blank" >edit</a>\]

-   Given given expression levels in the form of a time series from
    @DellaGatta:direct08.

``` {.python}
import numpy as np
import pods

data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=937)

x = data['X']
y = data['Y']

offset = y.mean()
scale = np.sqrt(y.var())

import matplotlib.pyplot as plt
import teaching_plots as plot
import mlai


xlim = (-20,260)
ylim = (5, 7.5)
yhat = (y-offset)/scale

fig, ax = plt.subplots(figsize=plot.big_wide_figsize)
_ = ax.plot(x, y, 'r.',markersize=10)
ax.set_xlabel('time/min', fontsize=20)
ax.set_ylabel('expression', fontsize=20)
ax.set_xlim(xlim)
ax.set_ylim(ylim)

mlai.write_figure(figure=fig, 
                  filename='../slides/diagrams/datasets/della-gatta-gene.svg', 
                  transparent=True, 
                  frameon=True)

import GPy

m_full = GPy.models.GPRegression(x,yhat)
m_full.kern.lengthscale=50
_ = m_full.optimize() # Optimize parameters of covariance function

xt = np.linspace(-20,260,200)[:,np.newaxis]
yt_mean, yt_var = m_full.predict(xt)
yt_sd=np.sqrt(yt_var)

import teaching_plots as plot

fig, ax = plt.subplots(figsize=plot.big_wide_figsize)
plot.model_output(m_full, scale=scale, offset=offset, ax=ax, xlabel='time/min', ylabel='expression', fontsize=20, portion=0.2)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.set_title('log likelihood: {ll:.3}'.format(ll=m_full.log_likelihood()), fontsize=20)
mlai.write_figure(figure=fig,
                  filename='../slides/diagrams/gp/della-gatta-gene-gp.svg', 
                  transparent=True, frameon=True)

m_full2 = GPy.models.GPRegression(x,yhat)
m_full2.kern.lengthscale=2000
_ = m_full2.optimize() # Optimize parameters of covariance function

import teaching_plots as plot

fig, ax = plt.subplots(figsize=plot.big_wide_figsize)
plot.model_output(m_full2, scale=scale, offset=offset, ax=ax, xlabel='time/min', ylabel='expression', fontsize=20, portion=0.2)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.set_title('log likelihood: {ll:.3}'.format(ll=m_full2.log_likelihood()), fontsize=20)
mlai.write_figure(figure=fig,
                  filename='../slides/diagrams/gp/della-gatta-gene-gp2.svg', 
                  transparent=True, frameon=True)

m_full3 = GPy.models.GPRegression(x,yhat)
m_full3.kern.lengthscale=20
m_full3.likelihood.variance=0.001
_ = m_full3.optimize() # Optimize parameters of covariance function

import teaching_plots as plot

fig, ax = plt.subplots(figsize=plot.big_wide_figsize)
plot.model_output(m_full3, scale=scale, offset=offset, ax=ax, xlabel='time/min', ylabel='expression', fontsize=20, portion=0.2)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.set_title('log likelihood: {ll:.3}'.format(ll=m_full3.log_likelihood()), fontsize=20)
mlai.write_figure(figure=fig,
                  filename='../slides/diagrams/gp/della-gatta-gene-gp3.svg', 
                  transparent=True, frameon=True)

%load -s basis_cov mlai.py

%load -s radial mlai.py

%load -s brownian_cov mlai.py

%load -s mlp_cov mlai.py

%load -s relu_cov mlai.py

%load -s sinc_cov mlai.py

%load -s lmc_cov mlai.py

from IPython.core.display import HTML

HTML(anim.to_jshtml())

plot.save_animation(anim, 
                    diagrams='../slides/diagrams/kern', 
                    filename='lmc_covariance.html')

%load -s icm_cov mlai.py

from IPython.core.display import HTML

HTML(anim.to_jshtml())

plot.save_animation(anim, 
                    diagrams='../slides/diagrams/kern', 
                    filename='icm_covariance.html')