# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 2
#     language: python
#     name: python2
# ---

# # Lab 3: Gaussian process regression
#
# ### Machine Learning 1, September 2015
#
# * The lab exercises should be made in groups of two, three or four people.
# * The deadline is October 25th (Sunday) 23:59.
# * Assignment should be sent to <NAME> (<EMAIL>). The subject line of your email should be "lab\#\_lastname1\_lastname2\_lastname3".
# * Put your and your teammates' names in the body of the email.
# * Attach the .IPYNB (IPython Notebook) file containing your code and answers. Naming of the file follows the same rule as the subject line. For example, if the subject line is "lab01\_Kingma\_Hu", the attached file should be "lab01\_Kingma\_Hu.ipynb". Only use underscores ("\_") to connect names, otherwise the files cannot be parsed.
#
# Notes on implementation:
#
# * You should write your code and answers in an IPython Notebook: http://ipython.org/notebook.html. If you have problems, please contact us.
# * Among the first lines of your notebook should be "%pylab inline". This imports all required modules, and your plots will appear inline.
# * NOTE: Make sure we can run your notebook / scripts!
# $\newcommand{\bx}{\mathbf{x}}$
# $\newcommand{\bxp}{\mathbf{x}^{'}}$
# $\newcommand{\bw}{\mathbf{w}}$
# $\newcommand{\bt}{\mathbf{t}}$
# $\newcommand{\by}{\mathbf{y}}$
# $\newcommand{\bm}{\mathbf{m}}$
# $\newcommand{\bb}{\mathbf{b}}$
# $\newcommand{\bS}{\mathbf{S}}$
# $\newcommand{\ba}{\mathbf{a}}$
# $\newcommand{\bz}{\mathbf{z}}$
# $\newcommand{\bv}{\mathbf{v}}$
# $\newcommand{\bq}{\mathbf{q}}$
# $\newcommand{\bp}{\mathbf{p}}$
# $\newcommand{\bh}{\mathbf{h}}$
# $\newcommand{\bI}{\mathbf{I}}$
# $\newcommand{\bX}{\mathbf{X}}$
# $\newcommand{\bT}{\mathbf{T}}$
# $\newcommand{\bPhi}{\mathbf{\Phi}}$
# $\newcommand{\bW}{\mathbf{W}}$
# $\newcommand{\bV}{\mathbf{V}}$
# $\newcommand{\xm}{\mathbf{x}_m}$
# $\newcommand{\xn}{\mathbf{x}_n}$
# $\newcommand{\y}{\mathbf{y}}$
# $\newcommand{\K}{\mathbf{K}}$
# $\newcommand{\zero}{\mathbf{0}}$
# $\newcommand{\yi}{\y_i}$
# $\newcommand{\thetav}{\mathbf{\theta}}$
# $\newcommand{\t}{\mathbf{t}}$
# $\newcommand{\x}{\mathbf{x}}$
# $\newcommand{\tN}{\mathbf{t}_N}$
# $\newcommand{\xN}{\mathbf{x}_N}$
# $\newcommand{\k}{\mathbf{k}}$
# $\newcommand{\C}{\mathbf{C}}$
# $\newcommand{\CN}{\mathbf{C}_N}$
# $\newcommand{\KN}{\mathbf{K}_N}$
# $\newcommand{\eyeN}{\mathbf{I}_N}$

# # Gaussian process regression
#
# For this Lab we will be refer to Bishop sections 6.4.2 and 6.4.3. You may also want to refer to Rasmussen's Gaussian Process text which is available online at http://www.gaussianprocess.org/gpml/chapters/ and especially to the project found at http://www.automaticstatistician.com/index.php by Ghahramani for some intuition in GP.  To understand Gaussian processes, it is highly recommended understand how marginal, partitioned Gaussian distributions can be converted into conditional Gaussian distributions.  This is covered in Bishop 2.3 and summarized in Eqns 2.94-2.98.
#
#
#

# ### Sinusoidal Data

# We will use the same data generating function that we used previously for regression.  You can change sigma/beta, but keep it reasonable.  Definitely play around once you have things working.  Make use of these functions as you wish.

# %pylab inline
import numpy as np
import pylab as pp
import matplotlib.pyplot as plt
import math
import numpy.matlib

sigma = 0.5
beta  = 1.0 / pow(sigma,2) # this is the beta used in Bishop Eqn. 6.59
N_test = 100
x_test = np.linspace(-1,1,N_test); 
mu_test = np.zeros( N_test )


# +
def true_mean_function( x ):
    return np.sin( 2*pi*(x+1) )

def add_noise( y, sigma ):
    return y + sigma*np.random.randn(len(y))

def generate_t( x, sigma ):
    return add_noise( true_mean_function( x), sigma )


# -

y_test = true_mean_function( x_test )
t_test = add_noise( y_test, sigma )
pp.plot( x_test, y_test, 'b-', lw=2)
pp.plot( x_test, t_test, 'go')


# ### 1. Sampling from the Gaussian process prior (30 points)

# We will implement Gaussian process regression using the kernel function in Bishop Eqn. 6.63.  
#
# #### 1.1 k_n_m( xn, xm, thetas ) (10 points)
# To start, implement function "k_n_m( xn, xm, thetas )" that takes scalars $\xn$ and $\xm$, and a vector of $4$ thetas, and computes the kernel function Bishop Eqn. 6.63 (10 points). 

def k_n_m(xn,xm,thetas):
    return thetas[0]*np.exp((-thetas[1]/2) * np.sum(pow(xn-xm,2))) + thetas[2]+thetas[3]*xn.transpose().dot(xm)


# #### 1.2 computeK( X1, X2, thetas ) (5 points)
# Eqn 6.60 is the marginal distribution of mean ouput of $N$ data vectors: $p(\y) = \mathcal{N}(\zero, \K)$.  Notice that the expected mean function is $0$ at all locations, and that the covariance is a $N$ by $N$ kernel matrix $\K$.  Write a function "computeK( X1, X2, thetas )" that computes the kernel matrix. Hint: use k_n_m as part of an innner loop (of course, there are more efficient ways of computing the kernel function making better use of vectorization, but that is not necessary) (5 points).  

def computeK( X1, X2, thetas ):

    K= np.empty([len(X1),len(X2)])
    for n in range(0,len(X1)):
     for m in range(0,len(X2)):
            K[n,m]=k_n_m(np.array(X1[n]),np.array(X2[m]),thetas)
    return K


# #### 1.3 Plot function samples (15 points)
# Now sample mean functions at the x_test locations for the theta values in Bishop Figure 6.5, make a figure with a 2 by 3 subplot and make sure the title reflects the theta values (make sure everything is legible).  In other words, sample $\yi \sim \mathcal{N}(\zero, \K_{\thetav})$.  Make use of numpy.random.multivariate_normal().  On your plots include the expected value of $\y$ with a dashed line and fill_between 2 standard deviations of the uncertainty due to $\K$ (the diagonal of $\K$ is the variance of the model uncertainty) (15 points).

# +
def plot_1():
  N_test = 100
  x_test = np.linspace(-1,1,N_test)
  fig, ax = plt.subplots(2,3)

  thetas = [[1.00,4.00,0.00,0.00],[9.00,4.00,0.00,0.00],[1.00,64.00,0.00,0.00],[1.00,0.25,0.00,0.00],[1.00,4.00,10.00,0.00],[1.00,4.00,0.00,5.00]]

  i=0
  for theta in thetas:
      col = i%3
      row = math.floor(i/3)
      i+=1


      K=computeK(x_test,x_test,theta)
      var =K.diagonal()
      std=np.array([math.sqrt(x) for x in var])

      ys=np.random.multivariate_normal(np.zeros(N_test),K,5)
      for y in ys:
         ax[row,col].plot(x_test,y,'-')
         ax[row,col].set_title('('+','.join(str(x) for x in theta)+")")
      #mean function
      ax[row,col].plot(x_test,np.zeros(N_test),'--')
      ax[row,col].fill_between(x_test,np.zeros(N_test)+2*std,np.zeros(N_test)-2*std,color='pink',alpha=0.1)
  #    ax.ylim((-1.5,1.5))
  plt.show()


plot_1()


# -

# ### 2. Predictive distribution (35 points)

# So far we have sampled mean functions from the prior.  We can draw actual data $\t$ two ways.  The first way is generatively, by first sampling $\y | \K$, then sampling $\t | \y, \beta$ (Eqns 6.60 followed by 6.59).  The second way is to integrate over $\y$ (the mean draw) and directly sample $\t | \K, \beta$ using Eqn 6.61.    This is the generative process for $\t$.  Note that we have not specified a distribution over inputs $\x$;  this is because Gaussian processes are conditional models.  Because of this we are free to generate locations $\x$ when playing around with the GP; obviously a dataset will give us input-output pairs.
#
# Once we have data, we are interested in the predictive distribution (note: the prior is the predictive distribution when there is no data).  Consider the joint distribution for $N+1$ targets, given by Eqn 6.64.  Its covariance matrix is composed of block components $\CN$, $\k$, and $c$.  The covariance matrix $CN$ for $\tN$ is $\CN = \KN + \eyeN / \beta$.  We have just made explicit the size $N$ of the matrix; $N$ is the number of training points.  The kernel vector $\k$ is a $N$ by $1$ vector of kernel function evaluations between the training input data and the test input vector.  The scalar $c$ is a kernel evaluation at the test input.
#
# #### 2.1 gp_predictive_distribution(...) (10 points)
# Write a function "gp_predictive_distribution(x_train, t_train, x_test, theta, beta, C = None)" that computes  Eqns 6.66 and 6.67, except allow for an arbitrary number of test points (not just one) and now the kernel matrix is for training data.  By having C as an optional parameter, we can avoid computing it more than once (for this problem it is unimportant, but for real problems this is an issue).  The function should compute $\C$, $\k$, and $c$ and the mean and noise functions.  Do not forget: the computeK function computes $\K$, not $\C$! (10 points)

def gp_predictive_distribution(x_train, t_train, x_test, theta, beta, C = None):

    # computing ingredients
    c=[]
    for i in range(0,len(x_test)):
        t= computeK([x_test[0]],[x_test[0]],theta)
        c.append(t[0][0]+(1/beta))

    k = computeK(x_train,x_test,theta)
    C= C if C!=None else computeK(x_train,x_train,theta)+np.matlib.identity(np.shape(x_train)[0])/beta

    k_C=k.transpose().dot(np.linalg.inv(C))
    #  mean computation
    mean = k_C.dot(t_train)
    #  var computation
    var = c-k_C.dot(k)
    return mean, var


# #### 2.2 gp_log_likelihood(...) (10 points)
# Later, to learn the hyperparameters, we will need to compute the log-likelihood of the of the training data.  Implicitly, this is conditioned on the value setting for $\thetav$.  Write a function "gp_log_likelihood( x_train, t_train, theta, C = None, invC = None )", where C and invC can be stored and reused.  (10 points)

# +
def gp_log_likelihood( x_train, t_train, theta, beta, C = None, invC = None ):
    C= C if C!=None else computeK(x_train,x_train,theta)+np.matlib.identity(np.shape(x_train)[0])/beta
    invC= invC if invC!=None else np.linalg.inv(C)
    N=np.shape(x_train)[0]
    res= (-np.log(np.linalg.det(C))-t_train.T.dot(invC).dot(t_train)-N*np.log(2*np.pi))/2
    return np.asarray(res)[0][0]

x_train = np.linspace(-1, 1, 5)
y_train = true_mean_function(x_train)
t_train = add_noise(y_train, sigma)
#beta = 1
x_test = np.linspace(-1, 1, 7)
thetas = np.array([(1, 4, 0, 0), (9, 4, 0, 0), (1, 64, 0, 0), (1, 0.25, 0, 0), (1, 4, 10, 0), (1, 4, 0, 5)])

print (gp_log_likelihood(x_train,t_train,thetas[0],beta))


# -

# #### 2.3 Plotting (10 points)
# Repeat the 6 plots above, but this time conditioned on the training points.  Use the sinuosoidal data generator to create 2 training points where x is sampled uniformly between $-1$ and $1$.  For these plots, feel free to use the provided function "gp_plot".  Make sure you put the parameters in the title and this time also the log-likelihood. (10 points)  Try to understand the two types of uncertainty!  If you do not use "gp_plot", please add a fill between for the model and target noise. 

def gp_plot( x_test, y_test, mu_test, var_test, x_train, t_train, theta, beta ):
    # x_test:   the test data
    # y_test:   the true function at x_test
    # mu_test:  predictive mean at x_test
    # var_test: predictive covariance at x_test 
    # t_train:  the training values
    # theta:    the kernel parameters
    # beta:     the precision (known)
    
    # the reason for the manipulation is to allow plots separating model and data stddevs.
    std_total = np.sqrt(np.diag(var_test))         # includes all uncertainty, model and target noise 
    std_model = np.sqrt( std_total**2 - 1.0/beta ) # remove data noise to get model uncertainty in stddev
    std_combo = std_model + np.sqrt( 1.0/beta )    # add stddev (note: not the same as full)
    
    pp.plot( x_test, y_test, 'b', lw=3)
    pp.plot( x_test, mu_test, 'k--', lw=2 )
    pp.fill_between( x_test, mu_test+2*std_combo,mu_test-2*std_combo, color='k', alpha=0.25 )
    pp.fill_between( x_test, mu_test+2*std_model,mu_test-2*std_model, color='r', alpha=0.25 )
    pp.plot( x_train, t_train, 'ro', ms=10 )


def plot_2_3(N_train,N_test):


  thetas = [[1.00,4.00,0.00,0.00],[9.00,4.00,0.00,0.00],[1.00,64.00,0.00,0.00],[1.00,0.25,0.00,0.00],[1.00,4.00,10.00,0.00],[1.00,4.00,0.00,5.00]]
  x_train = np.random.uniform(-1, 1,N_train)
  t_train = generate_t(x_train,sigma)

  x_test = np.linspace(-1, 1,N_test)
  t_test = generate_t(x_test,sigma)
  y_test =true_mean_function(x_test)

 # fig, ax = plt.subplots(2,3)


  for row in range(2):
        for col in range(3):
            idx = row * 3 + col
            plt.subplot(2, 3, idx + 1)
            mu,var = gp_predictive_distribution(x_train,t_train,x_test,thetas[idx],beta)
            mu=np.squeeze(np.asarray(mu.T))
            gp_plot(x_test,y_test,mu,var,x_train,t_train,thetas[idx],beta)
            plt.title(str(thetas[idx]))

  plt.show()

# +
N_train = 2
N_test=100

plot_2_3(N_train,N_test)
# -

# #### 2.4 More ploting (5 points)
# Repeat the 6 plots above, but this time conditioned a new set of 10 training points. (5 points)

# +

N_train = 10
N_test=100
plot_2_3(N_train,N_test)


# -

# ### 3. Learning the hyperparameters (45 points)

# Learning the values of the parameter $\thetav$ can be very tricky for Gaussian processes in general, but when the data is univariate like ours, we can visualize the fit and see how plausible it looks.
#
# #### 3.1 Derivatives (5 points)
# Maximum likelihood or MAP learning is the most common way of setting the parameters, though a fully Bayesian approach is possible too.  We will look at ML today.  For this, we start with the dervivative of the log-likelihood with respect to the parameters $\thetav$; this is Eqn 6.70.  This, in turn, requires the derivative of the kernel matrix $\CN$ wrt $\thetav$.  This is the matrix of element-wise derivatives of the kernel function.  Write the derivatives for $\theta_0$ to $\theta_3$ for our kernel function (5 points).  

# $\frac{\partial k(x_n,x_m)}{\partial \theta_0}= \exp(-\frac{\theta_1}{2} ||x_n-x_m||^2)$
#
# $\frac{\partial k(x_n,x_m)}{\partial \theta_1}= - \frac{\theta_0}{2} || x_n-x_m||^2 \exp(-\frac{\theta_1}{2} ||x_n-x_m||^2)$
#
# $\frac{\partial k(x_n,x_m)}{\partial \theta_2}= 1$
#
# $\frac{\partial k(x_n,x_m)}{\partial \theta_3}= x_n^Tx_m$
#

# #### 3.2 Questions (5 points)
# Which parameters in $\thetav$ are unconstrained, that is, where any positive/ negative values are valid? (5 points)

# We start by decomposing our kernal function into a sum of 3 kernal functions:
#
# $k(x_n,x_m)= \theta_0 k_1(x_n,x_m) + \theta_2  k_2(x_n,x_m) + \theta_3 k_3(x_n,x_m)$
#
# As we can see $\theta_0$, $\theta_2$, and $\theta_3$ are constants and therefore should be >0 for the kernal to be valid. On the other hand, $\theta_1$ plays a role of a bandwidth in $k_1$ and is not constrained.  

# #### 3.3 More derivatives (5 points)
# For parameters that are constrained to be positive, the usual approach is to use the exponential of the free-parameter in the kernel function, but perform gradient ascent on the unconstrained values.  Consider the case  $\theta_i = \exp( \phi_i)$, where $\phi_i$ is unconstrained.  Write the derivative for $\phi_i$ in terms of the derivatives you already computed (5 points).  Hint: use the chain rule and do not repeat the full derivation.
#

# ___answer___
#
# $\frac{\partial k(x_n,x_m)}{\partial \phi_0}= \exp(-\frac{\theta_1}{2} ||x_n-x_m||^2 + \phi_0)$
#
#
# $\frac{\partial k(x_n,x_m)}{\partial \phi_2}= \exp(\phi_2)$
#
# $\frac{\partial k(x_n,x_m)}{\partial \phi_3}= x_n^Tx_m \exp(\phi_3)$
#

# #### 3.4 Grid search (10 points)
# Grid-search: for the same training set you have above, perform a small grid search over $\thetav$ (try at least 20 combinations).  Have your grid-search loop or function print out rows of log-likelihood + $\thetav$ sorted by best to worst.  Use the log-likelihood to select the best $\thetav$ and the worst.  Plots both the same way as the subplots above (ie a 1 by 2 subplot of best and worst). (10 points)

# +
def grid_search(iter):

  sigma = 0.5
  beta  = 1.0 / pow(sigma,2) # this is the beta used in Bishop Eqn. 6.59
  N_train = 30
  N_test = 100

  x_train = np.random.uniform(-1, 1,N_train)
  t_train = generate_t(x_train,sigma)

  x_test = np.linspace(-1, 1,N_test)
  t_test = generate_t(x_test,sigma)
  y_test =true_mean_function(x_test)


  res ={}

  for i in range(0,iter):
    thetas=np.random.uniform(1,10,4)
    lh=gp_log_likelihood(x_train,t_train,thetas,beta)
    res[lh]=thetas
    print ('likelihood: '+str(lh)+' thetas: '+str(thetas))

  keys=sorted(res.keys())
  best_theta = np.array(res[keys[iter-1]])
  worst_theta = np.array(res[keys[0]])


  # best
  plt.subplot(1, 2, 1)
  mu,var = gp_predictive_distribution(x_train,t_train,x_test,best_theta,beta)
  mu=np.squeeze(np.asarray(mu.T))
  gp_plot(x_test,y_test,mu,var,x_train,t_train,best_theta,beta)
  plt.title("The best theta: "+str(best_theta))


  # worst
  plt.subplot(1, 2, 2)
  mu,var = gp_predictive_distribution(x_train,t_train,x_test,worst_theta,beta)
  mu=np.squeeze(np.asarray(mu.T))
  gp_plot(x_test,y_test,mu,var,x_train,t_train,worst_theta,beta)
  plt.title("The worst theta: "+ str(worst_theta))

  plt.show()



grid_search(300)
# -

# #### 3.5 Questions (10 points)
# Selecting kernel functions can be somewhat of an art.  There are charateristics of kernel functions that are useful for some data sets, but not others.  Complicating the matter is the ability to combine kernels with different characteristics (long term trends + seasonal fluctuations).  Describe the charactistics of the kernel function we are using in terms of (signal, scale, offsets, etc). You may want to play around with $\thetav$ and see what each parameter does/affects/etc.  (5 points)  Describe why the best parameters work well for the training data and explain why the bad parameter settings perform poorly (in terms of the first part of the question).  (5 points)

# ___answer___
#
# As we showed in 3.2 the kernal function can be decomposed as a combination of other kernel functions, we shall continue our discussion from that point.
# $k(x_n,x_m)= \theta_0 k_1(x_n,x_m) + \theta_2  k_2(x_n,x_m) + \theta_3 k_3(x_n,x_m)$
#
# $k_1$ kernel corresponds to the similarity measure between points (signal capturing). $k_2 $ and $k_3$ operate as offsets and it was observed that it’s possible to achieve a good fit to the true function without the offset. Theta parameters have also some interpretation. $\theta_1$ corresponds to the bandwidth of the unnormalised gaussian distribution, and $\theta_0$ corresponds to the scaling factor of the distribution. Since $k_2(x_n,x_m)=x_n^0x_m^0$, $\theta_2$ plays a role of the offset, and $\theta_3$ is the scaling factor of the linear kernel. 
#
#
# As we stated previously, it’s possible to achieve a good data fit without the offset (i.e. only with a similarity kernel and two thetas), the opposite is also true, i.e. it’s not possible to achieve a good data fit without  the similarity kernel and the corresponding thetas. The worst performance models have been observed to have low values for $\theta_0$ and $\theta_1$, and that implies that it can’t capture the signal. The key aspect in the kernel is obviously the similarity measure between two data points, and say, if both $\theta_0$ and $\theta_1$ are set close to 0, then the kernel loses its power, and the expectation on the plot becomes similar to a constant function.   

# #### 3.6 Bonus: Implementation (20 points)
# Implement gradient-ascent (or descent if you wish) using the combination of a) the log-likelihood objective function and b) the gradients you calculated above.  Run on the training data above and show the log-likehood curve as it learns and a plot of the final model.  Feel free to use available software (eg search for "minimize.py" which uses conjugate gradient descent, or something in scipy).  NB: log-likelihood should be monotonically increasing.  You are encouraged to also search and use "checkgrad".  (20 points)



StarcoderdataJupyterScriptsDedupFiltered