lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

updating

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2021-09-13 11:52:44 +02:00
parent a5a8a681f9
commit c30998b68a
18 changed files with 5355 additions and 8651 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/LectureNotes/chapter1.ipynb
+271 -85
View File
@@ -6,8 +6,6 @@
"source": [
"# Linear Regression\n",
"\n",
"[Video of Lecture](https://www.uio.no/studier/emner/matnat/fys/FYS-STK3155/h20/forelesningsvideoer/LectureAug21.mp4?vrtx=view-as-webpage)\n",
"\n",
"\n",
"## Introduction\n",
"\n",
@@ -15,7 +13,7 @@
"\n",
"\n",
"\n",
"Our emphasis throughout this series of lectures (small change) \n",
"Our emphasis throughout this series of lectures \n",
"is on understanding the mathematical aspects of\n",
"different algorithms used in the fields of data analysis and machine learning. \n",
"\n",
@@ -235,7 +233,9 @@
"\n",
"## Simple linear regression model using **scikit-learn**\n",
"\n",
"We start with perhaps our simplest possible example, using **Scikit-Learn** to perform linear regression analysis on a data set produced by us. \n",
"We start with perhaps our simplest possible example, using\n",
"**Scikit-Learn** to perform linear regression analysis on a data set\n",
"produced by us.\n",
"\n",
"What follows is a simple Python code where we have defined a function\n",
"$y$ in terms of the variable $x$. Both are defined as vectors with $100$ entries. \n",
@@ -437,7 +437,8 @@
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.linear_model import LinearRegression\n",
"\n",
"# Number of data points\n",
"n = 100\n",
"x = np.random.rand(100,1)\n",
"y = 5*x+0.01*np.random.randn(100,1)\n",
"linreg = LinearRegression()\n",
@@ -459,7 +460,8 @@
"Depending on the parameter in front of the normal distribution, we may\n",
"have a small or larger relative error. Try to play around with\n",
"different training data sets and study (graphically) the value of the\n",
"relative error.\n",
"relative error. Note also that **Scikit-Learn** requires a matrix as input for the input values $x$ and $y$. In the above code we have\n",
"solved this by declaring $x$ and $y$ as arrays of dimension $n\\times 1$.\n",
"\n",
"As mentioned above, **Scikit-Learn** has an impressive functionality.\n",
"We can for example extract the values of $\\alpha$ and $\\beta$ and\n",
@@ -629,8 +631,7 @@
"metadata": {},
"source": [
"$$\n",
"H_{\\delta}(\\boldsymbol{a})=\\left\\{\\begin{array}{cc}\\frac{1}{2} \\boldsymbol{a}^{2}& \\text{for }|\\boldsymbol{a}|\\leq \\delta\\\\ \\delta (|\\b\\\n",
"m{a}|-\\frac{1}{2}\\delta ),&\\text{otherwise}.\\end{array}\\right.\n",
"H_{\\delta}(\\boldsymbol{a})=\\left\\{\\begin{array}{cc}\\frac{1}{2} \\boldsymbol{a}^{2}& \\text{for }|\\boldsymbol{a}|\\leq \\delta\\\\ \\delta (|\\boldsymbol{a}|-\\frac{1}{2}\\delta ),&\\text{otherwise}.\\end{array}\\right.\n",
"$$"
]
},
@@ -641,6 +642,8 @@
"Here $\\boldsymbol{a}=\\boldsymbol{y} - \\boldsymbol{\\tilde{y}}$.\n",
"\n",
"\n",
"\n",
"\n",
"We will discuss in more\n",
"detail these and other functions in the various lectures. We conclude this part with another example. Instead of \n",
"a linear $x$-dependence we study now a cubic polynomial and use the polynomial regression analysis tools of scikit-learn."
@@ -1061,6 +1064,8 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Note well that we have made life simple here. We perform a fit in terms of the number of nucleons only. A more sophisticated fit can be done by including an explicit dependence on the number of protons and neutrons in the asymmetry and Coulomb terms.\n",
"\n",
"With **scikitlearn** we are now ready to use linear regression and fit our data."
]
},
@@ -1100,7 +1105,6 @@
"print('Variance score: %.2f' % r2_score(Energies, fity))\n",
"# Mean absolute error \n",
"print('Mean absolute error: %.2f' % mean_absolute_error(Energies, fity))\n",
"print(clf.coef_, clf.intercept_)\n",
"\n",
"Masses['Eapprox'] = fity\n",
"# Generate a plot comparing the experimental with the fitted values values.\n",
@@ -3046,6 +3050,119 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Splitting our Data in Training and Test data\n",
"\n",
"\n",
"It is normal in essentially all Machine Learning studies to split the\n",
"data in a training set and a test set (sometimes also an additional\n",
"validation set). **Scikit-Learn** has an own function for this. There\n",
"is no explicit recipe for how much data should be included as training\n",
"data and say test data. An accepted rule of thumb is to use\n",
"approximately $2/3$ to $4/5$ of the data as training data. We will\n",
"postpone a discussion of this splitting to the end of these notes and\n",
"our discussion of the so-called **bias-variance** tradeoff. Here we\n",
"limit ourselves to repeat the above equation of state fitting example\n",
"but now splitting the data into a training set and a test set.\n",
"\n",
"Let us study some examples. The first code here takes a simple\n",
"one-dimensional second-order polynomial and we fit it to a\n",
"second-order polynomial. Depending on the strength of the added noise,\n",
"the various measures like the $R2$ score or the mean-squared error,\n",
"the fit becomes better or worse."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import os\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import train_test_split\n",
"\n",
"\n",
"def R2(y_data, y_model):\n",
" return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n",
"def MSE(y_data,y_model):\n",
" n = np.size(y_model)\n",
" return np.sum((y_data-y_model)**2)/n\n",
"\n",
"x = np.random.rand(100)\n",
"y = 2.0+5*x*x+0.1*np.random.randn(100)\n",
"\n",
"\n",
"# The design matrix now as function of a given polynomial\n",
"X = np.zeros((len(x),3))\n",
"X[:,0] = 1.0\n",
"X[:,1] = x\n",
"X[:,2] = x**2\n",
"# We split the data in test and training data\n",
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"# matrix inversion to find beta\n",
"beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train\n",
"print(beta)\n",
"# and then make the prediction\n",
"ytilde = X_train @ beta\n",
"print(\"Training R2\")\n",
"print(R2(y_train,ytilde))\n",
"print(\"Training MSE\")\n",
"print(MSE(y_train,ytilde))\n",
"ypredict = X_test @ beta\n",
"print(\"Test R2\")\n",
"print(R2(y_test,ypredict))\n",
"print(\"Test MSE\")\n",
"print(MSE(y_test,ypredict))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Alternatively, you could write your own test-train splitting function as shown here."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"# equivalently in numpy\n",
"def train_test_split_numpy(inputs, labels, train_size, test_size):\n",
" n_inputs = len(inputs)\n",
" inputs_shuffled = inputs.copy()\n",
" labels_shuffled = labels.copy()\n",
"\n",
" np.random.shuffle(inputs_shuffled)\n",
" np.random.shuffle(labels_shuffled)\n",
"\n",
" train_end = int(n_inputs*train_size)\n",
" X_train, X_test = inputs_shuffled[:train_end], inputs_shuffled[train_end:]\n",
" Y_train, Y_test = labels_shuffled[:train_end], labels_shuffled[train_end:]\n",
"\n",
" return X_train, X_test, Y_train, Y_test"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But since **scikit-learn** has its own function for doing this and since\n",
"it interfaces easily with **tensorflow** and other libraries, we\n",
"normally recommend using the latter functionality.\n",
"\n",
"\n",
"\n",
"\n",
"## Reducing the number of degrees of freedom, overarching view\n",
"\n",
"Many Machine Learning problems involve thousands or even millions of\n",
@@ -3069,6 +3186,7 @@
"visualization.\n",
"\n",
"\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",
@@ -3076,6 +3194,15 @@
"to extreme values. Scaling the data renders our inputs much more\n",
"suitable for the algorithms we want to employ.\n",
"\n",
"For data sets gathered for real world applications, it is rather normal that\n",
"different features have very different units and\n",
"numerical scales. For example, a data set detailing health habits may include\n",
"features such as **age** in the range $0-80$, and **caloric intake** of order $2000$.\n",
"Many machine learning methods sensitive to the scales of the features and may perform poorly if they\n",
"are very different scales. Therefore, it is typical to scale\n",
"the features in a way to avoid such outlier values.\n",
"\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",
@@ -3105,7 +3232,38 @@
"techniques.\n",
"\n",
"\n",
"### Simple preprocessing examples, Franke function and regression"
"Many features are often scaled using standardization to improve\n",
"performance. In **Scikit-Learn** this is given by the **StandardScaler**\n",
"function as discussed above. It is easy however to write your own.\n",
"Mathematically, this involves subtracting the mean and divide by the\n",
"standard deviation over the data set, for each feature:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"x_j^{(i)} \\rightarrow \\frac{x_j^{(i)} - \\overline{x}_j}{\\sigma(x_j)},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $\\overline{x}_j$ and $\\sigma(x_j)$ are the mean and standard\n",
"deviation, respectively, of the feature $x_j$. This ensures that each\n",
"feature has zero mean and unit standard deviation. For data sets\n",
"where we do not have the standard deviation or don't wish to calculate\n",
"it, it is then common to simply set it to one.\n",
"\n",
"\n",
"\n",
"Let us consider the following vanilla example where we use both\n",
"**Scikit-Learn** and write our own function as well. We produce a\n",
"simple test design matrix with random numbers. Each column could then\n",
"represent a specific feature whose mean value is subracted."
]
},
{
@@ -3117,98 +3275,126 @@
},
"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",
"import numpy as np\n",
"import pandas as pd\n",
"from IPython.display import display\n",
"np.random.seed(100)\n",
"# setting up a 10 x 5 matrix\n",
"rows = 10\n",
"cols = 5\n",
"X = np.random.randn(rows,cols)\n",
"XPandas = pd.DataFrame(X)\n",
"display(XPandas)\n",
"print(XPandas.mean())\n",
"print(XPandas.std())\n",
"XPandas = (XPandas -XPandas.mean())\n",
"display(XPandas)\n",
"# This option does not include the standard deviation\n",
"scaler = StandardScaler(with_std=False)\n",
"scaler.fit(X)\n",
"Xscaled = scaler.transform(X)\n",
"display(XPandas-Xscaled)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Small exercise: perform the standard scaling by including the standard deviation and compare with what Scikit-Learn gives.\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",
"Another commonly used scaling method is min-max scaling. This is very\n",
"useful for when we want the features to lie in a certain interval. To\n",
"scale the feature $x_j$ to the interval $[a, b]$, we can apply the\n",
"transformation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"x_j^{(i)} \\rightarrow (b-a)\\frac{x_j^{(i)} - \\min(x_j)}{\\max(x_j) - \\min(x_j)} - a\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $\\min(x_j)$ and $\\max(x_j)$ return the minimum and maximum value of $x_j$ over the data set, respectively.\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",
"## Testing the Means Squared Error as function of Complexity\n",
"\n",
"\n",
"clf = skl.LinearRegression().fit(X_train, y_train)\n",
"Before we proceed with a more detailed analysis of the so-called\n",
"Bias-Variance tradeoff, we present here an example of the relation\n",
"between model complexity and the mean squared error for the triaining\n",
"data and the test data.\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",
"The results here tell us clearly that for the data not included in the\n",
"training, there is an optimal model as function of the complexity of\n",
"ourmodel (here in terms of the polynomial degree of the model).\n",
"\n",
"The results here will vary as function of model complexity and the amount od data used for training. \n",
"\n",
"\n",
"Our data is defined by $x\\in [-3,3]$ with a total of for example $100$ data points."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from sklearn.linear_model import LinearRegression, Ridge, Lasso\n",
"from sklearn.preprocessing import PolynomialFeatures\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.pipeline import make_pipeline\n",
"\n",
"\n",
"np.random.seed(2018)\n",
"n = 100\n",
"maxdegree = 14\n",
"# Make data set.\n",
"x = np.linspace(-3, 3, n).reshape(-1, 1)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)\n",
"TestError = np.zeros(maxdegree)\n",
"TrainError = np.zeros(maxdegree)\n",
"polydegree = np.zeros(maxdegree)\n",
"x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)\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",
"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",
"for degree in range(maxdegree):\n",
" model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False))\n",
" clf = model.fit(x_train_scaled,y_train)\n",
" y_fit = clf.predict(x_train_scaled)\n",
" y_pred = clf.predict(x_test_scaled) \n",
" polydegree[degree] = degree\n",
" TestError[degree] = np.mean( np.mean((y_test - y_pred)**2) )\n",
" TrainError[degree] = np.mean( np.mean((y_train - y_fit)**2) )\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)))"
"plt.plot(polydegree, TestError, label='Test Error')\n",
"plt.plot(polydegree, TrainError, label='Train Error')\n",
"plt.legend()\n",
"plt.show()"
]
},
{