diff --git a/doc/pub/Regression/ipynb/Regression.ipynb b/doc/pub/Regression/ipynb/Regression.ipynb index 452d9d097..e3e5908cd 100644 --- a/doc/pub/Regression/ipynb/Regression.ipynb +++ b/doc/pub/Regression/ipynb/Regression.ipynb @@ -10,7 +10,7 @@ " \n", "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n", "\n", - "Date: **Sep 7, 2020**\n", + "Date: **Sep 10, 2020**\n", "\n", "Copyright 1999-2020, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n", "\n", @@ -1895,6 +1895,182 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Reducing the number of degrees of freedom, overarching view\n", + "\n", + "Many Machine Learning problems involve thousands or even millions of\n", + "features for each training instance. Not only does this make training\n", + "extremely slow, it can also make it much harder to find a good\n", + "solution, as we will see. This problem is often referred to as the\n", + "curse of dimensionality. Fortunately, in real-world problems, it is\n", + "often possible to reduce the number of features considerably, turning\n", + "an intractable problem into a tractable one.\n", + "\n", + "Here we will discuss some of the most popular dimensionality reduction\n", + "techniques: the principal component analysis (PCA), Kernel PCA, and\n", + "Locally Linear Embedding (LLE). Furthermore, we will start by looking\n", + "at some simple preprocessing of the data which allow us to rescale the\n", + "data.\n", + "\n", + "Principal component analysis and its various variants deal with the\n", + "problem of fitting a low-dimensional [affine\n", + "subspace](https://en.wikipedia.org/wiki/Affine_space) to a set of of\n", + "data points in a high-dimensional space. With its family of methods it\n", + "is one of the most used tools in data modeling, compression and\n", + "visualization.\n", + "\n", + "\n", + "\n", + "\n", + "## Preprocessing our data\n", + "\n", + "Before we proceed however, we will discuss how to preprocess our\n", + "data. Till now and in connection with our previous examples we have\n", + "not met so many cases where we are too sensitive to the scaling of our\n", + "data. Normally the data may need a rescaling and/or may be sensitive\n", + "to extreme values. Scaling the data renders our inputs much more\n", + "suitable for the algorithms we want to employ.\n", + "\n", + "**Scikit-Learn** has several functions which allow us to rescale the\n", + "data, normally resulting in much better results in terms of various\n", + "accuracy scores. The **StandardScaler** function in **Scikit-Learn**\n", + "ensures that for each feature/predictor we study the mean value is\n", + "zero and the variance is one (every column in the design/feature\n", + "matrix). This scaling has the drawback that it does not ensure that\n", + "we have a particular maximum or minimum in our data set. Another\n", + "function included in **Scikit-Learn** is the **MinMaxScaler** which\n", + "ensures that all features are exactly between $0$ and $1$. The\n", + "\n", + "## More preprocessing\n", + "\n", + "\n", + "The **Normalizer** scales each data\n", + "point such that the feature vector has a euclidean length of one. In other words, it\n", + "projects a data point on the circle (or sphere in the case of higher dimensions) with a\n", + "radius of 1. This means every data point is scaled by a different number (by the\n", + "inverse of it’s length).\n", + "This normalization is often used when only the direction (or angle) of the data matters,\n", + "not the length of the feature vector.\n", + "\n", + "The **RobustScaler** works similarly to the StandardScaler in that it\n", + "ensures statistical properties for each feature that guarantee that\n", + "they are on the same scale. However, the RobustScaler uses the median\n", + "and quartiles, instead of mean and variance. This makes the\n", + "RobustScaler ignore data points that are very different from the rest\n", + "(like measurement errors). These odd data points are also called\n", + "outliers, and might often lead to trouble for other scaling\n", + "techniques.\n", + "\n", + "\n", + "\n", + "## Simple preprocessing examples, Franke function and regression" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Common imports\n", + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import sklearn.linear_model as skl\n", + "from sklearn.metrics import mean_squared_error\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.preprocessing import MinMaxScaler, StandardScaler, Normalizer\n", + "\n", + "# Where to save the figures and data files\n", + "PROJECT_ROOT_DIR = \"Results\"\n", + "FIGURE_ID = \"Results/FigureFiles\"\n", + "DATA_ID = \"DataFiles/\"\n", + "\n", + "if not os.path.exists(PROJECT_ROOT_DIR):\n", + " os.mkdir(PROJECT_ROOT_DIR)\n", + "\n", + "if not os.path.exists(FIGURE_ID):\n", + " os.makedirs(FIGURE_ID)\n", + "\n", + "if not os.path.exists(DATA_ID):\n", + " os.makedirs(DATA_ID)\n", + "\n", + "def image_path(fig_id):\n", + " return os.path.join(FIGURE_ID, fig_id)\n", + "\n", + "def data_path(dat_id):\n", + " return os.path.join(DATA_ID, dat_id)\n", + "\n", + "def save_fig(fig_id):\n", + " plt.savefig(image_path(fig_id) + \".png\", format='png')\n", + "\n", + "\n", + "def FrankeFunction(x,y):\n", + "\tterm1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))\n", + "\tterm2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))\n", + "\tterm3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))\n", + "\tterm4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)\n", + "\treturn term1 + term2 + term3 + term4\n", + "\n", + "\n", + "def create_X(x, y, n ):\n", + "\tif len(x.shape) > 1:\n", + "\t\tx = np.ravel(x)\n", + "\t\ty = np.ravel(y)\n", + "\n", + "\tN = len(x)\n", + "\tl = int((n+1)*(n+2)/2)\t\t# Number of elements in beta\n", + "\tX = np.ones((N,l))\n", + "\n", + "\tfor i in range(1,n+1):\n", + "\t\tq = int((i)*(i+1)/2)\n", + "\t\tfor k in range(i+1):\n", + "\t\t\tX[:,q+k] = (x**(i-k))*(y**k)\n", + "\n", + "\treturn X\n", + "\n", + "\n", + "# Making meshgrid of datapoints and compute Franke's function\n", + "n = 5\n", + "N = 1000\n", + "x = np.sort(np.random.uniform(0, 1, N))\n", + "y = np.sort(np.random.uniform(0, 1, N))\n", + "z = FrankeFunction(x, y)\n", + "X = create_X(x, y, n=n) \n", + "# split in training and test data\n", + "X_train, X_test, y_train, y_test = train_test_split(X,z,test_size=0.2)\n", + "\n", + "\n", + "clf = skl.LinearRegression().fit(X_train, y_train)\n", + "\n", + "# The mean squared error and R2 score\n", + "print(\"MSE before scaling: {:.2f}\".format(mean_squared_error(clf.predict(X_test), y_test)))\n", + "print(\"R2 score before scaling {:.2f}\".format(clf.score(X_test,y_test)))\n", + "\n", + "scaler = StandardScaler()\n", + "scaler.fit(X_train)\n", + "X_train_scaled = scaler.transform(X_train)\n", + "X_test_scaled = scaler.transform(X_test)\n", + "\n", + "print(\"Feature min values before scaling:\\n {}\".format(X_train.min(axis=0)))\n", + "print(\"Feature max values before scaling:\\n {}\".format(X_train.max(axis=0)))\n", + "\n", + "print(\"Feature min values after scaling:\\n {}\".format(X_train_scaled.min(axis=0)))\n", + "print(\"Feature max values after scaling:\\n {}\".format(X_train_scaled.max(axis=0)))\n", + "\n", + "clf = skl.LinearRegression().fit(X_train_scaled, y_train)\n", + "\n", + "\n", + "print(\"MSE after scaling: {:.2f}\".format(mean_squared_error(clf.predict(X_test_scaled), y_test)))\n", + "print(\"R2 score for scaled data: {:.2f}\".format(clf.score(X_test_scaled,y_test)))" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -2130,7 +2306,7 @@ "two orthogonal/unitary matrices. The [Singular Value Decompostion\n", "(SVD) theorem](https://en.wikipedia.org/wiki/Singular_value_decomposition)\n", "states that a general $m\\times n$ matrix $\\boldsymbol{X}$ can be written in\n", - "terms of a diagonal matrix $\\boldsymbol{\\Sigma}$ of dimensionality $n\\times n$\n", + "terms of a diagonal matrix $\\boldsymbol{\\Sigma}$ of dimensionality $m\\times n$\n", "and two orthognal matrices $\\boldsymbol{U}$ and $\\boldsymbol{V}$, where the first has\n", "dimensionality $m \\times m$ and the last dimensionality $n\\times n$.\n", "We have then" @@ -2664,7 +2840,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 23, "metadata": { "collapsed": false }, @@ -2757,9 +2933,618 @@ "\n", "The matrices $\\boldsymbol{U}$ and $\\boldsymbol{V}$ are unitary/orthonormal matrices, that is in case the matrices are real we have $\\boldsymbol{U}^T\\boldsymbol{U}=\\boldsymbol{U}\\boldsymbol{U}^T=\\boldsymbol{I}$ and $\\boldsymbol{V}^T\\boldsymbol{V}=\\boldsymbol{V}\\boldsymbol{V}^T=\\boldsymbol{I}$.\n", "\n", - "## Spectral Decomposition of the OLS\n", "\n", - "More material to be added here\n", + "\n", + "## Introducing the Covariance and Correlation functions\n", + "\n", + "Before we discuss the link between for example Ridge regression and the singular value decomposition, we need to remind ourselves about\n", + "the definition of the covariance and the correlation function. These are quantities \n", + "\n", + "Suppose we have defined two vectors\n", + "$\\hat{x}$ and $\\hat{y}$ with $n$ elements each. The covariance matrix $\\boldsymbol{C}$ is defined as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{C}[\\boldsymbol{x},\\boldsymbol{y}] = \\begin{bmatrix} \\mathrm{cov}[\\boldsymbol{x},\\boldsymbol{x}] & \\mathrm{cov}[\\boldsymbol{x},\\boldsymbol{y}] \\\\\n", + " \\mathrm{cov}[\\boldsymbol{y},\\boldsymbol{x}] & \\mathrm{cov}[\\boldsymbol{y},\\boldsymbol{y}] \\\\\n", + " \\end{bmatrix},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where for example" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\mathrm{cov}[\\boldsymbol{x},\\boldsymbol{y}] =\\frac{1}{n} \\sum_{i=0}^{n-1}(x_i- \\overline{x})(y_i- \\overline{y}).\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With this definition and recalling that the variance is defined as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\mathrm{var}[\\boldsymbol{x}]=\\frac{1}{n} \\sum_{i=0}^{n-1}(x_i- \\overline{x})^2,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "we can rewrite the covariance matrix as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{C}[\\boldsymbol{x},\\boldsymbol{y}] = \\begin{bmatrix} \\mathrm{var}[\\boldsymbol{x}] & \\mathrm{cov}[\\boldsymbol{x},\\boldsymbol{y}] \\\\\n", + " \\mathrm{cov}[\\boldsymbol{x},\\boldsymbol{y}] & \\mathrm{var}[\\boldsymbol{y}] \\\\\n", + " \\end{bmatrix}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The covariance takes values between zero and infinity and may thus\n", + "lead to problems with loss of numerical precision for particularly\n", + "large values. It is common to scale the covariance matrix by\n", + "introducing instead the correlation matrix defined via the so-called\n", + "correlation function" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\mathrm{corr}[\\boldsymbol{x},\\boldsymbol{y}]=\\frac{\\mathrm{cov}[\\boldsymbol{x},\\boldsymbol{y}]}{\\sqrt{\\mathrm{var}[\\boldsymbol{x}] \\mathrm{var}[\\boldsymbol{y}]}}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The correlation function is then given by values $\\mathrm{corr}[\\boldsymbol{x},\\boldsymbol{y}]\n", + "\\in [-1,1]$. This avoids eventual problems with too large values. We\n", + "can then define the correlation matrix for the two vectors $\\boldsymbol{x}$\n", + "and $\\boldsymbol{y}$ as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{K}[\\boldsymbol{x},\\boldsymbol{y}] = \\begin{bmatrix} 1 & \\mathrm{corr}[\\boldsymbol{x},\\boldsymbol{y}] \\\\\n", + " \\mathrm{corr}[\\boldsymbol{y},\\boldsymbol{x}] & 1 \\\\\n", + " \\end{bmatrix},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the above example this is the function we constructed using **pandas**.\n", + "\n", + "## Correlation Function and Design/Feature Matrix\n", + "\n", + "In our derivation of the various regression algorithms like **Ordinary Least Squares** or **Ridge regression**\n", + "we defined the design/feature matrix $\\boldsymbol{X}$ as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{X}=\\begin{bmatrix}\n", + "x_{0,0} & x_{0,1} & x_{0,2}& \\dots & \\dots x_{0,p-1}\\\\\n", + "x_{1,0} & x_{1,1} & x_{1,2}& \\dots & \\dots x_{1,p-1}\\\\\n", + "x_{2,0} & x_{2,1} & x_{2,2}& \\dots & \\dots x_{2,p-1}\\\\\n", + "\\dots & \\dots & \\dots & \\dots \\dots & \\dots \\\\\n", + "x_{n-2,0} & x_{n-2,1} & x_{n-2,2}& \\dots & \\dots x_{n-2,p-1}\\\\\n", + "x_{n-1,0} & x_{n-1,1} & x_{n-1,2}& \\dots & \\dots x_{n-1,p-1}\\\\\n", + "\\end{bmatrix},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "with $\\boldsymbol{X}\\in {\\mathbb{R}}^{n\\times p}$, with the predictors/features $p$ refering to the column numbers and the\n", + "entries $n$ being the row elements.\n", + "We can rewrite the design/feature matrix in terms of its column vectors as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{X}=\\begin{bmatrix} \\boldsymbol{x}_0 & \\boldsymbol{x}_1 & \\boldsymbol{x}_2 & \\dots & \\dots & \\boldsymbol{x}_{p-1}\\end{bmatrix},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "with a given vector" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{x}_i^T = \\begin{bmatrix}x_{0,i} & x_{1,i} & x_{2,i}& \\dots & \\dots x_{n-1,i}\\end{bmatrix}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With these definitions, we can now rewrite our $2\\times 2$\n", + "correaltion/covariance matrix in terms of a moe general design/feature\n", + "matrix $\\boldsymbol{X}\\in {\\mathbb{R}}^{n\\times p}$. This leads to a $p\\times p$\n", + "covariance matrix for the vectors $\\boldsymbol{x}_i$ with $i=0,1,\\dots,p-1$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{C}[\\boldsymbol{x}] = \\begin{bmatrix}\n", + "\\mathrm{var}[\\boldsymbol{x}_0] & \\mathrm{cov}[\\boldsymbol{x}_0,\\boldsymbol{x}_1] & \\mathrm{cov}[\\boldsymbol{x}_0,\\boldsymbol{x}_2] & \\dots & \\dots & \\mathrm{cov}[\\boldsymbol{x}_0,\\boldsymbol{x}_{p-1}]\\\\\n", + "\\mathrm{cov}[\\boldsymbol{x}_1,\\boldsymbol{x}_0] & \\mathrm{var}[\\boldsymbol{x}_1] & \\mathrm{cov}[\\boldsymbol{x}_1,\\boldsymbol{x}_2] & \\dots & \\dots & \\mathrm{cov}[\\boldsymbol{x}_1,\\boldsymbol{x}_{p-1}]\\\\\n", + "\\mathrm{cov}[\\boldsymbol{x}_2,\\boldsymbol{x}_0] & \\mathrm{cov}[\\boldsymbol{x}_2,\\boldsymbol{x}_1] & \\mathrm{var}[\\boldsymbol{x}_2] & \\dots & \\dots & \\mathrm{cov}[\\boldsymbol{x}_2,\\boldsymbol{x}_{p-1}]\\\\\n", + "\\dots & \\dots & \\dots & \\dots & \\dots & \\dots \\\\\n", + "\\dots & \\dots & \\dots & \\dots & \\dots & \\dots \\\\\n", + "\\mathrm{cov}[\\boldsymbol{x}_{p-1},\\boldsymbol{x}_0] & \\mathrm{cov}[\\boldsymbol{x}_{p-1},\\boldsymbol{x}_1] & \\mathrm{cov}[\\boldsymbol{x}_{p-1},\\boldsymbol{x}_{2}] & \\dots & \\dots & \\mathrm{var}[\\boldsymbol{x}_{p-1}]\\\\\n", + "\\end{bmatrix},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and the correlation matrix" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{K}[\\boldsymbol{x}] = \\begin{bmatrix}\n", + "1 & \\mathrm{corr}[\\boldsymbol{x}_0,\\boldsymbol{x}_1] & \\mathrm{corr}[\\boldsymbol{x}_0,\\boldsymbol{x}_2] & \\dots & \\dots & \\mathrm{corr}[\\boldsymbol{x}_0,\\boldsymbol{x}_{p-1}]\\\\\n", + "\\mathrm{corr}[\\boldsymbol{x}_1,\\boldsymbol{x}_0] & 1 & \\mathrm{corr}[\\boldsymbol{x}_1,\\boldsymbol{x}_2] & \\dots & \\dots & \\mathrm{corr}[\\boldsymbol{x}_1,\\boldsymbol{x}_{p-1}]\\\\\n", + "\\mathrm{corr}[\\boldsymbol{x}_2,\\boldsymbol{x}_0] & \\mathrm{corr}[\\boldsymbol{x}_2,\\boldsymbol{x}_1] & 1 & \\dots & \\dots & \\mathrm{corr}[\\boldsymbol{x}_2,\\boldsymbol{x}_{p-1}]\\\\\n", + "\\dots & \\dots & \\dots & \\dots & \\dots & \\dots \\\\\n", + "\\dots & \\dots & \\dots & \\dots & \\dots & \\dots \\\\\n", + "\\mathrm{corr}[\\boldsymbol{x}_{p-1},\\boldsymbol{x}_0] & \\mathrm{corr}[\\boldsymbol{x}_{p-1},\\boldsymbol{x}_1] & \\mathrm{corr}[\\boldsymbol{x}_{p-1},\\boldsymbol{x}_{2}] & \\dots & \\dots & 1\\\\\n", + "\\end{bmatrix},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Covariance Matrix Examples\n", + "\n", + "\n", + "The Numpy function **np.cov** calculates the covariance elements using\n", + "the factor $1/(n-1)$ instead of $1/n$ since it assumes we do not have\n", + "the exact mean values. The following simple function uses the\n", + "**np.vstack** function which takes each vector of dimension $1\\times n$\n", + "and produces a $2\\times n$ matrix $\\boldsymbol{W}$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{W} = \\begin{bmatrix} x_0 & y_0 \\\\\n", + " x_1 & y_1 \\\\\n", + " x_2 & y_2\\\\\n", + " \\dots & \\dots \\\\\n", + " x_{n-2} & y_{n-2}\\\\\n", + " x_{n-1} & y_{n-1} & \n", + " \\end{bmatrix},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which in turn is converted into into the $2\\times 2$ covariance matrix\n", + "$\\boldsymbol{C}$ via the Numpy function **np.cov()**. We note that we can also calculate\n", + "the mean value of each set of samples $\\boldsymbol{x}$ etc using the Numpy\n", + "function **np.mean(x)**. We can also extract the eigenvalues of the\n", + "covariance matrix through the **np.linalg.eig()** function." + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Importing various packages\n", + "import numpy as np\n", + "n = 100\n", + "x = np.random.normal(size=n)\n", + "print(np.mean(x))\n", + "y = 4+3*x+np.random.normal(size=n)\n", + "print(np.mean(y))\n", + "W = np.vstack((x, y))\n", + "C = np.cov(W)\n", + "print(C)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Correlation Matrix\n", + "\n", + "The previous example can be converted into the correlation matrix by\n", + "simply scaling the matrix elements with the variances. We should also\n", + "subtract the mean values for each column. This leads to the following\n", + "code which sets up the correlations matrix for the previous example in\n", + "a more brute force way. Here we scale the mean values for each column of the design matrix, calculate the relevant mean values and variances and then finally set up the $2\\times 2$ correlation matrix (since we have only two vectors)." + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "n = 100\n", + "# define two vectors \n", + "x = np.random.random(size=n)\n", + "y = 4+3*x+np.random.normal(size=n)\n", + "#scaling the x and y vectors \n", + "x = x - np.mean(x)\n", + "y = y - np.mean(y)\n", + "variance_x = np.sum(x@x)/n\n", + "variance_y = np.sum(y@y)/n\n", + "print(variance_x)\n", + "print(variance_y)\n", + "cov_xy = np.sum(x@y)/n\n", + "cov_xx = np.sum(x@x)/n\n", + "cov_yy = np.sum(y@y)/n\n", + "C = np.zeros((2,2))\n", + "C[0,0]= cov_xx/variance_x\n", + "C[1,1]= cov_yy/variance_y\n", + "C[0,1]= cov_xy/np.sqrt(variance_y*variance_x)\n", + "C[1,0]= C[0,1]\n", + "print(C)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that the matrix elements along the diagonal are one as they\n", + "should be and that the matrix is symmetric. Furthermore, diagonalizing\n", + "this matrix we easily see that it is a positive definite matrix.\n", + "\n", + "The above procedure with **numpy** can be made more compact if we use **pandas**.\n", + "\n", + "## Correlation Matrix with Pandas\n", + "\n", + "We whow here how we can set up the correlation matrix using **pandas**, as done in this simple code" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "n = 10\n", + "x = np.random.normal(size=n)\n", + "x = x - np.mean(x)\n", + "y = 4+3*x+np.random.normal(size=n)\n", + "y = y - np.mean(y)\n", + "X = (np.vstack((x, y))).T\n", + "print(X)\n", + "Xpd = pd.DataFrame(X)\n", + "print(Xpd)\n", + "correlation_matrix = Xpd.corr()\n", + "print(correlation_matrix)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We expand this model to the Franke function discussed above.\n", + "\n", + "## Correlation Matrix with Pandas and the Franke function" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Common imports\n", + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "\n", + "def FrankeFunction(x,y):\n", + "\tterm1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))\n", + "\tterm2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))\n", + "\tterm3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))\n", + "\tterm4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)\n", + "\treturn term1 + term2 + term3 + term4\n", + "\n", + "\n", + "def create_X(x, y, n ):\n", + "\tif len(x.shape) > 1:\n", + "\t\tx = np.ravel(x)\n", + "\t\ty = np.ravel(y)\n", + "\n", + "\tN = len(x)\n", + "\tl = int((n+1)*(n+2)/2)\t\t# Number of elements in beta\n", + "\tX = np.ones((N,l))\n", + "\n", + "\tfor i in range(1,n+1):\n", + "\t\tq = int((i)*(i+1)/2)\n", + "\t\tfor k in range(i+1):\n", + "\t\t\tX[:,q+k] = (x**(i-k))*(y**k)\n", + "\n", + "\treturn X\n", + "\n", + "\n", + "# Making meshgrid of datapoints and compute Franke's function\n", + "n = 4\n", + "N = 100\n", + "x = np.sort(np.random.uniform(0, 1, N))\n", + "y = np.sort(np.random.uniform(0, 1, N))\n", + "z = FrankeFunction(x, y)\n", + "X = create_X(x, y, n=n) \n", + "\n", + "Xpd = pd.DataFrame(X)\n", + "# subtract the mean values and set up the covariance matrix\n", + "Xpd = Xpd - Xpd.mean()\n", + "covariance_matrix = Xpd.cov()\n", + "print(covariance_matrix)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We note here that the covariance is zero for the first rows and\n", + "columns since all matrix elements in the design matrix were set to one\n", + "(we are fitting the function in terms of a polynomial of degree $n$).\n", + "\n", + "This means that the variance for these elements will be zero and will\n", + "cause problems when we set up the correlation matrix. We can simply\n", + "drop these elements and construct a correlation\n", + "matrix without these elements. \n", + "\n", + "\n", + "## Rewriting the Covariance and/or Correlation Matrix\n", + "\n", + "We can rewrite the covariance matrix in a more compact form in terms of the design/feature matrix $\\boldsymbol{X}$ as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{C}[\\boldsymbol{x}] = \\frac{1}{n}\\boldsymbol{X}\\boldsymbol{X}^T= \\mathbb{E}[\\boldsymbol{X}\\boldsymbol{X}^T].\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To see this let us simply look at a design matrix $\\boldsymbol{X}\\in {\\mathbb{R}}^{2\\times 2}$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{X}=\\begin{bmatrix}\n", + "x_{00} & x_{01}\\\\\n", + "x_{10} & x_{11}\\\\\n", + "\\end{bmatrix}=\\begin{bmatrix}\n", + "\\boldsymbol{x}_{0} & \\boldsymbol{x}_{1}\\\\\n", + "\\end{bmatrix}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we then compute the expectation value" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\mathbb{E}[\\boldsymbol{X}\\boldsymbol{X}^T] = \\frac{1}{n}\\boldsymbol{X}\\boldsymbol{X}^T=\\begin{bmatrix}\n", + "x_{00}^2+x_{01}^2 & x_{00}x_{10}+x_{01}x_{11}\\\\\n", + "x_{10}x_{00}+x_{11}x_{01} & x_{10}^2+x_{11}^2\\\\\n", + "\\end{bmatrix},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which is just" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{C}[\\boldsymbol{x}_0,\\boldsymbol{x}_1] = \\boldsymbol{C}[\\boldsymbol{x}]=\\begin{bmatrix} \\mathrm{var}[\\boldsymbol{x}_0] & \\mathrm{cov}[\\boldsymbol{x}_0,\\boldsymbol{x}_1] \\\\\n", + " \\mathrm{cov}[\\boldsymbol{x}_1,\\boldsymbol{x}_0] & \\mathrm{var}[\\boldsymbol{x}_1] \\\\\n", + " \\end{bmatrix},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where we wrote $$\\boldsymbol{C}[\\boldsymbol{x}_0,\\boldsymbol{x}_1] = \\boldsymbol{C}[\\boldsymbol{x}]$$ to indicate that this the covariance of the vectors $\\boldsymbol{x}$ of the design/feature matrix $\\boldsymbol{X}$.\n", + "\n", + "It is easy to generalize this to a matrix $\\boldsymbol{X}\\in {\\mathbb{R}}^{n\\times p}$.\n", + "\n", + "\n", + "## Linking with SVD\n", + "\n", + "We have that the covariance matrix (the correlation matrix involves a simple rescaling) is given as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{C}[\\boldsymbol{x}] = \\frac{1}{n}\\boldsymbol{X}\\boldsymbol{X}^T= \\mathbb{E}[\\boldsymbol{X}\\boldsymbol{X}^T].\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us now assume that we can perform a series of orthogonal transformations where we employ some orthogonal matrices $\\boldsymbol{S}$.\n", + "These matrices are defined as $\\boldsymbol{S}\\in {\\mathbb{R}}^{p\\times p}$ and obey the orthogonality requirements $\\boldsymbol{S}\\boldsymbol{S}^T=\\boldsymbol{S}^T\\boldsymbol{S}=\\boldsymbol{I}$. The matrix can be written out in terms of the column vectors $\\boldsymbol{s}_i$ as $\\boldsymbol{S}=[\\boldsymbol{s}_0,\\boldsymbol{s}_1,\\dots,\\boldsymbol{s}_{p-1}]$ and $\\boldsymbol{s}_i \\in {\\mathbb{R}}^{p}$.\n", + "\n", + "Assume also that there is a transformation $\\boldsymbol{S}\\boldsymbol{C}[\\boldsymbol{x}]\\boldsymbol{S}^T=\\boldsymbol{C}[\\boldsymbol{y}]$ such that the new matrix $\\boldsymbol{C}[\\boldsymbol{y}]$ is diagonal with elements $[\\lambda_0,\\lambda_1,\\lambda_2,\\dots,\\lambda_{p-1}]$. \n", + "\n", + "That is we have" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{C}[\\boldsymbol{y}] = \\mathbb{E}[\\boldsymbol{S}\\boldsymbol{X}\\boldsymbol{X}^T\\boldsymbol{S}^T]=\\boldsymbol{S}\\boldsymbol{C}[\\boldsymbol{x}]\\boldsymbol{S}^T,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "since the matrix $\\boldsymbol{S}$ is not a data dependent matrix. Multiplying with $\\boldsymbol{S}^T$ from the left we have" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{S}^T\\boldsymbol{C}[\\boldsymbol{y}] = \\boldsymbol{C}[\\boldsymbol{x}]\\boldsymbol{S}^T,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and since $\\boldsymbol{C}[\\boldsymbol{y}]$ is diagonal we have for a given eigenvalue $i$ of the covariance matrix that" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{S}^T_i\\lambda_i = \\boldsymbol{C}[\\boldsymbol{x}]\\boldsymbol{S}^T_i.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the derivation of the PCA theorem we will assume that the eigenvalues are ordered in descending order, that is\n", + "$\\lambda_0 > \\lambda_1 > \\dots > \\lambda_{p-1}$. \n", + "\n", + "\n", + "The eigenvalues tell us then how much we need to stretch the\n", + "corresponding eigenvectors. Dimensions with large eigenvalues have\n", + "thus large variations (large variance) and define therefore useful\n", + "dimensions. The data points are more spread out in the direction of\n", + "these eigenvectors. Smaller eigenvalues mean on the other hand that\n", + "the corresponding eigenvectors are shrunk accordingly and the data\n", + "points are tightly bunched together and there is not much variation in\n", + "these specific directions. Hopefully then we could leave it out\n", + "dimensions where the eigenvalues are very small. If $p$ is very large,\n", + "we could then aim at reducing $p$ to $l << p$ and handle only $l$\n", + "features/predictors.\n", + "\n", + "\n", + "\n", "\n", "## Where are we going?\n", "\n", @@ -3254,7 +4039,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 28, "metadata": { "collapsed": false }, @@ -3279,7 +4064,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 29, "metadata": { "collapsed": false }, @@ -3375,7 +4160,7 @@ "metadata": {}, "source": [ "1\n", - "1\n", + "3\n", "2\n", " \n", "<\n", @@ -4444,7 +5229,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 30, "metadata": { "collapsed": false }, @@ -4587,7 +5372,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 31, "metadata": { "collapsed": false }, @@ -4735,7 +5520,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 32, "metadata": { "collapsed": false }, @@ -4966,7 +5751,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 33, "metadata": { "collapsed": false }, @@ -5037,7 +5822,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 34, "metadata": { "collapsed": false }, @@ -5138,7 +5923,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 35, "metadata": { "collapsed": false }, @@ -5226,7 +6011,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 36, "metadata": { "collapsed": false }, @@ -5322,7 +6107,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 37, "metadata": { "collapsed": false }, @@ -5406,7 +6191,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 38, "metadata": { "collapsed": false }, @@ -5489,7 +6274,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 39, "metadata": { "collapsed": false }, @@ -5604,7 +6389,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 40, "metadata": { "collapsed": false }, @@ -5670,7 +6455,7 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 41, "metadata": { "collapsed": false }, @@ -5688,7 +6473,7 @@ }, { "cell_type": "code", - "execution_count": 37, + "execution_count": 42, "metadata": { "collapsed": false }, @@ -5775,7 +6560,7 @@ }, { "cell_type": "code", - "execution_count": 38, + "execution_count": 43, "metadata": { "collapsed": false }, @@ -5788,7 +6573,7 @@ }, { "cell_type": "code", - "execution_count": 39, + "execution_count": 44, "metadata": { "collapsed": false }, @@ -5806,7 +6591,7 @@ }, { "cell_type": "code", - "execution_count": 40, + "execution_count": 45, "metadata": { "collapsed": false }, @@ -5824,7 +6609,7 @@ }, { "cell_type": "code", - "execution_count": 41, + "execution_count": 46, "metadata": { "collapsed": false }, @@ -5892,7 +6677,7 @@ }, { "cell_type": "code", - "execution_count": 42, + "execution_count": 47, "metadata": { "collapsed": false }, @@ -6000,7 +6785,7 @@ }, { "cell_type": "code", - "execution_count": 43, + "execution_count": 48, "metadata": { "collapsed": false }, @@ -6032,7 +6817,7 @@ }, { "cell_type": "code", - "execution_count": 44, + "execution_count": 49, "metadata": { "collapsed": false }, @@ -6050,7 +6835,7 @@ }, { "cell_type": "code", - "execution_count": 45, + "execution_count": 50, "metadata": { "collapsed": false }, @@ -6068,7 +6853,7 @@ }, { "cell_type": "code", - "execution_count": 46, + "execution_count": 51, "metadata": { "collapsed": false }, @@ -6104,7 +6889,7 @@ "metadata": {}, "source": [ "1\n", - "7\n", + "9\n", "6\n", " \n", "<\n", @@ -6126,7 +6911,7 @@ }, { "cell_type": "code", - "execution_count": 47, + "execution_count": 52, "metadata": { "collapsed": false }, @@ -6179,7 +6964,7 @@ }, { "cell_type": "code", - "execution_count": 48, + "execution_count": 53, "metadata": { "collapsed": false }, @@ -6215,7 +7000,7 @@ }, { "cell_type": "code", - "execution_count": 49, + "execution_count": 54, "metadata": { "collapsed": false }, @@ -6280,7 +7065,7 @@ }, { "cell_type": "code", - "execution_count": 50, + "execution_count": 55, "metadata": { "collapsed": false },