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

adding codes

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2023-09-11 09:53:59 +02:00
parent 006b2407fe
commit 6f1eca1dc1
8 changed files with 1555 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/src/week37/programs/.ipynb_checkpoints/codeexamplesscaling-checkpoint.ipynb
+278
View File
File diff suppressed because one or more lines are too long
doc/src/week37/programs/LinearReg.py
+2 -2
View File
@@ -4,7 +4,7 @@ from sklearn.model_selection import train_test_split
class linregOwn:
"""
A class of linear regressions. Perform ordinarly least squares (OLS) and Ridge regression manually. Lasso
A class of linear regressions. Perform ordinary least squares (OLS) and Ridge regression manually. Lasso
is performed using scikit-learn functionality.
"""
def __init__(self, method = 'ols'):
@@ -267,4 +267,4 @@ class linregSKL:
if self.yHat is None :
self._sklPredict()
self._R2 = metrics.r2_score(self.y_test, self.yHat)
return self._R2
return self._R2
doc/src/week37/programs/codeexamplesscaling.do.txt
+370
View File
@@ -0,0 +1,370 @@
TITLE: Scaling examples with own code and the library Scikit-Learn
AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University
DATE: today
===== This note contains code examples with a simple scaling =====
The programs here use both ordinrary least squares and Ridge regression with one value only for
the hyperparameter $\lambda$. The first example has no scaling and includes the intercept as well and we are trying to fit a second-order
polynomial.
!bc pycod
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
def MSE(y_data,y_model):
n = np.size(y_model)
return np.sum((y_data-y_model)**2)/n
def OLS_fit_beta(X, y):
return np.linalg.pinv(X.T @ X) @ X.T @ y
def Ridge_fit_beta(X, y,L,d):
I = np.eye(d,d)
return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y
# Same random numbers for each test.
np.random.seed(2018)
n = 100
d = 3
Lambda = 0.01
true_beta = [2, 0.5, 3.7]
# Make data set.
x = np.linspace(-3, 3, n)
y_real = 2 + 0.5*x + 3.7*x**2
y = np.sum(
np.asarray([x ** p * b for p, b in enumerate(true_beta)]),
axis=0) + 0.1 * np.random.normal(size=len(x))
#Design matrix X includes the intercept and scaling is made
X = np.zeros((len(x), d))
for p in range(d):
X[:, p] = x ** (p)
#Split data, no scaling is used and we include the intercept
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
#Calculate beta, own code
beta_OLS = OLS_fit_beta(X_train, y_train)
beta_Ridge = Ridge_fit_beta(X_train, y_train,Lambda,d)
print(beta_OLS)
print(beta_Ridge)
#predict value
ytilde_test_OLS = X_test @ beta_OLS
ytilde_test_Ridge = X_test @ beta_Ridge
#Calculate MSE
print(" ")
print("test MSE of OLS:")
print(MSE(y_test,ytilde_test_OLS))
print(" ")
print("test MSE of Ridge")
print(MSE(y_test,ytilde_test_Ridge))
plt.scatter(x,y,label='Data')
plt.plot(x, X @ beta_OLS,'*', label="OLS_Fit")
plt.plot(x, X @ beta_Ridge, label="Ridge_Fit")
plt.grid()
plt.legend()
plt.show()
!ec
In this example we do not include the intercept and we scale the data by subtracting the mean values. This follows the discussion in the "lecture material":"https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/chapter3.html#more-on-rescaling-data".
see also the weekly slides "for week 36":"https://compphysics.github.io/MachineLearning/doc/pub/week36/html/._week36-bs029.html".
Before we discuss the code, we repeat some of the basic math from the slides of week 36.
Let us try to understand what this may imply mathematically when we
subtract the mean values, also known as *zero centering*. For
simplicity, we will focus on ordinary regression, as done in the above example.
The cost/loss function for regression is
!bt
\[
C(\beta_0, \beta_1, ... , \beta_{p-1}) = \frac{1}{n}\sum_{i=0}^{n} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij}\beta_j\right)^2,.
\]
!et
Recall also that we use the squared value. This expression can lead to an
increased penalty for higher differences between predicted and
output/target values.
What we have done is to single out the $\beta_0$ term in the
definition of the mean squared error (MSE). The design matrix $X$
does in this case not contain any intercept column. When we take the
derivative with respect to $\beta_0$, we want the derivative to obey
!bt
\[
\frac{\partial C}{\partial \beta_j} = 0,
\]
!et
for all $j$. For $\beta_0$ we have
!bt
\[
\frac{\partial C}{\partial \beta_0} = -\frac{2}{n}\sum_{i=0}^{n-1} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij} \beta_j\right).
\]
!et
Multiplying away the constant $2/n$, we obtain
!bt
\[
\sum_{i=0}^{n-1} \beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} \sum_{j=1}^{p-1} X_{ij} \beta_j.
\]
!et
Let us specialize first to the case where we have only two parameters $\beta_0$ and $\beta_1$.
Our result for $\beta_0$ simplifies then to
!bt
\[
n\beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \beta_1.
\]
!et
We obtain then
!bt
\[
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \beta_1\frac{1}{n}\sum_{i=0}^{n-1} X_{i1}.
\]
!et
If we define
!bt
\[
\mu_{\bm{x}_1}=\frac{1}{n}\sum_{i=0}^{n-1} X_{i1},
\]
!et
and the mean value of the outputs as
!bt
\[
\mu_y=\frac{1}{n}\sum_{i=0}^{n-1}y_i,
\]
!et
we have
!bt
\[
\beta_0 = \mu_y - \beta_1\mu_{\bm{x}_1}.
\]
!et
In the general case with more parameters than $\beta_0$ and $\beta_1$, we have
!bt
\[
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \frac{1}{n}\sum_{i=0}^{n-1}\sum_{j=1}^{p-1} X_{ij}\beta_j.
\]
!et
We can rewrite the latter equation as
!bt
\[
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \sum_{j=1}^{p-1} \mu_{\bm{x}_j}\beta_j,
\]
!et
where we have defined
!bt
\[
\mu_{\bm{x}_j}=\frac{1}{n}\sum_{i=0}^{n-1} X_{ij},
\]
!et
the mean value for all elements of the column vector $\bm{x}_j$.
Replacing $y_i$ with $y_i - y_i - \overline{\bm{y}}$ and centering also our design matrix results in a cost function (in vector-matrix disguise)
!bt
\[
C(\boldsymbol{\beta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta}).
\]
!et
If we minimize with respect to $\bm{\beta}$ we have then
!bt
\[
\hat{\bm{\beta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}},
\]
!et
where $\boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\bm{y}}$
and $\tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj}$.
For Ridge regression we need to add $\lambda \boldsymbol{\beta}^T\boldsymbol{\beta}$ to the cost function and get then
!bt
\[
\hat{\bm{\beta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}.
\]
!et
Now we try to implement this.
!bc pycod
np.random.seed(2018)
n = 100
d = 3
Lambda = 0.01
true_beta = [2, 0.5, 3.7]
# Make data set.
x = np.linspace(-3, 3, n)
y_real = 2 + 0.5*x + 3.7*x**2
y = np.sum(
np.asarray([x ** p * b for p, b in enumerate(true_beta)]),
axis=0) + 0.1 * np.random.normal(size=len(x))
#Design matrix X does not include the intercept.
X = np.zeros((len(x), d))
for p in range(d-1):
X[:, p] = x ** (p+1)
#Split data in train and test
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
# Scale data by subtracting mean value,own implementation
#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable
X_train_mean = np.mean(X_train,axis=0)
#Center by removing mean from each feature
X_train_scaled = X_train - X_train_mean
X_test_scaled = X_test - X_train_mean
#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note)
y_scaler = np.mean(y_train)
y_train_scaled = y_train - y_scaler
#Calculate beta
beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)
beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)
print(beta_OLS)
print(beta_Ridge)
interceptOLS = y_scaler - X_train_mean @ beta_OLS
interceptRidge = y_scaler - X_train_mean @ beta_Ridge
print(interceptOLS)
print(interceptRidge)
#predict value
ytilde_test_OLS = X_test_scaled @ beta_OLS+y_scaler
ytilde_test_Ridge = X_test_scaled @ beta_Ridge+y_scaler
#Calculate MSE
print(" ")
print("test MSE of OLS:")
print(MSE(y_test,ytilde_test_OLS))
print(" ")
print("test MSE of Ridge")
print(MSE(y_test,ytilde_test_Ridge))
plt.scatter(x,y,label='Data')
#plt.plot(x,y_real,label='no noise')
plt.plot(x, X @ beta_OLS+interceptOLS,'*', label="OLS_Fit")
plt.plot(x, X @ beta_Ridge+interceptRidge, label="Ridge_Fit")
plt.grid()
plt.legend()
plt.show()
!ec
We see that we get the same values for the parameters! As it should be. The MSE may however change (not the case here).
Finally, instead of using our own function we repeat the same example using the _standardscaler_ functionality of the library _Scikit-Learn_.
#!bc pycod
np.random.seed(2018)
n = 100
d = 3
Lambda = 0.01
true_beta = [2, 0.5, 3.7]
# Make data set.
x = np.linspace(-3, 3, n)
y_real = 2 + 0.5*x + 3.7*x**2
y = np.sum(
np.asarray([x ** p * b for p, b in enumerate(true_beta)]),
axis=0) + 0.1 * np.random.normal(size=len(x))
#Design matrix X does not include the intercept.
X = np.zeros((len(x), d))
for p in range(d-1):
X[:, p] = x ** (p+1)
#Split data in train and test
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
# Scale data by subtracting mean value,own implementation
#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable
X_train_mean = np.mean(X_train,axis=0)
#Center by removing mean from each feature
X_train_scaled = X_train - X_train_mean
X_test_scaled = X_test - X_train_mean
#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note)
y_scaler = np.mean(y_train)
y_train_scaled = y_train - y_scaler
#Calculate beta
beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)
beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)
print(beta_OLS)
print(beta_Ridge)
interceptOLS = y_scaler - X_train_mean @ beta_OLS
interceptRidge = y_scaler - X_train_mean @ beta_Ridge
print(interceptOLS)
print(interceptRidge)
#predict value
ytilde_test_OLS = X_test_scaled @ beta_OLS+y_scaler
ytilde_test_Ridge = X_test_scaled @ beta_Ridge+y_scaler
#Calculate MSE
print(" ")
print("test MSE of OLS:")
print(MSE(y_test,ytilde_test_OLS))
print(" ")
print("test MSE of Ridge")
print(MSE(y_test,ytilde_test_Ridge))
plt.scatter(x,y,label='Data')
#plt.plot(x,y_real,label='no noise')
plt.plot(x, X @ beta_OLS+interceptOLS,'*', label="OLS_Fit")
plt.plot(x, X @ beta_Ridge+interceptRidge, label="Ridge_Fit")
plt.grid()
plt.legend()
plt.show()
#!ec
doc/src/week37/programs/codeexamplesscaling.ipynb
+649
View File
@@ -0,0 +1,649 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "718b0cee",
"metadata": {
"editable": true
},
"source": [
"<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)\n",
"doconce format html codeexamplesscaling.do.txt -->\n",
"<!-- dom:TITLE: Scaling examples with own code and the library Scikit-Learn -->"
]
},
{
"cell_type": "markdown",
"id": "cdce7555",
"metadata": {
"editable": true
},
"source": [
"# Scaling examples with own code and the library Scikit-Learn\n",
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University\n",
"\n",
"Date: **Sep 11, 2023**\n",
"\n",
"Copyright 1999-2023, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license"
]
},
{
"cell_type": "markdown",
"id": "67634fc9",
"metadata": {
"editable": true
},
"source": [
"## This note contains code examples with a simple scaling\n",
"\n",
"The programs here use both ordinrary least squares and Ridge regression with one value only for\n",
"the hyperparameter $\\lambda$. The first example has no scaling and includes the intercept as well and we are trying to fit a second-order\n",
"polynomial."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "e75a3307",
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from sklearn.linear_model import LinearRegression\n",
"from sklearn.preprocessing import PolynomialFeatures\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.preprocessing import StandardScaler\n",
"\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",
"def OLS_fit_beta(X, y):\n",
" return np.linalg.pinv(X.T @ X) @ X.T @ y\n",
"\n",
"def Ridge_fit_beta(X, y,L,d):\n",
" I = np.eye(d,d)\n",
" return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y\n",
"\n",
"# Same random numbers for each test.\n",
"np.random.seed(2018)\n",
"n = 100\n",
"d = 3\n",
"Lambda = 0.01\n",
"true_beta = [2, 0.5, 3.7]\n",
"\n",
"# Make data set.\n",
"x = np.linspace(-3, 3, n)\n",
"y_real = 2 + 0.5*x + 3.7*x**2\n",
"\n",
"y = np.sum(\n",
" np.asarray([x ** p * b for p, b in enumerate(true_beta)]), \n",
" axis=0) + 0.1 * np.random.normal(size=len(x))\n",
"\n",
"\n",
"#Design matrix X includes the intercept and scaling is made\n",
"X = np.zeros((len(x), d))\n",
"for p in range(d): \n",
" X[:, p] = x ** (p) \n",
"\n",
"\n",
"#Split data, no scaling is used and we include the intercept\n",
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"\n",
"\n",
"#Calculate beta, own code\n",
"beta_OLS = OLS_fit_beta(X_train, y_train)\n",
"beta_Ridge = Ridge_fit_beta(X_train, y_train,Lambda,d)\n",
"print(beta_OLS)\n",
"print(beta_Ridge)\n",
"#predict value\n",
"ytilde_test_OLS = X_test @ beta_OLS\n",
"ytilde_test_Ridge = X_test @ beta_Ridge\n",
"\n",
"#Calculate MSE\n",
"print(\" \")\n",
"print(\"test MSE of OLS:\")\n",
"print(MSE(y_test,ytilde_test_OLS))\n",
"print(\" \")\n",
"print(\"test MSE of Ridge\")\n",
"print(MSE(y_test,ytilde_test_Ridge))\n",
"\n",
"plt.scatter(x,y,label='Data')\n",
"plt.plot(x, X @ beta_OLS,'*', label=\"OLS_Fit\")\n",
"plt.plot(x, X @ beta_Ridge, label=\"Ridge_Fit\")\n",
"plt.grid()\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "25d58cba",
"metadata": {
"editable": true
},
"source": [
"In this example we do not include the intercept and we scale the data by subtracting the mean values. This follows the discussion in the [lecture material](https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/chapter3.html#more-on-rescaling-data).\n",
"see also the weekly slides [for week 36](https://compphysics.github.io/MachineLearning/doc/pub/week36/html/._week36-bs029.html).\n",
"\n",
"Before we discuss the code, we repeat some of the basic math from the slides of week 36.\n",
"\n",
"Let us try to understand what this may imply mathematically when we\n",
"subtract the mean values, also known as *zero centering*. For\n",
"simplicity, we will focus on ordinary regression, as done in the above example.\n",
"\n",
"The cost/loss function for regression is"
]
},
{
"cell_type": "markdown",
"id": "5b820b6d",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"C(\\beta_0, \\beta_1, ... , \\beta_{p-1}) = \\frac{1}{n}\\sum_{i=0}^{n} \\left(y_i - \\beta_0 - \\sum_{j=1}^{p-1} X_{ij}\\beta_j\\right)^2,.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "9d3946fb",
"metadata": {
"editable": true
},
"source": [
"Recall also that we use the squared value. This expression can lead to an\n",
"increased penalty for higher differences between predicted and\n",
"output/target values.\n",
"\n",
"What we have done is to single out the $\\beta_0$ term in the\n",
"definition of the mean squared error (MSE). The design matrix $X$\n",
"does in this case not contain any intercept column. When we take the\n",
"derivative with respect to $\\beta_0$, we want the derivative to obey"
]
},
{
"cell_type": "markdown",
"id": "11228e0d",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\frac{\\partial C}{\\partial \\beta_j} = 0,\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "6dd5cd47",
"metadata": {
"editable": true
},
"source": [
"for all $j$. For $\\beta_0$ we have"
]
},
{
"cell_type": "markdown",
"id": "4bf8184e",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\frac{\\partial C}{\\partial \\beta_0} = -\\frac{2}{n}\\sum_{i=0}^{n-1} \\left(y_i - \\beta_0 - \\sum_{j=1}^{p-1} X_{ij} \\beta_j\\right).\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "c8c7dc9d",
"metadata": {
"editable": true
},
"source": [
"Multiplying away the constant $2/n$, we obtain"
]
},
{
"cell_type": "markdown",
"id": "3c6a746e",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\sum_{i=0}^{n-1} \\beta_0 = \\sum_{i=0}^{n-1}y_i - \\sum_{i=0}^{n-1} \\sum_{j=1}^{p-1} X_{ij} \\beta_j.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "00b74b33",
"metadata": {
"editable": true
},
"source": [
"Let us specialize first to the case where we have only two parameters $\\beta_0$ and $\\beta_1$.\n",
"Our result for $\\beta_0$ simplifies then to"
]
},
{
"cell_type": "markdown",
"id": "20d80ba8",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"n\\beta_0 = \\sum_{i=0}^{n-1}y_i - \\sum_{i=0}^{n-1} X_{i1} \\beta_1.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "6cf3aa5d",
"metadata": {
"editable": true
},
"source": [
"We obtain then"
]
},
{
"cell_type": "markdown",
"id": "3492889a",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\beta_0 = \\frac{1}{n}\\sum_{i=0}^{n-1}y_i - \\beta_1\\frac{1}{n}\\sum_{i=0}^{n-1} X_{i1}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "6f6526ea",
"metadata": {
"editable": true
},
"source": [
"If we define"
]
},
{
"cell_type": "markdown",
"id": "a34541ec",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\mu_{\\boldsymbol{x}_1}=\\frac{1}{n}\\sum_{i=0}^{n-1} X_{i1},\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "8d08709c",
"metadata": {
"editable": true
},
"source": [
"and the mean value of the outputs as"
]
},
{
"cell_type": "markdown",
"id": "0efce920",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\mu_y=\\frac{1}{n}\\sum_{i=0}^{n-1}y_i,\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "5a0488c6",
"metadata": {
"editable": true
},
"source": [
"we have"
]
},
{
"cell_type": "markdown",
"id": "9727879c",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\beta_0 = \\mu_y - \\beta_1\\mu_{\\boldsymbol{x}_1}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "d0297a4d",
"metadata": {
"editable": true
},
"source": [
"In the general case with more parameters than $\\beta_0$ and $\\beta_1$, we have"
]
},
{
"cell_type": "markdown",
"id": "5b4f7606",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\beta_0 = \\frac{1}{n}\\sum_{i=0}^{n-1}y_i - \\frac{1}{n}\\sum_{i=0}^{n-1}\\sum_{j=1}^{p-1} X_{ij}\\beta_j.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "cd6fa027",
"metadata": {
"editable": true
},
"source": [
"We can rewrite the latter equation as"
]
},
{
"cell_type": "markdown",
"id": "73827116",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\beta_0 = \\frac{1}{n}\\sum_{i=0}^{n-1}y_i - \\sum_{j=1}^{p-1} \\mu_{\\boldsymbol{x}_j}\\beta_j,\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "83069a4e",
"metadata": {
"editable": true
},
"source": [
"where we have defined"
]
},
{
"cell_type": "markdown",
"id": "1fd7a5df",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\mu_{\\boldsymbol{x}_j}=\\frac{1}{n}\\sum_{i=0}^{n-1} X_{ij},\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "e6e597a6",
"metadata": {
"editable": true
},
"source": [
"the mean value for all elements of the column vector $\\boldsymbol{x}_j$.\n",
"\n",
"Replacing $y_i$ with $y_i - y_i - \\overline{\\boldsymbol{y}}$ and centering also our design matrix results in a cost function (in vector-matrix disguise)"
]
},
{
"cell_type": "markdown",
"id": "6ef1dd83",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"C(\\boldsymbol{\\beta}) = (\\boldsymbol{\\tilde{y}} - \\tilde{X}\\boldsymbol{\\beta})^T(\\boldsymbol{\\tilde{y}} - \\tilde{X}\\boldsymbol{\\beta}).\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "190376f5",
"metadata": {
"editable": true
},
"source": [
"If we minimize with respect to $\\boldsymbol{\\beta}$ we have then"
]
},
{
"cell_type": "markdown",
"id": "7bfde500",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\hat{\\boldsymbol{\\beta}} = (\\tilde{X}^T\\tilde{X})^{-1}\\tilde{X}^T\\boldsymbol{\\tilde{y}},\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "9fdcab64",
"metadata": {
"editable": true
},
"source": [
"where $\\boldsymbol{\\tilde{y}} = \\boldsymbol{y} - \\overline{\\boldsymbol{y}}$\n",
"and $\\tilde{X}_{ij} = X_{ij} - \\frac{1}{n}\\sum_{k=0}^{n-1}X_{kj}$.\n",
"\n",
"For Ridge regression we need to add $\\lambda \\boldsymbol{\\beta}^T\\boldsymbol{\\beta}$ to the cost function and get then"
]
},
{
"cell_type": "markdown",
"id": "1ce3f20e",
"metadata": {
"editable": true
},
"source": [
"$$\n",
"\\hat{\\boldsymbol{\\beta}} = (\\tilde{X}^T\\tilde{X} + \\lambda I)^{-1}\\tilde{X}^T\\boldsymbol{\\tilde{y}}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "0431630c",
"metadata": {
"editable": true
},
"source": [
"Now we try to implement this."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "9da6e338",
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"\n",
"np.random.seed(2018)\n",
"n = 100\n",
"d = 3\n",
"Lambda = 0.01\n",
"true_beta = [2, 0.5, 3.7]\n",
"\n",
"# Make data set.\n",
"x = np.linspace(-3, 3, n)\n",
"y_real = 2 + 0.5*x + 3.7*x**2\n",
"\n",
"y = np.sum(\n",
" np.asarray([x ** p * b for p, b in enumerate(true_beta)]), \n",
" axis=0) + 0.1 * np.random.normal(size=len(x))\n",
"\n",
"\n",
"#Design matrix X does not include the intercept. \n",
"X = np.zeros((len(x), d))\n",
"for p in range(d-1): \n",
" X[:, p] = x ** (p+1)\n",
"\n",
"\n",
"#Split data in train and test\n",
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"\n",
"# Scale data by subtracting mean value,own implementation\n",
"#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable\n",
"X_train_mean = np.mean(X_train,axis=0)\n",
"#Center by removing mean from each feature\n",
"X_train_scaled = X_train - X_train_mean\n",
"X_test_scaled = X_test - X_train_mean\n",
"#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note)\n",
"y_scaler = np.mean(y_train)\n",
"y_train_scaled = y_train - y_scaler\n",
"\n",
"\n",
"#Calculate beta\n",
"beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)\n",
"beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)\n",
"print(beta_OLS)\n",
"print(beta_Ridge)\n",
"\n",
"interceptOLS = y_scaler - X_train_mean @ beta_OLS\n",
"interceptRidge = y_scaler - X_train_mean @ beta_Ridge\n",
"print(interceptOLS)\n",
"print(interceptRidge)\n",
"#predict value\n",
"ytilde_test_OLS = X_test_scaled @ beta_OLS+y_scaler\n",
"ytilde_test_Ridge = X_test_scaled @ beta_Ridge+y_scaler\n",
"\n",
"\n",
"#Calculate MSE\n",
"\n",
"print(\" \")\n",
"print(\"test MSE of OLS:\")\n",
"print(MSE(y_test,ytilde_test_OLS))\n",
"print(\" \")\n",
"print(\"test MSE of Ridge\")\n",
"print(MSE(y_test,ytilde_test_Ridge))\n",
"\n",
"\n",
"plt.scatter(x,y,label='Data')\n",
"#plt.plot(x,y_real,label='no noise')\n",
"plt.plot(x, X @ beta_OLS+interceptOLS,'*', label=\"OLS_Fit\")\n",
"plt.plot(x, X @ beta_Ridge+interceptRidge, label=\"Ridge_Fit\")\n",
"plt.grid()\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "c2271c6b",
"metadata": {
"editable": true
},
"source": [
"We see that we get the same values for the parameters! As it should be. The MSE may however change (not the case here).\n",
"\n",
"Finally, instead of using our own function we repeat the same example using the **standardscaler** functionality of the library **Scikit-Learn**.\n",
"\n",
"<!-- !bc pycod -->\n",
"np.random.seed(2018)\n",
"n = 100\n",
"d = 3\n",
"Lambda = 0.01\n",
"true_beta = [2, 0.5, 3.7]\n",
"\n",
"<!-- Make data set. -->\n",
"x = np.linspace(-3, 3, n)\n",
"y_real = 2 + 0.5*x + 3.7*x**2\n",
"\n",
"y = np.sum(\n",
" np.asarray([x ** p * b for p, b in enumerate(true_beta)]), \n",
" axis=0) + 0.1 * np.random.normal(size=len(x))\n",
"\n",
"<!-- Design matrix X does not include the intercept. -->\n",
"X = np.zeros((len(x), d))\n",
"for p in range(d-1): \n",
" X[:, p] = x ** (p+1)\n",
"\n",
"<!-- Split data in train and test -->\n",
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"\n",
"<!-- Scale data by subtracting mean value,own implementation -->\n",
"<!-- For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable -->\n",
"X_train_mean = np.mean(X_train,axis=0)\n",
"<!-- Center by removing mean from each feature -->\n",
"X_train_scaled = X_train - X_train_mean\n",
"X_test_scaled = X_test - X_train_mean\n",
"<!-- The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note) -->\n",
"y_scaler = np.mean(y_train)\n",
"y_train_scaled = y_train - y_scaler\n",
"\n",
"<!-- Calculate beta -->\n",
"beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)\n",
"beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)\n",
"print(beta_OLS)\n",
"print(beta_Ridge)\n",
"\n",
"interceptOLS = y_scaler - X_train_mean @ beta_OLS\n",
"interceptRidge = y_scaler - X_train_mean @ beta_Ridge\n",
"print(interceptOLS)\n",
"print(interceptRidge)\n",
"<!-- predict value -->\n",
"ytilde_test_OLS = X_test_scaled @ beta_OLS+y_scaler\n",
"ytilde_test_Ridge = X_test_scaled @ beta_Ridge+y_scaler\n",
"\n",
"<!-- Calculate MSE -->\n",
"\n",
"print(\" \")\n",
"print(\"test MSE of OLS:\")\n",
"print(MSE(y_test,ytilde_test_OLS))\n",
"print(\" \")\n",
"print(\"test MSE of Ridge\")\n",
"print(MSE(y_test,ytilde_test_Ridge))\n",
"\n",
"plt.scatter(x,y,label='Data')\n",
"<!-- plt.plot(x,y_real,label='no noise') -->\n",
"plt.plot(x, X @ beta_OLS+interceptOLS,'*', label=\"OLS_Fit\")\n",
"plt.plot(x, X @ beta_Ridge+interceptRidge, label=\"Ridge_Fit\")\n",
"plt.grid()\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"<!-- !ec -->"
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 5
}
doc/src/week37/programs/noscaling.py
+72
View File
@@ -0,0 +1,72 @@
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
def MSE(y_data,y_model):
n = np.size(y_model)
return np.sum((y_data-y_model)**2)/n
def OLS_fit_beta(X, y):
return np.linalg.pinv(X.T @ X) @ X.T @ y
def Ridge_fit_beta(X, y,L,d):
I = np.eye(d,d)
return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y
np.random.seed(2018)
n = 100
d = 3
L = 0.001
true_beta = [2, 0.5, 3.7]
# Make data set.
x = np.linspace(-3, 3, n)
y_real = 2 + 0.5*x + 3.7*x**2
y = np.sum(
np.asarray([x ** p * b for p, b in enumerate(true_beta)]),
axis=0) + 0.1 * np.random.normal(size=len(x))
#Design matrix X including the intercept
X = np.zeros((len(x), d))
for p in range(d): # (d-1)
X[:, p] = x ** (p) # (p+1 if not intercept included)
#Split datamatrix
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
#Calculate beta, own code
beta_OLS = OLS_fit_beta(X_train, y_train)
beta_Ridge = Ridge_fit_beta(X_train, y_train,L,d)
print(beta_OLS)
print(beta_Ridge)
#predict value
ytilde_test_OLS = X_test @ beta_OLS
ytilde_test_Ridge = X_test @ beta_Ridge
#Calculate MSE
print(" ")
print("test MSE of OLS:")
print(MSE(y_test,ytilde_test_OLS))
print(" ")
print("test MSE of Ridge")
print(MSE(y_test,ytilde_test_Ridge))
plt.scatter(x,y,label='Data')
#plt.plot(x,y_real,label='no noise')
plt.plot(x, X @ beta_OLS,'*', label="OLS_Fit")
plt.plot(x, X @ beta_Ridge, label="Ridge_Fit")
plt.grid()
plt.legend()
plt.show()
doc/src/week37/programs/pfaffian.py
+10
View File
@@ -0,0 +1,10 @@
from pfapack import pfaffian as pf
import numpy.matlib
A = numpy.matlib.rand(100, 100)
A = A - A.T
pfa1 = pf.pfaffian(A)
pfa2 = pf.pfaffian(A, method="H")
pfa3 = pf.pfaffian_schur(A)
print(pfa1, pfa2, pfa3)
doc/src/week37/programs/scaling.py
+86
View File
@@ -0,0 +1,86 @@
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
def MSE(y_data,y_model):
n = np.size(y_model)
return np.sum((y_data-y_model)**2)/n
def OLS_fit_beta(X, y):
return np.linalg.pinv(X.T @ X) @ X.T @ y
def Ridge_fit_beta(X, y,L,d):
I = np.eye(d,d)
return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y
np.random.seed(2018)
n = 1000
d = 3
L = 0.001
true_beta = [2, 0.5, 3.7]
# Make data set.
x = np.linspace(-3, 3, n)
y_real = 2 + 0.5*x + 3.7*x**2
y = np.sum(
np.asarray([x ** p * b for p, b in enumerate(true_beta)]),
axis=0) + 0.1 * np.random.normal(size=len(x))
#Design matrix X including the intercept
X = np.zeros((len(x), d))
for p in range(d): # (d-1)
X[:, p] = x ** (p+1) # (p+1 if not intercept included)
#Split datamatrix
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
# Scale data by subtracting mean value,own function
#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable
X_train_mean = np.mean(X_train,axis=0)
#Center by removing mean from each feature
X_train_scaled = X_train - X_train_mean
X_test_scaled = X_test - X_train_mean
#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note)
y_scaler = np.mean(y_train)
y_train_scaled = y_train - y_scaler
#Calculate beta
beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)
beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,L,d)
print(beta_OLS)
print(beta_Ridge)
interceptOLS = y_scaler - X_train_mean @ beta_OLS
interceptRidge = y_scaler - X_train_mean @ beta_Ridge
print(interceptOLS)
print(interceptRidge)
#predict value
ytilde_test_OLS = X_test_scaled @ beta_OLS+y_scaler
ytilde_test_Ridge = X_test_scaled @ beta_Ridge+y_scaler
#Calculate MSE
print(" ")
print("test MSE of OLS:")
print(MSE(y_test,ytilde_test_OLS))
print(" ")
print("test MSE of Ridge")
print(MSE(y_test,ytilde_test_Ridge))
plt.scatter(x,y,label='Data')
#plt.plot(x,y_real,label='no noise')
plt.plot(x, X @ beta_OLS+interceptOLS,'*', label="OLS_Fit")
plt.plot(x, X @ beta_Ridge+interceptRidge, label="Ridge_Fit")
plt.grid()
plt.legend()
plt.show()
doc/src/week37/programs/sklearnscaling.py
+88
View File
@@ -0,0 +1,88 @@
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
def MSE(y_data,y_model):
n = np.size(y_model)
return np.sum((y_data-y_model)**2)/n
def OLS_fit_beta(X, y):
return np.linalg.pinv(X.T @ X) @ X.T @ y
def Ridge_fit_beta(X, y,L,d):
I = np.eye(d,d)
return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y
np.random.seed(2018)
n = 1000
d = 3
L = 0.001
true_beta = [2, 0.5, 3.7]
# Make data set.
x = np.linspace(-3, 3, n)
y_real = 2 + 0.5*x + 3.7*x**2
y = np.sum(
np.asarray([x ** p * b for p, b in enumerate(true_beta)]),
axis=0) + 0.1 * np.random.normal(size=len(x))
#Design matrix X does include the intercept
X = np.zeros((len(x), d))
for p in range(d): # (d-1)
X[:, p] = x ** (p+1) # (p+1 if not intercept included)
#Split datamatrix
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
scaler = StandardScaler()
yscaler = StandardScaler()
scaler.fit(X_train)
yscaler.fit(y_train)
X_train_scaled = scaler.transform(X_train)
X_test_scaled = scaler.transform(X_test)
y_train_scaled = yscaler.transform(y_train)
y_test_scaled = yscaler.transform(y_test)
#Calculate beta
beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)
beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,L,d)
print(beta_OLS)
print(beta_Ridge)
"""
interceptOLS = y_scaler - X_train_mean @ beta_OLS
interceptRidge = y_scaler - X_train_mean @ beta_Ridge
print(interceptOLS)
print(interceptRidge)
"""
#predict value
ytilde_test_OLS = X_test_scaled @ beta_OLS
ytilde_test_Ridge = X_test_scaled @ beta_Ridge
#Calculate MSE
print(" ")
print("test MSE of OLS:")
print(MSE(y_test,ytilde_test_OLS))
print(" ")
print("test MSE of Ridge")
print(MSE(y_test,ytilde_test_Ridge))
plt.scatter(x,y,label='Data')
#plt.plot(x,y_real,label='no noise')
plt.plot(x, X @ beta_OLS,'*', label="OLS_Fit")
plt.plot(x, X @ beta_Ridge, label="Ridge_Fit")
plt.grid()
plt.legend()
plt.show()