<jupyter_start><jupyter_text>Parametric ML and Bayesian regression    Notebook version: 1.2 (Sep 28, 2018)    Authors: <NAME>             <NAME> (<EMAIL>)             <NAME> (<EMAIL>)     Changes: v.1.0 - First version. Python version             v.1.1 - Python 3 compatibility. ML section.             v.1.2 - Revised content. 2D visualization removed.        Pending changes:<jupyter_code># Import some libraries that will be necessary for working with data and displaying plots

# To visualize plots in the notebook
%matplotlib inline 

import matplotlib
import matplotlib.pyplot as plt
import matplotlib.cm as cm

import numpy as np
import scipy.io       # To read matlab files
from scipy import spatial
import pylab
pylab.rcParams['figure.figsize'] = 8, 5<jupyter_output><empty_output><jupyter_text>1. IntroductionIn this exercise the student will review several key concepts of Maximum Likelihood and Bayesian regression. To do so, we will assume the regression model$$s = f({\bf x}) + \varepsilon$$where $s$ is the output corresponding to input ${\bf x}$, $f({\bf x})$ is an unobservable latent function, and $\varepsilon$ is white zero-mean Gaussian noise, i.e., $$\varepsilon \sim {\cal N}(0,\sigma_\varepsilon^2).$$In addition, we will assume that the latent function is *linear in the parameters*$$f({\bf x}) = {\bf w}^\top {\bf z}$$where ${\bf z} = T({\bf x})$ is a possibly non-linear transformation of the input. Along this notebook, we will explore different types of transformations.Also, we will assume an a priori distribution for ${\bf w}$ given by$${\bf w} \sim {\cal N}({\bf 0}, \sigma_p^2~{\bf I})$$  Practical considerations   - Though sometimes unavoidable, it is recommended not to use explicit matrix inversion whenever possible. For instance, if an operation like ${\mathbf A}^{-1} {\mathbf b}$ must be performed, it is preferable to code it using python $\mbox{numpy.linalg.lstsq}$ function (see http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html), which provides the LS solution to the overdetermined system ${\mathbf A} {\mathbf w} = {\mathbf b}$.      - Sometimes, the computation of $\log|{\mathbf A}|$ (where ${\mathbf A}$ is a positive definite matrix) can overflow available precision, producing incorrect results. A numerically more stable alternative, providing the same result is $2\sum_i \log([{\mathbf L}]_{ii})$, where $\mathbf L$ is the Cholesky decomposition of $\mathbf A$ (i.e., ${\mathbf A} = {\mathbf L}^\top {\mathbf L}$), and $[{\mathbf L}]_{ii}$ is the $i$th element of the diagonal of ${\mathbf L}$.      - Non-degenerate covariance matrices, such as the ones in this exercise, are always positive definite. It may happen, as a consequence of chained rounding errors, that a matrix which was mathematically expected to be positive definite, turns out not to be so. This implies its Cholesky decomposition will not be available. A quick way to palliate this problem is by adding a small number (such as $10^{-6}$) to the diagonal of such matrix.  Reproducibility of computationsTo guarantee the exact reproducibility of the experiments, it may be useful to start your code initializing the seed of the random numbers generator, so that you can compare your results with the ones given in this notebook.<jupyter_code>np.random.seed(3)<jupyter_output><empty_output><jupyter_text>2. Data generation with a linear modelDuring this section, we will assume affine transformation$${\bf z} = T({\bf x}) = (1, {\bf x}^\top)^\top$$.The a priori distribution of ${\bf w}$ is assumed to be$${\bf w} \sim {\cal N}({\bf 0}, \sigma_p^2~{\bf I})$$  2.1. Synthetic data generationFirst, we are going to generate synthetic data (so that we have the ground-truth model) and use them to make sure everything works correctly and our estimations are sensible.* [1] Set parameters $\sigma_p^2 = 2$ and $\sigma_{\varepsilon}^2 = 0.2$. To do so, define variables `sigma_p` and `sigma_eps` containing the respectiv standard deviations.<jupyter_code># Parameter settings
# sigma_p = <FILL IN>
# sigma_eps = <FILL IN><jupyter_output><empty_output><jupyter_text>* [2] Generate a weight vector $\mbox{true_w}$ with two elements from the a priori distribution of the weights. This vector determines the regression line that we want to find (i.e., the optimum unknown solution).<jupyter_code># Data dimension:
dim_x = 2

# Generate a parameter vector taking a random sample from the prior distributions
# (the np.random module may be usefull for this purpose)
# true_w = <FILL IN>

print('The true parameter vector is:')
print(true_w)<jupyter_output><empty_output><jupyter_text>* [3] Generate an input matrix ${\bf X}$ (in this case, a single column) containing 20 samples with equally spaced values between 0 and 2 in the second column (method `linspace` from numpy can be useful for this)<jupyter_code># <SOL>
# </SOL><jupyter_output><empty_output><jupyter_text>* [4] Finally, generate the output vector ${\mbox s}$ as the product $\mbox{X} \ast \mbox{true_w}$ plus Gaussian noise of pdf ${\cal N}(0,\sigma_\varepsilon^2)$ at each element.<jupyter_code># Expand input matrix with an all-ones column
col_1 = np.ones((n_points, 1))
# Z = <FILL IN>

# Generate values of the target variable
# s = <FILL IN><jupyter_output><empty_output><jupyter_text>2.2. Data visualization * Plot the generated data. You will notice a linear behavior, but the presence of noise makes it hard to estimate precisely the original straight line that generated them (which is stored in $\mbox{true_w}$).<jupyter_code># <SOL> 
# </SOL><jupyter_output><empty_output><jupyter_text>3. Maximum Likelihood (ML) regression 3.1. Likelihood function * [1] Define a function `predict(we, Z)` that computes the linear predictions for all inputs in data matrix `Z` (a 2-D numpy arry), for a given parameter vector `we` (a 1-D numpy array). The output should be a 1-D array. Test your function with the given dataset and `we = [0.4, 0.7]`<jupyter_code># <SOL>
# </SOL>

# Print predictions
print(p)<jupyter_output><empty_output><jupyter_text>* [2] Define a function `sse(we, Z, s)` that computes the sum of squared errors (SSE) for the linear prediction with parameters `we ` (1D numpy array), inputs `Z `  (2D numpy array) and targets `s ` (1D numpy array). Using this function, compute the SSE of the true parameter vector in `true_w`.<jupyter_code># <SOL>
# </SOL>

print(" The SSE is: {0}".format(SSE))<jupyter_output><empty_output><jupyter_text>* [3] Define a function `likelihood(we, Z, s, sigma_eps)` that computes the likelihood of parameter vector `we` for a given dataset in matrix `Z` and vector `s`, assuming Gaussian noise with varianze $\sigma_\epsilon^2$. Note that this function can use the `sse` function defined above. Using this function, compute the likelihood of the true parameter vector in `true_w`.<jupyter_code># <SOL>
# </SOL>

print("The likelihood of the true parameter vector is {0}".format(L_w_true))<jupyter_output><empty_output><jupyter_text>* [4] Define a function `LL(we, Xe, s)` that computes the log-likelihood of parameter vector `we` for a given dataset in matrix `Z` and vector `s`. Note that this function can use the `likelihood` function defined above. However, for a highe numerical precission, implemening a direct expression for the log-likelihood is recommended.    Using this function, compute the likelihood of the true parameter vector in `true_w`.<jupyter_code># <SOL>
# </SOL>

print("The log-likelihood of the true parameter vector is {0}".format(LL_w_true))<jupyter_output><empty_output><jupyter_text>3.2. ML estimate* [1] Compute the ML estimate of $w_e$ given the data.<jupyter_code># <SOL>
# </SOL>

print(w_ML)<jupyter_output><empty_output><jupyter_text>* [2] Compute the maximum likelihood, and the maximum log-likelihood.<jupyter_code># <SOL>
# </SOL>

print('Maximum likelihood: {0}'.format(L_w_ML))
print('Maximum log-likelihood: {0}'.format(LL_w_ML))<jupyter_output><empty_output><jupyter_text>Just as an illustration, the code below generates a set of points in a two dimensional grid going from $(-\sigma_p, -\sigma_p)$ to $(\sigma_p, \sigma_p)$, computes the log-likelihood for all these points and visualize them using a 2-dimensional plot. You can see the difference between the true value of the parameter ${\bf w}$ (black) and the ML estimate (red).<jupyter_code># First construct a grid of (theta0, theta1) parameter pairs and their
# corresponding cost function values.
N = 200    # Number of points along each dimension.
w0_grid = np.linspace(-2.5*sigma_p, 2.5*sigma_p, N)
w1_grid = np.linspace(-2.5*sigma_p, 2.5*sigma_p, N)

Lw = np.zeros((N,N))
# Fill Lw with the likelihood values
for i, w0i in enumerate(w0_grid):
    for j, w1j in enumerate(w1_grid):
        we = np.array((w0i, w1j)) 
        Lw[i, j] = LL(we, Z, s, sigma_eps)

WW0, WW1 = np.meshgrid(w0_grid, w1_grid, indexing='ij')
contours = plt.contour(WW0, WW1, Lw, 20)

plt.figure
plt.clabel(contours)
plt.scatter([true_w[0]]*2, [true_w[1]]*2, s=[50,10], color=['k','w'])
plt.scatter([w_ML[0]]*2, [w_ML[1]]*2, s=[50,10], color=['r','w'])
plt.xlabel('$w_0$')
plt.ylabel('$w_1$')
plt.show()<jupyter_output><empty_output><jupyter_text>3.3. Convergence of the ML estimate for the true modelNote that the likelihood of the true parameter vector is, in general, smaller than that of the ML estimate. However, as the sample size increasis, both should converge to the same value.* [1] Generate a longer dataset, with $K_\text{max}=2^{16}$ samples, uniformly spaces between 0 and 2. Store it in the 2D-array `X2` and the 1D-array `s2`<jupyter_code># Parameter settings
x_min = 0
x_max = 2
n_points = 2**16

# <SOL>
# </SOL><jupyter_output><empty_output><jupyter_text>* [2] Compute the ML estimate based on the first $2^k$ samples, for $k=2,3,\ldots, 16$. For each value of $k$ compute the squared euclidean distance between the true parameter vector and the ML estimate. Represent it graphically (using a logarithmic scale in the y-axis).<jupyter_code># <SOL>
# </SOL><jupyter_output><empty_output><jupyter_text>4. ML estimation with real data. The stocks dataset.Once our code has been tested on synthetic data, we will use it with real data.   4.1. Dataset * [1] Load data corresponding to the evolution of the stocks of 10 airline companies. This data set is an adaptation of the Stock dataset from http://www.dcc.fc.up.pt/~ltorgo/Regression/DataSets.html, which in turn was taken from the StatLib Repository, http://lib.stat.cmu.edu/<jupyter_code># <SOL>
# </SOL><jupyter_output><empty_output><jupyter_text>* [2] Normalize the data so all training sample components have zero mean and unit standard deviation. Store the normalized training and test samples in 2D numpy arrays `Xtrain` and `Xtest`, respectively.<jupyter_code># <SOL>
# </SOL><jupyter_output><empty_output><jupyter_text>4.2. Polynomial ML regression with a single variableIn this first part, we will work with the first component of the input only. * [1] Take the first column of `Xtrain` and `Xtest` into arrays `X0train` and `X0test`, respectively.  * [2] Visualize, in a single scatter plot, the target variable (in the vertical axes) versus the input variable) * [3] Since the data have been taken from a real scenario, we do not have any *true* mathematical model of the process that generated the data. Thus, we will explore different models trying to take the one that fits better the training data.  Assume a polinomial model given by  $$  {\bf z} = T({\bf x}) = (1, x_0, x_0^2, \ldots, x_0^{g-1})^\top.  $$  Built a method       `    Ztrain, Ztest = T_poly(Xtrain, Xtest, g)`  that, for a given value of $g$, computes normailzed data matrices `Ztrain` and `Ztest` that result from applying the polinomial transformation to the inputs in `X0train` and `X0test` for an arbitrary value of $g$.    Note that, despite `X0train` and `X0test` where normalized, you will need to re-normalize the transformed variables.<jupyter_code># <SOL>
# </SOL><jupyter_output><empty_output><jupyter_text>* [4] Fit a polynomial model with degree $g$ for $g$ ranging from 0 to 10. Store the weights of all models in a list of weight vectors, named `models`, such that `models[g]` returns the parameters estimated for the polynomial model with degree $g$.   We will use these models in the following sections. * [5] Plot the polynomial models with degrees 1, 3 and 10, superimposed over a scatter plot of the training data (in blue) and the test data (in red). * [6] Show, in the same plot:  - The log-likelihood function corresponding to each model, as a function of $g$, computed over the training set.  - The log-likelihood function corresponding to each model, as a function of $g$, computed over the test set. * [7] [OPTIONAL] You may have seen the the likelihood function grows with the degree of the polynomial. However, large values of $g$ produce a strong data overfitting. For this reasong, $g$ cannot be selected with the same data used to fit the model.   This kind of parameters, like $g$ are usually called *hyperparameters* and need to be selected by cross validation.    Another hyperparameter is $\sigma_\epsilon^2$. Plot the log-likelihood function corresponding to the polynomial model with degree 3 for different values of $\sigma_\epsilon^2$, for the training set and the test set. What would be the optimal value of this hyperparameters according to the training set?    In any case, not that  the model coefficients do not depend on $\sigma_eps^2$. Therefore, we do not need to estimat its value for ML regression. * [8] Select the optimal value of $g$ by cross-validation. To do so, the cross validation methods provided by sklearn will simplify this task. * [9] For the selected model:    - Plot the regresion function over the scater plot of the data.    - Compute the log-likelihood and the SSE over the test set.  5. Bayesian regression. The stock dataset.In this section we will keep using the first component of the data from the stock dataset, assuming the same kind of plolynomial model. We will explore the potential advantages of using a Bayesian model. To do so, we will asume that the a priori distribution of ${\bf w}$ is$${\bf w} \sim {\cal N}({\bf 0}, \sigma_p^2~{\bf I})$$  5.1. Posterior pdf of the weight vectorIn this section we will visualize prior and the posterior distribution functions. First, we will restore the dataset at the begining of this notebook: * [1] Define a function `posterior_stats(Z, s, sigma_eps, sigma_p)` that computes the parameters of the posterior coefficient distribution given the dataset in matrix `Z` and vector `s`, for given values of the hyperparameters.This function should return the posterior mean, the covariance matrix and the precision matrix (the inverse of the covariance matrix). Test the function to the given dataset, for $g=3$.<jupyter_code># <SOL>
# </SOL>

mean_w, Cov_w, iCov_w = posterior_stats(Z, s, sigma_eps, sigma_p)

print('true_w = {0}'.format(true_w))
print('mean_w = {0}'.format(mean_w))
print('Cov_w = {0}'.format(Cov_w))
print('iCov_w = {0}'.format(iCov_w))<jupyter_output><empty_output><jupyter_text>* [2] Define a function `gauss_pdf(we, mean_w, iCov_w)` that computes the Gaussian pdf with mean `mean_w` and precision matrix `iCov_w`. Use this function to compute and compare the posterior pdf value of the true coefficients, the ML estimate and the MSE estimate, given the dataset.<jupyter_code># <SOL>
# </SOL>

print('p(true_w | s)  = {0}'.format(gauss_pdf(true_w, mean_w, iCov_w)))
print('p(w_ML | s)  = {0}'.format(gauss_pdf(w_ML, mean_w, iCov_w)))
print('p(w_MSE | s)  = {0}'.format(gauss_pdf(mean_w, mean_w, iCov_w)))<jupyter_output><empty_output><jupyter_text>* [3] Define a function `log_gauss_pdf(we, mean_w, iCov_w)` that computes the log of the Gaussian pdf with mean `mean_w` and precision matrix `iCov_w`. Use this function to compute and compare the log of the posterior pdf value of the true coefficients, the ML estimate and the MSE estimate, given the dataset.  5.2. Hyperparameter selectionSince the values $\sigma_p$ and $\sigma_\varepsilon$ are no longer known, we have to select them in some way. To see their influence, assume $g=3$ and plot the regression function for different values of $\sigma_p$To this end, we will adjust them using the LS solution to the regression problem:   - $\sigma_p^2$ will be taken as the average of the square values of ${\hat {\bf w}}_{LS}$   - $\sigma_\varepsilon^2$ will be taken as two times the average of the square of the residuals when using ${\hat {\bf w}}_{LS}$  5.1. Hyperparameter selectionSince the values $\sigma_p$ and $\sigma_\varepsilon$ are no longer known, a first rough estimation is needed (we will soon see how to estimate these values in a principled way).To this end, we will adjust them using the LS solution to the regression problem:   - $\sigma_p^2$ will be taken as the average of the square values of ${\hat {\bf w}}_{LS}$   - $\sigma_\varepsilon^2$ will be taken as two times the average of the square of the residuals when using ${\hat {\bf w}}_{LS}$<jupyter_code># w_LS, residuals, rank, s = <FILL IN>
# sigma_p = <FILL IN>
# sigma_eps = <FILL IN>

print(sigma_eps)<jupyter_output><empty_output><jupyter_text>5.2. Sampling regression curves from the posteriorIn this section we will plot the functions corresponding to different samples drawn from the posterior distribution of the weight vector. To this end, we will first generate an input dataset of equally spaced samples. We will compute the functions at these points<jupyter_code># Definition of the interval for representation purposes 
x2_min = -1
x2_max = 3
n_points = 100   # Only two points are needed to plot a straigh line

# Build the input data matrix:
# Input values for representation of the regression curves
X2 = np.linspace(x2_min, x2_max, n_points)
col_1 = np.ones((n_points,))
X2e = np.vstack((col_1, X2)).T<jupyter_output><empty_output><jupyter_text>Generate random vectors ${\bf w}_l$ with $l = 1,\dots, 50$, from the posterior density of the weights, $p({\bf w}\mid{\bf s})$, and use them to generate 50 polinomial regression functions, $f({\bf x}^\ast) = {{\bf z}^\ast}^\top {\bf w}_l$, with ${\bf x}^\ast$ between $-1.2$ and $1.2$, with step $0.1$.Plot the line corresponding to the model with the posterior mean parameters, along with the $50$ generated straight lines and the original samples, all in the same plot. As you can check, the Bayesian model is not providing a single answer, but instead a density over them, from which we have extracted 50 options.<jupyter_code># Drawing weights from the posterior
# First, compute the cholesky decomposition of the covariance matrix
# L = <FILL IN>

for l in range(50):
    # Generate a random sample from the posterior distribution
    # w_l = <FILL IN>

    # Compute predictions for the inputs in the data matrix
    # p_l = <FILL IN>

    # Plot prediction function
    # plt.plot(<FILL IN>, 'c:');

# Compute predictions for the inputs in the data matrix and using the true model
# p_truew = <FILL IN>

# Plot the true model
plt.plot(X2, p_truew, 'b', label='True model', linewidth=2);

# Plot the training points
plt.plot(X,s,'r.',markersize=12);

plt.xlim((x2_min,x2_max));
plt.legend(loc='best')
plt.xlabel('$x$',fontsize=14);
plt.ylabel('$s$',fontsize=14);<jupyter_output><empty_output><jupyter_text>5.3. Plotting the confidence intervalsOn top of the previous figure (copy here your code from the previous section), plot functions$${\mathbb E}\left\{f({\bf x}^\ast)\mid{\bf s}\right\}$$and$${\mathbb E}\left\{f({\bf x}^\ast)\mid{\bf s}\right\} \pm 2 \sqrt{{\mathbb V}\left\{f({\bf x}^\ast)\mid{\bf s}\right\}}$$(i.e., the posterior mean of $f({\bf x}^\ast)$, as well as two standard deviations above and below).It is possible to show analytically that this region comprises $95.45\%$ probability of the posterior probability $p(f({\bf x}^\ast)\mid {\bf s})$ at each ${\bf x}^\ast$.<jupyter_code># Note that you can re-use code from sect. 4.2 to solve this exercise

# Plot sample functions from the posterior, and the training points
# <SOL>
# </SOL>    

# Plot the posterior mean.
# mean_ast = <FILL IN>
plt.plot(X2, mean_ast, 'm', label='Predictive mean', linewidth=2);

# Plot the posterior mean \pm 2 std
# std_ast = <FILL IN>
# plt.plot(<FILL IN>, 'm--', label='Predictive mean $\pm$ 2std', linewidth=2);
# plt.plot(<FILL IN>, 'm--', linewidth=3);


plt.legend(loc='best')
plt.xlabel('$x$',fontsize=14);
plt.ylabel('$s$',fontsize=14);<jupyter_output><empty_output><jupyter_text>Plot now ${\mathbb E}\left\{s({\bf x}^\ast)\mid{\bf s}\right\} \pm 2 \sqrt{{\mathbb V}\left\{s({\bf x}^\ast)\mid{\bf s}\right\}}$ (note that the posterior means of $f({\bf x}^\ast)$ and $s({\bf x}^\ast)$ are the same, so there is no need to plot it again). Notice that $95.45\%$ of observed data lie now within the newly designated region. These new limits establish a confidence range for our predictions. See how the uncertainty grows as we move away from the interpolation region to the extrapolation areas.<jupyter_code># Plot sample functions confidence intervals and sampling points
# Note that you can simply copy and paste most of the code used in the cell above.

# <SOL>
# </SOL>    

# Compute the standad deviations for s and plot the confidence intervals
# <SOL>
# </SOL>

plt.legend(loc='best')
plt.xlabel('$x$',fontsize=14);
plt.ylabel('$s$',fontsize=14);<jupyter_output><empty_output><jupyter_text>5.3. Model assessment[OPTIONAL. You can skip this section]In order to verify the performance of the resulting model, compute the posterior mean and variance of each of the test outputs from the posterior over ${\bf w}$. I.e, compute ${\mathbb E}\left\{s({\bf x}^\ast)\mid{\bf s}\right\}$ and $\sqrt{{\mathbb V}\left\{s({\bf x}^\ast)\mid{\bf s}\right\}}$ for each test sample ${\bf x}^\ast$ contained in each row of `Xtest`. Be sure not to use the outputs `Ytest` at any point during this process.Store the predictive mean and variance of all test samples in two column vectors called `m_y` and `v_y`, respectively.<jupyter_code># <SOL>
# </SOL><jupyter_output><empty_output><jupyter_text>Compute now the mean square error (MSE) and the negative log-predictive density (NLPD) with the following code:<jupyter_code># <SOL>
# </SOL><jupyter_output><empty_output><jupyter_text>Results should be:<jupyter_code>print('MSE = {0}'.format(MSE))
print('NLPD = {0}'.format(NLPD))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11313347608850821, 'fraction_numerical': 0.013160242683797287, 'mean_word_length': 4.77608863694924, 'pattern_counts': {'":': 0, '<': 179, '<?xml version=': 0, '>': 179, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '29', 'dir': 'jupyter-structured-clean-dedup', 'id': '44971', 'license': 'MIT', 'n_tokens_mistral': 7011, 'n_tokens_neox': 6687, 'n_words': 2933, 'path': 'R_lab1_ML_Bay_Regresion/old/Pract_regression_student.ipynb', 'repo_name': 'ML4DS/ML4all', 'size': '22480'}