<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 A Brief Guide of Using Hybrid Monte Carlo for Gaussian Process Models
Zhenwen Dai and Mu Niu, November 2014
Step1:  Example 1  HMC Inference for GP Regression
Let's first generate some synthetic data.
Step2:  Let us Make a GP Regression model and give some general prior distributions to model parameters.
Step3:  Sampling the posterior distribution of model parameters
Draw 1000 samples from the GP model
Step4:  Plot the samples
Step5:  Plot the posterior marginal distribution of model parameters
Step6:  Plot the model parameters (lengthscale, variance and noise variance) against each other
Step7:  By setting the model parameters to the posterior mean, we can visualize the model fit
Step8:  Sample the posterior distribution of X given some new Y
Given some new observations, inferring the posterior distribution of the corresponding inputs is difficult, because it can lead to multi-modal distributions. 
Assume we have a new observation $1.5$, and try to infer its input distribution.
Step9:  Generate the inference model for the new observations. X_new are the MAP estimations by optimizing the log likelihood. As plotted with a red dot, the MAP estimation corresponds to only one of the modes.
Step10:  Draw 10,000 samples from the inference model
Step11:  Plot the samples
Step12:  Plot the marginal distribution of inferred inputs. The two modes of inputs are clearly visible from the sampled posterior distribution.
Step13:  Example 2 HMC for lengthscale and variance with marathon data
we set prior for lengthscale and variance of kernel. The mean of the prior is close to the result of GP optimisation. we then allow a big variance. In the case below,we set gamma prior to lengthscale and variance. E(lengthscale) = 120, Var(lengthscale)=2000, E(variance) = 25, Var(variance) = 150
Step14:  Run HMC
we plot the full length of hmc iteration. The first 5000 could be burn in stage and can be ignored.
Step15:  Seaborn
Using Seaborn for plotting distributions over Hyperparameters
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import GPy
%matplotlib inline
#%config InlineBackend.figure_format = 'svg'
from pylab import *
Explanation: A Brief Guide of Using Hybrid Monte Carlo for Gaussian Process Models
Zhenwen Dai and Mu Niu, November 2014
End of explanation
# Let's make some synthetic data
x = np.linspace(0.,2*np.pi,100)[:,None]
y = -cos(x)+np.random.randn(*x.shape)*0.3+1
_ = plot(x,y,'.')
Explanation: Example 1  HMC Inference for GP Regression
Let's first generate some synthetic data.
End of explanation
# Make a GP regression model
m = GPy.models.GPRegression(x,y)
# Give some general prior distributions for model parameters
m.kern.lengthscale.set_prior(GPy.priors.Gamma.from_EV(1.,10.))
m.kern.variance.set_prior(GPy.priors.Gamma.from_EV(1.,10.))
m.likelihood.variance.set_prior(GPy.priors.Gamma.from_EV(1.,10.))
_=m.plot()
Explanation: Let us Make a GP Regression model and give some general prior distributions to model parameters.
End of explanation
hmc = GPy.inference.mcmc.HMC(m,stepsize=5e-2)
s = hmc.sample(num_samples=1000) # Burnin
s = hmc.sample(num_samples=1000)
Explanation: Sampling the posterior distribution of model parameters
Draw 1000 samples from the GP model
End of explanation
plot(s)
Explanation: Plot the samples:
End of explanation
labels = ['kern variance', 'kern lengthscale','noise variance']
samples = s[300:] # cut out the burn-in period
from scipy import stats
xmin = samples.min()
xmax = samples.max()
xs = np.linspace(xmin,xmax,100)
for i in xrange(samples.shape[1]):
    kernel = stats.gaussian_kde(samples[:,i])
    plot(xs,kernel(xs),label=labels[i])
_ = legend()
Explanation: Plot the posterior marginal distribution of model parameters:
End of explanation
fig = figure(figsize=(14,4))
ax = fig.add_subplot(131)
_=ax.plot(samples[:,0],samples[:,1],'.')
ax.set_xlabel(labels[0]); ax.set_ylabel(labels[1])
ax = fig.add_subplot(132)
_=ax.plot(samples[:,1],samples[:,2],'.')
ax.set_xlabel(labels[1]); ax.set_ylabel(labels[2])
ax = fig.add_subplot(133)
_=ax.plot(samples[:,0],samples[:,2],'.')
ax.set_xlabel(labels[0]); ax.set_ylabel(labels[2])
Explanation: Plot the model parameters (lengthscale, variance and noise variance) against each other:
End of explanation
# Set the model parameters as the posterior mean
m.kern.variance[:] = samples[:,0].mean()
m.kern.lengthscale[:] = samples[:,1].mean()
m.likelihood.variance[:] = samples[:,2].mean()
print m
_=m.plot()
Explanation: By setting the model parameters to the posterior mean, we can visualize the model fit:
End of explanation
y_new = np.array([1.5])[:,None]
Explanation: Sample the posterior distribution of X given some new Y
Given some new observations, inferring the posterior distribution of the corresponding inputs is difficult, because it can lead to multi-modal distributions. 
Assume we have a new observation $1.5$, and try to infer its input distribution.
End of explanation
x_new,mi = m.infer_newX(y_new)
print mi
m.plot()
plot(x_new,y_new,'or')
Explanation: Generate the inference model for the new observations. X_new are the MAP estimations by optimizing the log likelihood. As plotted with a red dot, the MAP estimation corresponds to only one of the modes.
End of explanation
hmc_new = GPy.inference.mcmc.HMC(mi,stepsize=2e-1)
s_new = hmc_new.sample(num_samples=10000,hmc_iters=10)
Explanation: Draw 10,000 samples from the inference model:
End of explanation
_ = plot(s_new[:,:])
Explanation: Plot the samples:
End of explanation
from scipy import stats
samples_new = s_new[:]
xmin = samples_new.min()
xmax = samples_new.max()
xs = np.linspace(xmin,xmax,100)
for i in xrange(samples_new.shape[1]):
    kernel = stats.gaussian_kde(samples_new[:,i])
    plot(xs,kernel(xs))
Explanation: Plot the marginal distribution of inferred inputs. The two modes of inputs are clearly visible from the sampled posterior distribution.
End of explanation
m = GPy.examples.regression.olympic_marathon_men()
#
#set prior for lengthscale and variance.
m.kern.variance.set_prior(GPy.priors.Gamma.from_EV(25.,150.))
m.kern.lengthscale.set_prior(GPy.priors.Gamma.from_EV(120.,2000.))
print m
Explanation: Example 2 HMC for lengthscale and variance with marathon data
we set prior for lengthscale and variance of kernel. The mean of the prior is close to the result of GP optimisation. we then allow a big variance. In the case below,we set gamma prior to lengthscale and variance. E(lengthscale) = 120, Var(lengthscale)=2000, E(variance) = 25, Var(variance) = 150
End of explanation
# initialise hmc
hmc = GPy.inference.mcmc.HMC(m,stepsize=2e-1)
# run hmc
t = hmc.sample(num_samples=20000,hmc_iters=20)
# Sample parameters 
#hmc = GPy.inference.optimization.HMC(m, stepsize=5e-1)
#t = hmc.sample(m_iters=50000,hmc_iters=20)
_=plot(t)
print t.mean(axis=0)
print t.std(axis=0)
_=hist(t[:,:2],50)
Explanation: Run HMC
we plot the full length of hmc iteration. The first 5000 could be burn in stage and can be ignored.
End of explanation
import seaborn as sns, pandas as pd
plt.rcParams['text.usetex'] = False
df = pd.DataFrame(t, columns=m.parameter_names_flat())
ax = sns.kdeplot(df['rbf.variance'],
                 color="b", shade=True, shade_lowest=False)
ax = sns.kdeplot(df['rbf.lengthscale'],
                 color="r", shade=True, shade_lowest=False)
sns.set(style="white", color_codes=True)
_ = sns.jointplot(data=df, x='rbf.variance', y='rbf.lengthscale', kind="hex", 
                  marginal_kws=dict(kde=True, hist=True, kde_kws=dict(shade=False)),
                 stat_func=None
                 )
df
Explanation: Seaborn
Using Seaborn for plotting distributions over Hyperparameters:
End of explanation
<END_TASK>