<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Spike and Slab GP-LVM for Factor Analysis
Model Definition
GP-LVM (Lawrence 2005) with a linear kernel
Step1:  The obersved data $Y$ is generated by projecting the samples of the 3 sine waves onto a 100D space with a randomly generated weight matrix $W$. The first two dimensions of $Y$ is shown in the figure.
Step2:  Applying Spike and Slab GP-LVM
We create a Spike and Slab GPLVM model object with a linear kernel
Step3:  Apply Bayesian GP-LVM
As a comparison, we also apply the synthetic data to the original Bayesian GPLVM model.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import GPy
%matplotlib inline
#%config InlineBackend.figure_format = 'svg'
from pylab import *
np.random.seed(113321)
# Prepare the data
N,D,Q = 500, 100, 3
pi = 0.2
# sample from 3 random sine waves
X = np.sin(2*np.pi*(np.random.rand(Q)[None,:]+.5)*(np.linspace(0.,3.,N)[:,None]-np.random.rand(Q)[None,:]))
# set 80% of the values to zeros
X = X*((np.random.rand(N,Q)<pi)*1)
# plot the samples
_ = plot(X)
Explanation: Spike and Slab GP-LVM for Factor Analysis
Model Definition
GP-LVM (Lawrence 2005) with a linear kernel:
\begin{equation}
p(Y|X) = \prod_{d=1}^D p(y_{:,d}| X), \quad p(y_{:,d}|X) = \mathcal{N}(y_{:,d}; 0, XX^T+\beta^{-1}I)
\end{equation}
It is related to probabilistic PCA, while giving a prior distribution over the weight matrix $W$ instead of $X$ and marginalizing $W$:
\begin{equation}
p(Y|X) = \prod_{d=1}^D \mathcal{N}(y_{:,d};Xw_d, \beta^{-1}I), \quad
p(W) = \prod_{d=1}^D \mathcal{N}(w_d; 0, I).
\end{equation}
GP-LVM has been extended by giving a prior over $X$, known as Bayesian GP-LVM (Titsias & Lawrence 2010).
\begin{equation}
p(X) = \prod_{n=1}^N \mathcal{N}(x_n; 0, I)
\end{equation}
For factor analysis, it is often useful to assume $X$ to be sparse. Therefore, we extend Bayesian GP-LVM by assuming a Spike and Slab prior over $X$. Then, the GP model definition becomes:
\begin{align}
p(X)=\prod_{n=1}^N \mathcal{N}(x_n; 0, I), &\quad p(B) = \prod_{n=1}^N \prod_{q=1}^Q \pi^{b_{nq}} (1-\pi)^{1-b_{nq}},\
p(Y|X) = \prod_{d=1}^D p(y_{:,d}| X), &\quad p(y_{:,d}|X) = \mathcal{N}(y_{:,d}; 0, K+\beta^{-1}I), 
\end{align}
where $K$ is calculated according to the corresponding linear kernel function:
\begin{equation}
k(x,b,x',b') = \sum_{q=1}^Q \sigma^2_q x_q b_q x'_q b'_q
\end{equation}
An Experiment with Synthetic Data
We generate some artifical data and mixing them with a weight matrix.
End of explanation
noise_sigma = 0.1
# Generate a weight matrix
W = np.random.randn(D,Q)
# Project to the observed space plus some Gaussian noise
y = np.dot(X,W.T)+np.random.randn(N,D)*noise_sigma
y = (y-y.mean())/y.std()
# Plot the 
_ = plot(y[:,:2])
Explanation: The obersved data $Y$ is generated by projecting the samples of the 3 sine waves onto a 100D space with a randomly generated weight matrix $W$. The first two dimensions of $Y$ is shown in the figure.
End of explanation
m = GPy.models.SSGPLVM(y, Q, kernel=GPy.kern.Linear(Q,ARD=True),num_inducing=Q)
# Show the model parameters
m.likelihood.variance = .01
m.X.variance = .1
print m
# Optimize the model parameters
m.optimize('bfgs',messages=1,max_iters=10000)
# Show the model parameters after optimization
print m
print m.variational_prior.pi
# Plot the posterior distribution of X
# The left figure shows the mean and variance (the thickness) of 
# the slab part (the gray color shows the other dimensions for comparison)
# The right figure shows the posterior probability of the binary variable for having the slab
a = m.X.plot()
# In order to compare the recovered X, we flipped the sign of the 2nd, 3rd dimensions, and rescale the values.
# The figure shows the absolute difference between the recovered X and the original X
Xcom = m.X.mean.copy()*np.array([1.,-1.,-1.])
Xcom = Xcom/np.abs(Xcom).max(axis=0)
_ = plot(np.abs(Xcom-X))
# Plot the distribution of the recovered X (only 1st and 3rd dimension)
# As shown in the figure, it follows a spike and slab distribution
_ = m.plot_latent()
Explanation: Applying Spike and Slab GP-LVM
We create a Spike and Slab GPLVM model object with a linear kernel:
End of explanation
m = GPy.models.BayesianGPLVM(y, Q, kernel=GPy.kern.Linear(Q,ARD=True),num_inducing=Q)
# Show the model parameters
print m
# Optimize the model parameters
m.optimize('bfgs',messages=1,max_iters=10000)
# Show the model parameters after optimization
print m
# Plot the posterior distribution of X
a = m.X.plot()
# In order to compare the recovered X, we swapped the 2nd and 3rd dimension, 
# flipped the sign of the 2nd dimension, and rescale the values.
# The figure shows the absolute difference between the recovered X and the original X
Xcom = np.dot(m.X.mean.copy(),np.array([[1.,0.,0.],[0.,0.,1.],[0.,-1.,0.]]))
Xcom = Xcom/np.abs(Xcom).max(axis=0)
_ = plot(np.abs(Xcom-X))
# Plot the distribution of the recovered X (only 1st and 2nd dimension)
# As shown in the figure, it follows a spike and slab distribution
_ = m.plot_latent()
Explanation: Apply Bayesian GP-LVM
As a comparison, we also apply the synthetic data to the original Bayesian GPLVM model.
End of explanation
<END_TASK>