Morten Hjorth-Jensen
3008 lines
228 KiB
HTML
3008 lines
228 KiB
HTML
|
||
<!DOCTYPE html>
|
||
|
||
<html>
|
||
<head>
|
||
<meta charset="utf-8" />
|
||
<meta name="viewport" content="width=device-width, initial-scale=1.0" />
|
||
<title>Week 36: Linear Regression and Statistical interpretations — Applied Data Analysis and Machine Learning</title>
|
||
|
||
<link href="_static/css/theme.css" rel="stylesheet">
|
||
<link href="_static/css/index.ff1ffe594081f20da1ef19478df9384b.css" rel="stylesheet">
|
||
|
||
|
||
<link rel="stylesheet"
|
||
href="_static/vendor/fontawesome/5.13.0/css/all.min.css">
|
||
<link rel="preload" as="font" type="font/woff2" crossorigin
|
||
href="_static/vendor/fontawesome/5.13.0/webfonts/fa-solid-900.woff2">
|
||
<link rel="preload" as="font" type="font/woff2" crossorigin
|
||
href="_static/vendor/fontawesome/5.13.0/webfonts/fa-brands-400.woff2">
|
||
|
||
|
||
|
||
|
||
|
||
<link rel="stylesheet" type="text/css" href="_static/pygments.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/sphinx-book-theme.css?digest=c3fdc42140077d1ad13ad2f1588a4309" />
|
||
<link rel="stylesheet" type="text/css" href="_static/togglebutton.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/copybutton.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/mystnb.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/sphinx-thebe.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/panels-main.c949a650a448cc0ae9fd3441c0e17fb0.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/panels-variables.06eb56fa6e07937060861dad626602ad.css" />
|
||
|
||
<link rel="preload" as="script" href="_static/js/index.be7d3bbb2ef33a8344ce.js">
|
||
|
||
<script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
|
||
<script src="_static/jquery.js"></script>
|
||
<script src="_static/underscore.js"></script>
|
||
<script src="_static/doctools.js"></script>
|
||
<script src="_static/clipboard.min.js"></script>
|
||
<script src="_static/copybutton.js"></script>
|
||
<script>let toggleHintShow = 'Click to show';</script>
|
||
<script>let toggleHintHide = 'Click to hide';</script>
|
||
<script>let toggleOpenOnPrint = 'true';</script>
|
||
<script src="_static/togglebutton.js"></script>
|
||
<script>var togglebuttonSelector = '.toggle, .admonition.dropdown, .tag_hide_input div.cell_input, .tag_hide-input div.cell_input, .tag_hide_output div.cell_output, .tag_hide-output div.cell_output, .tag_hide_cell.cell, .tag_hide-cell.cell';</script>
|
||
<script src="_static/sphinx-book-theme.d59cb220de22ca1c485ebbdc042f0030.js"></script>
|
||
<script>const THEBE_JS_URL = "https://unpkg.com/thebe@0.8.2/lib/index.js"
|
||
const thebe_selector = ".thebe,.cell"
|
||
const thebe_selector_input = "pre"
|
||
const thebe_selector_output = ".output, .cell_output"
|
||
</script>
|
||
<script async="async" src="_static/sphinx-thebe.js"></script>
|
||
<script>window.MathJax = {"options": {"processHtmlClass": "tex2jax_process|mathjax_process|math|output_area"}}</script>
|
||
<script defer="defer" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
|
||
<link rel="index" title="Index" href="genindex.html" />
|
||
<link rel="search" title="Search" href="search.html" />
|
||
<link rel="next" title="Exercises week 37" href="exercisesweek37.html" />
|
||
<link rel="prev" title="Exercises week 36" href="exercisesweek36.html" />
|
||
<meta name="viewport" content="width=device-width, initial-scale=1" />
|
||
<meta name="docsearch:language" content="None">
|
||
|
||
|
||
<!-- Google Analytics -->
|
||
|
||
</head>
|
||
<body data-spy="scroll" data-target="#bd-toc-nav" data-offset="80">
|
||
|
||
<div class="container-fluid" id="banner"></div>
|
||
|
||
|
||
|
||
<div class="container-xl">
|
||
<div class="row">
|
||
|
||
<div class="col-12 col-md-3 bd-sidebar site-navigation show" id="site-navigation">
|
||
|
||
<div class="navbar-brand-box">
|
||
<a class="navbar-brand text-wrap" href="index.html">
|
||
|
||
<!-- `logo` is deprecated in Sphinx 4.0, so remove this when we stop supporting 3 -->
|
||
|
||
|
||
|
||
<img src="_static/logo.png" class="logo" alt="logo">
|
||
|
||
|
||
<h1 class="site-logo" id="site-title">Applied Data Analysis and Machine Learning</h1>
|
||
|
||
</a>
|
||
</div><form class="bd-search d-flex align-items-center" action="search.html" method="get">
|
||
<i class="icon fas fa-search"></i>
|
||
<input type="search" class="form-control" name="q" id="search-input" placeholder="Search this book..." aria-label="Search this book..." autocomplete="off" >
|
||
</form><nav class="bd-links" id="bd-docs-nav" aria-label="Main">
|
||
<div class="bd-toc-item active">
|
||
<ul class="nav bd-sidenav">
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="intro.html">
|
||
Applied Data Analysis and Machine Learning
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
<p aria-level="2" class="caption" role="heading">
|
||
<span class="caption-text">
|
||
About the course
|
||
</span>
|
||
</p>
|
||
<ul class="nav bd-sidenav">
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="schedule.html">
|
||
Teaching schedule with links to material
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="teachers.html">
|
||
Teachers and Grading
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="textbooks.html">
|
||
Textbooks
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
<p aria-level="2" class="caption" role="heading">
|
||
<span class="caption-text">
|
||
Review of Statistics with Resampling Techniques and Linear Algebra
|
||
</span>
|
||
</p>
|
||
<ul class="nav bd-sidenav">
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="statistics.html">
|
||
1. Elements of Probability Theory and Statistical Data Analysis
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="linalg.html">
|
||
2. Linear Algebra, Handling of Arrays and more Python Features
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
<p aria-level="2" class="caption" role="heading">
|
||
<span class="caption-text">
|
||
From Regression to Support Vector Machines
|
||
</span>
|
||
</p>
|
||
<ul class="nav bd-sidenav">
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter1.html">
|
||
3. Linear Regression
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter2.html">
|
||
4. Ridge and Lasso Regression
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter3.html">
|
||
5. Resampling Methods
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter4.html">
|
||
6. Logistic Regression
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapteroptimization.html">
|
||
7. Optimization, the central part of any Machine Learning algortithm
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter5.html">
|
||
8. Support Vector Machines, overarching aims
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
<p aria-level="2" class="caption" role="heading">
|
||
<span class="caption-text">
|
||
Decision Trees, Ensemble Methods and Boosting
|
||
</span>
|
||
</p>
|
||
<ul class="nav bd-sidenav">
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter6.html">
|
||
9. Decision trees, overarching aims
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter7.html">
|
||
10. Ensemble Methods: From a Single Tree to Many Trees and Extreme Boosting, Meet the Jungle of Methods
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
<p aria-level="2" class="caption" role="heading">
|
||
<span class="caption-text">
|
||
Dimensionality Reduction
|
||
</span>
|
||
</p>
|
||
<ul class="nav bd-sidenav">
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter8.html">
|
||
11. Basic ideas of the Principal Component Analysis (PCA)
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="clustering.html">
|
||
12. Clustering and Unsupervised Learning
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
<p aria-level="2" class="caption" role="heading">
|
||
<span class="caption-text">
|
||
Deep Learning Methods
|
||
</span>
|
||
</p>
|
||
<ul class="nav bd-sidenav">
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter9.html">
|
||
13. Neural networks
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter10.html">
|
||
14. Building a Feed Forward Neural Network
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter11.html">
|
||
15. Solving Differential Equations with Deep Learning
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter12.html">
|
||
16. Convolutional Neural Networks
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="chapter13.html">
|
||
17. Recurrent neural networks: Overarching view
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
<p aria-level="2" class="caption" role="heading">
|
||
<span class="caption-text">
|
||
Weekly material, notes and exercises
|
||
</span>
|
||
</p>
|
||
<ul class="current nav bd-sidenav">
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="exercisesweek34.html">
|
||
Exercises week 34
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="week34.html">
|
||
Week 34: Introduction to the course, Logistics and Practicalities
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="exercisesweek35.html">
|
||
Exercises week 35
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="week35.html">
|
||
Week 35: From Ordinary Linear Regression to Ridge and Lasso Regression
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="exercisesweek36.html">
|
||
Exercises week 36
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1 current active">
|
||
<a class="current reference internal" href="#">
|
||
Week 36: Linear Regression and Statistical interpretations
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="exercisesweek37.html">
|
||
Exercises week 37
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="week37.html">
|
||
Week 37: Statistical interpretations and Resampling Methods
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="exercisesweek38.html">
|
||
Exercises week 38
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="week38.html">
|
||
Week 38: Logistic Regression and Optimization
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="exercisesweek39.html">
|
||
Exercises week 39
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="week39.html">
|
||
Week 39: Optimization and Gradient Methods
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="week40.html">
|
||
Week 40: Gradient descent methods (continued) and start Neural networks
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="exercisesweek41.html">
|
||
Exercises week 41
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="week41.html">
|
||
Week 41 Neural networks and constructing a neural network code
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="exercisesweek42.html">
|
||
Exercises week 42
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="week42.html">
|
||
Week 42 Constructing a Neural Network code with examples
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
<p aria-level="2" class="caption" role="heading">
|
||
<span class="caption-text">
|
||
Projects
|
||
</span>
|
||
</p>
|
||
<ul class="nav bd-sidenav">
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="project1.html">
|
||
Project 1 on Machine Learning, deadline October 7 (midnight), 2024
|
||
</a>
|
||
</li>
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="project2.html">
|
||
Project 2 on Machine Learning, deadline November 4 (Midnight)
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
|
||
</div>
|
||
</nav> <!-- To handle the deprecated key -->
|
||
|
||
<div class="navbar_extra_footer">
|
||
Powered by <a href="https://jupyterbook.org">Jupyter Book</a>
|
||
</div>
|
||
|
||
</div>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<main class="col py-md-3 pl-md-4 bd-content overflow-auto" role="main">
|
||
|
||
<div class="topbar container-xl fixed-top">
|
||
<div class="topbar-contents row">
|
||
<div class="col-12 col-md-3 bd-topbar-whitespace site-navigation show"></div>
|
||
<div class="col pl-md-4 topbar-main">
|
||
|
||
<button id="navbar-toggler" class="navbar-toggler ml-0" type="button" data-toggle="collapse"
|
||
data-toggle="tooltip" data-placement="bottom" data-target=".site-navigation" aria-controls="navbar-menu"
|
||
aria-expanded="true" aria-label="Toggle navigation" aria-controls="site-navigation"
|
||
title="Toggle navigation" data-toggle="tooltip" data-placement="left">
|
||
<i class="fas fa-bars"></i>
|
||
<i class="fas fa-arrow-left"></i>
|
||
<i class="fas fa-arrow-up"></i>
|
||
</button>
|
||
|
||
|
||
<div class="dropdown-buttons-trigger">
|
||
<button id="dropdown-buttons-trigger" class="btn btn-secondary topbarbtn" aria-label="Download this page"><i
|
||
class="fas fa-download"></i></button>
|
||
|
||
<div class="dropdown-buttons">
|
||
<!-- ipynb file if we had a myst markdown file -->
|
||
|
||
<!-- Download raw file -->
|
||
<a class="dropdown-buttons" href="_sources/week36.ipynb"><button type="button"
|
||
class="btn btn-secondary topbarbtn" title="Download source file" data-toggle="tooltip"
|
||
data-placement="left">.ipynb</button></a>
|
||
<!-- Download PDF via print -->
|
||
<button type="button" id="download-print" class="btn btn-secondary topbarbtn" title="Print to PDF"
|
||
onclick="printPdf(this)" data-toggle="tooltip" data-placement="left">.pdf</button>
|
||
</div>
|
||
</div>
|
||
|
||
<!-- Source interaction buttons -->
|
||
|
||
<!-- Full screen (wrap in <a> to have style consistency -->
|
||
|
||
<a class="full-screen-button"><button type="button" class="btn btn-secondary topbarbtn" data-toggle="tooltip"
|
||
data-placement="bottom" onclick="toggleFullScreen()" aria-label="Fullscreen mode"
|
||
title="Fullscreen mode"><i
|
||
class="fas fa-expand"></i></button></a>
|
||
|
||
<!-- Launch buttons -->
|
||
|
||
</div>
|
||
|
||
<!-- Table of contents -->
|
||
<div class="d-none d-md-block col-md-2 bd-toc show noprint">
|
||
|
||
<div class="tocsection onthispage pt-5 pb-3">
|
||
<i class="fas fa-list"></i> Contents
|
||
</div>
|
||
<nav id="bd-toc-nav" aria-label="Page">
|
||
<ul class="visible nav section-nav flex-column">
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#plans-for-week-36">
|
||
Plans for week 36
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#material-for-lecture-monday-september-2">
|
||
Material for lecture Monday September 2
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#important-technicalities-more-on-rescaling-data">
|
||
Important technicalities: More on Rescaling data
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#test-function-for-what-happens-with-ols-ridge-and-lasso">
|
||
Test Function for what happens with OLS, Ridge and Lasso
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#linking-the-regression-analysis-with-a-statistical-interpretation">
|
||
Linking the regression analysis with a statistical interpretation
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#assumptions-made">
|
||
Assumptions made
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#expectation-value-and-variance">
|
||
Expectation value and variance
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#expectation-value-and-variance-for-boldsymbol-beta">
|
||
Expectation value and variance for
|
||
<span class="math notranslate nohighlight">
|
||
\(\boldsymbol{\beta}\)
|
||
</span>
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#deriving-ols-from-a-probability-distribution">
|
||
Deriving OLS from a probability distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#independent-and-identically-distributed-iid">
|
||
Independent and Identically Distributed (iid)
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#maximum-likelihood-estimation-mle">
|
||
Maximum Likelihood Estimation (MLE)
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#a-new-cost-function">
|
||
A new Cost Function
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#more-basic-statistics-and-bayes-theorem">
|
||
More basic Statistics and Bayes’ theorem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#marginal-probability">
|
||
Marginal Probability
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#conditional-probability">
|
||
Conditional Probability
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#bayes-theorem">
|
||
Bayes’ Theorem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#interpretations-of-bayes-theorem">
|
||
Interpretations of Bayes’ Theorem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#example-of-usage-of-bayes-theorem">
|
||
Example of Usage of Bayes’ theorem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#doing-it-correctly">
|
||
Doing it correctly
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#bayes-theorem-and-ridge-and-lasso-regression">
|
||
Bayes’ Theorem and Ridge and Lasso Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#ridge-and-bayes">
|
||
Ridge and Bayes
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#lasso-and-bayes">
|
||
Lasso and Bayes
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#material-for-the-active-learning-sessions-tuesday-and-wednesday">
|
||
Material for the active learning sessions Tuesday and Wednesday
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#linear-regression-and-the-svd">
|
||
Linear Regression and the SVD
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#what-does-it-mean">
|
||
What does it mean?
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#and-finally-boldsymbol-x-boldsymbol-x-t">
|
||
And finally
|
||
<span class="math notranslate nohighlight">
|
||
\(\boldsymbol{X}\boldsymbol{X}^T\)
|
||
</span>
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#code-for-svd-and-inversion-of-matrices">
|
||
Code for SVD and Inversion of Matrices
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#inverse-of-rectangular-matrix">
|
||
Inverse of Rectangular Matrix
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#ridge-and-lasso-regression">
|
||
Ridge and LASSO Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#from-ols-to-ridge-and-lasso">
|
||
From OLS to Ridge and Lasso
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#deriving-the-ridge-regression-equations">
|
||
Deriving the Ridge Regression Equations
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#note-on-scikit-learn">
|
||
Note on Scikit-Learn
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#comparison-with-ols">
|
||
Comparison with OLS
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#svd-analysis">
|
||
SVD analysis
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#interpreting-the-ridge-results">
|
||
Interpreting the Ridge results
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#more-interpretations">
|
||
More interpretations
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#deriving-the-lasso-regression-equations">
|
||
Deriving the Lasso Regression Equations
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#simple-example-to-illustrate-ordinary-least-squares-ridge-and-lasso-regression">
|
||
Simple example to illustrate Ordinary Least Squares, Ridge and Lasso Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#ridge-regression">
|
||
Ridge Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#lasso-regression">
|
||
Lasso Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#yet-another-example">
|
||
Yet another Example
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#the-ols-case">
|
||
The OLS case
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#the-ridge-case">
|
||
The Ridge case
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#writing-the-cost-function">
|
||
Writing the Cost Function
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#lasso-case">
|
||
Lasso case
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#the-first-case">
|
||
The first Case
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#simple-code-for-solving-the-above-problem">
|
||
Simple code for solving the above problem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#with-lasso-regression">
|
||
With Lasso Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#another-example-now-with-a-polynomial-fit">
|
||
Another Example, now with a polynomial fit
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
|
||
</nav>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<div id="main-content" class="row">
|
||
<div class="col-12 col-md-9 pl-md-3 pr-md-0">
|
||
<!-- Table of contents that is only displayed when printing the page -->
|
||
<div id="jb-print-docs-body" class="onlyprint">
|
||
<h1>Week 36: Linear Regression and Statistical interpretations</h1>
|
||
<!-- Table of contents -->
|
||
<div id="print-main-content">
|
||
<div id="jb-print-toc">
|
||
|
||
<div>
|
||
<h2> Contents </h2>
|
||
</div>
|
||
<nav aria-label="Page">
|
||
<ul class="visible nav section-nav flex-column">
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#plans-for-week-36">
|
||
Plans for week 36
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#material-for-lecture-monday-september-2">
|
||
Material for lecture Monday September 2
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#important-technicalities-more-on-rescaling-data">
|
||
Important technicalities: More on Rescaling data
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#test-function-for-what-happens-with-ols-ridge-and-lasso">
|
||
Test Function for what happens with OLS, Ridge and Lasso
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#linking-the-regression-analysis-with-a-statistical-interpretation">
|
||
Linking the regression analysis with a statistical interpretation
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#assumptions-made">
|
||
Assumptions made
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#expectation-value-and-variance">
|
||
Expectation value and variance
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#expectation-value-and-variance-for-boldsymbol-beta">
|
||
Expectation value and variance for
|
||
<span class="math notranslate nohighlight">
|
||
\(\boldsymbol{\beta}\)
|
||
</span>
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#deriving-ols-from-a-probability-distribution">
|
||
Deriving OLS from a probability distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#independent-and-identically-distributed-iid">
|
||
Independent and Identically Distributed (iid)
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#maximum-likelihood-estimation-mle">
|
||
Maximum Likelihood Estimation (MLE)
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#a-new-cost-function">
|
||
A new Cost Function
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#more-basic-statistics-and-bayes-theorem">
|
||
More basic Statistics and Bayes’ theorem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#marginal-probability">
|
||
Marginal Probability
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#conditional-probability">
|
||
Conditional Probability
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#bayes-theorem">
|
||
Bayes’ Theorem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#interpretations-of-bayes-theorem">
|
||
Interpretations of Bayes’ Theorem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#example-of-usage-of-bayes-theorem">
|
||
Example of Usage of Bayes’ theorem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#doing-it-correctly">
|
||
Doing it correctly
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#bayes-theorem-and-ridge-and-lasso-regression">
|
||
Bayes’ Theorem and Ridge and Lasso Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#ridge-and-bayes">
|
||
Ridge and Bayes
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#lasso-and-bayes">
|
||
Lasso and Bayes
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#material-for-the-active-learning-sessions-tuesday-and-wednesday">
|
||
Material for the active learning sessions Tuesday and Wednesday
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#linear-regression-and-the-svd">
|
||
Linear Regression and the SVD
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#what-does-it-mean">
|
||
What does it mean?
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#and-finally-boldsymbol-x-boldsymbol-x-t">
|
||
And finally
|
||
<span class="math notranslate nohighlight">
|
||
\(\boldsymbol{X}\boldsymbol{X}^T\)
|
||
</span>
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#code-for-svd-and-inversion-of-matrices">
|
||
Code for SVD and Inversion of Matrices
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#inverse-of-rectangular-matrix">
|
||
Inverse of Rectangular Matrix
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#ridge-and-lasso-regression">
|
||
Ridge and LASSO Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#from-ols-to-ridge-and-lasso">
|
||
From OLS to Ridge and Lasso
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#deriving-the-ridge-regression-equations">
|
||
Deriving the Ridge Regression Equations
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#note-on-scikit-learn">
|
||
Note on Scikit-Learn
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#comparison-with-ols">
|
||
Comparison with OLS
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#svd-analysis">
|
||
SVD analysis
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#interpreting-the-ridge-results">
|
||
Interpreting the Ridge results
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#more-interpretations">
|
||
More interpretations
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#deriving-the-lasso-regression-equations">
|
||
Deriving the Lasso Regression Equations
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#simple-example-to-illustrate-ordinary-least-squares-ridge-and-lasso-regression">
|
||
Simple example to illustrate Ordinary Least Squares, Ridge and Lasso Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#ridge-regression">
|
||
Ridge Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#lasso-regression">
|
||
Lasso Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#yet-another-example">
|
||
Yet another Example
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#the-ols-case">
|
||
The OLS case
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#the-ridge-case">
|
||
The Ridge case
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#writing-the-cost-function">
|
||
Writing the Cost Function
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#lasso-case">
|
||
Lasso case
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#the-first-case">
|
||
The first Case
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#simple-code-for-solving-the-above-problem">
|
||
Simple code for solving the above problem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#with-lasso-regression">
|
||
With Lasso Regression
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#another-example-now-with-a-polynomial-fit">
|
||
Another Example, now with a polynomial fit
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
|
||
</nav>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<div>
|
||
|
||
<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)
|
||
doconce format html week36.do.txt --no_mako -->
|
||
<!-- dom:TITLE: Week 36: Linear Regression and Statistical interpretations --><div class="tex2jax_ignore mathjax_ignore section" id="week-36-linear-regression-and-statistical-interpretations">
|
||
<h1>Week 36: Linear Regression and Statistical interpretations<a class="headerlink" href="#week-36-linear-regression-and-statistical-interpretations" title="Permalink to this headline">¶</a></h1>
|
||
<p><strong>Morten Hjorth-Jensen</strong>, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</p>
|
||
<p>Date: <strong>September 2-6, 2024</strong></p>
|
||
<div class="section" id="plans-for-week-36">
|
||
<h2>Plans for week 36<a class="headerlink" href="#plans-for-week-36" title="Permalink to this headline">¶</a></h2>
|
||
<ul class="simple">
|
||
<li><p>Material for the lecture on Monday September 2</p>
|
||
<ul>
|
||
<li><p>Technicalities related to scaling and other issues with data handling</p></li>
|
||
<li><p>Linear Regression, ordinary least squares, Ridge and Lasso and links with Statistics</p></li>
|
||
<li><p>Recommended Reading: Goodfellow et al chapter 3 (till 3.11) on probability theory, see <a class="reference external" href="https://www.deeplearningbook.org/">https://www.deeplearningbook.org/</a></p></li>
|
||
<li><p>Raschka et al, chapter 4 pages 105-134 and chapter 6 pages 171-185. Chapter 4 and 6 contain many useful hints which will also be relevant for the various projects as well.</p></li>
|
||
<li><p><a class="reference external" href="https://youtu.be/oHjqjUB36KE">Video of lecture</a></p></li>
|
||
<li><p><a class="reference external" href="https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2024/NotesSeptember2.pdf">Whiteboard notes</a></p></li>
|
||
</ul>
|
||
</li>
|
||
<li><p>Material for the active learning sessions on Tuesday and Wednesday</p>
|
||
<ul>
|
||
<li><p>Summary from last week on discussion of SVD, Ridge and Lasso linear regression.</p></li>
|
||
<li><p>Raschka et al, chapter 4 pages 105-134 and chapter 6 pages 171-185. Chapter 4 and 6, as for the lecture discussion</p></li>
|
||
<li><p>Presentation and discussion of first project</p></li>
|
||
</ul>
|
||
</li>
|
||
</ul>
|
||
</div>
|
||
<div class="section" id="material-for-lecture-monday-september-2">
|
||
<h2>Material for lecture Monday September 2<a class="headerlink" href="#material-for-lecture-monday-september-2" title="Permalink to this headline">¶</a></h2>
|
||
</div>
|
||
<div class="section" id="important-technicalities-more-on-rescaling-data">
|
||
<h2>Important technicalities: More on Rescaling data<a class="headerlink" href="#important-technicalities-more-on-rescaling-data" title="Permalink to this headline">¶</a></h2>
|
||
<p>When you are comparing your own code with for example <strong>Scikit-Learn</strong>’s
|
||
library, there are some technicalities to keep in mind. The examples
|
||
here demonstrate some of these aspects with potential pitfalls.</p>
|
||
<p>The discussion here focuses on the role of the intercept, how we can
|
||
set up the design matrix, what scaling we should use and other topics
|
||
which tend confuse us.</p>
|
||
<p>The intercept can be interpreted as the expected value of our
|
||
target/output variables when all other predictors are set to zero.
|
||
Thus, if we cannot assume that the expected outputs/targets are zero
|
||
when all predictors are zero (the columns in the design matrix), it
|
||
may be a bad idea to implement a model which penalizes the intercept.
|
||
Furthermore, in for example Ridge and Lasso regression, the default solutions
|
||
from the library <strong>Scikit-Learn</strong> (when not shrinking <span class="math notranslate nohighlight">\(\beta_0\)</span>) for the unknown parameters
|
||
<span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>, are derived under the assumption that both <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> and
|
||
<span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> are zero centered, that is we subtract the mean values.</p>
|
||
<p>If our predictors represent different scales, then it is important to
|
||
standardize the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> by subtracting the mean of each
|
||
column from the corresponding column and dividing the column with its
|
||
standard deviation. Most machine learning libraries do this as a default. This means that if you compare your code with the results from a given library,
|
||
the results may differ.</p>
|
||
<p>The
|
||
<a class="reference external" href="https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html">Standardscaler</a>
|
||
function in <strong>Scikit-Learn</strong> does this for us. For the data sets we
|
||
have been studying in our various examples, the data are in many cases
|
||
already scaled and there is no need to scale them. You as a user of different machine learning algorithms, should always perform a
|
||
survey of your data, with a critical assessment of them in case you need to scale the data.</p>
|
||
<p>If you need to scale the data, not doing so will give an <em>unfair</em>
|
||
penalization of the parameters since their magnitude depends on the
|
||
scale of their corresponding predictor.</p>
|
||
<p>The <strong>Scikit-Learn</strong> site <a class="reference external" href="https://scikit-learn.org/stable/auto_examples/preprocessing/plot_all_scaling.html#plot-all-scaling-standard-scaler-section">https://scikit-learn.org/stable/auto_examples/preprocessing/plot_all_scaling.html#plot-all-scaling-standard-scaler-section</a> has a good discussion of different ways of preprocessing data.</p>
|
||
<p>Suppose as an example that you
|
||
you have an input variable given by the heights of different persons.
|
||
Human height might be measured in inches or meters or
|
||
kilometers. If measured in kilometers, a standard linear regression
|
||
model with this predictor would probably give a much bigger
|
||
coefficient term, than if measured in millimeters.
|
||
This can clearly lead to problems in evaluating the cost/loss functions.</p>
|
||
<p>Keep in mind that when you transform your data set before training a model, the same transformation needs to be done
|
||
on your eventual new data set before making a prediction. If we translate this into a Python code, it would could be implemented as</p>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="sd">"""</span>
|
||
<span class="sd">#Model training, we compute the mean value of y and X</span>
|
||
<span class="sd">y_train_mean = np.mean(y_train)</span>
|
||
<span class="sd">X_train_mean = np.mean(X_train,axis=0)</span>
|
||
<span class="sd">X_train = X_train - X_train_mean</span>
|
||
<span class="sd">y_train = y_train - y_train_mean</span>
|
||
|
||
<span class="sd"># The we fit our model with the training data</span>
|
||
<span class="sd">trained_model = some_model.fit(X_train,y_train)</span>
|
||
|
||
|
||
<span class="sd">#Model prediction, we need also to transform our data set used for the prediction.</span>
|
||
<span class="sd">X_test = X_test - X_train_mean #Use mean from training data</span>
|
||
<span class="sd">y_pred = trained_model(X_test)</span>
|
||
<span class="sd">y_pred = y_pred + y_train_mean</span>
|
||
<span class="sd">"""</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<div class="output text_plain highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>'\n#Model training, we compute the mean value of y and X\ny_train_mean = np.mean(y_train)\nX_train_mean = np.mean(X_train,axis=0)\nX_train = X_train - X_train_mean\ny_train = y_train - y_train_mean\n\n# The we fit our model with the training data\ntrained_model = some_model.fit(X_train,y_train)\n\n\n#Model prediction, we need also to transform our data set used for the prediction.\nX_test = X_test - X_train_mean #Use mean from training data\ny_pred = trained_model(X_test)\ny_pred = y_pred + y_train_mean\n'
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<p>Let us try to understand what this may imply mathematically when we
|
||
subtract the mean values, also known as <em>zero centering</em>. For
|
||
simplicity, we will focus on ordinary regression, as done in the above example.</p>
|
||
<p>The cost/loss function for regression is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
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,.
|
||
\]</div>
|
||
<p>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.</p>
|
||
<p>What we have done is to single out the <span class="math notranslate nohighlight">\(\beta_0\)</span> term in the
|
||
definition of the mean squared error (MSE). The design matrix <span class="math notranslate nohighlight">\(X\)</span>
|
||
does in this case not contain any intercept column. When we take the
|
||
derivative with respect to <span class="math notranslate nohighlight">\(\beta_0\)</span>, we want the derivative to obey</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\frac{\partial C}{\partial \beta_j} = 0,
|
||
\]</div>
|
||
<p>for all <span class="math notranslate nohighlight">\(j\)</span>. For <span class="math notranslate nohighlight">\(\beta_0\)</span> we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\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).
|
||
\]</div>
|
||
<p>Multiplying away the constant <span class="math notranslate nohighlight">\(2/n\)</span>, we obtain</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\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.
|
||
\]</div>
|
||
<p>Let us specialize first to the case where we have only two parameters <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span>.
|
||
Our result for <span class="math notranslate nohighlight">\(\beta_0\)</span> simplifies then to</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
n\beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \beta_1.
|
||
\]</div>
|
||
<p>We obtain then</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \beta_1\frac{1}{n}\sum_{i=0}^{n-1} X_{i1}.
|
||
\]</div>
|
||
<p>If we define</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\mu_{\boldsymbol{x}_1}=\frac{1}{n}\sum_{i=0}^{n-1} X_{i1},
|
||
\]</div>
|
||
<p>and the mean value of the outputs as</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\mu_y=\frac{1}{n}\sum_{i=0}^{n-1}y_i,
|
||
\]</div>
|
||
<p>we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\beta_0 = \mu_y - \beta_1\mu_{\boldsymbol{x}_1}.
|
||
\]</div>
|
||
<p>In the general case with more parameters than <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span>, we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\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.
|
||
\]</div>
|
||
<p>We can rewrite the latter equation as</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \sum_{j=1}^{p-1} \mu_{\boldsymbol{x}_j}\beta_j,
|
||
\]</div>
|
||
<p>where we have defined</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\mu_{\boldsymbol{x}_j}=\frac{1}{n}\sum_{i=0}^{n-1} X_{ij},
|
||
\]</div>
|
||
<p>the mean value for all elements of the column vector <span class="math notranslate nohighlight">\(\boldsymbol{x}_j\)</span>.</p>
|
||
<p>Replacing <span class="math notranslate nohighlight">\(y_i\)</span> with <span class="math notranslate nohighlight">\(y_i - y_i - \overline{\boldsymbol{y}}\)</span> and centering also our design matrix results in a cost function (in vector-matrix disguise)</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta}).
|
||
\]</div>
|
||
<p>If we minimize with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> we have then</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\hat{\boldsymbol{\beta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}},
|
||
\]</div>
|
||
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\boldsymbol{y}}\)</span>
|
||
and <span class="math notranslate nohighlight">\(\tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj}\)</span>.</p>
|
||
<p>For Ridge regression we need to add <span class="math notranslate nohighlight">\(\lambda \boldsymbol{\beta}^T\boldsymbol{\beta}\)</span> to the cost function and get then</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\hat{\boldsymbol{\beta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}.
|
||
\]</div>
|
||
<p>What does this mean? And why do we insist on all this? Let us look at some examples.</p>
|
||
<p>This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (<em>code example thanks to Øyvind Sigmundson Schøyen</em>). Here our scaling of the data is done by subtracting the mean values only.
|
||
Note also that we do not split the data into training and test.</p>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="o">%</span><span class="k">matplotlib</span> inline
|
||
|
||
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
|
||
<span class="kn">from</span> <span class="nn">sklearn.linear_model</span> <span class="kn">import</span> <span class="n">LinearRegression</span>
|
||
|
||
|
||
<span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">seed</span><span class="p">(</span><span class="mi">2021</span><span class="p">)</span>
|
||
|
||
<span class="k">def</span> <span class="nf">MSE</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span><span class="n">y_model</span><span class="p">):</span>
|
||
<span class="n">n</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">size</span><span class="p">(</span><span class="n">y_model</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span><span class="o">-</span><span class="n">y_model</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">/</span><span class="n">n</span>
|
||
|
||
|
||
<span class="k">def</span> <span class="nf">fit_beta</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">):</span>
|
||
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">pinv</span><span class="p">(</span><span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X</span><span class="p">)</span> <span class="o">@</span> <span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y</span>
|
||
|
||
|
||
<span class="n">true_beta</span> <span class="o">=</span> <span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mf">0.5</span><span class="p">,</span> <span class="mf">3.7</span><span class="p">]</span>
|
||
|
||
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">11</span><span class="p">)</span>
|
||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span>
|
||
<span class="n">np</span><span class="o">.</span><span class="n">asarray</span><span class="p">([</span><span class="n">x</span> <span class="o">**</span> <span class="n">p</span> <span class="o">*</span> <span class="n">b</span> <span class="k">for</span> <span class="n">p</span><span class="p">,</span> <span class="n">b</span> <span class="ow">in</span> <span class="nb">enumerate</span><span class="p">(</span><span class="n">true_beta</span><span class="p">)]),</span> <span class="n">axis</span><span class="o">=</span><span class="mi">0</span>
|
||
<span class="p">)</span> <span class="o">+</span> <span class="mf">0.1</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">))</span>
|
||
|
||
<span class="n">degree</span> <span class="o">=</span> <span class="mi">3</span>
|
||
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">),</span> <span class="n">degree</span><span class="p">))</span>
|
||
|
||
<span class="c1"># Include the intercept in the design matrix</span>
|
||
<span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">degree</span><span class="p">):</span>
|
||
<span class="n">X</span><span class="p">[:,</span> <span class="n">p</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span> <span class="o">**</span> <span class="n">p</span>
|
||
|
||
<span class="n">beta</span> <span class="o">=</span> <span class="n">fit_beta</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">)</span>
|
||
|
||
<span class="c1"># Intercept is included in the design matrix</span>
|
||
<span class="n">skl</span> <span class="o">=</span> <span class="n">LinearRegression</span><span class="p">(</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">)</span>
|
||
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"True beta: </span><span class="si">{</span><span class="n">true_beta</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Fitted beta: </span><span class="si">{</span><span class="n">beta</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Sklearn fitted beta: </span><span class="si">{</span><span class="n">skl</span><span class="o">.</span><span class="n">coef_</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||
<span class="n">ypredictOwn</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta</span>
|
||
<span class="n">ypredictSKL</span> <span class="o">=</span> <span class="n">skl</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"MSE with intercept column"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">ypredictOwn</span><span class="p">))</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"MSE with intercept column from SKL"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">ypredictSKL</span><span class="p">))</span>
|
||
|
||
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">scatter</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">"Data"</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">"Fit"</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">skl</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X</span><span class="p">),</span> <span class="n">label</span><span class="o">=</span><span class="s2">"Sklearn (fit_intercept=False)"</span><span class="p">)</span>
|
||
|
||
|
||
<span class="c1"># Do not include the intercept in the design matrix</span>
|
||
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">),</span> <span class="n">degree</span> <span class="o">-</span> <span class="mi">1</span><span class="p">))</span>
|
||
|
||
<span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">degree</span> <span class="o">-</span> <span class="mi">1</span><span class="p">):</span>
|
||
<span class="n">X</span><span class="p">[:,</span> <span class="n">p</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span> <span class="o">**</span> <span class="p">(</span><span class="n">p</span> <span class="o">+</span> <span class="mi">1</span><span class="p">)</span>
|
||
|
||
<span class="c1"># Intercept is not included in the design matrix</span>
|
||
<span class="n">skl</span> <span class="o">=</span> <span class="n">LinearRegression</span><span class="p">(</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">)</span>
|
||
|
||
<span class="c1"># Use centered values for X and y when computing coefficients</span>
|
||
<span class="n">y_offset</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">average</span><span class="p">(</span><span class="n">y</span><span class="p">,</span> <span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
|
||
<span class="n">X_offset</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">average</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
|
||
|
||
<span class="n">beta</span> <span class="o">=</span> <span class="n">fit_beta</span><span class="p">(</span><span class="n">X</span> <span class="o">-</span> <span class="n">X_offset</span><span class="p">,</span> <span class="n">y</span> <span class="o">-</span> <span class="n">y_offset</span><span class="p">)</span>
|
||
<span class="n">intercept</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_offset</span> <span class="o">-</span> <span class="n">X_offset</span> <span class="o">@</span> <span class="n">beta</span><span class="p">)</span>
|
||
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Manual intercept: </span><span class="si">{</span><span class="n">intercept</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Fitted beta (without intercept): </span><span class="si">{</span><span class="n">beta</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Sklearn intercept: </span><span class="si">{</span><span class="n">skl</span><span class="o">.</span><span class="n">intercept_</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"Sklearn fitted beta (without intercept): </span><span class="si">{</span><span class="n">skl</span><span class="o">.</span><span class="n">coef_</span><span class="si">}</span><span class="s2">"</span><span class="p">)</span>
|
||
<span class="n">ypredictOwn</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta</span>
|
||
<span class="n">ypredictSKL</span> <span class="o">=</span> <span class="n">skl</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"MSE with Manual intercept"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">ypredictOwn</span><span class="o">+</span><span class="n">intercept</span><span class="p">))</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">"MSE with Sklearn intercept"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">ypredictSKL</span><span class="p">))</span>
|
||
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta</span> <span class="o">+</span> <span class="n">intercept</span><span class="p">,</span> <span class="s2">"--"</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">"Fit (manual intercept)"</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">skl</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X</span><span class="p">),</span> <span class="s2">"--"</span><span class="p">,</span> <span class="n">label</span><span class="o">=</span><span class="s2">"Sklearn (fit_intercept=True)"</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">grid</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
|
||
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>True beta: [2, 0.5, 3.7]
|
||
Fitted beta: [2.08376632 0.19569961 3.97898392]
|
||
Sklearn fitted beta: [2.08376632 0.19569961 3.97898392]
|
||
MSE with intercept column
|
||
0.00411363461744314
|
||
MSE with intercept column from SKL
|
||
0.004113634617443147
|
||
Manual intercept: 2.083766322923899
|
||
Fitted beta (without intercept): [0.19569961 3.97898392]
|
||
Sklearn intercept: 2.0837663229239043
|
||
Sklearn fitted beta (without intercept): [0.19569961 3.97898392]
|
||
MSE with Manual intercept
|
||
0.00411363461744314
|
||
MSE with Sklearn intercept
|
||
0.004113634617443131
|
||
</pre></div>
|
||
</div>
|
||
<img alt="_images/week36_37_1.png" src="_images/week36_37_1.png" />
|
||
</div>
|
||
</div>
|
||
<p>The intercept is the value of our output/target variable
|
||
when all our features are zero and our function crosses the <span class="math notranslate nohighlight">\(y\)</span>-axis (for a one-dimensional case).</p>
|
||
<p>Printing the MSE, we see first that both methods give the same MSE, as
|
||
they should. However, when we move to for example Ridge regression,
|
||
the way we treat the intercept may give a larger or smaller MSE,
|
||
meaning that the MSE can be penalized by the value of the
|
||
intercept. Not including the intercept in the fit, means that the
|
||
regularization term does not include <span class="math notranslate nohighlight">\(\beta_0\)</span>. For different values
|
||
of <span class="math notranslate nohighlight">\(\lambda\)</span>, this may lead to different MSE values.</p>
|
||
<p>To remind the reader, the regularization term, with the intercept in Ridge regression, is given by</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\lambda \vert\vert \boldsymbol{\beta} \vert\vert_2^2 = \lambda \sum_{j=0}^{p-1}\beta_j^2,
|
||
\]</div>
|
||
<p>but when we take out the intercept, this equation becomes</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\lambda \vert\vert \boldsymbol{\beta} \vert\vert_2^2 = \lambda \sum_{j=1}^{p-1}\beta_j^2.
|
||
\]</div>
|
||
<p>For Lasso regression we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\lambda \vert\vert \boldsymbol{\beta} \vert\vert_1 = \lambda \sum_{j=1}^{p-1}\vert\beta_j\vert.
|
||
\]</div>
|
||
<p>It means that, when scaling the design matrix and the outputs/targets,
|
||
by subtracting the mean values, we have an optimization problem which
|
||
is not penalized by the intercept. The MSE value can then be smaller
|
||
since it focuses only on the remaining quantities. If we however bring
|
||
back the intercept, we will get a MSE which then contains the
|
||
intercept.</p>
|
||
<p>Armed with this wisdom, we attempt first to simply set the intercept equal to <strong>False</strong> in our implementation of Ridge regression for our well-known vanilla data set.</p>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">train_test_split</span>
|
||
<span class="kn">from</span> <span class="nn">sklearn</span> <span class="kn">import</span> <span class="n">linear_model</span>
|
||
|
||
<span class="k">def</span> <span class="nf">MSE</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span><span class="n">y_model</span><span class="p">):</span>
|
||
<span class="n">n</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">size</span><span class="p">(</span><span class="n">y_model</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span><span class="o">-</span><span class="n">y_model</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">/</span><span class="n">n</span>
|
||
|
||
|
||
<span class="c1"># A seed just to ensure that the random numbers are the same for every run.</span>
|
||
<span class="c1"># Useful for eventual debugging.</span>
|
||
<span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">seed</span><span class="p">(</span><span class="mi">3155</span><span class="p">)</span>
|
||
|
||
<span class="n">n</span> <span class="o">=</span> <span class="mi">100</span>
|
||
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">rand</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">+</span> <span class="mf">1.5</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="mi">2</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
|
||
|
||
<span class="n">Maxpolydegree</span> <span class="o">=</span> <span class="mi">20</span>
|
||
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="n">Maxpolydegree</span><span class="p">))</span>
|
||
<span class="c1">#We include explicitely the intercept column</span>
|
||
<span class="k">for</span> <span class="n">degree</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">Maxpolydegree</span><span class="p">):</span>
|
||
<span class="n">X</span><span class="p">[:,</span><span class="n">degree</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="n">degree</span>
|
||
<span class="c1"># We split the data in test and training data</span>
|
||
<span class="n">X_train</span><span class="p">,</span> <span class="n">X_test</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span> <span class="n">y_test</span> <span class="o">=</span> <span class="n">train_test_split</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">test_size</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
|
||
|
||
<span class="n">p</span> <span class="o">=</span> <span class="n">Maxpolydegree</span>
|
||
<span class="n">I</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="n">p</span><span class="p">,</span><span class="n">p</span><span class="p">)</span>
|
||
<span class="c1"># Decide which values of lambda to use</span>
|
||
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">6</span>
|
||
<span class="n">MSEOwnRidgePredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">MSERidgePredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
|
||
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
|
||
<span class="n">OwnRidgeBeta</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">pinv</span><span class="p">(</span><span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">+</span><span class="n">lmb</span><span class="o">*</span><span class="n">I</span><span class="p">)</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y_train</span>
|
||
<span class="c1"># Note: we include the intercept column and no scaling</span>
|
||
<span class="n">RegRidge</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Ridge</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
|
||
<span class="n">RegRidge</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span><span class="n">y_train</span><span class="p">)</span>
|
||
<span class="c1"># and then make the prediction</span>
|
||
<span class="n">ytildeOwnRidge</span> <span class="o">=</span> <span class="n">X_train</span> <span class="o">@</span> <span class="n">OwnRidgeBeta</span>
|
||
<span class="n">ypredictOwnRidge</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">@</span> <span class="n">OwnRidgeBeta</span>
|
||
<span class="n">ytildeRidge</span> <span class="o">=</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_train</span><span class="p">)</span>
|
||
<span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test</span><span class="p">)</span>
|
||
<span class="n">MSEOwnRidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictOwnRidge</span><span class="p">)</span>
|
||
<span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictRidge</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Beta values for own Ridge implementation"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">OwnRidgeBeta</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Beta values for Scikit-Learn Ridge implementation"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">RegRidge</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"MSE values for own Ridge implementation"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSEOwnRidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"MSE values for Scikit-Learn Ridge implementation"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
|
||
|
||
<span class="c1"># Now plot the results</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSEOwnRidgePredict</span><span class="p">,</span> <span class="s1">'r'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE own Ridge Test'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSERidgePredict</span><span class="p">,</span> <span class="s1">'g'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Ridge Test'</span><span class="p">)</span>
|
||
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'log10(lambda)'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'MSE'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Beta values for own Ridge implementation
|
||
[ 1.03032441e+00 6.28336218e-02 -6.24175744e-01 5.21169159e-02
|
||
2.80847477e-01 2.12552073e-01 8.13220608e-02 -1.69634577e-02
|
||
-6.50846111e-02 -7.38962192e-02 -5.94226022e-02 -3.50227564e-02
|
||
-9.80609616e-03 1.08299273e-02 2.41882037e-02 2.93492130e-02
|
||
2.64742912e-02 1.63249532e-02 -5.01831050e-05 -2.15098090e-02]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[ 1.03032441e+00 6.28336218e-02 -6.24175744e-01 5.21169159e-02
|
||
2.80847477e-01 2.12552073e-01 8.13220608e-02 -1.69634577e-02
|
||
-6.50846112e-02 -7.38962192e-02 -5.94226022e-02 -3.50227564e-02
|
||
-9.80609615e-03 1.08299273e-02 2.41882037e-02 2.93492130e-02
|
||
2.64742912e-02 1.63249532e-02 -5.01831190e-05 -2.15098090e-02]
|
||
MSE values for own Ridge implementation
|
||
4.3632959111950474e-07
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
4.363295916366933e-07
|
||
Beta values for own Ridge implementation
|
||
[ 1.03630548 -0.01963611 -0.37900111 -0.07062318 0.12182967 0.16343471
|
||
0.13003291 0.07490892 0.02365049 -0.01449782 -0.03814292 -0.04909093
|
||
-0.05009826 -0.04389027 -0.03279636 -0.01866537 -0.00289724 0.01348565
|
||
0.02976145 0.04543942]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[ 1.03630548 -0.01963611 -0.37900111 -0.07062318 0.12182967 0.16343471
|
||
0.13003291 0.07490892 0.02365049 -0.01449782 -0.03814292 -0.04909093
|
||
-0.05009826 -0.04389027 -0.03279636 -0.01866537 -0.00289724 0.01348565
|
||
0.02976145 0.04543942]
|
||
MSE values for own Ridge implementation
|
||
5.194042826649355e-06
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
5.194042826815211e-06
|
||
Beta values for own Ridge implementation
|
||
[ 1.04220758 -0.10931453 -0.17641709 -0.06020587 0.02208512 0.05789007
|
||
0.06491736 0.05785343 0.04537385 0.03196357 0.01969145 0.00934499
|
||
0.00107405 -0.00526348 -0.00992331 -0.01318643 -0.01531845 -0.01655318
|
||
-0.01708852 -0.01708781]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[ 1.04220758 -0.10931453 -0.17641709 -0.06020587 0.02208512 0.05789007
|
||
0.06491736 0.05785343 0.04537385 0.03196357 0.01969145 0.00934499
|
||
0.00107405 -0.00526348 -0.00992331 -0.01318643 -0.01531845 -0.01655318
|
||
-0.01708852 -0.01708781]
|
||
MSE values for own Ridge implementation
|
||
2.0940821989652176e-05
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
2.0940821989627646e-05
|
||
Beta values for own Ridge implementation
|
||
[ 1.01219292 -0.06043581 -0.10391807 -0.05651951 -0.01898855 0.00312361
|
||
0.01463049 0.01975848 0.02123176 0.02068067 0.01905883 0.01691985
|
||
0.01458337 0.01223198 0.00996754 0.00784393 0.00588657 0.00410387
|
||
0.00249435 0.00105081]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[ 1.01219292 -0.06043581 -0.10391807 -0.05651951 -0.01898855 0.00312361
|
||
0.01463049 0.01975848 0.02123176 0.02068067 0.01905883 0.01691985
|
||
0.01458337 0.01223198 0.00996754 0.00784393 0.00588657 0.00410387
|
||
0.00249435 0.00105081]
|
||
MSE values for own Ridge implementation
|
||
0.00031535148309577417
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
0.0003153514830958095
|
||
Beta values for own Ridge implementation
|
||
[ 8.38916861e-01 1.31276579e-01 8.97497404e-03 -1.72271878e-02
|
||
-2.11744554e-02 -1.91492986e-02 -1.57201944e-02 -1.23002365e-02
|
||
-9.30466214e-03 -6.81048318e-03 -4.78184120e-03 -3.15130074e-03
|
||
-1.84923989e-03 -8.13661243e-04 7.46984697e-06 6.56636616e-04
|
||
1.16805821e-03 1.56912044e-03 1.88168312e-03 2.12318726e-03]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[ 8.38916861e-01 1.31276579e-01 8.97497404e-03 -1.72271878e-02
|
||
-2.11744554e-02 -1.91492986e-02 -1.57201944e-02 -1.23002365e-02
|
||
-9.30466214e-03 -6.81048318e-03 -4.78184120e-03 -3.15130074e-03
|
||
-1.84923989e-03 -8.13661243e-04 7.46984697e-06 6.56636616e-04
|
||
1.16805821e-03 1.56912044e-03 1.88168312e-03 2.12318726e-03]
|
||
MSE values for own Ridge implementation
|
||
0.01507238889517717
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
0.01507238889517706
|
||
Beta values for own Ridge implementation
|
||
[0.37396662 0.14174745 0.0764924 0.04892055 0.03447512 0.02586427
|
||
0.02024962 0.01633913 0.01347916 0.0113104 0.0096208 0.00827728
|
||
0.00719176 0.00630331 0.00556826 0.0049544 0.00443743 0.0039987
|
||
0.0036237 0.003301 ]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[0.37396662 0.14174745 0.0764924 0.04892055 0.03447512 0.02586427
|
||
0.02024962 0.01633913 0.01347916 0.0113104 0.0096208 0.00827728
|
||
0.00719176 0.00630331 0.00556826 0.0049544 0.00443743 0.0039987
|
||
0.0036237 0.003301 ]
|
||
MSE values for own Ridge implementation
|
||
0.2640931530791004
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
0.26409315307910025
|
||
</pre></div>
|
||
</div>
|
||
<img alt="_images/week36_45_1.png" src="_images/week36_45_1.png" />
|
||
</div>
|
||
</div>
|
||
<p>The results here agree when we force <strong>Scikit-Learn</strong>’s Ridge function to include the first column in our design matrix.
|
||
We see that the results agree very well. Here we have thus explicitely included the intercept column in the design matrix.
|
||
What happens if we do not include the intercept in our fit?
|
||
Let us see how we can change this code by zero centering.</p>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">train_test_split</span>
|
||
<span class="kn">from</span> <span class="nn">sklearn</span> <span class="kn">import</span> <span class="n">linear_model</span>
|
||
<span class="kn">from</span> <span class="nn">sklearn.preprocessing</span> <span class="kn">import</span> <span class="n">StandardScaler</span>
|
||
|
||
<span class="k">def</span> <span class="nf">MSE</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span><span class="n">y_model</span><span class="p">):</span>
|
||
<span class="n">n</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">size</span><span class="p">(</span><span class="n">y_model</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span><span class="o">-</span><span class="n">y_model</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">/</span><span class="n">n</span>
|
||
<span class="c1"># A seed just to ensure that the random numbers are the same for every run.</span>
|
||
<span class="c1"># Useful for eventual debugging.</span>
|
||
<span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">seed</span><span class="p">(</span><span class="mi">315</span><span class="p">)</span>
|
||
|
||
<span class="n">n</span> <span class="o">=</span> <span class="mi">100</span>
|
||
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">rand</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">+</span> <span class="mf">1.5</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="mi">2</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
|
||
|
||
<span class="n">Maxpolydegree</span> <span class="o">=</span> <span class="mi">20</span>
|
||
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="n">Maxpolydegree</span><span class="o">-</span><span class="mi">1</span><span class="p">))</span>
|
||
|
||
<span class="k">for</span> <span class="n">degree</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="n">Maxpolydegree</span><span class="p">):</span> <span class="c1">#No intercept column</span>
|
||
<span class="n">X</span><span class="p">[:,</span><span class="n">degree</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="p">(</span><span class="n">degree</span><span class="p">)</span>
|
||
|
||
<span class="c1"># We split the data in test and training data</span>
|
||
<span class="n">X_train</span><span class="p">,</span> <span class="n">X_test</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span> <span class="n">y_test</span> <span class="o">=</span> <span class="n">train_test_split</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">test_size</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
|
||
|
||
<span class="c1">#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable</span>
|
||
<span class="n">X_train_mean</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
|
||
<span class="c1">#Center by removing mean from each feature</span>
|
||
<span class="n">X_train_scaled</span> <span class="o">=</span> <span class="n">X_train</span> <span class="o">-</span> <span class="n">X_train_mean</span>
|
||
<span class="n">X_test_scaled</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">-</span> <span class="n">X_train_mean</span>
|
||
<span class="c1">#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered)</span>
|
||
<span class="c1">#Remove the intercept from the training data.</span>
|
||
<span class="n">y_scaler</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_train</span><span class="p">)</span>
|
||
<span class="n">y_train_scaled</span> <span class="o">=</span> <span class="n">y_train</span> <span class="o">-</span> <span class="n">y_scaler</span>
|
||
|
||
<span class="n">p</span> <span class="o">=</span> <span class="n">Maxpolydegree</span><span class="o">-</span><span class="mi">1</span>
|
||
<span class="n">I</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="n">p</span><span class="p">,</span><span class="n">p</span><span class="p">)</span>
|
||
<span class="c1"># Decide which values of lambda to use</span>
|
||
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">6</span>
|
||
<span class="n">MSEOwnRidgePredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">MSERidgePredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
|
||
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
|
||
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
|
||
<span class="n">OwnRidgeBeta</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">pinv</span><span class="p">(</span><span class="n">X_train_scaled</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X_train_scaled</span><span class="o">+</span><span class="n">lmb</span><span class="o">*</span><span class="n">I</span><span class="p">)</span> <span class="o">@</span> <span class="n">X_train_scaled</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="p">(</span><span class="n">y_train_scaled</span><span class="p">)</span>
|
||
<span class="n">intercept_</span> <span class="o">=</span> <span class="n">y_scaler</span> <span class="o">-</span> <span class="n">X_train_mean</span><span class="nd">@OwnRidgeBeta</span> <span class="c1">#The intercept can be shifted so the model can predict on uncentered data</span>
|
||
<span class="c1">#Add intercept to prediction</span>
|
||
<span class="n">ypredictOwnRidge</span> <span class="o">=</span> <span class="n">X_test_scaled</span> <span class="o">@</span> <span class="n">OwnRidgeBeta</span> <span class="o">+</span> <span class="n">y_scaler</span>
|
||
<span class="n">RegRidge</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Ridge</span><span class="p">(</span><span class="n">lmb</span><span class="p">)</span>
|
||
<span class="n">RegRidge</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span><span class="n">y_train</span><span class="p">)</span>
|
||
<span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test</span><span class="p">)</span>
|
||
<span class="n">MSEOwnRidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictOwnRidge</span><span class="p">)</span>
|
||
<span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictRidge</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Beta values for own Ridge implementation"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">OwnRidgeBeta</span><span class="p">)</span> <span class="c1">#Intercept is given by mean of target variable</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Beta values for Scikit-Learn Ridge implementation"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">RegRidge</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s1">'Intercept from own implementation:'</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">intercept_</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s1">'Intercept from Scikit-Learn Ridge implementation'</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">RegRidge</span><span class="o">.</span><span class="n">intercept_</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"MSE values for own Ridge implementation"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSEOwnRidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"MSE values for Scikit-Learn Ridge implementation"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
|
||
|
||
|
||
<span class="c1"># Now plot the results</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSEOwnRidgePredict</span><span class="p">,</span> <span class="s1">'b--'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE own Ridge Test'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSERidgePredict</span><span class="p">,</span> <span class="s1">'g--'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE SL Ridge Test'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'log10(lambda)'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'MSE'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Beta values for own Ridge implementation
|
||
[ 3.43579948e-02 -5.43330971e-01 -3.10141413e-03 2.47116868e-01
|
||
2.18613217e-01 1.02054837e-01 -4.25617658e-04 -5.90475506e-02
|
||
-7.68534263e-02 -6.68929213e-02 -4.24906604e-02 -1.40927184e-02
|
||
1.11482289e-02 2.88529063e-02 3.67047975e-02 3.38135733e-02
|
||
2.02198703e-02 -3.46383926e-03 -3.63025821e-02]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[ 3.43579948e-02 -5.43330971e-01 -3.10141413e-03 2.47116868e-01
|
||
2.18613217e-01 1.02054837e-01 -4.25617655e-04 -5.90475506e-02
|
||
-7.68534263e-02 -6.68929213e-02 -4.24906604e-02 -1.40927184e-02
|
||
1.11482289e-02 2.88529063e-02 3.67047975e-02 3.38135733e-02
|
||
2.02198702e-02 -3.46383925e-03 -3.63025821e-02]
|
||
Intercept from own implementation:
|
||
1.0330308045187757
|
||
Intercept from Scikit-Learn Ridge implementation
|
||
1.0330308045183194
|
||
MSE values for own Ridge implementation
|
||
3.139255958997547e-06
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
3.1392559585020426e-06
|
||
Beta values for own Ridge implementation
|
||
[-0.05807125 -0.29822833 -0.08551306 0.08156108 0.13679863 0.12333649
|
||
0.08251519 0.03815288 0.00111756 -0.02498832 -0.04010697 -0.04566964
|
||
-0.04355837 -0.03562355 -0.02348765 -0.00848904 0.00831018 0.0260906
|
||
0.04423486]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[-0.05807125 -0.29822833 -0.08551306 0.08156108 0.13679863 0.12333649
|
||
0.08251519 0.03815288 0.00111756 -0.02498832 -0.04010697 -0.04566964
|
||
-0.04355837 -0.03562355 -0.02348765 -0.00848904 0.00831018 0.0260906
|
||
0.04423486]
|
||
Intercept from own implementation:
|
||
1.0411487294305088
|
||
Intercept from Scikit-Learn Ridge implementation
|
||
1.0411487294305226
|
||
MSE values for own Ridge implementation
|
||
1.9601304850035702e-05
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
1.9601304850073734e-05
|
||
Beta values for own Ridge implementation
|
||
[-0.1416398 -0.14021063 -0.05383795 0.01367553 0.04784395 0.05796251
|
||
0.05447415 0.044613 0.03267527 0.02098261 0.01066519 0.00217499
|
||
-0.00440346 -0.00917248 -0.01231917 -0.01405935 -0.0146081 -0.01416528
|
||
-0.01290947]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[-0.1416398 -0.14021063 -0.05383795 0.01367553 0.04784395 0.05796251
|
||
0.05447415 0.044613 0.03267527 0.02098261 0.01066519 0.00217499
|
||
-0.00440346 -0.00917248 -0.01231917 -0.01405935 -0.0146081 -0.01416528
|
||
-0.01290947]
|
||
Intercept from own implementation:
|
||
1.049556996627824
|
||
Intercept from Scikit-Learn Ridge implementation
|
||
1.0495569966278269
|
||
MSE values for own Ridge implementation
|
||
5.4959161509357395e-05
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
5.4959161509366685e-05
|
||
Beta values for own Ridge implementation
|
||
[-0.13535942 -0.08593216 -0.03568439 -0.0036367 0.01397146 0.02229529
|
||
0.02503753 0.0245528 0.02228115 0.01908936 0.01549377 0.01179792
|
||
0.00817631 0.00472512 0.00149311 -0.00149956 -0.00424967 -0.00676387
|
||
-0.00905423]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[-0.13535942 -0.08593216 -0.03568439 -0.0036367 0.01397146 0.02229529
|
||
0.02503753 0.0245528 0.02228115 0.01908936 0.01549377 0.01179792
|
||
0.00817631 0.00472512 0.00149311 -0.00149956 -0.00424967 -0.00676387
|
||
-0.00905423]
|
||
Intercept from own implementation:
|
||
1.039967668952797
|
||
Intercept from Scikit-Learn Ridge implementation
|
||
1.0399676689527975
|
||
MSE values for own Ridge implementation
|
||
7.571105947979344e-05
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
7.571105947979412e-05
|
||
Beta values for own Ridge implementation
|
||
[-0.05100875 -0.04063602 -0.02723445 -0.01713366 -0.0100706 -0.00517114
|
||
-0.00174276 0.00068734 0.00243186 0.00369758 0.00462287 0.0053018
|
||
0.00579953 0.006162 0.00642221 0.00660427 0.00672607 0.0068011
|
||
0.00683964]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[-0.05100875 -0.04063602 -0.02723445 -0.01713366 -0.0100706 -0.00517114
|
||
-0.00174276 0.00068734 0.00243186 0.00369758 0.00462287 0.0053018
|
||
0.00579953 0.006162 0.00642221 0.00660427 0.00672607 0.0068011
|
||
0.00683964]
|
||
Intercept from own implementation:
|
||
0.999955585168597
|
||
Intercept from Scikit-Learn Ridge implementation
|
||
0.999955585168597
|
||
MSE values for own Ridge implementation
|
||
0.0007698473260556343
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
0.0007698473260556325
|
||
Beta values for own Ridge implementation
|
||
[-0.00834567 -0.00803064 -0.00673407 -0.00554552 -0.00458878 -0.0038335
|
||
-0.00323332 -0.00274989 -0.0023548 -0.00202756 -0.00175331 -0.00152117
|
||
-0.001323 -0.0011526 -0.00100519 -0.00087697 -0.00076495 -0.00066668
|
||
-0.00058016]
|
||
Beta values for Scikit-Learn Ridge implementation
|
||
[-0.00834567 -0.00803064 -0.00673407 -0.00554552 -0.00458878 -0.0038335
|
||
-0.00323332 -0.00274989 -0.0023548 -0.00202756 -0.00175331 -0.00152117
|
||
-0.001323 -0.0011526 -0.00100519 -0.00087697 -0.00076495 -0.00066668
|
||
-0.00058016]
|
||
Intercept from own implementation:
|
||
0.9637117593816477
|
||
Intercept from Scikit-Learn Ridge implementation
|
||
0.9637117593816477
|
||
MSE values for own Ridge implementation
|
||
0.0023813163025848865
|
||
MSE values for Scikit-Learn Ridge implementation
|
||
0.002381316302584885
|
||
</pre></div>
|
||
</div>
|
||
<img alt="_images/week36_47_1.png" src="_images/week36_47_1.png" />
|
||
</div>
|
||
</div>
|
||
<p>We see here, when compared to the code which includes explicitely the
|
||
intercept column, that our MSE value is actually smaller. This is
|
||
because the regularization term does not include the intercept value
|
||
<span class="math notranslate nohighlight">\(\beta_0\)</span> in the fitting. This applies to Lasso regularization as
|
||
well. It means that our optimization is now done only with the
|
||
centered matrix and/or vector that enter the fitting procedure.</p>
|
||
</div>
|
||
<div class="section" id="test-function-for-what-happens-with-ols-ridge-and-lasso">
|
||
<h2>Test Function for what happens with OLS, Ridge and Lasso<a class="headerlink" href="#test-function-for-what-happens-with-ols-ridge-and-lasso" title="Permalink to this headline">¶</a></h2>
|
||
<p>Hitherto we have discussed Ridge and Lasso regression in terms of a
|
||
linear analysis. This may to many of you feel rather technical and
|
||
perhaps not that intuitive. The question is whether we can develop a
|
||
more intuitive way of understanding what Ridge and Lasso express.</p>
|
||
<p>Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.</p>
|
||
<p>We will play around with a study of the values for the optimal
|
||
parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> using OLS, Ridge and Lasso regression. For
|
||
OLS, you will notice as function of the noise and polynomial degree,
|
||
that the parameters <span class="math notranslate nohighlight">\(\beta\)</span> will fluctuate from order to order in the
|
||
polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS.</p>
|
||
<p>For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.</p>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">train_test_split</span>
|
||
<span class="kn">from</span> <span class="nn">sklearn</span> <span class="kn">import</span> <span class="n">linear_model</span>
|
||
|
||
<span class="k">def</span> <span class="nf">R2</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span> <span class="n">y_model</span><span class="p">):</span>
|
||
<span class="k">return</span> <span class="mi">1</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span> <span class="o">-</span> <span class="n">y_model</span><span class="p">)</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_data</span><span class="p">))</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span>
|
||
<span class="k">def</span> <span class="nf">MSE</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span><span class="n">y_model</span><span class="p">):</span>
|
||
<span class="n">n</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">size</span><span class="p">(</span><span class="n">y_model</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span><span class="o">-</span><span class="n">y_model</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">/</span><span class="n">n</span>
|
||
|
||
<span class="c1"># Make data set.</span>
|
||
<span class="n">n</span> <span class="o">=</span> <span class="mi">10000</span>
|
||
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">rand</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">+</span> <span class="mf">1.5</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="mi">2</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">+</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">randn</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
|
||
|
||
<span class="n">Maxpolydegree</span> <span class="o">=</span> <span class="mi">5</span>
|
||
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">),</span><span class="n">Maxpolydegree</span><span class="p">))</span>
|
||
<span class="n">X</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>
|
||
|
||
|
||
<span class="k">for</span> <span class="n">polydegree</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="n">Maxpolydegree</span><span class="p">):</span>
|
||
<span class="n">X</span><span class="p">[:,</span><span class="n">polydegree</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="p">(</span><span class="n">polydegree</span><span class="p">)</span>
|
||
|
||
<span class="c1"># We split the data in test and training data</span>
|
||
<span class="n">X_train</span><span class="p">,</span> <span class="n">X_test</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span> <span class="n">y_test</span> <span class="o">=</span> <span class="n">train_test_split</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">test_size</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
|
||
|
||
<span class="c1"># matrix inversion to find beta</span>
|
||
<span class="n">OLSbeta</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">pinv</span><span class="p">(</span><span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X_train</span><span class="p">)</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y_train</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">OLSbeta</span><span class="p">)</span>
|
||
<span class="n">ypredictOLS</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">@</span> <span class="n">OLSbeta</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Test MSE OLS"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictOLS</span><span class="p">))</span>
|
||
<span class="c1"># Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn</span>
|
||
<span class="c1"># Decide which values of lambda to use</span>
|
||
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">4</span>
|
||
<span class="n">MSERidgePredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">MSELassoPredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
|
||
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
|
||
<span class="c1"># Make the fit using Ridge and Lasso</span>
|
||
<span class="n">RegRidge</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Ridge</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
|
||
<span class="n">RegRidge</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span><span class="n">y_train</span><span class="p">)</span>
|
||
<span class="n">RegLasso</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Lasso</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
|
||
<span class="n">RegLasso</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span><span class="n">y_train</span><span class="p">)</span>
|
||
<span class="c1"># and then make the prediction</span>
|
||
<span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test</span><span class="p">)</span>
|
||
<span class="n">ypredictLasso</span> <span class="o">=</span> <span class="n">RegLasso</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test</span><span class="p">)</span>
|
||
<span class="c1"># Compute the MSE and print it</span>
|
||
<span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictRidge</span><span class="p">)</span>
|
||
<span class="n">MSELassoPredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictLasso</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">RegRidge</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">RegLasso</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||
<span class="c1"># Now plot the results</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSERidgePredict</span><span class="p">,</span> <span class="s1">'b'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Ridge Test'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSELassoPredict</span><span class="p">,</span> <span class="s1">'r'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Lasso Test'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'log10(lambda)'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'MSE'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[ 0.99762142 0.46713267 -2.86273371 4.1801709 -1.87068171]
|
||
Test MSE OLS
|
||
0.9990927157696581
|
||
0.001 [ 0.99884986 0.44395301 -2.76190591 4.02693858 -1.7953443 ]
|
||
0.001 [ 1.00497837 -0.11830757 -0. 0. 0.0064472 ]
|
||
0.021544346900318832 [ 1.01267007 0.18307821 -1.62725226 2.30283188 -0.94781064]
|
||
0.021544346900318832 [ 0.92756256 -0. -0. -0. -0. ]
|
||
0.46415888336127775 [ 1.02895918 -0.12952339 -0.27065633 0.25185716 0.05519999]
|
||
0.46415888336127775 [0.48494803 0. 0. 0. 0. ]
|
||
10.0 [ 1.02077431 -0.14568914 -0.10041812 0.03846169 0.11382896]
|
||
10.0 [0. 0. 0. 0. 0.]
|
||
</pre></div>
|
||
</div>
|
||
<img alt="_images/week36_50_1.png" src="_images/week36_50_1.png" />
|
||
</div>
|
||
</div>
|
||
<p>How can we understand this?</p>
|
||
</div>
|
||
<div class="section" id="linking-the-regression-analysis-with-a-statistical-interpretation">
|
||
<h2>Linking the regression analysis with a statistical interpretation<a class="headerlink" href="#linking-the-regression-analysis-with-a-statistical-interpretation" title="Permalink to this headline">¶</a></h2>
|
||
<p>We will now couple the discussions of ordinary least squares, Ridge
|
||
and Lasso regression with a statistical interpretation, that is we
|
||
move from a linear algebra analysis to a statistical analysis. In
|
||
particular, we will focus on what the regularization terms can result
|
||
in. We will amongst other things show that the regularization
|
||
parameter can reduce considerably the variance of the parameters
|
||
<span class="math notranslate nohighlight">\(\beta\)</span>.</p>
|
||
<p>The
|
||
advantage of doing linear regression is that we actually end up with
|
||
analytical expressions for several statistical quantities.<br />
|
||
Standard least squares and Ridge regression allow us to
|
||
derive quantities like the variance and other expectation values in a
|
||
rather straightforward way.</p>
|
||
<p>It is assumed that <span class="math notranslate nohighlight">\(\varepsilon_i
|
||
\sim \mathcal{N}(0, \sigma^2)\)</span> and the <span class="math notranslate nohighlight">\(\varepsilon_{i}\)</span> are
|
||
independent, i.e.:</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\begin{align*}
|
||
\mbox{Cov}(\varepsilon_{i_1},
|
||
\varepsilon_{i_2}) & = \left\{ \begin{array}{lcc} \sigma^2 & \mbox{if}
|
||
& i_1 = i_2, \\ 0 & \mbox{if} & i_1 \not= i_2. \end{array} \right.
|
||
\end{align*}
|
||
\end{split}\]</div>
|
||
<p>The randomness of <span class="math notranslate nohighlight">\(\varepsilon_i\)</span> implies that
|
||
<span class="math notranslate nohighlight">\(\mathbf{y}_i\)</span> is also a random variable. In particular,
|
||
<span class="math notranslate nohighlight">\(\mathbf{y}_i\)</span> is normally distributed, because <span class="math notranslate nohighlight">\(\varepsilon_i \sim
|
||
\mathcal{N}(0, \sigma^2)\)</span> and <span class="math notranslate nohighlight">\(\mathbf{X}_{i,\ast} \, \boldsymbol{\beta}\)</span> is a
|
||
non-random scalar. To specify the parameters of the distribution of
|
||
<span class="math notranslate nohighlight">\(\mathbf{y}_i\)</span> we need to calculate its first two moments.</p>
|
||
<p>Recall that <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is a matrix of dimensionality <span class="math notranslate nohighlight">\(n\times p\)</span>. The
|
||
notation above <span class="math notranslate nohighlight">\(\mathbf{X}_{i,\ast}\)</span> means that we are looking at the
|
||
row number <span class="math notranslate nohighlight">\(i\)</span> and perform a sum over all values <span class="math notranslate nohighlight">\(p\)</span>.</p>
|
||
</div>
|
||
<div class="section" id="assumptions-made">
|
||
<h2>Assumptions made<a class="headerlink" href="#assumptions-made" title="Permalink to this headline">¶</a></h2>
|
||
<p>The assumption we have made here can be summarized as (and this is going to be useful when we discuss the bias-variance trade off)
|
||
that there exists a function <span class="math notranslate nohighlight">\(f(\boldsymbol{x})\)</span> and a normal distributed error <span class="math notranslate nohighlight">\(\boldsymbol{\varepsilon}\sim \mathcal{N}(0, \sigma^2)\)</span>
|
||
which describe our data</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{y} = f(\boldsymbol{x})+\boldsymbol{\varepsilon}
|
||
\]</div>
|
||
<p>We approximate this function with our model from the solution of the linear regression equations, that is our
|
||
function <span class="math notranslate nohighlight">\(f\)</span> is approximated by <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}}\)</span> where we want to minimize <span class="math notranslate nohighlight">\((\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\)</span>, our MSE, with</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{\tilde{y}} = \boldsymbol{X}\boldsymbol{\beta}.
|
||
\]</div>
|
||
</div>
|
||
<div class="section" id="expectation-value-and-variance">
|
||
<h2>Expectation value and variance<a class="headerlink" href="#expectation-value-and-variance" title="Permalink to this headline">¶</a></h2>
|
||
<p>We can calculate the expectation value of <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> for a given element <span class="math notranslate nohighlight">\(i\)</span></p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\begin{align*}
|
||
\mathbb{E}(y_i) & =
|
||
\mathbb{E}(\mathbf{X}_{i, \ast} \, \boldsymbol{\beta}) + \mathbb{E}(\varepsilon_i)
|
||
\, \, \, = \, \, \, \mathbf{X}_{i, \ast} \, \beta,
|
||
\end{align*}
|
||
\]</div>
|
||
<p>while
|
||
its variance is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\begin{align*} \mbox{Var}(y_i) & = \mathbb{E} \{ [y_i
|
||
- \mathbb{E}(y_i)]^2 \} \, \, \, = \, \, \, \mathbb{E} ( y_i^2 ) -
|
||
[\mathbb{E}(y_i)]^2 \\ & = \mathbb{E} [ ( \mathbf{X}_{i, \ast} \,
|
||
\beta + \varepsilon_i )^2] - ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 \\ &
|
||
= \mathbb{E} [ ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 + 2 \varepsilon_i
|
||
\mathbf{X}_{i, \ast} \, \boldsymbol{\beta} + \varepsilon_i^2 ] - ( \mathbf{X}_{i,
|
||
\ast} \, \beta)^2 \\ & = ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 + 2
|
||
\mathbb{E}(\varepsilon_i) \mathbf{X}_{i, \ast} \, \boldsymbol{\beta} +
|
||
\mathbb{E}(\varepsilon_i^2 ) - ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2
|
||
\\ & = \mathbb{E}(\varepsilon_i^2 ) \, \, \, = \, \, \,
|
||
\mbox{Var}(\varepsilon_i) \, \, \, = \, \, \, \sigma^2.
|
||
\end{align*}
|
||
\end{split}\]</div>
|
||
<p>Hence, <span class="math notranslate nohighlight">\(y_i \sim \mathcal{N}( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta}, \sigma^2)\)</span>, that is <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> follows a normal distribution with
|
||
mean value <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{\beta}\)</span> and variance <span class="math notranslate nohighlight">\(\sigma^2\)</span> (not be confused with the singular values of the SVD).</p>
|
||
</div>
|
||
<div class="section" id="expectation-value-and-variance-for-boldsymbol-beta">
|
||
<h2>Expectation value and variance for <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span><a class="headerlink" href="#expectation-value-and-variance-for-boldsymbol-beta" title="Permalink to this headline">¶</a></h2>
|
||
<p>With the OLS expressions for the optimal parameters <span class="math notranslate nohighlight">\(\boldsymbol{\hat{\beta}}\)</span> we can evaluate the expectation value</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\mathbb{E}(\boldsymbol{\hat{\beta}}) = \mathbb{E}[ (\mathbf{X}^{\top} \mathbf{X})^{-1}\mathbf{X}^{T} \mathbf{Y}]=(\mathbf{X}^{T} \mathbf{X})^{-1}\mathbf{X}^{T} \mathbb{E}[ \mathbf{Y}]=(\mathbf{X}^{T} \mathbf{X})^{-1} \mathbf{X}^{T}\mathbf{X}\boldsymbol{\beta}=\boldsymbol{\beta}.
|
||
\]</div>
|
||
<p>This means that the estimator of the regression parameters is unbiased.</p>
|
||
<p>We can also calculate the variance</p>
|
||
<p>The variance of the optimal value <span class="math notranslate nohighlight">\(\boldsymbol{\hat{\beta}}\)</span> is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\begin{eqnarray*}
|
||
\mbox{Var}(\boldsymbol{\hat{\beta}}) & = & \mathbb{E} \{ [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})] [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})]^{T} \}
|
||
\\
|
||
& = & \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}]^{T} \}
|
||
\\
|
||
% & = & \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}]^{T} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
|
||
% \\
|
||
% & = & \mathbb{E} \{ (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} \, \mathbf{Y}^{T} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
|
||
% \\
|
||
& = & (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \mathbb{E} \{ \mathbf{Y} \, \mathbf{Y}^{T} \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
|
||
\\
|
||
& = & (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \{ \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} + \sigma^2 \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
|
||
% \\
|
||
% & = & (\mathbf{X}^T \mathbf{X})^{-1} \, \mathbf{X}^T \, \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^T \, \mathbf{X}^T \, \mathbf{X} \, (\mathbf{X}^T % \mathbf{X})^{-1}
|
||
% \\
|
||
% & & + \, \, \sigma^2 \, (\mathbf{X}^T \mathbf{X})^{-1} \, \mathbf{X}^T \, \mathbf{X} \, (\mathbf{X}^T \mathbf{X})^{-1} - \boldsymbol{\beta} \boldsymbol{\beta}^T
|
||
\\
|
||
& = & \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} + \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
|
||
\, \, \, = \, \, \, \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1},
|
||
\end{eqnarray*}
|
||
\end{split}\]</div>
|
||
<p>where we have used that <span class="math notranslate nohighlight">\(\mathbb{E} (\mathbf{Y} \mathbf{Y}^{T}) =
|
||
\mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} +
|
||
\sigma^2 \, \mathbf{I}_{nn}\)</span>. From <span class="math notranslate nohighlight">\(\mbox{Var}(\boldsymbol{\beta}) = \sigma^2
|
||
\, (\mathbf{X}^{T} \mathbf{X})^{-1}\)</span>, one obtains an estimate of the
|
||
variance of the estimate of the <span class="math notranslate nohighlight">\(j\)</span>-th regression coefficient:
|
||
<span class="math notranslate nohighlight">\(\boldsymbol{\sigma}^2 (\boldsymbol{\beta}_j ) = \boldsymbol{\sigma}^2 [(\mathbf{X}^{T} \mathbf{X})^{-1}]_{jj} \)</span>. This may be used to
|
||
construct a confidence interval for the estimates.</p>
|
||
<p>In a similar way, we can obtain analytical expressions for say the
|
||
expectation values of the parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> and their variance
|
||
when we employ Ridge regression, allowing us again to define a confidence interval.</p>
|
||
<p>It is rather straightforward to show that</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\boldsymbol{\beta}^{\mathrm{OLS}}.
|
||
\]</div>
|
||
<p>We see clearly that
|
||
<span class="math notranslate nohighlight">\(\mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big] \not= \boldsymbol{\beta}^{\mathrm{OLS}}\)</span> for any <span class="math notranslate nohighlight">\(\lambda > 0\)</span>. We say then that the ridge estimator is biased.</p>
|
||
<p>We can also compute the variance as</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\mbox{Var}[\boldsymbol{\beta}^{\mathrm{Ridge}}]=\sigma^2[ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1} \mathbf{X}^{T} \mathbf{X} \{ [ \mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T},
|
||
\]</div>
|
||
<p>and it is easy to see that if the parameter <span class="math notranslate nohighlight">\(\lambda\)</span> goes to infinity then the variance of Ridge parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> goes to zero.</p>
|
||
<p>With this, we can compute the difference</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\mbox{Var}[\boldsymbol{\beta}^{\mathrm{OLS}}]-\mbox{Var}(\boldsymbol{\beta}^{\mathrm{Ridge}})=\sigma^2 [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}[ 2\lambda\mathbf{I} + \lambda^2 (\mathbf{X}^{T} \mathbf{X})^{-1} ] \{ [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}.
|
||
\]</div>
|
||
<p>The difference is non-negative definite since each component of the
|
||
matrix product is non-negative definite.
|
||
This means the variance we obtain with the standard OLS will always for <span class="math notranslate nohighlight">\(\lambda > 0\)</span> be larger than the variance of <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> obtained with the Ridge estimator. This has interesting consequences when we discuss the so-called bias-variance trade-off below.</p>
|
||
</div>
|
||
<div class="section" id="deriving-ols-from-a-probability-distribution">
|
||
<h2>Deriving OLS from a probability distribution<a class="headerlink" href="#deriving-ols-from-a-probability-distribution" title="Permalink to this headline">¶</a></h2>
|
||
<p>Our basic assumption when we derived the OLS equations was to assume
|
||
that our output is determined by a given continuous function
|
||
<span class="math notranslate nohighlight">\(f(\boldsymbol{x})\)</span> and a random noise <span class="math notranslate nohighlight">\(\boldsymbol{\epsilon}\)</span> given by the normal
|
||
distribution with zero mean value and an undetermined variance
|
||
<span class="math notranslate nohighlight">\(\sigma^2\)</span>.</p>
|
||
<p>We found above that the outputs <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> have a mean value given by
|
||
<span class="math notranslate nohighlight">\(\boldsymbol{X}\hat{\boldsymbol{\beta}}\)</span> and variance <span class="math notranslate nohighlight">\(\sigma^2\)</span>. Since the entries to
|
||
the design matrix are not stochastic variables, we can assume that the
|
||
probability distribution of our targets is also a normal distribution
|
||
but now with mean value <span class="math notranslate nohighlight">\(\boldsymbol{X}\hat{\boldsymbol{\beta}}\)</span>. This means that a
|
||
single output <span class="math notranslate nohighlight">\(y_i\)</span> is given by the Gaussian distribution</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
y_i\sim \mathcal{N}(\boldsymbol{X}_{i,*}\boldsymbol{\beta}, \sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}.
|
||
\]</div>
|
||
</div>
|
||
<div class="section" id="independent-and-identically-distributed-iid">
|
||
<h2>Independent and Identically Distributed (iid)<a class="headerlink" href="#independent-and-identically-distributed-iid" title="Permalink to this headline">¶</a></h2>
|
||
<p>We assume now that the various <span class="math notranslate nohighlight">\(y_i\)</span> values are stochastically distributed according to the above Gaussian distribution.
|
||
We define this distribution as</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(y_i, \boldsymbol{X}\vert\boldsymbol{\beta})=\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]},
|
||
\]</div>
|
||
<p>which reads as finding the likelihood of an event <span class="math notranslate nohighlight">\(y_i\)</span> with the input variables <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> given the parameters (to be determined) <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>.</p>
|
||
<p>Since these events are assumed to be independent and identicall distributed we can build the probability distribution function (PDF) for all possible event <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> as the product of the single events, that is we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(\boldsymbol{y},\boldsymbol{X}\vert\boldsymbol{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}=\prod_{i=0}^{n-1}p(y_i,\boldsymbol{X}\vert\boldsymbol{\beta}).
|
||
\]</div>
|
||
<p>We will write this in a more compact form reserving <span class="math notranslate nohighlight">\(\boldsymbol{D}\)</span> for the domain of events, including the ouputs (targets) and the inputs. That is
|
||
in case we have a simple one-dimensional input and output case</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\dots, (x_{n-1},y_{n-1})].
|
||
\]</div>
|
||
<p>In the more general case the various inputs should be replaced by the possible features represented by the input data set <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.
|
||
We can now rewrite the above probability as</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(\boldsymbol{D}\vert\boldsymbol{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}.
|
||
\]</div>
|
||
<p>It is a conditional probability (see below) and reads as the likelihood of a domain of events <span class="math notranslate nohighlight">\(\boldsymbol{D}\)</span> given a set of parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>.</p>
|
||
</div>
|
||
<div class="section" id="maximum-likelihood-estimation-mle">
|
||
<h2>Maximum Likelihood Estimation (MLE)<a class="headerlink" href="#maximum-likelihood-estimation-mle" title="Permalink to this headline">¶</a></h2>
|
||
<p>In statistics, maximum likelihood estimation (MLE) is a method of
|
||
estimating the parameters of an assumed probability distribution,
|
||
given some observed data. This is achieved by maximizing a likelihood
|
||
function so that, under the assumed statistical model, the observed
|
||
data is the most probable.</p>
|
||
<p>We will assume here that our events are given by the above Gaussian
|
||
distribution and we will determine the optimal parameters <span class="math notranslate nohighlight">\(\beta\)</span> by
|
||
maximizing the above PDF. However, computing the derivatives of a
|
||
product function is cumbersome and can easily lead to overflow and/or
|
||
underflowproblems, with potentials for loss of numerical precision.</p>
|
||
<p>In practice, it is more convenient to maximize the logarithm of the
|
||
PDF because it is a monotonically increasing function of the argument.
|
||
Alternatively, and this will be our option, we will minimize the
|
||
negative of the logarithm since this is a monotonically decreasing
|
||
function.</p>
|
||
<p>Note also that maximization/minimization of the logarithm of the PDF
|
||
is equivalent to the maximization/minimization of the function itself.</p>
|
||
</div>
|
||
<div class="section" id="a-new-cost-function">
|
||
<h2>A new Cost Function<a class="headerlink" href="#a-new-cost-function" title="Permalink to this headline">¶</a></h2>
|
||
<p>We could now define a new cost function to minimize, namely the negative logarithm of the above PDF</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta}=-\log{\prod_{i=0}^{n-1}p(y_i,\boldsymbol{X}\vert\boldsymbol{\beta})}=-\sum_{i=0}^{n-1}\log{p(y_i,\boldsymbol{X}\vert\boldsymbol{\beta})},
|
||
\]</div>
|
||
<p>which becomes</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta}=\frac{n}{2}\log{2\pi\sigma^2}+\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}.
|
||
\]</div>
|
||
<p>Taking the derivative of the <em>new</em> cost function with respect to the parameters <span class="math notranslate nohighlight">\(\beta\)</span> we recognize our familiar OLS equation, namely</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{X}^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right) =0,
|
||
\]</div>
|
||
<p>which leads to the well-known OLS equation for the optimal paramters <span class="math notranslate nohighlight">\(\beta\)</span></p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\hat{\boldsymbol{\beta}}^{\mathrm{OLS}}=\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}!
|
||
\]</div>
|
||
<p>Before we make a similar analysis for Ridge and Lasso regression, we need a short reminder on statistics.</p>
|
||
</div>
|
||
<div class="section" id="more-basic-statistics-and-bayes-theorem">
|
||
<h2>More basic Statistics and Bayes’ theorem<a class="headerlink" href="#more-basic-statistics-and-bayes-theorem" title="Permalink to this headline">¶</a></h2>
|
||
<p>A central theorem in statistics is Bayes’ theorem. This theorem plays a similar role as the good old Pythagoras’ theorem in geometry.
|
||
Bayes’ theorem is extremely simple to derive. But to do so we need some basic axioms from statistics.</p>
|
||
<p>Assume we have two domains of events <span class="math notranslate nohighlight">\(X=[x_0,x_1,\dots,x_{n-1}]\)</span> and <span class="math notranslate nohighlight">\(Y=[y_0,y_1,\dots,y_{n-1}]\)</span>.</p>
|
||
<p>We define also the likelihood for <span class="math notranslate nohighlight">\(X\)</span> and <span class="math notranslate nohighlight">\(Y\)</span> as <span class="math notranslate nohighlight">\(p(X)\)</span> and <span class="math notranslate nohighlight">\(p(Y)\)</span> respectively.
|
||
The likelihood of a specific event <span class="math notranslate nohighlight">\(x_i\)</span> (or <span class="math notranslate nohighlight">\(y_i\)</span>) is then written as <span class="math notranslate nohighlight">\(p(X=x_i)\)</span> or just <span class="math notranslate nohighlight">\(p(x_i)=p_i\)</span>.</p>
|
||
<p><strong>Union of events is given by.</strong></p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(X \cup Y)= p(X)+p(Y)-p(X \cap Y).
|
||
\]</div>
|
||
<p><strong>The product rule (aka joint probability) is given by.</strong></p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(X \cup Y)= p(X,Y)= p(X\vert Y)p(Y)=p(Y\vert X)p(X),
|
||
\]</div>
|
||
<p>where we read <span class="math notranslate nohighlight">\(p(X\vert Y)\)</span> as the likelihood of obtaining <span class="math notranslate nohighlight">\(X\)</span> given <span class="math notranslate nohighlight">\(Y\)</span>.</p>
|
||
<p>If we have independent events then <span class="math notranslate nohighlight">\(p(X,Y)=p(X)p(Y)\)</span>.</p>
|
||
</div>
|
||
<div class="section" id="marginal-probability">
|
||
<h2>Marginal Probability<a class="headerlink" href="#marginal-probability" title="Permalink to this headline">¶</a></h2>
|
||
<p>The marginal probability is defined in terms of only one of the set of variables <span class="math notranslate nohighlight">\(X,Y\)</span>. For a discrete probability we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(X)=\sum_{i=0}^{n-1}p(X,Y=y_i)=\sum_{i=0}^{n-1}p(X\vert Y=y_i)p(Y=y_i)=\sum_{i=0}^{n-1}p(X\vert y_i)p(y_i).
|
||
\]</div>
|
||
</div>
|
||
<div class="section" id="conditional-probability">
|
||
<h2>Conditional Probability<a class="headerlink" href="#conditional-probability" title="Permalink to this headline">¶</a></h2>
|
||
<p>The conditional probability, if <span class="math notranslate nohighlight">\(p(Y) > 0\)</span>, is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(X\vert Y)= \frac{p(X,Y)}{p(Y)}=\frac{p(X,Y)}{\sum_{i=0}^{n-1}p(Y\vert X=x_i)p(x_i)}.
|
||
\]</div>
|
||
</div>
|
||
<div class="section" id="bayes-theorem">
|
||
<h2>Bayes’ Theorem<a class="headerlink" href="#bayes-theorem" title="Permalink to this headline">¶</a></h2>
|
||
<p>If we combine the conditional probability with the marginal probability and the standard product rule, we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(X\vert Y)= \frac{p(X,Y)}{p(Y)},
|
||
\]</div>
|
||
<p>which we can rewrite as</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(X\vert Y)= \frac{p(X,Y)}{\sum_{i=0}^{n-1}p(Y\vert X=x_i)p(x_i)}=\frac{p(Y\vert X)p(X)}{\sum_{i=0}^{n-1}p(Y\vert X=x_i)p(x_i)},
|
||
\]</div>
|
||
<p>which is Bayes’ theorem. It allows us to evaluate the uncertainty in in <span class="math notranslate nohighlight">\(X\)</span> after we have observed <span class="math notranslate nohighlight">\(Y\)</span>. We can easily interchange <span class="math notranslate nohighlight">\(X\)</span> with <span class="math notranslate nohighlight">\(Y\)</span>.</p>
|
||
</div>
|
||
<div class="section" id="interpretations-of-bayes-theorem">
|
||
<h2>Interpretations of Bayes’ Theorem<a class="headerlink" href="#interpretations-of-bayes-theorem" title="Permalink to this headline">¶</a></h2>
|
||
<p>The quantity <span class="math notranslate nohighlight">\(p(Y\vert X)\)</span> on the right-hand side of the theorem is
|
||
evaluated for the observed data <span class="math notranslate nohighlight">\(Y\)</span> and can be viewed as a function of
|
||
the parameter space represented by <span class="math notranslate nohighlight">\(X\)</span>. This function is not
|
||
necesseraly normalized and is normally called the likelihood function.</p>
|
||
<p>The function <span class="math notranslate nohighlight">\(p(X)\)</span> on the right hand side is called the prior while the function on the left hand side is the called the posterior probability. The denominator on the right hand side serves as a normalization factor for the posterior distribution.</p>
|
||
<p>Let us try to illustrate Bayes’ theorem through an example.</p>
|
||
</div>
|
||
<div class="section" id="example-of-usage-of-bayes-theorem">
|
||
<h2>Example of Usage of Bayes’ theorem<a class="headerlink" href="#example-of-usage-of-bayes-theorem" title="Permalink to this headline">¶</a></h2>
|
||
<p>Let us suppose that you are undergoing a series of mammography scans in
|
||
order to rule out possible breast cancer cases. We define the
|
||
sensitivity for a positive event by the variable <span class="math notranslate nohighlight">\(X\)</span>. It takes binary
|
||
values with <span class="math notranslate nohighlight">\(X=1\)</span> representing a positive event and <span class="math notranslate nohighlight">\(X=0\)</span> being a
|
||
negative event. We reserve <span class="math notranslate nohighlight">\(Y\)</span> as a classification parameter for
|
||
either a negative or a positive breast cancer confirmation. (Short note on wordings: positive here means having breast cancer, although none of us would consider this being a positive thing).</p>
|
||
<p>We let <span class="math notranslate nohighlight">\(Y=1\)</span> represent the the case of having breast cancer and <span class="math notranslate nohighlight">\(Y=0\)</span> as not.</p>
|
||
<p>Let us assume that if you have breast cancer, the test will be positive with a probability of <span class="math notranslate nohighlight">\(0.8\)</span>, that is we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(X=1\vert Y=1) =0.8.
|
||
\]</div>
|
||
<p>This obviously sounds scary since many would conclude that if the test is positive, there is a likelihood of <span class="math notranslate nohighlight">\(80\%\)</span> for having cancer.
|
||
It is however not correct, as the following Bayesian analysis shows.</p>
|
||
</div>
|
||
<div class="section" id="doing-it-correctly">
|
||
<h2>Doing it correctly<a class="headerlink" href="#doing-it-correctly" title="Permalink to this headline">¶</a></h2>
|
||
<p>If we look at various national surveys on breast cancer, the general likelihood of developing breast cancer is a very small number.
|
||
Let us assume that the prior probability in the population as a whole is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(Y=1) =0.004.
|
||
\]</div>
|
||
<p>We need also to account for the fact that the test may produce a false positive result (false alarm). Let us here assume that we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(X=1\vert Y=0) =0.1.
|
||
\]</div>
|
||
<p>Using Bayes’ theorem we can then find the posterior probability that the person has breast cancer in case of a positive test, that is we can compute</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(Y=1\vert X=1)=\frac{p(X=1\vert Y=1)p(Y=1)}{p(X=1\vert Y=1)p(Y=1)+p(X=1\vert Y=0)p(Y=0)}=\frac{0.8\times 0.004}{0.8\times 0.004+0.1\times 0.996}=0.031.
|
||
\]</div>
|
||
<p>That is, in case of a positive test, there is only a <span class="math notranslate nohighlight">\(3\%\)</span> chance of having breast cancer!</p>
|
||
</div>
|
||
<div class="section" id="bayes-theorem-and-ridge-and-lasso-regression">
|
||
<h2>Bayes’ Theorem and Ridge and Lasso Regression<a class="headerlink" href="#bayes-theorem-and-ridge-and-lasso-regression" title="Permalink to this headline">¶</a></h2>
|
||
<p>Using Bayes’ theorem we can gain a better intuition about Ridge and Lasso regression.</p>
|
||
<p>For ordinary least squares we postulated that the maximum likelihood for the doamin of events <span class="math notranslate nohighlight">\(\boldsymbol{D}\)</span> (one-dimensional case)</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\dots, (x_{n-1},y_{n-1})],
|
||
\]</div>
|
||
<p>is given by</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(\boldsymbol{D}\vert\boldsymbol{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}.
|
||
\]</div>
|
||
<p>In Bayes’ theorem this function plays the role of the so-called likelihood. We could now ask the question what is the posterior probability of a parameter set <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> given a domain of events <span class="math notranslate nohighlight">\(\boldsymbol{D}\)</span>? That is, how can we define the posterior probability</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(\boldsymbol{\beta}\vert\boldsymbol{D}).
|
||
\]</div>
|
||
<p>Bayes’ theorem comes to our rescue here since (omitting the normalization constant)</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(\boldsymbol{\beta}\vert\boldsymbol{D})\propto p(\boldsymbol{D}\vert\boldsymbol{\beta})p(\boldsymbol{\beta}).
|
||
\]</div>
|
||
<p>We have a model for <span class="math notranslate nohighlight">\(p(\boldsymbol{D}\vert\boldsymbol{\beta})\)</span> but need one for the <strong>prior</strong> <span class="math notranslate nohighlight">\(p(\boldsymbol{\beta}\)</span>!</p>
|
||
</div>
|
||
<div class="section" id="ridge-and-bayes">
|
||
<h2>Ridge and Bayes<a class="headerlink" href="#ridge-and-bayes" title="Permalink to this headline">¶</a></h2>
|
||
<p>With the posterior probability defined by a likelihood which we have
|
||
already modeled and an unknown prior, we are now ready to make
|
||
additional models for the prior.</p>
|
||
<p>We can, based on our discussions of the variance of <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> and the mean value, assume that the prior for the values <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> is given by a Gaussian with mean value zero and variance <span class="math notranslate nohighlight">\(\tau^2\)</span>, that is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(\boldsymbol{\beta})=\prod_{j=0}^{p-1}\exp{\left(-\frac{\beta_j^2}{2\tau^2}\right)}.
|
||
\]</div>
|
||
<p>Our posterior probability becomes then (omitting the normalization factor which is just a constant)</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(\boldsymbol{\beta\vert\boldsymbol{D})}=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}\prod_{j=0}^{p-1}\exp{\left(-\frac{\beta_j^2}{2\tau^2}\right)}.
|
||
\]</div>
|
||
<p>We can now optimize this quantity with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>. As we
|
||
did for OLS, this is most conveniently done by taking the negative
|
||
logarithm of the posterior probability. Doing so and leaving out the
|
||
constants terms that do not depend on <span class="math notranslate nohighlight">\(\beta\)</span>, we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta})=\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}+\frac{1}{2\tau^2}\vert\vert\boldsymbol{\beta}\vert\vert_2^2,
|
||
\]</div>
|
||
<p>and replacing <span class="math notranslate nohighlight">\(1/2\tau^2\)</span> with <span class="math notranslate nohighlight">\(\lambda\)</span> we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta})=\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}+\lambda\vert\vert\boldsymbol{\beta}\vert\vert_2^2,
|
||
\]</div>
|
||
<p>which is our Ridge cost function! Nice, isn’t it?</p>
|
||
</div>
|
||
<div class="section" id="lasso-and-bayes">
|
||
<h2>Lasso and Bayes<a class="headerlink" href="#lasso-and-bayes" title="Permalink to this headline">¶</a></h2>
|
||
<p>To derive the Lasso cost function, we simply replace the Gaussian prior with an exponential distribution (<a class="reference external" href="https://en.wikipedia.org/wiki/Laplace_distribution">Laplace in this case</a>) with zero mean value, that is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(\boldsymbol{\beta})=\prod_{j=0}^{p-1}\exp{\left(-\frac{\vert\beta_j\vert}{\tau}\right)}.
|
||
\]</div>
|
||
<p>Our posterior probability becomes then (omitting the normalization factor which is just a constant)</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(\boldsymbol{\beta}\vert\boldsymbol{D})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}\prod_{j=0}^{p-1}\exp{\left(-\frac{\vert\beta_j\vert}{\tau}\right)}.
|
||
\]</div>
|
||
<p>Taking the negative
|
||
logarithm of the posterior probability and leaving out the
|
||
constants terms that do not depend on <span class="math notranslate nohighlight">\(\beta\)</span>, we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta}=\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}+\frac{1}{\tau}\vert\vert\boldsymbol{\beta}\vert\vert_1,
|
||
\]</div>
|
||
<p>and replacing <span class="math notranslate nohighlight">\(1/\tau\)</span> with <span class="math notranslate nohighlight">\(\lambda\)</span> we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta}=\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}+\lambda\vert\vert\boldsymbol{\beta}\vert\vert_1,
|
||
\]</div>
|
||
<p>which is our Lasso cost function!</p>
|
||
</div>
|
||
<div class="section" id="material-for-the-active-learning-sessions-tuesday-and-wednesday">
|
||
<h2>Material for the active learning sessions Tuesday and Wednesday<a class="headerlink" href="#material-for-the-active-learning-sessions-tuesday-and-wednesday" title="Permalink to this headline">¶</a></h2>
|
||
<p>The material here contains a summary from last week and discussion of SVD, Ridge and Lasso regression with examples</p>
|
||
</div>
|
||
<div class="section" id="linear-regression-and-the-svd">
|
||
<h2>Linear Regression and the SVD<a class="headerlink" href="#linear-regression-and-the-svd" title="Permalink to this headline">¶</a></h2>
|
||
<p>We used the SVD to analyse the matrix to invert in ordinary lineat regression</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{X}^T\boldsymbol{X}=\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T=\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}\boldsymbol{V}^T.
|
||
\]</div>
|
||
<p>Since the matrices here have dimension <span class="math notranslate nohighlight">\(p\times p\)</span>, with <span class="math notranslate nohighlight">\(p\)</span> corresponding to the singular values, we defined last week the matrix</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\boldsymbol{\Sigma}^T\boldsymbol{\Sigma} = \begin{bmatrix} \tilde{\boldsymbol{\Sigma}} & \boldsymbol{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} \\ \boldsymbol{0}\end{bmatrix},
|
||
\end{split}\]</div>
|
||
<p>where the tilde-matrix <span class="math notranslate nohighlight">\(\tilde{\boldsymbol{\Sigma}}\)</span> is a matrix of dimension <span class="math notranslate nohighlight">\(p\times p\)</span> containing only the singular values <span class="math notranslate nohighlight">\(\sigma_i\)</span>, that is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\tilde{\boldsymbol{\Sigma}}=\begin{bmatrix} \sigma_0 & 0 & 0 & \dots & 0 & 0 \\
|
||
0 & \sigma_1 & 0 & \dots & 0 & 0 \\
|
||
0 & 0 & \sigma_2 & \dots & 0 & 0 \\
|
||
0 & 0 & 0 & \dots & \sigma_{p-2} & 0 \\
|
||
0 & 0 & 0 & \dots & 0 & \sigma_{p-1} \\
|
||
\end{bmatrix},
|
||
\end{split}\]</div>
|
||
<p>meaning we can write</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{X}^T\boldsymbol{X}=\boldsymbol{V}\tilde{\boldsymbol{\Sigma}}^2\boldsymbol{V}^T.
|
||
\]</div>
|
||
<p>Multiplying from the right with <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span> (using the orthogonality of <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>) we get</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\left(\boldsymbol{X}^T\boldsymbol{X}\right)\boldsymbol{V}=\boldsymbol{V}\tilde{\boldsymbol{\Sigma}}^2.
|
||
\]</div>
|
||
</div>
|
||
<div class="section" id="what-does-it-mean">
|
||
<h2>What does it mean?<a class="headerlink" href="#what-does-it-mean" title="Permalink to this headline">¶</a></h2>
|
||
<p>This means the vectors <span class="math notranslate nohighlight">\(\boldsymbol{v}_i\)</span> of the orthogonal matrix <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>
|
||
are the eigenvectors of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> with eigenvalues
|
||
given by the singular values squared, that is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\left(\boldsymbol{X}^T\boldsymbol{X}\right)\boldsymbol{v}_i=\boldsymbol{v}_i\sigma_i^2.
|
||
\]</div>
|
||
<p>In other words, each non-zero singular value of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is a positive
|
||
square root of an eigenvalue of <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. It means also that
|
||
the columns of <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span> are the eigenvectors of
|
||
<span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. Since we have ordered the singular values of
|
||
<span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> in a descending order, it means that the column vectors
|
||
<span class="math notranslate nohighlight">\(\boldsymbol{v}_i\)</span> are hierarchically ordered by how much correlation they
|
||
encode from the columns of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
|
||
<p>Note that these are also the eigenvectors and eigenvalues of the
|
||
Hessian matrix.</p>
|
||
<p>If we now recall the definition of the covariance matrix (not using
|
||
Bessel’s correction) we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{C}[\boldsymbol{X}]=\frac{1}{n}\boldsymbol{X}^T\boldsymbol{X},
|
||
\]</div>
|
||
<p>meaning that every squared non-singular value of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> divided by <span class="math notranslate nohighlight">\(n\)</span> (
|
||
the number of samples) are the eigenvalues of the covariance
|
||
matrix. Every singular value of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is thus a positive square
|
||
root of an eigenvalue of <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. If the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is
|
||
self-adjoint, the singular values of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> are equal to the
|
||
absolute value of the eigenvalues of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
|
||
</div>
|
||
<div class="section" id="and-finally-boldsymbol-x-boldsymbol-x-t">
|
||
<h2>And finally <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^T\)</span><a class="headerlink" href="#and-finally-boldsymbol-x-boldsymbol-x-t" title="Permalink to this headline">¶</a></h2>
|
||
<p>For <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^T\)</span> we found</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{X}\boldsymbol{X}^T=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T=\boldsymbol{U}\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}\boldsymbol{U}^T.
|
||
\]</div>
|
||
<p>Since the matrices here have dimension <span class="math notranslate nohighlight">\(n\times n\)</span>, we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\boldsymbol{\Sigma}\boldsymbol{\Sigma}^T = \begin{bmatrix} \tilde{\boldsymbol{\Sigma}} \\ \boldsymbol{0}\\ \end{bmatrix}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} \boldsymbol{0}\\ \end{bmatrix}=\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0}\\ \end{bmatrix},
|
||
\end{split}\]</div>
|
||
<p>leading to</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\boldsymbol{X}\boldsymbol{X}^T=\boldsymbol{U}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0}\\ \end{bmatrix}\boldsymbol{U}^T.
|
||
\end{split}\]</div>
|
||
<p>Multiplying with <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> from the right gives us the eigenvalue problem</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
(\boldsymbol{X}\boldsymbol{X}^T)\boldsymbol{U}=\boldsymbol{U}\begin{bmatrix} \tilde{\boldsymbol{\Sigma}} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0}\\ \end{bmatrix}.
|
||
\end{split}\]</div>
|
||
<p>It means that the eigenvalues of <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^T\)</span> are again given by
|
||
the non-zero singular values plus now a series of zeros. The column
|
||
vectors of <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> are the eigenvectors of <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{X}^T\)</span> and
|
||
measure how much correlations are contained in the rows of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
|
||
<p>Since we will mainly be interested in the correlations among the features
|
||
of our data (the columns of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>, the quantity of interest for us are the non-zero singular
|
||
values and the column vectors of <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>.</p>
|
||
</div>
|
||
<div class="section" id="code-for-svd-and-inversion-of-matrices">
|
||
<h2>Code for SVD and Inversion of Matrices<a class="headerlink" href="#code-for-svd-and-inversion-of-matrices" title="Permalink to this headline">¶</a></h2>
|
||
<p>How do we use the SVD to invert a matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^\boldsymbol{X}\)</span> which is singular or near singular?
|
||
The simple answer is to use the linear algebra function for the computation of the pseudoinverse of a given matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>, that is</p>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span> <span class="p">[</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">],[</span><span class="mi">2</span><span class="p">,</span><span class="mi">4</span><span class="p">,</span><span class="mi">5</span><span class="p">],[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="p">,</span><span class="mi">6</span><span class="p">]])</span>
|
||
<span class="n">Xinv</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linlag</span><span class="o">.</span><span class="n">pinv</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<div class="output traceback highlight-ipythontb notranslate"><div class="highlight"><pre><span></span><span class="gt">---------------------------------------------------------------------------</span>
|
||
<span class="ne">AttributeError</span><span class="g g-Whitespace"> </span>Traceback (most recent call last)
|
||
<span class="n">Cell</span> <span class="n">In</span><span class="p">[</span><span class="mi">6</span><span class="p">],</span> <span class="n">line</span> <span class="mi">3</span>
|
||
<span class="g g-Whitespace"> </span><span class="mi">1</span> <span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="g g-Whitespace"> </span><span class="mi">2</span> <span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span> <span class="p">[</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">],[</span><span class="mi">2</span><span class="p">,</span><span class="mi">4</span><span class="p">,</span><span class="mi">5</span><span class="p">],[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="p">,</span><span class="mi">6</span><span class="p">]])</span>
|
||
<span class="ne">----> </span><span class="mi">3</span> <span class="n">Xinv</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linlag</span><span class="o">.</span><span class="n">pinv</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
|
||
|
||
<span class="nn">File ~/miniforge3/envs/myenv/lib/python3.9/site-packages/numpy/__init__.py:315,</span> in <span class="ni">__getattr__</span><span class="nt">(attr)</span>
|
||
<span class="g g-Whitespace"> </span><span class="mi">312</span> <span class="kn">from</span> <span class="nn">.testing</span> <span class="kn">import</span> <span class="n">Tester</span>
|
||
<span class="g g-Whitespace"> </span><span class="mi">313</span> <span class="k">return</span> <span class="n">Tester</span>
|
||
<span class="ne">--> </span><span class="mi">315</span> <span class="k">raise</span> <span class="ne">AttributeError</span><span class="p">(</span><span class="s2">"module </span><span class="si">{!r}</span><span class="s2"> has no attribute "</span>
|
||
<span class="g g-Whitespace"> </span><span class="mi">316</span> <span class="s2">"</span><span class="si">{!r}</span><span class="s2">"</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="vm">__name__</span><span class="p">,</span> <span class="n">attr</span><span class="p">))</span>
|
||
|
||
<span class="ne">AttributeError</span>: module 'numpy' has no attribute 'linlag'
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<p>Let us first look at a matrix which does not causes problems and write our own function where we just use the SVD.</p>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="c1"># SVD inversion</span>
|
||
<span class="k">def</span> <span class="nf">SVDinv</span><span class="p">(</span><span class="n">A</span><span class="p">):</span>
|
||
<span class="w"> </span><span class="sd">''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD).</span>
|
||
<span class="sd"> SVD is numerically more stable than the inversion algorithms provided by</span>
|
||
<span class="sd"> numpy and scipy.linalg at the cost of being slower.</span>
|
||
<span class="sd"> '''</span>
|
||
<span class="n">U</span><span class="p">,</span> <span class="n">s</span><span class="p">,</span> <span class="n">VT</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">svd</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s1">'test U'</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">transpose</span><span class="p">(</span><span class="n">U</span><span class="p">)</span> <span class="o">@</span> <span class="n">U</span> <span class="o">-</span> <span class="n">U</span> <span class="nd">@np</span><span class="o">.</span><span class="n">transpose</span><span class="p">(</span><span class="n">U</span><span class="p">)))</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s1">'test VT'</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">transpose</span><span class="p">(</span><span class="n">VT</span><span class="p">)</span> <span class="o">@</span> <span class="n">VT</span> <span class="o">-</span> <span class="n">VT</span> <span class="nd">@np</span><span class="o">.</span><span class="n">transpose</span><span class="p">(</span><span class="n">VT</span><span class="p">)))</span>
|
||
|
||
|
||
<span class="n">D</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">U</span><span class="p">),</span><span class="nb">len</span><span class="p">(</span><span class="n">VT</span><span class="p">)))</span>
|
||
<span class="n">D</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">s</span><span class="p">)</span>
|
||
<span class="n">UT</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">transpose</span><span class="p">(</span><span class="n">U</span><span class="p">);</span> <span class="n">V</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">transpose</span><span class="p">(</span><span class="n">VT</span><span class="p">);</span> <span class="n">invD</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">D</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">matmul</span><span class="p">(</span><span class="n">V</span><span class="p">,</span><span class="n">np</span><span class="o">.</span><span class="n">matmul</span><span class="p">(</span><span class="n">invD</span><span class="p">,</span><span class="n">UT</span><span class="p">))</span>
|
||
|
||
|
||
<span class="c1">#X = np.array([ [1.0, -1.0, 2.0], [1.0, 0.0, 1.0], [1.0, 2.0, -1.0], [1.0, 1.0, 0.0] ])</span>
|
||
<span class="c1"># Non-singular square matrix</span>
|
||
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span> <span class="p">[</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">],[</span><span class="mi">2</span><span class="p">,</span><span class="mi">4</span><span class="p">,</span><span class="mi">5</span><span class="p">],[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="p">,</span><span class="mi">6</span><span class="p">]])</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
|
||
<span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">transpose</span><span class="p">(</span><span class="n">X</span><span class="p">)</span> <span class="o">@</span> <span class="n">X</span>
|
||
<span class="c1"># Brute force inversion</span>
|
||
<span class="n">B</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">A</span><span class="p">)</span> <span class="c1"># here we could use np.linalg.pinv(A)</span>
|
||
<span class="n">C</span> <span class="o">=</span> <span class="n">SVDinv</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">B</span><span class="o">-</span><span class="n">C</span><span class="p">))</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="inverse-of-rectangular-matrix">
|
||
<h2>Inverse of Rectangular Matrix<a class="headerlink" href="#inverse-of-rectangular-matrix" title="Permalink to this headline">¶</a></h2>
|
||
<p>Although our matrix to invert <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> is a square matrix, our matrix may be singular.</p>
|
||
<p>The pseudoinverse is the generalization of the matrix inverse for square matrices to
|
||
rectangular matrices where the number of rows and columns are not equal.</p>
|
||
<p>It is also called the the Moore-Penrose Inverse after two independent discoverers of the method or the Generalized Inverse.
|
||
It is used for the calculation of the inverse for singular or near singular matrices and for rectangular matrices.</p>
|
||
<p>Using the SVD we can obtain the pseudoinverse of a matrix <span class="math notranslate nohighlight">\(\boldsymbol{A}\)</span> (labeled here as <span class="math notranslate nohighlight">\(\boldsymbol{A}_{\mathrm{PI}}\)</span>)</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{A}_{\mathrm{PI}}= \boldsymbol{V}\boldsymbol{D}_{\mathrm{PI}}\boldsymbol{U}^T,
|
||
\]</div>
|
||
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{D}_{\mathrm{PI}}\)</span> can be calculated by creating a diagonal matrix from <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span> where we only keep the singular values (the non-zero values). The following code computes the pseudoinvers of the matrix based on the SVD.</p>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="c1"># SVD inversion</span>
|
||
<span class="k">def</span> <span class="nf">SVDinv</span><span class="p">(</span><span class="n">A</span><span class="p">):</span>
|
||
<span class="n">U</span><span class="p">,</span> <span class="n">s</span><span class="p">,</span> <span class="n">VT</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">svd</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
|
||
<span class="c1"># reciprocals of singular values of s</span>
|
||
<span class="n">d</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="o">/</span> <span class="n">s</span>
|
||
<span class="c1"># create m x n D matrix</span>
|
||
<span class="n">D</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">A</span><span class="o">.</span><span class="n">shape</span><span class="p">)</span>
|
||
<span class="c1"># populate D with n x n diagonal matrix</span>
|
||
<span class="n">D</span><span class="p">[:</span><span class="n">A</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="p">:</span><span class="n">A</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">]]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">d</span><span class="p">)</span>
|
||
<span class="n">UT</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">transpose</span><span class="p">(</span><span class="n">U</span><span class="p">)</span>
|
||
<span class="n">V</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">transpose</span><span class="p">(</span><span class="n">VT</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">matmul</span><span class="p">(</span><span class="n">V</span><span class="p">,</span><span class="n">np</span><span class="o">.</span><span class="n">matmul</span><span class="p">(</span><span class="n">D</span><span class="o">.</span><span class="n">T</span><span class="p">,</span><span class="n">UT</span><span class="p">))</span>
|
||
|
||
|
||
<span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span> <span class="p">[</span><span class="mf">0.3</span><span class="p">,</span> <span class="mf">0.4</span><span class="p">],</span> <span class="p">[</span><span class="mf">0.5</span><span class="p">,</span> <span class="mf">0.6</span><span class="p">],</span> <span class="p">[</span><span class="mf">0.7</span><span class="p">,</span> <span class="mf">0.8</span><span class="p">],[</span><span class="mf">0.9</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">]])</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
|
||
<span class="c1"># Brute force inversion of super-collinear matrix</span>
|
||
<span class="n">B</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">pinv</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">B</span><span class="p">)</span>
|
||
<span class="c1"># Compare our own algorithm with pinv</span>
|
||
<span class="n">C</span> <span class="o">=</span> <span class="n">SVDinv</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">C</span><span class="o">-</span><span class="n">B</span><span class="p">))</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<p>As you can see from this example, our own decomposition based on the SVD agrees with the pseudoinverse algorithm provided by <strong>Numpy</strong>.</p>
|
||
</div>
|
||
<div class="section" id="ridge-and-lasso-regression">
|
||
<h2>Ridge and LASSO Regression<a class="headerlink" href="#ridge-and-lasso-regression" title="Permalink to this headline">¶</a></h2>
|
||
<p>Let us remind ourselves about the expression for the standard Mean Squared Error (MSE) which we used to define our cost function and the equations for the ordinary least squares (OLS) method, that is
|
||
our optimization problem is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
{\displaystyle \min_{\boldsymbol{\beta}\in {\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}.
|
||
\]</div>
|
||
<p>or we can state it as</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
{\displaystyle \min_{\boldsymbol{\beta}\in
|
||
{\mathbb{R}}^{p}}}\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\vert\vert_2^2,
|
||
\]</div>
|
||
<p>where we have used the definition of a norm-2 vector, that is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\vert\vert \boldsymbol{x}\vert\vert_2 = \sqrt{\sum_i x_i^2}.
|
||
\]</div>
|
||
</div>
|
||
<div class="section" id="from-ols-to-ridge-and-lasso">
|
||
<h2>From OLS to Ridge and Lasso<a class="headerlink" href="#from-ols-to-ridge-and-lasso" title="Permalink to this headline">¶</a></h2>
|
||
<p>By minimizing the above equation with respect to the parameters
|
||
<span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> we could then obtain an analytical expression for the
|
||
parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>. We can add a regularization parameter <span class="math notranslate nohighlight">\(\lambda\)</span> by
|
||
defining a new cost function to be optimized, that is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
{\displaystyle \min_{\boldsymbol{\beta}\in
|
||
{\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\vert\vert_2^2+\lambda\vert\vert \boldsymbol{\beta}\vert\vert_2^2
|
||
\]</div>
|
||
<p>which leads to the Ridge regression minimization problem where we
|
||
require that <span class="math notranslate nohighlight">\(\vert\vert \boldsymbol{\beta}\vert\vert_2^2\le t\)</span>, where <span class="math notranslate nohighlight">\(t\)</span> is
|
||
a finite number larger than zero. We do not include such a constraints in the discussions here.</p>
|
||
<p>By defining</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{X},\boldsymbol{\beta})=\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\vert\vert_2^2+\lambda\vert\vert \boldsymbol{\beta}\vert\vert_1,
|
||
\]</div>
|
||
<p>we have a new optimization equation</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
{\displaystyle \min_{\boldsymbol{\beta}\in
|
||
{\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\vert\vert_2^2+\lambda\vert\vert \boldsymbol{\beta}\vert\vert_1
|
||
\]</div>
|
||
<p>which leads to Lasso regression. Lasso stands for least absolute shrinkage and selection operator.</p>
|
||
<p>Here we have defined the norm-1 as</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\vert\vert \boldsymbol{x}\vert\vert_1 = \sum_i \vert x_i\vert.
|
||
\]</div>
|
||
</div>
|
||
<div class="section" id="deriving-the-ridge-regression-equations">
|
||
<h2>Deriving the Ridge Regression Equations<a class="headerlink" href="#deriving-the-ridge-regression-equations" title="Permalink to this headline">¶</a></h2>
|
||
<p>Using the matrix-vector expression for Ridge regression and dropping the parameter <span class="math notranslate nohighlight">\(1/n\)</span> in front of the standard means squared error equation, we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{X},\boldsymbol{\beta})=\left\{(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\right\}+\lambda\boldsymbol{\beta}^T\boldsymbol{\beta},
|
||
\]</div>
|
||
<p>and
|
||
taking the derivatives with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> we obtain then
|
||
a slightly modified matrix inversion problem which for finite values
|
||
of <span class="math notranslate nohighlight">\(\lambda\)</span> does not suffer from singularity problems. We obtain
|
||
the optimal parameters</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\hat{\boldsymbol{\beta}}_{\mathrm{Ridge}} = \left(\boldsymbol{X}^T\boldsymbol{X}+\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y},
|
||
\]</div>
|
||
<p>with <span class="math notranslate nohighlight">\(\boldsymbol{I}\)</span> being a <span class="math notranslate nohighlight">\(p\times p\)</span> identity matrix with the constraint that</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\sum_{i=0}^{p-1} \beta_i^2 \leq t,
|
||
\]</div>
|
||
<p>with <span class="math notranslate nohighlight">\(t\)</span> a finite positive number.</p>
|
||
</div>
|
||
<div class="section" id="note-on-scikit-learn">
|
||
<h2>Note on Scikit-Learn<a class="headerlink" href="#note-on-scikit-learn" title="Permalink to this headline">¶</a></h2>
|
||
<p>Note well that a library like <strong>Scikit-Learn</strong> does not include the <span class="math notranslate nohighlight">\(1/n\)</span> factor in the expression for the mean-squared error. If you include it, the optimal parameter <span class="math notranslate nohighlight">\(\beta\)</span> becomes</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\hat{\boldsymbol{\beta}}_{\mathrm{Ridge}} = \left(\boldsymbol{X}^T\boldsymbol{X}+n\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
|
||
\]</div>
|
||
<p>In our codes where we compare our own codes with <strong>Scikit-Learn</strong>, we do thus not include the <span class="math notranslate nohighlight">\(1/n\)</span> factor in the cost function.</p>
|
||
</div>
|
||
<div class="section" id="comparison-with-ols">
|
||
<h2>Comparison with OLS<a class="headerlink" href="#comparison-with-ols" title="Permalink to this headline">¶</a></h2>
|
||
<p>When we compare this with the ordinary least squares result we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\hat{\boldsymbol{\beta}}_{\mathrm{OLS}} = \left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y},
|
||
\]</div>
|
||
<p>which can lead to singular matrices. However, with the SVD, we can always compute the inverse of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>.</p>
|
||
<p>We see that Ridge regression is nothing but the standard OLS with a
|
||
modified diagonal term added to <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span>. The consequences, in
|
||
particular for our discussion of the bias-variance tradeoff are rather
|
||
interesting. We will see that for specific values of <span class="math notranslate nohighlight">\(\lambda\)</span>, we may
|
||
even reduce the variance of the optimal parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>. These topics and other related ones, will be discussed after the more linear algebra oriented analysis here.</p>
|
||
</div>
|
||
<div class="section" id="svd-analysis">
|
||
<h2>SVD analysis<a class="headerlink" href="#svd-analysis" title="Permalink to this headline">¶</a></h2>
|
||
<p>Using our insights about the SVD of the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>
|
||
We have already analyzed the OLS solutions in terms of the eigenvectors (the columns) of the right singular value matrix <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> as</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\tilde{\boldsymbol{y}}_{\mathrm{OLS}}=\boldsymbol{X}\boldsymbol{\beta} =\boldsymbol{U}\boldsymbol{U}^T\boldsymbol{y}.
|
||
\]</div>
|
||
<p>For Ridge regression this becomes</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\tilde{\boldsymbol{y}}_{\mathrm{Ridge}}=\boldsymbol{X}\boldsymbol{\beta}_{\mathrm{Ridge}} = \boldsymbol{U\Sigma V^T}\left(\boldsymbol{V}\boldsymbol{\Sigma}^2\boldsymbol{V}^T+\lambda\boldsymbol{I} \right)^{-1}(\boldsymbol{U\Sigma V^T})^T\boldsymbol{y}=\sum_{j=0}^{p-1}\boldsymbol{u}_j\boldsymbol{u}_j^T\frac{\sigma_j^2}{\sigma_j^2+\lambda}\boldsymbol{y},
|
||
\]</div>
|
||
<p>with the vectors <span class="math notranslate nohighlight">\(\boldsymbol{u}_j\)</span> being the columns of <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> from the SVD of the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
|
||
</div>
|
||
<div class="section" id="interpreting-the-ridge-results">
|
||
<h2>Interpreting the Ridge results<a class="headerlink" href="#interpreting-the-ridge-results" title="Permalink to this headline">¶</a></h2>
|
||
<p>Since <span class="math notranslate nohighlight">\(\lambda \geq 0\)</span>, it means that compared to OLS, we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\frac{\sigma_j^2}{\sigma_j^2+\lambda} \leq 1.
|
||
\]</div>
|
||
<p>Ridge regression finds the coordinates of <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> with respect to the
|
||
orthonormal basis <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span>, it then shrinks the coordinates by
|
||
<span class="math notranslate nohighlight">\(\frac{\sigma_j^2}{\sigma_j^2+\lambda}\)</span>. Recall that the SVD has
|
||
eigenvalues ordered in a descending way, that is <span class="math notranslate nohighlight">\(\sigma_i \geq
|
||
\sigma_{i+1}\)</span>.</p>
|
||
<p>For small eigenvalues <span class="math notranslate nohighlight">\(\sigma_i\)</span> it means that their contributions become less important, a fact which can be used to reduce the number of degrees of freedom. More about this when we have covered the material on a statistical interpretation of various linear regression methods.</p>
|
||
</div>
|
||
<div class="section" id="more-interpretations">
|
||
<h2>More interpretations<a class="headerlink" href="#more-interpretations" title="Permalink to this headline">¶</a></h2>
|
||
<p>For the sake of simplicity, let us assume that the design matrix is orthonormal, that is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{X}^T\boldsymbol{X}=(\boldsymbol{X}^T\boldsymbol{X})^{-1} =\boldsymbol{I}.
|
||
\]</div>
|
||
<p>In this case the standard OLS results in</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{\beta}^{\mathrm{OLS}} = \boldsymbol{X}^T\boldsymbol{y}=\sum_{i=0}^{n-1}\boldsymbol{u}_i\boldsymbol{u}_i^T\boldsymbol{y},
|
||
\]</div>
|
||
<p>and</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{\beta}^{\mathrm{Ridge}} = \left(\boldsymbol{I}+\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}=\left(1+\lambda\right)^{-1}\boldsymbol{\beta}^{\mathrm{OLS}},
|
||
\]</div>
|
||
<p>that is the Ridge estimator scales the OLS estimator by the inverse of a factor <span class="math notranslate nohighlight">\(1+\lambda\)</span>, and
|
||
the Ridge estimator converges to zero when the hyperparameter goes to
|
||
infinity.</p>
|
||
<p>We will come back to more interpreations after we have gone through some of the statistical analysis part.</p>
|
||
<p>For more discussions of Ridge and Lasso regression, <a class="reference external" href="https://arxiv.org/abs/1509.09169">Wessel van Wieringen’s</a> article is highly recommended.
|
||
Similarly, <a class="reference external" href="https://arxiv.org/abs/1803.08823">Mehta et al’s article</a> is also recommended.</p>
|
||
</div>
|
||
<div class="section" id="deriving-the-lasso-regression-equations">
|
||
<h2>Deriving the Lasso Regression Equations<a class="headerlink" href="#deriving-the-lasso-regression-equations" title="Permalink to this headline">¶</a></h2>
|
||
<p>Using the matrix-vector expression for Lasso regression, we have the following <strong>cost</strong> function</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{X},\boldsymbol{\beta})=\frac{1}{n}\left\{(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\right\}+\lambda\vert\vert\boldsymbol{\beta}\vert\vert_1,
|
||
\]</div>
|
||
<p>Taking the derivative with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> and recalling that the derivative of the absolute value is (we drop the boldfaced vector symbol for simplicity)</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\frac{d \vert \beta\vert}{d \beta}=\mathrm{sgn}(\beta)=\left\{\begin{array}{cc} 1 & \beta > 0 \\-1 & \beta < 0, \end{array}\right.
|
||
\end{split}\]</div>
|
||
<p>we have that the derivative of the cost function is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\frac{\partial C(\boldsymbol{X},\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}=-\frac{2}{n}\boldsymbol{X}^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})+\lambda sgn(\boldsymbol{\beta})=0,
|
||
\]</div>
|
||
<p>and reordering we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta}+\lambda sgn(\boldsymbol{\beta})=\boldsymbol{X}^T\boldsymbol{y}.
|
||
\]</div>
|
||
<p>This equation does not lead to a nice analytical equation as in Ridge regression or ordinary least squares. We have absorbed the factor <span class="math notranslate nohighlight">\(2/n\)</span> in a redefinition of the parameter <span class="math notranslate nohighlight">\(\lambda\)</span>. We will solve this type of problems using libraries like <strong>scikit-learn</strong>.</p>
|
||
</div>
|
||
<div class="section" id="simple-example-to-illustrate-ordinary-least-squares-ridge-and-lasso-regression">
|
||
<h2>Simple example to illustrate Ordinary Least Squares, Ridge and Lasso Regression<a class="headerlink" href="#simple-example-to-illustrate-ordinary-least-squares-ridge-and-lasso-regression" title="Permalink to this headline">¶</a></h2>
|
||
<p>Let us assume that our design matrix is given by unit (identity) matrix, that is a square diagonal matrix with ones only along the
|
||
diagonal. In this case we have an equal number of rows and columns <span class="math notranslate nohighlight">\(n=p\)</span>.</p>
|
||
<p>Our model approximation is just <span class="math notranslate nohighlight">\(\tilde{\boldsymbol{y}}=\boldsymbol{\beta}\)</span> and the mean squared error and thereby the cost function for ordinary least sqquares (OLS) is then (we drop the term <span class="math notranslate nohighlight">\(1/n\)</span>)</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta})=\sum_{i=0}^{p-1}(y_i-\beta_i)^2,
|
||
\]</div>
|
||
<p>and minimizing we have that</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\hat{\beta}_i^{\mathrm{OLS}} = y_i.
|
||
\]</div>
|
||
</div>
|
||
<div class="section" id="ridge-regression">
|
||
<h2>Ridge Regression<a class="headerlink" href="#ridge-regression" title="Permalink to this headline">¶</a></h2>
|
||
<p>For Ridge regression our cost function is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta})=\sum_{i=0}^{p-1}(y_i-\beta_i)^2+\lambda\sum_{i=0}^{p-1}\beta_i^2,
|
||
\]</div>
|
||
<p>and minimizing we have that</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\hat{\beta}_i^{\mathrm{Ridge}} = \frac{y_i}{1+\lambda}.
|
||
\]</div>
|
||
</div>
|
||
<div class="section" id="lasso-regression">
|
||
<h2>Lasso Regression<a class="headerlink" href="#lasso-regression" title="Permalink to this headline">¶</a></h2>
|
||
<p>For Lasso regression our cost function is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta})=\sum_{i=0}^{p-1}(y_i-\beta_i)^2+\lambda\sum_{i=0}^{p-1}\vert\beta_i\vert=\sum_{i=0}^{p-1}(y_i-\beta_i)^2+\lambda\sum_{i=0}^{p-1}\sqrt{\beta_i^2},
|
||
\]</div>
|
||
<p>and minimizing we have that</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
-2\sum_{i=0}^{p-1}(y_i-\beta_i)+\lambda \sum_{i=0}^{p-1}\frac{(\beta_i)}{\vert\beta_i\vert}=0,
|
||
\]</div>
|
||
<p>which leads to</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\hat{\boldsymbol{\beta}}_i^{\mathrm{Lasso}} = \left\{\begin{array}{ccc}y_i-\frac{\lambda}{2} &\mathrm{if} & y_i> \frac{\lambda}{2}\\
|
||
y_i+\frac{\lambda}{2} &\mathrm{if} & y_i< -\frac{\lambda}{2}\\
|
||
0 &\mathrm{if} & \vert y_i\vert\le \frac{\lambda}{2}\end{array}\right.\\.
|
||
\end{split}\]</div>
|
||
<p>Plotting these results shows clearly that Lasso regression suppresses (sets to zero) values of <span class="math notranslate nohighlight">\(\beta_i\)</span> for specific values of <span class="math notranslate nohighlight">\(\lambda\)</span>. Ridge regression reduces on the other hand the values of <span class="math notranslate nohighlight">\(\beta_i\)</span> as function of <span class="math notranslate nohighlight">\(\lambda\)</span>.</p>
|
||
</div>
|
||
<div class="section" id="yet-another-example">
|
||
<h2>Yet another Example<a class="headerlink" href="#yet-another-example" title="Permalink to this headline">¶</a></h2>
|
||
<p>Let us assume we have a data set with outputs/targets given by the vector</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\boldsymbol{y}=\begin{bmatrix}4 \\ 2 \\3\end{bmatrix},
|
||
\end{split}\]</div>
|
||
<p>and our inputs as a <span class="math notranslate nohighlight">\(3\times 2\)</span> design matrix</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\boldsymbol{X}=\begin{bmatrix}2 & 0\\ 0 & 1 \\ 0 & 0\end{bmatrix},
|
||
\end{split}\]</div>
|
||
<p>meaning that we have two features and two unknown parameters <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span> to be determined either by ordinary least squares, Ridge or Lasso regression.</p>
|
||
</div>
|
||
<div class="section" id="the-ols-case">
|
||
<h2>The OLS case<a class="headerlink" href="#the-ols-case" title="Permalink to this headline">¶</a></h2>
|
||
<p>For ordinary least squares (OLS) we know that the optimal solution is</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\hat{\boldsymbol{\beta}}^{\mathrm{OLS}}=\left( \boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
|
||
\]</div>
|
||
<p>Inserting the above values we obtain that</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\hat{\boldsymbol{\beta}}^{\mathrm{OLS}}=\begin{bmatrix}2 \\ 2\end{bmatrix},
|
||
\end{split}\]</div>
|
||
<p>The code which implements this simpler case is presented after the discussion of Ridge and Lasso.</p>
|
||
</div>
|
||
<div class="section" id="the-ridge-case">
|
||
<h2>The Ridge case<a class="headerlink" href="#the-ridge-case" title="Permalink to this headline">¶</a></h2>
|
||
<p>For Ridge regression we have</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\hat{\boldsymbol{\beta}}^{\mathrm{Ridge}}=\left( \boldsymbol{X}^T\boldsymbol{X}+\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}.
|
||
\]</div>
|
||
<p>Inserting the above values we obtain that</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\hat{\boldsymbol{\beta}}^{\mathrm{Ridge}}=\begin{bmatrix}\frac{8}{4+\lambda} \\ \frac{2}{1+\lambda}\end{bmatrix},
|
||
\end{split}\]</div>
|
||
<p>There is normally a constraint on the value of <span class="math notranslate nohighlight">\(\vert\vert \boldsymbol{\beta}\vert\vert_2\)</span> via the parameter <span class="math notranslate nohighlight">\(\lambda\)</span>.
|
||
Let us for simplicity assume that <span class="math notranslate nohighlight">\(\beta_0^2+\beta_1^2=1\)</span> as constraint. This will allow us to find an expression for the optimal values of <span class="math notranslate nohighlight">\(\beta\)</span> and <span class="math notranslate nohighlight">\(\lambda\)</span>.</p>
|
||
<p>To see this, let us write the cost function for Ridge regression.</p>
|
||
</div>
|
||
<div class="section" id="writing-the-cost-function">
|
||
<h2>Writing the Cost Function<a class="headerlink" href="#writing-the-cost-function" title="Permalink to this headline">¶</a></h2>
|
||
<p>We define the MSE without the <span class="math notranslate nohighlight">\(1/n\)</span> factor and have then, using that</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\boldsymbol{X}\boldsymbol{\beta}=\begin{bmatrix} 2\beta_0 \\ \beta_1 \\0 \end{bmatrix},
|
||
\end{split}\]</div>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta})=(4-2\beta_0)^2+(2-\beta_1)^2+\lambda(\beta_0^2+\beta_1^2),
|
||
\]</div>
|
||
<p>and taking the derivative with respect to <span class="math notranslate nohighlight">\(\beta_0\)</span> we get</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\beta_0=\frac{8}{4+\lambda},
|
||
\]</div>
|
||
<p>and for <span class="math notranslate nohighlight">\(\beta_1\)</span> we obtain</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\beta_1=\frac{2}{1+\lambda},
|
||
\]</div>
|
||
<p>Using the constraint for <span class="math notranslate nohighlight">\(\beta_0^2+\beta_1^2=1\)</span> we can constrain <span class="math notranslate nohighlight">\(\lambda\)</span> by solving</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\left(\frac{8}{4+\lambda}\right)^2+\left(\frac{2}{1+\lambda}\right)^2=1,
|
||
\]</div>
|
||
<p>which gives <span class="math notranslate nohighlight">\(\lambda=4.571\)</span> and <span class="math notranslate nohighlight">\(\beta_0=0.933\)</span> and <span class="math notranslate nohighlight">\(\beta_1=0.359\)</span>.</p>
|
||
</div>
|
||
<div class="section" id="lasso-case">
|
||
<h2>Lasso case<a class="headerlink" href="#lasso-case" title="Permalink to this headline">¶</a></h2>
|
||
<p>For Lasso we need now, keeping a constraint on <span class="math notranslate nohighlight">\(\vert\beta_0\vert+\vert\beta_1\vert=1\)</span>, to take the derivative of the absolute values of <span class="math notranslate nohighlight">\(\beta_0\)</span>
|
||
and <span class="math notranslate nohighlight">\(\beta_1\)</span>. This gives us the following derivatives of the cost function</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
C(\boldsymbol{\beta})=(4-2\beta_0)^2+(2-\beta_1)^2+\lambda(\vert\beta_0\vert+\vert\beta_1\vert),
|
||
\]</div>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\frac{\partial C(\boldsymbol{\beta})}{\partial \beta_0}=-4(4-2\beta_0)+\lambda\mathrm{sgn}(\beta_0)=0,
|
||
\]</div>
|
||
<p>and</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\frac{\partial C(\boldsymbol{\beta})}{\partial \beta_1}=-2(2-\beta_1)+\lambda\mathrm{sgn}(\beta_1)=0.
|
||
\]</div>
|
||
<p>We have now four cases to solve besides the trivial cases <span class="math notranslate nohighlight">\(\beta_0\)</span> and/or <span class="math notranslate nohighlight">\(\beta_1\)</span> are zero, namely</p>
|
||
<ol class="simple">
|
||
<li><p><span class="math notranslate nohighlight">\(\beta_0 > 0\)</span> and <span class="math notranslate nohighlight">\(\beta_1 > 0\)</span>,</p></li>
|
||
<li><p><span class="math notranslate nohighlight">\(\beta_0 > 0\)</span> and <span class="math notranslate nohighlight">\(\beta_1 < 0\)</span>,</p></li>
|
||
<li><p><span class="math notranslate nohighlight">\(\beta_0 < 0\)</span> and <span class="math notranslate nohighlight">\(\beta_1 > 0\)</span>,</p></li>
|
||
<li><p><span class="math notranslate nohighlight">\(\beta_0 < 0\)</span> and <span class="math notranslate nohighlight">\(\beta_1 < 0\)</span>.</p></li>
|
||
</ol>
|
||
</div>
|
||
<div class="section" id="the-first-case">
|
||
<h2>The first Case<a class="headerlink" href="#the-first-case" title="Permalink to this headline">¶</a></h2>
|
||
<p>If we consider the first case, we have then</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
-4(4-2\beta_0)+\lambda=0,
|
||
\]</div>
|
||
<p>and</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
-2(2-\beta_1)+\lambda=0.
|
||
\]</div>
|
||
<p>which yields</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\beta_0=\frac{16+\lambda}{8},
|
||
\]</div>
|
||
<p>and</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\beta_1=\frac{4+\lambda}{2}.
|
||
\]</div>
|
||
<p>Using the constraint on <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span> we can then find the optimal value of <span class="math notranslate nohighlight">\(\lambda\)</span> for the different cases. We leave this as an exercise to you.</p>
|
||
</div>
|
||
<div class="section" id="simple-code-for-solving-the-above-problem">
|
||
<h2>Simple code for solving the above problem<a class="headerlink" href="#simple-code-for-solving-the-above-problem" title="Permalink to this headline">¶</a></h2>
|
||
<p>Here we set up the OLS, Ridge and Lasso functionality in order to study the above example. Note that here we have opted for a set of values of <span class="math notranslate nohighlight">\(\lambda\)</span>, meaning that we need to perform a search in order to find the optimal values.</p>
|
||
<p>First we study and compare the OLS and Ridge results. The next code compares all three methods.</p>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">os</span>
|
||
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
|
||
<span class="k">def</span> <span class="nf">R2</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span> <span class="n">y_model</span><span class="p">):</span>
|
||
<span class="k">return</span> <span class="mi">1</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span> <span class="o">-</span> <span class="n">y_model</span><span class="p">)</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_data</span><span class="p">))</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span>
|
||
<span class="k">def</span> <span class="nf">MSE</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span><span class="n">y_model</span><span class="p">):</span>
|
||
<span class="n">n</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">size</span><span class="p">(</span><span class="n">y_model</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span><span class="o">-</span><span class="n">y_model</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">/</span><span class="n">n</span>
|
||
|
||
|
||
<span class="c1"># A seed just to ensure that the random numbers are the same for every run.</span>
|
||
<span class="c1"># Useful for eventual debugging.</span>
|
||
|
||
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span> <span class="p">[</span> <span class="p">[</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]])</span>
|
||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span> <span class="p">[</span><span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">])</span>
|
||
|
||
|
||
<span class="c1"># matrix inversion to find beta</span>
|
||
<span class="n">OLSbeta</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X</span><span class="p">)</span> <span class="o">@</span> <span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">OLSbeta</span><span class="p">)</span>
|
||
<span class="c1"># and then make the prediction</span>
|
||
<span class="n">ytildeOLS</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">OLSbeta</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Training MSE for OLS"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">ytildeOLS</span><span class="p">))</span>
|
||
<span class="n">ypredictOLS</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">OLSbeta</span>
|
||
|
||
<span class="c1"># Repeat now for Ridge regression and various values of the regularization parameter</span>
|
||
<span class="n">I</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">)</span>
|
||
<span class="c1"># Decide which values of lambda to use</span>
|
||
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">100</span>
|
||
<span class="n">MSEPredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
|
||
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
|
||
<span class="n">Ridgebeta</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X</span><span class="o">+</span><span class="n">lmb</span><span class="o">*</span><span class="n">I</span><span class="p">)</span> <span class="o">@</span> <span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y</span>
|
||
<span class="c1"># print(Ridgebeta)</span>
|
||
<span class="c1"># and then make the prediction</span>
|
||
<span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">Ridgebeta</span>
|
||
<span class="n">MSEPredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">ypredictRidge</span><span class="p">)</span>
|
||
<span class="c1"># print(MSEPredict[i])</span>
|
||
<span class="c1"># Now plot the results</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSEPredict</span><span class="p">,</span> <span class="s1">'r--'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Ridge Train'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'log10(lambda)'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'MSE'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<p>We see here that we reach a plateau. What is actually happening?</p>
|
||
</div>
|
||
<div class="section" id="with-lasso-regression">
|
||
<h2>With Lasso Regression<a class="headerlink" href="#with-lasso-regression" title="Permalink to this headline">¶</a></h2>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">os</span>
|
||
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
<span class="kn">from</span> <span class="nn">sklearn</span> <span class="kn">import</span> <span class="n">linear_model</span>
|
||
|
||
<span class="k">def</span> <span class="nf">R2</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span> <span class="n">y_model</span><span class="p">):</span>
|
||
<span class="k">return</span> <span class="mi">1</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span> <span class="o">-</span> <span class="n">y_model</span><span class="p">)</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_data</span><span class="p">))</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span>
|
||
<span class="k">def</span> <span class="nf">MSE</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span><span class="n">y_model</span><span class="p">):</span>
|
||
<span class="n">n</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">size</span><span class="p">(</span><span class="n">y_model</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span><span class="o">-</span><span class="n">y_model</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">/</span><span class="n">n</span>
|
||
|
||
|
||
<span class="c1"># A seed just to ensure that the random numbers are the same for every run.</span>
|
||
<span class="c1"># Useful for eventual debugging.</span>
|
||
|
||
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span> <span class="p">[</span> <span class="p">[</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]])</span>
|
||
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span> <span class="p">[</span><span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">])</span>
|
||
|
||
|
||
<span class="c1"># matrix inversion to find beta</span>
|
||
<span class="n">OLSbeta</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X</span><span class="p">)</span> <span class="o">@</span> <span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">OLSbeta</span><span class="p">)</span>
|
||
<span class="c1"># and then make the prediction</span>
|
||
<span class="n">ytildeOLS</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">OLSbeta</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Training MSE for OLS"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">ytildeOLS</span><span class="p">))</span>
|
||
<span class="n">ypredictOLS</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">OLSbeta</span>
|
||
|
||
<span class="c1"># Repeat now for Ridge regression and various values of the regularization parameter</span>
|
||
<span class="n">I</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">)</span>
|
||
<span class="c1"># Decide which values of lambda to use</span>
|
||
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">100</span>
|
||
<span class="n">MSERidgePredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">MSELassoPredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
|
||
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
|
||
<span class="n">Ridgebeta</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X</span><span class="o">+</span><span class="n">lmb</span><span class="o">*</span><span class="n">I</span><span class="p">)</span> <span class="o">@</span> <span class="n">X</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">Ridgebeta</span><span class="p">)</span>
|
||
<span class="c1"># and then make the prediction</span>
|
||
<span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">Ridgebeta</span>
|
||
<span class="n">MSERidgePredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">ypredictRidge</span><span class="p">)</span>
|
||
<span class="n">RegLasso</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Lasso</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
|
||
<span class="n">RegLasso</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X</span><span class="p">,</span><span class="n">y</span><span class="p">)</span>
|
||
<span class="n">ypredictLasso</span> <span class="o">=</span> <span class="n">RegLasso</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">RegLasso</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
|
||
<span class="n">MSELassoPredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">ypredictLasso</span><span class="p">)</span>
|
||
<span class="c1"># Now plot the results</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSERidgePredict</span><span class="p">,</span> <span class="s1">'r--'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Ridge Train'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSELassoPredict</span><span class="p">,</span> <span class="s1">'r--'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Lasso Train'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'log10(lambda)'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'MSE'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="another-example-now-with-a-polynomial-fit">
|
||
<h2>Another Example, now with a polynomial fit<a class="headerlink" href="#another-example-now-with-a-polynomial-fit" title="Permalink to this headline">¶</a></h2>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">os</span>
|
||
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">train_test_split</span>
|
||
<span class="kn">from</span> <span class="nn">sklearn</span> <span class="kn">import</span> <span class="n">linear_model</span>
|
||
|
||
<span class="k">def</span> <span class="nf">R2</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span> <span class="n">y_model</span><span class="p">):</span>
|
||
<span class="k">return</span> <span class="mi">1</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span> <span class="o">-</span> <span class="n">y_model</span><span class="p">)</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_data</span><span class="p">))</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span>
|
||
<span class="k">def</span> <span class="nf">MSE</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span><span class="n">y_model</span><span class="p">):</span>
|
||
<span class="n">n</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">size</span><span class="p">(</span><span class="n">y_model</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span><span class="o">-</span><span class="n">y_model</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">/</span><span class="n">n</span>
|
||
|
||
|
||
<span class="c1"># A seed just to ensure that the random numbers are the same for every run.</span>
|
||
<span class="c1"># Useful for eventual debugging.</span>
|
||
<span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">seed</span><span class="p">(</span><span class="mi">3155</span><span class="p">)</span>
|
||
|
||
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">rand</span><span class="p">(</span><span class="mi">100</span><span class="p">)</span>
|
||
<span class="n">y</span> <span class="o">=</span> <span class="mf">2.0</span><span class="o">+</span><span class="mi">5</span><span class="o">*</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="o">+</span><span class="mf">0.1</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">randn</span><span class="p">(</span><span class="mi">100</span><span class="p">)</span>
|
||
|
||
<span class="c1"># number of features p (here degree of polynomial</span>
|
||
<span class="n">p</span> <span class="o">=</span> <span class="mi">3</span>
|
||
<span class="c1"># The design matrix now as function of a given polynomial</span>
|
||
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">),</span><span class="n">p</span><span class="p">))</span>
|
||
<span class="n">X</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>
|
||
<span class="n">X</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span>
|
||
<span class="n">X</span><span class="p">[:,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="o">*</span><span class="n">x</span>
|
||
<span class="c1"># We split the data in test and training data</span>
|
||
<span class="n">X_train</span><span class="p">,</span> <span class="n">X_test</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span> <span class="n">y_test</span> <span class="o">=</span> <span class="n">train_test_split</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">test_size</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
|
||
|
||
<span class="c1"># matrix inversion to find beta</span>
|
||
<span class="n">OLSbeta</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X_train</span><span class="p">)</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y_train</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">OLSbeta</span><span class="p">)</span>
|
||
<span class="c1"># and then make the prediction</span>
|
||
<span class="n">ytildeOLS</span> <span class="o">=</span> <span class="n">X_train</span> <span class="o">@</span> <span class="n">OLSbeta</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Training MSE for OLS"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_train</span><span class="p">,</span><span class="n">ytildeOLS</span><span class="p">))</span>
|
||
<span class="n">ypredictOLS</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">@</span> <span class="n">OLSbeta</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="s2">"Test MSE OLS"</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictOLS</span><span class="p">))</span>
|
||
|
||
<span class="c1"># Repeat now for Lasso and Ridge regression and various values of the regularization parameter</span>
|
||
<span class="n">I</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="n">p</span><span class="p">,</span><span class="n">p</span><span class="p">)</span>
|
||
<span class="c1"># Decide which values of lambda to use</span>
|
||
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">100</span>
|
||
<span class="n">MSEPredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">MSETrain</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">MSELassoPredict</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">MSELassoTrain</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">4</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
|
||
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
|
||
<span class="n">Ridgebeta</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">+</span><span class="n">lmb</span><span class="o">*</span><span class="n">I</span><span class="p">)</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y_train</span>
|
||
<span class="c1"># include lasso using Scikit-Learn</span>
|
||
<span class="n">RegLasso</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Lasso</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
|
||
<span class="n">RegLasso</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span><span class="n">y_train</span><span class="p">)</span>
|
||
<span class="c1"># and then make the prediction</span>
|
||
<span class="n">ytildeRidge</span> <span class="o">=</span> <span class="n">X_train</span> <span class="o">@</span> <span class="n">Ridgebeta</span>
|
||
<span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">@</span> <span class="n">Ridgebeta</span>
|
||
<span class="n">ytildeLasso</span> <span class="o">=</span> <span class="n">RegLasso</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_train</span><span class="p">)</span>
|
||
<span class="n">ypredictLasso</span> <span class="o">=</span> <span class="n">RegLasso</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test</span><span class="p">)</span>
|
||
<span class="n">MSEPredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictRidge</span><span class="p">)</span>
|
||
<span class="n">MSETrain</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_train</span><span class="p">,</span><span class="n">ytildeRidge</span><span class="p">)</span>
|
||
<span class="n">MSELassoPredict</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredictLasso</span><span class="p">)</span>
|
||
<span class="n">MSELassoTrain</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">MSE</span><span class="p">(</span><span class="n">y_train</span><span class="p">,</span><span class="n">ytildeLasso</span><span class="p">)</span>
|
||
|
||
<span class="c1"># Now plot the results</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSETrain</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Ridge train'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSEPredict</span><span class="p">,</span> <span class="s1">'r--'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Ridge Test'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSELassoTrain</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Lasso train'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">lambdas</span><span class="p">),</span> <span class="n">MSELassoPredict</span><span class="p">,</span> <span class="s1">'r--'</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">'MSE Lasso Test'</span><span class="p">)</span>
|
||
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'log10(lambda)'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'MSE'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">()</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<script type="text/x-thebe-config">
|
||
{
|
||
requestKernel: true,
|
||
binderOptions: {
|
||
repo: "binder-examples/jupyter-stacks-datascience",
|
||
ref: "master",
|
||
},
|
||
codeMirrorConfig: {
|
||
theme: "abcdef",
|
||
mode: "python"
|
||
},
|
||
kernelOptions: {
|
||
kernelName: "python3",
|
||
path: "./."
|
||
},
|
||
predefinedOutput: true
|
||
}
|
||
</script>
|
||
<script>kernelName = 'python3'</script>
|
||
|
||
</div>
|
||
|
||
|
||
<!-- Previous / next buttons -->
|
||
<div class='prev-next-area'>
|
||
<a class='left-prev' id="prev-link" href="exercisesweek36.html" title="previous page">
|
||
<i class="fas fa-angle-left"></i>
|
||
<div class="prev-next-info">
|
||
<p class="prev-next-subtitle">previous</p>
|
||
<p class="prev-next-title">Exercises week 36</p>
|
||
</div>
|
||
</a>
|
||
<a class='right-next' id="next-link" href="exercisesweek37.html" title="next page">
|
||
<div class="prev-next-info">
|
||
<p class="prev-next-subtitle">next</p>
|
||
<p class="prev-next-title">Exercises week 37</p>
|
||
</div>
|
||
<i class="fas fa-angle-right"></i>
|
||
</a>
|
||
</div>
|
||
|
||
</div>
|
||
</div>
|
||
<footer class="footer">
|
||
<p>
|
||
|
||
By Morten Hjorth-Jensen<br/>
|
||
|
||
© Copyright 2021.<br/>
|
||
</p>
|
||
</footer>
|
||
</main>
|
||
|
||
|
||
</div>
|
||
</div>
|
||
|
||
<script src="_static/js/index.be7d3bbb2ef33a8344ce.js"></script>
|
||
|
||
</body>
|
||
</html> |