<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Overlapping Mixtures of Gaussian Processses
Valentine Svensson 2015 <br> (with small edits by James Hensman November 2015)
This illustrates use of the OMGP model described in
Overlapping Mixtures of Gaussian Processes for the data association problem
Miguel Lázaro-Gredilla, Steven Van Vaerenbergh, Neil D. Lawrence
Pattern Recognition 2012
The GPclust implementation makes use of the collapsed variational mixture model for GP assignment.
Step1:  Diverging trend seperation
One application of the OMGP model could be to find diverging trends among populations over time. Imagine for example two species evolving from a common ancestor over time.
We load some pre-generated data which diverge over time.
Step2:  We define a model assuming K = 2 trends. By default the model will be populated by K RBF kernels. The OMGP implementation is compatible with most kernels in GPy, so that you for example can encode periodicity in the model.
Step3:  A simple plot function is included which illustrates the asignment probability for each data point, it also shows the posterior mean and confidence intervals for each Gaussian Process.
Step4:  There is also a function for plotting the assignment probability for a given GP directly. Since we haven't optimized the mixture parameters yet the assignment probability is just a random draw from the prior.
Step5:  We can first performa a quick optimization to find the rough trends.
Step6:  The model identifies the branches of the time series, and in particular the non-branched region have ambigous GP assignment. In this region the two trends share information for prediction.
Like any GPy model the hyper parameters can be inspected.
Step7:  We continue by letting the model optimize some more, and also allow it to optimize the hyper parameters. The hyper parameter optimization works best if the mixture parameters have converged or are close to converging.
Step8:  Separating signal from noise
An interesting application of the OMGP model pointed out in the original publication is the use for robust GP regression.
Let's illustrate this by creating sinusoidal test data with background noise.
Step9:  First we make a model with only one mixture component / kernel. This is equivalent to normal GP regression.
Step10:  Now we in stead view this is a mixture problem, and consider two different kinds of kernels for the different GP components. One encoding white noise, and another which can encode a trend over time (an RBF kernel in this case).
Step11:  The trend over time is much more noticable, and the confidence intervals are smaller.
Noisy points will have high assignment probability to the 'noise GP', while the assignment of the sinusoidal points is ambiguous. We can use this to seperate the points which are more likely to be noise from the remaining points.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import GPy
from GPclust import OMGP
import matplotlib
matplotlib.rcParams['figure.figsize'] = (12,6)
from matplotlib import pyplot as plt
Explanation: Overlapping Mixtures of Gaussian Processses
Valentine Svensson 2015 <br> (with small edits by James Hensman November 2015)
This illustrates use of the OMGP model described in
Overlapping Mixtures of Gaussian Processes for the data association problem
Miguel Lázaro-Gredilla, Steven Van Vaerenbergh, Neil D. Lawrence
Pattern Recognition 2012
The GPclust implementation makes use of the collapsed variational mixture model for GP assignment.
End of explanation
XY = np.loadtxt('../data/split_data_test.csv', delimiter=',', skiprows=1, usecols=[1, 2])
X = XY[:, 0, None]
Y = XY[:, 1, None]
plt.scatter(X, Y);
Explanation: Diverging trend seperation
One application of the OMGP model could be to find diverging trends among populations over time. Imagine for example two species evolving from a common ancestor over time.
We load some pre-generated data which diverge over time.
End of explanation
m = OMGP(X, Y, K=2, variance=0.01, prior_Z='DP')
m.log_likelihood()
Explanation: We define a model assuming K = 2 trends. By default the model will be populated by K RBF kernels. The OMGP implementation is compatible with most kernels in GPy, so that you for example can encode periodicity in the model.
End of explanation
m.plot()
Explanation: A simple plot function is included which illustrates the asignment probability for each data point, it also shows the posterior mean and confidence intervals for each Gaussian Process.
End of explanation
m.plot_probs(gp_num=0)
Explanation: There is also a function for plotting the assignment probability for a given GP directly. Since we haven't optimized the mixture parameters yet the assignment probability is just a random draw from the prior.
End of explanation
m.optimize(step_length=0.01, maxiter=20)
m.plot()
m.plot_probs()
Explanation: We can first performa a quick optimization to find the rough trends.
End of explanation
m
Explanation: The model identifies the branches of the time series, and in particular the non-branched region have ambigous GP assignment. In this region the two trends share information for prediction.
Like any GPy model the hyper parameters can be inspected.
End of explanation
m.optimize(step_length=0.01, maxiter=200)
m
m.plot()
m.plot_probs()
Explanation: We continue by letting the model optimize some more, and also allow it to optimize the hyper parameters. The hyper parameter optimization works best if the mixture parameters have converged or are close to converging.
End of explanation
x1 = np.random.uniform(0, 10, (100, 1))
x2 = np.random.uniform(0, 10, (100, 1))
y1 = 4 * np.random.randn(*x1.shape)
y2 = 3 * np.sin(x2) + 0.5 * np.random.randn(*x2.shape)
x = np.vstack((x1, x2))
y = np.vstack((y1, y2))
plt.scatter(x, y);
Explanation: Separating signal from noise
An interesting application of the OMGP model pointed out in the original publication is the use for robust GP regression.
Let's illustrate this by creating sinusoidal test data with background noise.
End of explanation
kernels = [GPy.kern.RBF(1)]
m = OMGP(x, y, K=1, prior_Z='DP', kernels=kernels)
m.variance = 3
m.hyperparam_interval = 100
m.rbf.lengthscale = 2
m.optimize(verbose=False)
m.plot()
Explanation: First we make a model with only one mixture component / kernel. This is equivalent to normal GP regression.
End of explanation
kernels = [GPy.kern.White(1, name='Noise'), GPy.kern.RBF(1, name='Signal')]
m = OMGP(x, y, K=2, prior_Z='DP', kernels=kernels)
m.variance = 3
m.hyperparam_interval = 250
m.Signal.lengthscale = 2
m.plot(0)
m.optimize(step_length=0.01, verbose=False)
m
m.plot()
Explanation: Now we in stead view this is a mixture problem, and consider two different kinds of kernels for the different GP components. One encoding white noise, and another which can encode a trend over time (an RBF kernel in this case).
End of explanation
m.plot_probs(0)
plt.axhline(0.75);
thr = 0.75
idx = np.where(m.phi[:,0] < thr)[0]
nidx = np.where(m.phi[:,0] >= thr)[0]
plt.figure(figsize=(12,10))
plt.subplot(211)
plt.scatter(x[idx], y[idx]);
plt.title('Signal')
plt.subplot(212, sharey=plt.gca())
plt.scatter(x[nidx], y[nidx]);
plt.title('Noise');
Explanation: The trend over time is much more noticable, and the confidence intervals are smaller.
Noisy points will have high assignment probability to the 'noise GP', while the assignment of the sinusoidal points is ambiguous. We can use this to seperate the points which are more likely to be noise from the remaining points.
End of explanation
<END_TASK>