lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
de6a6d8bcb2b9a6cbdc410ef9c4c3f5a293af83d
FYS-STK4155/doc/LectureNotes/_build/html/week38.html
T
Morten Hjorth-Jensen de6a6d8bcb update
2024-10-21 13:45:48 +02:00

2777 lines
214 KiB
HTML
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" />
<title>Week 38: Logistic Regression and Optimization &#8212; 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 39" href="exercisesweek39.html" />
<link rel="prev" title="Exercises week 38" href="exercisesweek38.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">
<a class="reference internal" href="week36.html">
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 current active">
<a class="current reference internal" href="#">
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>
<li class="toctree-l1">
<a class="reference internal" href="additionweek42.html">
Exercises Week 42: Logistic Regression and Optimization, reminders from week 38 and week 40
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="week43.html">
Week 43: Deep Learning: Constructing a Neural Network code and solving differential equations
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="exercisesweek43.html">
Exercises week 43
</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/week38.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-38-lecture-monday-september-16">
Plans for week 38, lecture Monday September 16
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#suggested-reading-and-videos">
Suggested reading and videos
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#plans-for-the-lab-sessions">
Plans for the lab sessions
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#material-for-lecture-monday-september-16">
Material for lecture Monday September 16
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#logistic-regression">
Logistic Regression
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#classification-problems">
Classification problems
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#optimization-and-deep-learning">
Optimization and Deep learning
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#basics">
Basics
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#linear-classifier">
Linear classifier
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#some-selected-properties">
Some selected properties
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#simple-example">
Simple example
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#plotting-the-mean-value-for-each-group">
Plotting the mean value for each group
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-logistic-function">
The logistic function
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#examples-of-likelihood-functions-used-in-logistic-regression-and-nueral-networks">
Examples of likelihood functions used in logistic regression and nueral networks
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#two-parameters">
Two parameters
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#maximum-likelihood">
Maximum likelihood
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-cost-function-rewritten">
The cost function rewritten
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#minimizing-the-cross-entropy">
Minimizing the cross entropy
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#a-more-compact-expression">
A more compact expression
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#extending-to-more-predictors">
Extending to more predictors
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#including-more-classes">
Including more classes
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#more-classes">
More classes
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#searching-for-optimal-regularization-parameters-lambda">
Searching for Optimal Regularization Parameters
<span class="math notranslate nohighlight">
\(\lambda\)
</span>
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#grid-search">
Grid Search
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#randomized-grid-search">
Randomized Grid Search
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#wisconsin-cancer-data">
Wisconsin Cancer Data
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#using-the-correlation-matrix">
Using the correlation matrix
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#discussing-the-correlation-data">
Discussing the correlation data
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#other-measures-in-classification-studies-cancer-data-again">
Other measures in classification studies: Cancer Data again
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#optimization-the-central-part-of-any-machine-learning-algortithm">
Optimization, the central part of any Machine Learning algortithm
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#revisiting-our-logistic-regression-case">
Revisiting our Logistic Regression case
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-equations-to-solve">
The equations to solve
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#solving-using-newton-raphson-s-method">
Solving using Newton-Raphson’s method
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#brief-reminder-on-newton-raphson-s-method">
Brief reminder on Newton-Raphson’s method
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-equations">
The equations
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#simple-geometric-interpretation">
Simple geometric interpretation
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#extending-to-more-than-one-variable">
Extending to more than one variable
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#steepest-descent">
Steepest descent
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#more-on-steepest-descent">
More on Steepest descent
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-ideal">
The ideal
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-sensitiveness-of-the-gradient-descent">
The sensitiveness of the gradient descent
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#convex-functions">
Convex functions
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#convex-function">
Convex function
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#conditions-on-convex-functions">
Conditions on convex functions
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#more-on-convex-functions">
More on convex functions
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#some-simple-problems">
Some simple problems
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#revisiting-our-first-homework">
Revisiting our first homework
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#gradient-descent-example">
Gradient descent example
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-derivative-of-the-cost-loss-function">
The derivative of the cost/loss function
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-hessian-matrix">
The Hessian matrix
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#simple-program">
Simple program
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#id1">
Gradient Descent Example
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#and-a-corresponding-example-using-scikit-learn">
And a corresponding example using
<strong>
scikit-learn
</strong>
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#gradient-descent-and-ridge">
Gradient descent and Ridge
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-hessian-matrix-for-ridge-regression">
The Hessian matrix for Ridge Regression
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#program-example-for-gradient-descent-with-ridge-regression">
Program example for gradient descent with Ridge Regression
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#using-gradient-descent-methods-limitations">
Using gradient descent methods, limitations
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#challenge-yourself-the-coming-weekend">
Challenge yourself the coming weekend
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#lab-session-material-from-last-week-and-relevant-for-the-first-project">
Lab session: Material from last week and relevant for the first project
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#various-steps-in-cross-validation">
Various steps in cross-validation
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#how-to-set-up-the-cross-validation-for-ridge-and-or-lasso">
How to set up the cross-validation for Ridge and/or Lasso
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#cross-validation-in-brief">
Cross-validation in brief
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#code-example-for-cross-validation-and-k-fold-cross-validation">
Code Example for Cross-validation and
<span class="math notranslate nohighlight">
\(k\)
</span>
-fold Cross-validation
</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 38: Logistic Regression and Optimization</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-38-lecture-monday-september-16">
Plans for week 38, lecture Monday September 16
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#suggested-reading-and-videos">
Suggested reading and videos
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#plans-for-the-lab-sessions">
Plans for the lab sessions
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#material-for-lecture-monday-september-16">
Material for lecture Monday September 16
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#logistic-regression">
Logistic Regression
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#classification-problems">
Classification problems
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#optimization-and-deep-learning">
Optimization and Deep learning
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#basics">
Basics
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#linear-classifier">
Linear classifier
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#some-selected-properties">
Some selected properties
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#simple-example">
Simple example
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#plotting-the-mean-value-for-each-group">
Plotting the mean value for each group
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-logistic-function">
The logistic function
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#examples-of-likelihood-functions-used-in-logistic-regression-and-nueral-networks">
Examples of likelihood functions used in logistic regression and nueral networks
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#two-parameters">
Two parameters
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#maximum-likelihood">
Maximum likelihood
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-cost-function-rewritten">
The cost function rewritten
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#minimizing-the-cross-entropy">
Minimizing the cross entropy
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#a-more-compact-expression">
A more compact expression
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#extending-to-more-predictors">
Extending to more predictors
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#including-more-classes">
Including more classes
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#more-classes">
More classes
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#searching-for-optimal-regularization-parameters-lambda">
Searching for Optimal Regularization Parameters
<span class="math notranslate nohighlight">
\(\lambda\)
</span>
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#grid-search">
Grid Search
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#randomized-grid-search">
Randomized Grid Search
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#wisconsin-cancer-data">
Wisconsin Cancer Data
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#using-the-correlation-matrix">
Using the correlation matrix
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#discussing-the-correlation-data">
Discussing the correlation data
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#other-measures-in-classification-studies-cancer-data-again">
Other measures in classification studies: Cancer Data again
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#optimization-the-central-part-of-any-machine-learning-algortithm">
Optimization, the central part of any Machine Learning algortithm
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#revisiting-our-logistic-regression-case">
Revisiting our Logistic Regression case
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-equations-to-solve">
The equations to solve
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#solving-using-newton-raphson-s-method">
Solving using Newton-Raphson’s method
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#brief-reminder-on-newton-raphson-s-method">
Brief reminder on Newton-Raphson’s method
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-equations">
The equations
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#simple-geometric-interpretation">
Simple geometric interpretation
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#extending-to-more-than-one-variable">
Extending to more than one variable
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#steepest-descent">
Steepest descent
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#more-on-steepest-descent">
More on Steepest descent
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-ideal">
The ideal
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-sensitiveness-of-the-gradient-descent">
The sensitiveness of the gradient descent
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#convex-functions">
Convex functions
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#convex-function">
Convex function
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#conditions-on-convex-functions">
Conditions on convex functions
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#more-on-convex-functions">
More on convex functions
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#some-simple-problems">
Some simple problems
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#revisiting-our-first-homework">
Revisiting our first homework
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#gradient-descent-example">
Gradient descent example
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-derivative-of-the-cost-loss-function">
The derivative of the cost/loss function
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-hessian-matrix">
The Hessian matrix
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#simple-program">
Simple program
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#id1">
Gradient Descent Example
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#and-a-corresponding-example-using-scikit-learn">
And a corresponding example using
<strong>
scikit-learn
</strong>
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#gradient-descent-and-ridge">
Gradient descent and Ridge
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#the-hessian-matrix-for-ridge-regression">
The Hessian matrix for Ridge Regression
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#program-example-for-gradient-descent-with-ridge-regression">
Program example for gradient descent with Ridge Regression
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#using-gradient-descent-methods-limitations">
Using gradient descent methods, limitations
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#challenge-yourself-the-coming-weekend">
Challenge yourself the coming weekend
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#lab-session-material-from-last-week-and-relevant-for-the-first-project">
Lab session: Material from last week and relevant for the first project
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#various-steps-in-cross-validation">
Various steps in cross-validation
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#how-to-set-up-the-cross-validation-for-ridge-and-or-lasso">
How to set up the cross-validation for Ridge and/or Lasso
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#cross-validation-in-brief">
Cross-validation in brief
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#code-example-for-cross-validation-and-k-fold-cross-validation">
Code Example for Cross-validation and
<span class="math notranslate nohighlight">
\(k\)
</span>
-fold Cross-validation
</a>
</li>
</ul>
</nav>
</div>
</div>
</div>
<div>
<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)
doconce format html week38.do.txt --no_mako -->
<!-- dom:TITLE: Week 38: Logistic Regression and Optimization --><div class="tex2jax_ignore mathjax_ignore section" id="week-38-logistic-regression-and-optimization">
<h1>Week 38: Logistic Regression and Optimization<a class="headerlink" href="#week-38-logistic-regression-and-optimization" title="Permalink to this headline"></a></h1>
<p><strong>Morten Hjorth-Jensen</strong>, Department of Physics and Center for Computing in Science Education, University of Oslo and Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University</p>
<p>Date: <strong>September 16-20, 2024</strong></p>
<div class="section" id="plans-for-week-38-lecture-monday-september-16">
<h2>Plans for week 38, lecture Monday September 16<a class="headerlink" href="#plans-for-week-38-lecture-monday-september-16" title="Permalink to this headline"></a></h2>
<p><strong>Material for the lecture on Monday September 16.</strong></p>
<ul class="simple">
<li><p>Logistic regression as our first encounter of classification methods. From binary cases to several categories.</p></li>
<li><p>Start gradient and optimization methods</p></li>
<li><p><a class="reference external" href="https://youtu.be/c9DIfNHy2ks">Video of lecture</a></p></li>
<li><p>Whiteboard notes at <a class="reference external" href="https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2024/NotesSeptember16.pdf">https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2024/NotesSeptember16.pdf</a></p></li>
</ul>
</div>
<div class="section" id="suggested-reading-and-videos">
<h2>Suggested reading and videos<a class="headerlink" href="#suggested-reading-and-videos" title="Permalink to this headline"></a></h2>
<ul class="simple">
<li><p>Readings and Videos:</p>
<ul>
<li><p>Hastie et al 4.1, 4.2 and 4.3 on logistic regression</p></li>
<li><p>Raschka et al, pages 53-76 on Logistic regression and pages 37-52 on gradient optimization</p></li>
<li><p>For a good discussion on gradient methods, see Goodfellow et al section 4.3-4.5 and chapter 8. We will come back to the latter chapter in our discussion of Neural networks as well.</p></li>
<li><p><a class="reference external" href="https://www.youtube.com/watch?v=C5268D9t9Ak">Video on Logistic regression</a></p></li>
<li><p><a class="reference external" href="https://www.youtube.com/watch?v=yIYKR4sgzI8">Yet another video on logistic regression</a></p></li>
<li><p><a class="reference external" href="https://www.youtube.com/watch?v=sDv4f4s2SB8">Video on gradient descent</a></p></li>
</ul>
</li>
</ul>
</div>
<div class="section" id="plans-for-the-lab-sessions">
<h2>Plans for the lab sessions<a class="headerlink" href="#plans-for-the-lab-sessions" title="Permalink to this headline"></a></h2>
<p><strong>Material for the active learning sessions on Tuesday and Wednesday.</strong></p>
<ul class="simple">
<li><p>Repetition from last week on the bias-variance tradeoff</p></li>
<li><p>Resampling techniques, cross-validation examples included here, see also the lectures from last week on the bootstrap method</p></li>
<li><p>Exercise for week 38 on the bias-variance tradeoff, see also the video from the lab session from week 37 at <a class="reference external" href="https://youtu.be/omLmp_kkie0">https://youtu.be/omLmp_kkie0</a></p></li>
<li><p>Work on project 1, in particular resampling methods like cross-validation and bootstrap.</p></li>
<li><p><a class="reference external" href="https://youtu.be/T9jjWsmsd1o">Video on cross-validation from exercise session</a></p></li>
</ul>
</div>
<div class="section" id="material-for-lecture-monday-september-16">
<h2>Material for lecture Monday September 16<a class="headerlink" href="#material-for-lecture-monday-september-16" title="Permalink to this headline"></a></h2>
</div>
<div class="section" id="logistic-regression">
<h2>Logistic Regression<a class="headerlink" href="#logistic-regression" title="Permalink to this headline"></a></h2>
<p>In linear regression our main interest was centered on learning the
coefficients of a functional fit (say a polynomial) in order to be
able to predict the response of a continuous variable on some unseen
data. The fit to the continuous variable <span class="math notranslate nohighlight">\(y_i\)</span> is based on some
independent variables <span class="math notranslate nohighlight">\(\boldsymbol{x}_i\)</span>. Linear regression resulted in
analytical expressions for standard ordinary Least Squares or Ridge
regression (in terms of matrices to invert) for several quantities,
ranging from the variance and thereby the confidence intervals of the
parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> to the mean squared error. If we can invert
the product of the design matrices, linear regression gives then a
simple recipe for fitting our data.</p>
</div>
<div class="section" id="classification-problems">
<h2>Classification problems<a class="headerlink" href="#classification-problems" title="Permalink to this headline"></a></h2>
<p>Classification problems, however, are concerned with outcomes taking
the form of discrete variables (i.e. categories). We may for example,
on the basis of DNA sequencing for a number of patients, like to find
out which mutations are important for a certain disease; or based on
scans of various patients’ brains, figure out if there is a tumor or
not; or given a specific physical system, we’d like to identify its
state, say whether it is an ordered or disordered system (typical
situation in solid state physics); or classify the status of a
patient, whether she/he has a stroke or not and many other similar
situations.</p>
<p>The most common situation we encounter when we apply logistic
regression is that of two possible outcomes, normally denoted as a
binary outcome, true or false, positive or negative, success or
failure etc.</p>
</div>
<div class="section" id="optimization-and-deep-learning">
<h2>Optimization and Deep learning<a class="headerlink" href="#optimization-and-deep-learning" title="Permalink to this headline"></a></h2>
<p>Logistic regression will also serve as our stepping stone towards
neural network algorithms and supervised deep learning. For logistic
learning, the minimization of the cost function leads to a non-linear
equation in the parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>. The optimization of the
problem calls therefore for minimization algorithms. This forms the
bottle neck of all machine learning algorithms, namely how to find
reliable minima of a multi-variable function. This leads us to the
family of gradient descent methods. The latter are the working horses
of basically all modern machine learning algorithms.</p>
<p>We note also that many of the topics discussed here on logistic
regression are also commonly used in modern supervised Deep Learning
models, as we will see later.</p>
</div>
<div class="section" id="basics">
<h2>Basics<a class="headerlink" href="#basics" title="Permalink to this headline"></a></h2>
<p>We consider the case where the dependent variables, also called the
responses or the outcomes, <span class="math notranslate nohighlight">\(y_i\)</span> are discrete and only take values
from <span class="math notranslate nohighlight">\(k=0,\dots,K-1\)</span> (i.e. <span class="math notranslate nohighlight">\(K\)</span> classes).</p>
<p>The goal is to predict the
output classes from the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\in\mathbb{R}^{n\times p}\)</span>
made of <span class="math notranslate nohighlight">\(n\)</span> samples, each of which carries <span class="math notranslate nohighlight">\(p\)</span> features or predictors. The
primary goal is to identify the classes to which new unseen samples
belong.</p>
<p>Let us specialize to the case of two classes only, with outputs
<span class="math notranslate nohighlight">\(y_i=0\)</span> and <span class="math notranslate nohighlight">\(y_i=1\)</span>. Our outcomes could represent the status of a
credit card user that could default or not on her/his credit card
debt. That is</p>
<div class="math notranslate nohighlight">
\[\begin{split}
y_i = \begin{bmatrix} 0 &amp; \mathrm{no}\\ 1 &amp; \mathrm{yes} \end{bmatrix}.
\end{split}\]</div>
</div>
<div class="section" id="linear-classifier">
<h2>Linear classifier<a class="headerlink" href="#linear-classifier" title="Permalink to this headline"></a></h2>
<p>Before moving to the logistic model, let us try to use our linear
regression model to classify these two outcomes. We could for example
fit a linear model to the default case if <span class="math notranslate nohighlight">\(y_i &gt; 0.5\)</span> and the no
default case <span class="math notranslate nohighlight">\(y_i \leq 0.5\)</span>.</p>
<p>We would then have our
weighted linear combination, namely</p>
<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\boldsymbol{y} = \boldsymbol{X}^T\boldsymbol{\beta} + \boldsymbol{\epsilon},
\label{_auto1} \tag{1}
\end{equation}
\]</div>
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> is a vector representing the possible outcomes, <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is our
<span class="math notranslate nohighlight">\(n\times p\)</span> design matrix and <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> represents our estimators/predictors.</p>
</div>
<div class="section" id="some-selected-properties">
<h2>Some selected properties<a class="headerlink" href="#some-selected-properties" title="Permalink to this headline"></a></h2>
<p>The main problem with our function is that it takes values on the
entire real axis. In the case of logistic regression, however, the
labels <span class="math notranslate nohighlight">\(y_i\)</span> are discrete variables. A typical example is the credit
card data discussed below here, where we can set the state of
defaulting the debt to <span class="math notranslate nohighlight">\(y_i=1\)</span> and not to <span class="math notranslate nohighlight">\(y_i=0\)</span> for one the persons
in the data set (see the full example below).</p>
<p>One simple way to get a discrete output is to have sign
functions that map the output of a linear regressor to values <span class="math notranslate nohighlight">\(\{0,1\}\)</span>,
<span class="math notranslate nohighlight">\(f(s_i)=sign(s_i)=1\)</span> if <span class="math notranslate nohighlight">\(s_i\ge 0\)</span> and 0 if otherwise.
We will encounter this model in our first demonstration of neural networks.</p>
<p>Historically it is called the <strong>perceptron</strong> model in the machine learning
literature. This model is extremely simple. However, in many cases it is more
favorable to use a ``soft” classifier that outputs
the probability of a given category. This leads us to the logistic function.</p>
</div>
<div class="section" id="simple-example">
<h2>Simple example<a class="headerlink" href="#simple-example" title="Permalink to this headline"></a></h2>
<p>The following example on data for coronary heart disease (CHD) as function of age may serve as an illustration. In the code here we read and plot whether a person has had CHD (output = 1) or not (output = 0). This ouput is plotted the person’s against age. Clearly, the figure shows that attempting to make a standard linear regression fit may not be very meaningful.</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="c1"># Common imports</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.linear_model</span> <span class="kn">import</span> <span class="n">LinearRegression</span><span class="p">,</span> <span class="n">Ridge</span><span class="p">,</span> <span class="n">Lasso</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.utils</span> <span class="kn">import</span> <span class="n">resample</span>
<span class="kn">from</span> <span class="nn">sklearn.metrics</span> <span class="kn">import</span> <span class="n">mean_squared_error</span>
<span class="kn">from</span> <span class="nn">IPython.display</span> <span class="kn">import</span> <span class="n">display</span>
<span class="kn">from</span> <span class="nn">pylab</span> <span class="kn">import</span> <span class="n">plt</span><span class="p">,</span> <span class="n">mpl</span>
<span class="n">plt</span><span class="o">.</span><span class="n">style</span><span class="o">.</span><span class="n">use</span><span class="p">(</span><span class="s1">&#39;seaborn&#39;</span><span class="p">)</span>
<span class="n">mpl</span><span class="o">.</span><span class="n">rcParams</span><span class="p">[</span><span class="s1">&#39;font.family&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="s1">&#39;serif&#39;</span>
<span class="c1"># Where to save the figures and data files</span>
<span class="n">PROJECT_ROOT_DIR</span> <span class="o">=</span> <span class="s2">&quot;Results&quot;</span>
<span class="n">FIGURE_ID</span> <span class="o">=</span> <span class="s2">&quot;Results/FigureFiles&quot;</span>
<span class="n">DATA_ID</span> <span class="o">=</span> <span class="s2">&quot;DataFiles/&quot;</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">exists</span><span class="p">(</span><span class="n">PROJECT_ROOT_DIR</span><span class="p">):</span>
<span class="n">os</span><span class="o">.</span><span class="n">mkdir</span><span class="p">(</span><span class="n">PROJECT_ROOT_DIR</span><span class="p">)</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">exists</span><span class="p">(</span><span class="n">FIGURE_ID</span><span class="p">):</span>
<span class="n">os</span><span class="o">.</span><span class="n">makedirs</span><span class="p">(</span><span class="n">FIGURE_ID</span><span class="p">)</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">exists</span><span class="p">(</span><span class="n">DATA_ID</span><span class="p">):</span>
<span class="n">os</span><span class="o">.</span><span class="n">makedirs</span><span class="p">(</span><span class="n">DATA_ID</span><span class="p">)</span>
<span class="k">def</span> <span class="nf">image_path</span><span class="p">(</span><span class="n">fig_id</span><span class="p">):</span>
<span class="k">return</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">join</span><span class="p">(</span><span class="n">FIGURE_ID</span><span class="p">,</span> <span class="n">fig_id</span><span class="p">)</span>
<span class="k">def</span> <span class="nf">data_path</span><span class="p">(</span><span class="n">dat_id</span><span class="p">):</span>
<span class="k">return</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">join</span><span class="p">(</span><span class="n">DATA_ID</span><span class="p">,</span> <span class="n">dat_id</span><span class="p">)</span>
<span class="k">def</span> <span class="nf">save_fig</span><span class="p">(</span><span class="n">fig_id</span><span class="p">):</span>
<span class="n">plt</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="n">image_path</span><span class="p">(</span><span class="n">fig_id</span><span class="p">)</span> <span class="o">+</span> <span class="s2">&quot;.png&quot;</span><span class="p">,</span> <span class="nb">format</span><span class="o">=</span><span class="s1">&#39;png&#39;</span><span class="p">)</span>
<span class="n">infile</span> <span class="o">=</span> <span class="nb">open</span><span class="p">(</span><span class="n">data_path</span><span class="p">(</span><span class="s2">&quot;chddata.csv&quot;</span><span class="p">),</span><span class="s1">&#39;r&#39;</span><span class="p">)</span>
<span class="c1"># Read the chd data as csv file and organize the data into arrays with age group, age, and chd</span>
<span class="n">chd</span> <span class="o">=</span> <span class="n">pd</span><span class="o">.</span><span class="n">read_csv</span><span class="p">(</span><span class="n">infile</span><span class="p">,</span> <span class="n">names</span><span class="o">=</span><span class="p">(</span><span class="s1">&#39;ID&#39;</span><span class="p">,</span> <span class="s1">&#39;Age&#39;</span><span class="p">,</span> <span class="s1">&#39;Agegroup&#39;</span><span class="p">,</span> <span class="s1">&#39;CHD&#39;</span><span class="p">))</span>
<span class="n">chd</span><span class="o">.</span><span class="n">columns</span> <span class="o">=</span> <span class="p">[</span><span class="s1">&#39;ID&#39;</span><span class="p">,</span> <span class="s1">&#39;Age&#39;</span><span class="p">,</span> <span class="s1">&#39;Agegroup&#39;</span><span class="p">,</span> <span class="s1">&#39;CHD&#39;</span><span class="p">]</span>
<span class="n">output</span> <span class="o">=</span> <span class="n">chd</span><span class="p">[</span><span class="s1">&#39;CHD&#39;</span><span class="p">]</span>
<span class="n">age</span> <span class="o">=</span> <span class="n">chd</span><span class="p">[</span><span class="s1">&#39;Age&#39;</span><span class="p">]</span>
<span class="n">agegroup</span> <span class="o">=</span> <span class="n">chd</span><span class="p">[</span><span class="s1">&#39;Agegroup&#39;</span><span class="p">]</span>
<span class="n">numberID</span> <span class="o">=</span> <span class="n">chd</span><span class="p">[</span><span class="s1">&#39;ID&#39;</span><span class="p">]</span>
<span class="n">display</span><span class="p">(</span><span class="n">chd</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">age</span><span class="p">,</span> <span class="n">output</span><span class="p">,</span> <span class="n">marker</span><span class="o">=</span><span class="s1">&#39;o&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">axis</span><span class="p">([</span><span class="mi">18</span><span class="p">,</span><span class="mf">70.0</span><span class="p">,</span><span class="o">-</span><span class="mf">0.1</span><span class="p">,</span> <span class="mf">1.2</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="sa">r</span><span class="s1">&#39;Age&#39;</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="sa">r</span><span class="s1">&#39;CHD&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="sa">r</span><span class="s1">&#39;Age distribution and Coronary heart disease&#39;</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 traceback highlight-ipythontb notranslate"><div class="highlight"><pre><span></span><span class="gt">---------------------------------------------------------------------------</span>
<span class="ne">FileNotFoundError</span><span class="g g-Whitespace"> </span>Traceback (most recent call last)
<span class="nn">File ~/miniforge3/envs/myenv/lib/python3.9/site-packages/matplotlib/style/core.py:137,</span> in <span class="ni">use</span><span class="nt">(style)</span>
<span class="g g-Whitespace"> </span><span class="mi">136</span> <span class="k">try</span><span class="p">:</span>
<span class="ne">--&gt; </span><span class="mi">137</span> <span class="n">style</span> <span class="o">=</span> <span class="n">_rc_params_in_file</span><span class="p">(</span><span class="n">style</span><span class="p">)</span>
<span class="g g-Whitespace"> </span><span class="mi">138</span> <span class="k">except</span> <span class="ne">OSError</span> <span class="k">as</span> <span class="n">err</span><span class="p">:</span>
<span class="nn">File ~/miniforge3/envs/myenv/lib/python3.9/site-packages/matplotlib/__init__.py:870,</span> in <span class="ni">_rc_params_in_file</span><span class="nt">(fname, transform, fail_on_error)</span>
<span class="g g-Whitespace"> </span><span class="mi">869</span> <span class="n">rc_temp</span> <span class="o">=</span> <span class="p">{}</span>
<span class="ne">--&gt; </span><span class="mi">870</span> <span class="k">with</span> <span class="n">_open_file_or_url</span><span class="p">(</span><span class="n">fname</span><span class="p">)</span> <span class="k">as</span> <span class="n">fd</span><span class="p">:</span>
<span class="g g-Whitespace"> </span><span class="mi">871</span> <span class="k">try</span><span class="p">:</span>
<span class="nn">File ~/miniforge3/envs/myenv/lib/python3.9/contextlib.py:119,</span> in <span class="ni">_GeneratorContextManager.__enter__</span><span class="nt">(self)</span>
<span class="g g-Whitespace"> </span><span class="mi">118</span> <span class="k">try</span><span class="p">:</span>
<span class="ne">--&gt; </span><span class="mi">119</span> <span class="k">return</span> <span class="nb">next</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">gen</span><span class="p">)</span>
<span class="g g-Whitespace"> </span><span class="mi">120</span> <span class="k">except</span> <span class="ne">StopIteration</span><span class="p">:</span>
<span class="nn">File ~/miniforge3/envs/myenv/lib/python3.9/site-packages/matplotlib/__init__.py:847,</span> in <span class="ni">_open_file_or_url</span><span class="nt">(fname)</span>
<span class="g g-Whitespace"> </span><span class="mi">846</span> <span class="n">fname</span> <span class="o">=</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">expanduser</span><span class="p">(</span><span class="n">fname</span><span class="p">)</span>
<span class="ne">--&gt; </span><span class="mi">847</span> <span class="k">with</span> <span class="nb">open</span><span class="p">(</span><span class="n">fname</span><span class="p">,</span> <span class="n">encoding</span><span class="o">=</span><span class="s1">&#39;utf-8&#39;</span><span class="p">)</span> <span class="k">as</span> <span class="n">f</span><span class="p">:</span>
<span class="g g-Whitespace"> </span><span class="mi">848</span> <span class="k">yield</span> <span class="n">f</span>
<span class="ne">FileNotFoundError</span>: [Errno 2] No such file or directory: &#39;seaborn&#39;
<span class="n">The</span> <span class="n">above</span> <span class="n">exception</span> <span class="n">was</span> <span class="n">the</span> <span class="n">direct</span> <span class="n">cause</span> <span class="n">of</span> <span class="n">the</span> <span class="n">following</span> <span class="n">exception</span><span class="p">:</span>
<span class="ne">OSError</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">1</span><span class="p">],</span> <span class="n">line</span> <span class="mi">14</span>
<span class="g g-Whitespace"> </span><span class="mi">12</span> <span class="kn">from</span> <span class="nn">IPython.display</span> <span class="kn">import</span> <span class="n">display</span>
<span class="g g-Whitespace"> </span><span class="mi">13</span> <span class="kn">from</span> <span class="nn">pylab</span> <span class="kn">import</span> <span class="n">plt</span><span class="p">,</span> <span class="n">mpl</span>
<span class="ne">---&gt; </span><span class="mi">14</span> <span class="n">plt</span><span class="o">.</span><span class="n">style</span><span class="o">.</span><span class="n">use</span><span class="p">(</span><span class="s1">&#39;seaborn&#39;</span><span class="p">)</span>
<span class="g g-Whitespace"> </span><span class="mi">15</span> <span class="n">mpl</span><span class="o">.</span><span class="n">rcParams</span><span class="p">[</span><span class="s1">&#39;font.family&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="s1">&#39;serif&#39;</span>
<span class="g g-Whitespace"> </span><span class="mi">17</span> <span class="c1"># Where to save the figures and data files</span>
<span class="nn">File ~/miniforge3/envs/myenv/lib/python3.9/site-packages/matplotlib/style/core.py:139,</span> in <span class="ni">use</span><span class="nt">(style)</span>
<span class="g g-Whitespace"> </span><span class="mi">137</span> <span class="n">style</span> <span class="o">=</span> <span class="n">_rc_params_in_file</span><span class="p">(</span><span class="n">style</span><span class="p">)</span>
<span class="g g-Whitespace"> </span><span class="mi">138</span> <span class="k">except</span> <span class="ne">OSError</span> <span class="k">as</span> <span class="n">err</span><span class="p">:</span>
<span class="ne">--&gt; </span><span class="mi">139</span> <span class="k">raise</span> <span class="ne">OSError</span><span class="p">(</span>
<span class="g g-Whitespace"> </span><span class="mi">140</span> <span class="sa">f</span><span class="s2">&quot;</span><span class="si">{</span><span class="n">style</span><span class="si">!r}</span><span class="s2"> is not a valid package style, path of style &quot;</span>
<span class="g g-Whitespace"> </span><span class="mi">141</span> <span class="sa">f</span><span class="s2">&quot;file, URL of style file, or library style name (library &quot;</span>
<span class="g g-Whitespace"> </span><span class="mi">142</span> <span class="sa">f</span><span class="s2">&quot;styles are listed in `style.available`)&quot;</span><span class="p">)</span> <span class="kn">from</span> <span class="nn">err</span>
<span class="g g-Whitespace"> </span><span class="mi">143</span> <span class="n">filtered</span> <span class="o">=</span> <span class="p">{}</span>
<span class="nn"> 144 for k</span> in <span class="ni">style: # don&#39;t trigger RcParams.__getitem__</span><span class="nt">(&#39;backend&#39;)</span>
<span class="ne">OSError</span>: &#39;seaborn&#39; is not a valid package style, path of style file, URL of style file, or library style name (library styles are listed in `style.available`)
</pre></div>
</div>
</div>
</div>
</div>
<div class="section" id="plotting-the-mean-value-for-each-group">
<h2>Plotting the mean value for each group<a class="headerlink" href="#plotting-the-mean-value-for-each-group" title="Permalink to this headline"></a></h2>
<p>What we could attempt however is to plot the mean value for each group.</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="n">agegroupmean</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="mf">0.1</span><span class="p">,</span> <span class="mf">0.133</span><span class="p">,</span> <span class="mf">0.250</span><span class="p">,</span> <span class="mf">0.333</span><span class="p">,</span> <span class="mf">0.462</span><span class="p">,</span> <span class="mf">0.625</span><span class="p">,</span> <span class="mf">0.765</span><span class="p">,</span> <span class="mf">0.800</span><span class="p">])</span>
<span class="n">group</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="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">4</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="mi">7</span><span class="p">,</span> <span class="mi">8</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">group</span><span class="p">,</span> <span class="n">agegroupmean</span><span class="p">,</span> <span class="s2">&quot;r-&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">axis</span><span class="p">([</span><span class="mi">0</span><span class="p">,</span><span class="mi">9</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span> <span class="mf">1.0</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="sa">r</span><span class="s1">&#39;Age group&#39;</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="sa">r</span><span class="s1">&#39;CHD mean values&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="sa">r</span><span class="s1">&#39;Mean values for each age group&#39;</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 are now trying to find a function <span class="math notranslate nohighlight">\(f(y\vert x)\)</span>, that is a function which gives us an expected value for the output <span class="math notranslate nohighlight">\(y\)</span> with a given input <span class="math notranslate nohighlight">\(x\)</span>.
In standard linear regression with a linear dependence on <span class="math notranslate nohighlight">\(x\)</span>, we would write this in terms of our model</p>
<div class="math notranslate nohighlight">
\[
f(y_i\vert x_i)=\beta_0+\beta_1 x_i.
\]</div>
<p>This expression implies however that <span class="math notranslate nohighlight">\(f(y_i\vert x_i)\)</span> could take any
value from minus infinity to plus infinity. If we however let
<span class="math notranslate nohighlight">\(f(y\vert y)\)</span> be represented by the mean value, the above example
shows us that we can constrain the function to take values between
zero and one, that is we have <span class="math notranslate nohighlight">\(0 \le f(y_i\vert x_i) \le 1\)</span>. Looking
at our last curve we see also that it has an S-shaped form. This leads
us to a very popular model for the function <span class="math notranslate nohighlight">\(f\)</span>, namely the so-called
Sigmoid function or logistic model. We will consider this function as
representing the probability for finding a value of <span class="math notranslate nohighlight">\(y_i\)</span> with a given
<span class="math notranslate nohighlight">\(x_i\)</span>.</p>
</div>
<div class="section" id="the-logistic-function">
<h2>The logistic function<a class="headerlink" href="#the-logistic-function" title="Permalink to this headline"></a></h2>
<p>Another widely studied model, is the so-called
perceptron model, which is an example of a “hard classification” model. We
will encounter this model when we discuss neural networks as
well. Each datapoint is deterministically assigned to a category (i.e
<span class="math notranslate nohighlight">\(y_i=0\)</span> or <span class="math notranslate nohighlight">\(y_i=1\)</span>). In many cases, and the coronary heart disease data forms one of many such examples, it is favorable to have a “soft”
classifier that outputs the probability of a given category rather
than a single value. For example, given <span class="math notranslate nohighlight">\(x_i\)</span>, the classifier
outputs the probability of being in a category <span class="math notranslate nohighlight">\(k\)</span>. Logistic regression
is the most common example of a so-called soft classifier. In logistic
regression, the probability that a data point <span class="math notranslate nohighlight">\(x_i\)</span>
belongs to a category <span class="math notranslate nohighlight">\(y_i=\{0,1\}\)</span> is given by the so-called logit function (or Sigmoid) which is meant to represent the likelihood for a given event,</p>
<div class="math notranslate nohighlight">
\[
p(t) = \frac{1}{1+\mathrm \exp{-t}}=\frac{\exp{t}}{1+\mathrm \exp{t}}.
\]</div>
<p>Note that <span class="math notranslate nohighlight">\(1-p(t)= p(-t)\)</span>.</p>
</div>
<div class="section" id="examples-of-likelihood-functions-used-in-logistic-regression-and-nueral-networks">
<h2>Examples of likelihood functions used in logistic regression and nueral networks<a class="headerlink" href="#examples-of-likelihood-functions-used-in-logistic-regression-and-nueral-networks" title="Permalink to this headline"></a></h2>
<p>The following code plots the logistic function, the step function and other functions we will encounter from here and on.</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">&quot;&quot;&quot;The sigmoid function (or the logistic curve) is a</span>
<span class="sd">function that takes any real number, z, and outputs a number (0,1).</span>
<span class="sd">It is useful in neural networks for assigning weights on a relative scale.</span>
<span class="sd">The value z is the weighted sum of parameters involved in the learning algorithm.&quot;&quot;&quot;</span>
<span class="kn">import</span> <span class="nn">numpy</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">import</span> <span class="nn">math</span> <span class="k">as</span> <span class="nn">mt</span>
<span class="n">z</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="o">-</span><span class="mi">5</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="mf">.1</span><span class="p">)</span>
<span class="n">sigma_fn</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">vectorize</span><span class="p">(</span><span class="k">lambda</span> <span class="n">z</span><span class="p">:</span> <span class="mi">1</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">numpy</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">z</span><span class="p">)))</span>
<span class="n">sigma</span> <span class="o">=</span> <span class="n">sigma_fn</span><span class="p">(</span><span class="n">z</span><span class="p">)</span>
<span class="n">fig</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
<span class="n">ax</span> <span class="o">=</span> <span class="n">fig</span><span class="o">.</span><span class="n">add_subplot</span><span class="p">(</span><span class="mi">111</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">sigma</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_ylim</span><span class="p">([</span><span class="o">-</span><span class="mf">0.1</span><span class="p">,</span> <span class="mf">1.1</span><span class="p">])</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_xlim</span><span class="p">([</span><span class="o">-</span><span class="mi">5</span><span class="p">,</span><span class="mi">5</span><span class="p">])</span>
<span class="n">ax</span><span class="o">.</span><span class="n">grid</span><span class="p">(</span><span class="kc">True</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;z&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_title</span><span class="p">(</span><span class="s1">&#39;sigmoid function&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
<span class="sd">&quot;&quot;&quot;Step Function&quot;&quot;&quot;</span>
<span class="n">z</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="o">-</span><span class="mi">5</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="mf">.02</span><span class="p">)</span>
<span class="n">step_fn</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">vectorize</span><span class="p">(</span><span class="k">lambda</span> <span class="n">z</span><span class="p">:</span> <span class="mf">1.0</span> <span class="k">if</span> <span class="n">z</span> <span class="o">&gt;=</span> <span class="mf">0.0</span> <span class="k">else</span> <span class="mf">0.0</span><span class="p">)</span>
<span class="n">step</span> <span class="o">=</span> <span class="n">step_fn</span><span class="p">(</span><span class="n">z</span><span class="p">)</span>
<span class="n">fig</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
<span class="n">ax</span> <span class="o">=</span> <span class="n">fig</span><span class="o">.</span><span class="n">add_subplot</span><span class="p">(</span><span class="mi">111</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">step</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_ylim</span><span class="p">([</span><span class="o">-</span><span class="mf">0.5</span><span class="p">,</span> <span class="mf">1.5</span><span class="p">])</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_xlim</span><span class="p">([</span><span class="o">-</span><span class="mi">5</span><span class="p">,</span><span class="mi">5</span><span class="p">])</span>
<span class="n">ax</span><span class="o">.</span><span class="n">grid</span><span class="p">(</span><span class="kc">True</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;z&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_title</span><span class="p">(</span><span class="s1">&#39;step function&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
<span class="sd">&quot;&quot;&quot;tanh Function&quot;&quot;&quot;</span>
<span class="n">z</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">mt</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">mt</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span> <span class="mf">0.1</span><span class="p">)</span>
<span class="n">t</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">tanh</span><span class="p">(</span><span class="n">z</span><span class="p">)</span>
<span class="n">fig</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
<span class="n">ax</span> <span class="o">=</span> <span class="n">fig</span><span class="o">.</span><span class="n">add_subplot</span><span class="p">(</span><span class="mi">111</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">t</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_ylim</span><span class="p">([</span><span class="o">-</span><span class="mf">1.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">])</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_xlim</span><span class="p">([</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">mt</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span><span class="mi">2</span><span class="o">*</span><span class="n">mt</span><span class="o">.</span><span class="n">pi</span><span class="p">])</span>
<span class="n">ax</span><span class="o">.</span><span class="n">grid</span><span class="p">(</span><span class="kc">True</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;z&#39;</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">set_title</span><span class="p">(</span><span class="s1">&#39;tanh function&#39;</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="two-parameters">
<h2>Two parameters<a class="headerlink" href="#two-parameters" title="Permalink to this headline"></a></h2>
<p>We assume now that we have two classes with <span class="math notranslate nohighlight">\(y_i\)</span> either <span class="math notranslate nohighlight">\(0\)</span> or <span class="math notranslate nohighlight">\(1\)</span>. Furthermore we assume also that we have only two parameters <span class="math notranslate nohighlight">\(\beta\)</span> in our fitting of the Sigmoid function, that is we define probabilities</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{align*}
p(y_i=1|x_i,\boldsymbol{\beta}) &amp;= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\
p(y_i=0|x_i,\boldsymbol{\beta}) &amp;= 1 - p(y_i=1|x_i,\boldsymbol{\beta}),
\end{align*}
\end{split}\]</div>
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> are the weights we wish to extract from data, in our case <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span>.</p>
<p>Note that we used</p>
<div class="math notranslate nohighlight">
\[
p(y_i=0\vert x_i, \boldsymbol{\beta}) = 1-p(y_i=1\vert x_i, \boldsymbol{\beta}).
\]</div>
</div>
<div class="section" id="maximum-likelihood">
<h2>Maximum likelihood<a class="headerlink" href="#maximum-likelihood" title="Permalink to this headline"></a></h2>
<p>In order to define the total likelihood for all possible outcomes from a<br />
dataset <span class="math notranslate nohighlight">\(\mathcal{D}=\{(y_i,x_i)\}\)</span>, with the binary labels
<span class="math notranslate nohighlight">\(y_i\in\{0,1\}\)</span> and where the data points are drawn independently, we use the so-called <a class="reference external" href="https://en.wikipedia.org/wiki/Maximum_likelihood_estimation">Maximum Likelihood Estimation</a> (MLE) principle.
We aim thus at maximizing
the probability of seeing the observed data. We can then approximate the
likelihood in terms of the product of the individual probabilities of a specific outcome <span class="math notranslate nohighlight">\(y_i\)</span>, that is</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{align*}
P(\mathcal{D}|\boldsymbol{\beta})&amp; = \prod_{i=1}^n \left[p(y_i=1|x_i,\boldsymbol{\beta})\right]^{y_i}\left[1-p(y_i=1|x_i,\boldsymbol{\beta}))\right]^{1-y_i}\nonumber \\
\end{align*}
\end{split}\]</div>
<p>from which we obtain the log-likelihood and our <strong>cost/loss</strong> function</p>
<div class="math notranslate nohighlight">
\[
\mathcal{C}(\boldsymbol{\beta}) = \sum_{i=1}^n \left( y_i\log{p(y_i=1|x_i,\boldsymbol{\beta})} + (1-y_i)\log\left[1-p(y_i=1|x_i,\boldsymbol{\beta}))\right]\right).
\]</div>
</div>
<div class="section" id="the-cost-function-rewritten">
<h2>The cost function rewritten<a class="headerlink" href="#the-cost-function-rewritten" title="Permalink to this headline"></a></h2>
<p>Reordering the logarithms, we can rewrite the <strong>cost/loss</strong> function as</p>
<div class="math notranslate nohighlight">
\[
\mathcal{C}(\boldsymbol{\beta}) = \sum_{i=1}^n \left(y_i(\beta_0+\beta_1x_i) -\log{(1+\exp{(\beta_0+\beta_1x_i)})}\right).
\]</div>
<p>The maximum likelihood estimator is defined as the set of parameters that maximize the log-likelihood where we maximize with respect to <span class="math notranslate nohighlight">\(\beta\)</span>.
Since the cost (error) function is just the negative log-likelihood, for logistic regression we have that</p>
<div class="math notranslate nohighlight">
\[
\mathcal{C}(\boldsymbol{\beta})=-\sum_{i=1}^n \left(y_i(\beta_0+\beta_1x_i) -\log{(1+\exp{(\beta_0+\beta_1x_i)})}\right).
\]</div>
<p>This equation is known in statistics as the <strong>cross entropy</strong>. Finally, we note that just as in linear regression,
in practice we often supplement the cross-entropy with additional regularization terms, usually <span class="math notranslate nohighlight">\(L_1\)</span> and <span class="math notranslate nohighlight">\(L_2\)</span> regularization as we did for Ridge and Lasso regression.</p>
</div>
<div class="section" id="minimizing-the-cross-entropy">
<h2>Minimizing the cross entropy<a class="headerlink" href="#minimizing-the-cross-entropy" title="Permalink to this headline"></a></h2>
<p>The cross entropy is a convex function of the weights <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> and,
therefore, any local minimizer is a global minimizer.</p>
<p>Minimizing this
cost function with respect to the two parameters <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span> we obtain</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial \mathcal{C}(\boldsymbol{\beta})}{\partial \beta_0} = -\sum_{i=1}^n \left(y_i -\frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}}\right),
\]</div>
<p>and</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial \mathcal{C}(\boldsymbol{\beta})}{\partial \beta_1} = -\sum_{i=1}^n \left(y_ix_i -x_i\frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}}\right).
\]</div>
</div>
<div class="section" id="a-more-compact-expression">
<h2>A more compact expression<a class="headerlink" href="#a-more-compact-expression" title="Permalink to this headline"></a></h2>
<p>Let us now define a vector <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> with <span class="math notranslate nohighlight">\(n\)</span> elements <span class="math notranslate nohighlight">\(y_i\)</span>, an
<span class="math notranslate nohighlight">\(n\times p\)</span> matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> which contains the <span class="math notranslate nohighlight">\(x_i\)</span> values and a
vector <span class="math notranslate nohighlight">\(\boldsymbol{p}\)</span> of fitted probabilities <span class="math notranslate nohighlight">\(p(y_i\vert x_i,\boldsymbol{\beta})\)</span>. We can rewrite in a more compact form the first
derivative of cost function as</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial \mathcal{C}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -\boldsymbol{X}^T\left(\boldsymbol{y}-\boldsymbol{p}\right).
\]</div>
<p>If we in addition define a diagonal matrix <span class="math notranslate nohighlight">\(\boldsymbol{W}\)</span> with elements
<span class="math notranslate nohighlight">\(p(y_i\vert x_i,\boldsymbol{\beta})(1-p(y_i\vert x_i,\boldsymbol{\beta})\)</span>, we can obtain a compact expression of the second derivative as</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial^2 \mathcal{C}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}\partial \boldsymbol{\beta}^T} = \boldsymbol{X}^T\boldsymbol{W}\boldsymbol{X}.
\]</div>
</div>
<div class="section" id="extending-to-more-predictors">
<h2>Extending to more predictors<a class="headerlink" href="#extending-to-more-predictors" title="Permalink to this headline"></a></h2>
<p>Within a binary classification problem, we can easily expand our model to include multiple predictors. Our ratio between likelihoods is then with <span class="math notranslate nohighlight">\(p\)</span> predictors</p>
<div class="math notranslate nohighlight">
\[
\log{ \frac{p(\boldsymbol{\beta}\boldsymbol{x})}{1-p(\boldsymbol{\beta}\boldsymbol{x})}} = \beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p.
\]</div>
<p>Here we defined <span class="math notranslate nohighlight">\(\boldsymbol{x}=[1,x_1,x_2,\dots,x_p]\)</span> and <span class="math notranslate nohighlight">\(\boldsymbol{\beta}=[\beta_0, \beta_1, \dots, \beta_p]\)</span> leading to</p>
<div class="math notranslate nohighlight">
\[
p(\boldsymbol{\beta}\boldsymbol{x})=\frac{ \exp{(\beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p)}}{1+\exp{(\beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p)}}.
\]</div>
</div>
<div class="section" id="including-more-classes">
<h2>Including more classes<a class="headerlink" href="#including-more-classes" title="Permalink to this headline"></a></h2>
<p>Till now we have mainly focused on two classes, the so-called binary
system. Suppose we wish to extend to <span class="math notranslate nohighlight">\(K\)</span> classes. Let us for the sake
of simplicity assume we have only two predictors. We have then following model</p>
<div class="math notranslate nohighlight">
\[
\log{\frac{p(C=1\vert x)}{p(K\vert x)}} = \beta_{10}+\beta_{11}x_1,
\]</div>
<p>and</p>
<div class="math notranslate nohighlight">
\[
\log{\frac{p(C=2\vert x)}{p(K\vert x)}} = \beta_{20}+\beta_{21}x_1,
\]</div>
<p>and so on till the class <span class="math notranslate nohighlight">\(C=K-1\)</span> class</p>
<div class="math notranslate nohighlight">
\[
\log{\frac{p(C=K-1\vert x)}{p(K\vert x)}} = \beta_{(K-1)0}+\beta_{(K-1)1}x_1,
\]</div>
<p>and the model is specified in term of <span class="math notranslate nohighlight">\(K-1\)</span> so-called log-odds or
<strong>logit</strong> transformations.</p>
</div>
<div class="section" id="more-classes">
<h2>More classes<a class="headerlink" href="#more-classes" title="Permalink to this headline"></a></h2>
<p>In our discussion of neural networks we will encounter the above again
in terms of a slightly modified function, the so-called <strong>Softmax</strong> function.</p>
<p>The softmax function is used in various multiclass classification
methods, such as multinomial logistic regression (also known as
softmax regression), multiclass linear discriminant analysis, naive
Bayes classifiers, and artificial neural networks. Specifically, in
multinomial logistic regression and linear discriminant analysis, the
input to the function is the result of <span class="math notranslate nohighlight">\(K\)</span> distinct linear functions,
and the predicted probability for the <span class="math notranslate nohighlight">\(k\)</span>-th class given a sample
vector <span class="math notranslate nohighlight">\(\boldsymbol{x}\)</span> and a weighting vector <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> is (with two
predictors):</p>
<div class="math notranslate nohighlight">
\[
p(C=k\vert \mathbf {x} )=\frac{\exp{(\beta_{k0}+\beta_{k1}x_1)}}{1+\sum_{l=1}^{K-1}\exp{(\beta_{l0}+\beta_{l1}x_1)}}.
\]</div>
<p>It is easy to extend to more predictors. The final class is</p>
<div class="math notranslate nohighlight">
\[
p(C=K\vert \mathbf {x} )=\frac{1}{1+\sum_{l=1}^{K-1}\exp{(\beta_{l0}+\beta_{l1}x_1)}},
\]</div>
<p>and they sum to one. Our earlier discussions were all specialized to
the case with two classes only. It is easy to see from the above that
what we derived earlier is compatible with these equations.</p>
<p>To find the optimal parameters we would typically use a gradient
descent method. Newton’s method and gradient descent methods are
discussed in the material on <a class="reference external" href="https://compphysics.github.io/MachineLearning/doc/pub/Splines/html/Splines-bs.html">optimization
methods</a>.</p>
</div>
<div class="section" id="searching-for-optimal-regularization-parameters-lambda">
<h2>Searching for Optimal Regularization Parameters <span class="math notranslate nohighlight">\(\lambda\)</span><a class="headerlink" href="#searching-for-optimal-regularization-parameters-lambda" title="Permalink to this headline"></a></h2>
<p>In project 1, when using Ridge and Lasso regression, we end up
searching for the optimal parameter <span class="math notranslate nohighlight">\(\lambda\)</span> which minimizes our
selected scores (MSE or <span class="math notranslate nohighlight">\(R2\)</span> values for example). The brute force
approach, as discussed in the code here for Ridge regression, consists
in evaluating the MSE as function of different <span class="math notranslate nohighlight">\(\lambda\)</span> values.
Based on these calculations, one tries then to determine the value of the hyperparameter <span class="math notranslate nohighlight">\(\lambda\)</span>
which results in optimal scores (for example the smallest MSE or an <span class="math notranslate nohighlight">\(R2=1\)</span>).</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">2021</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="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="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"># Decide which values of lambda to use</span>
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">500</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">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">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="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">&#39;g--&#39;</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">&#39;MSE SL Ridge Test&#39;</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">&#39;log10(lambda)&#39;</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">&#39;MSE&#39;</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>Here we have performed a rather data greedy calculation as function of the regularization parameter <span class="math notranslate nohighlight">\(\lambda\)</span>. There is no resampling here. The latter can easily be added by employing the function <strong>RidgeCV</strong> instead of just calling the <strong>Ridge</strong> function. For <strong>RidgeCV</strong> we need to pass the array of <span class="math notranslate nohighlight">\(\lambda\)</span> values.
By inspecting the figure we can in turn determine which is the optimal regularization parameter.
This becomes however less functional in the long run.</p>
</div>
<div class="section" id="grid-search">
<h2>Grid Search<a class="headerlink" href="#grid-search" title="Permalink to this headline"></a></h2>
<p>An alternative is to use the so-called grid search functionality
included with the library <strong>Scikit-Learn</strong>, as demonstrated for the same
example here.</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">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.linear_model</span> <span class="kn">import</span> <span class="n">Ridge</span>
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">GridSearchCV</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">2021</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="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="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"># Decide which values of lambda to use</span>
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">10</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="c1"># create and fit a ridge regression model, testing each alpha</span>
<span class="n">model</span> <span class="o">=</span> <span class="n">Ridge</span><span class="p">()</span>
<span class="n">gridsearch</span> <span class="o">=</span> <span class="n">GridSearchCV</span><span class="p">(</span><span class="n">estimator</span><span class="o">=</span><span class="n">model</span><span class="p">,</span> <span class="n">param_grid</span><span class="o">=</span><span class="nb">dict</span><span class="p">(</span><span class="n">alpha</span><span class="o">=</span><span class="n">lambdas</span><span class="p">))</span>
<span class="n">gridsearch</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="nb">print</span><span class="p">(</span><span class="n">gridsearch</span><span class="p">)</span>
<span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">gridsearch</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"># summarize the results of the grid search</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;Best estimated lambda-value: </span><span class="si">{</span><span class="n">gridsearch</span><span class="o">.</span><span class="n">best_estimator_</span><span class="o">.</span><span class="n">alpha</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;MSE score: </span><span class="si">{</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="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;R2 score: </span><span class="si">{</span><span class="n">R2</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="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>
<p>By default the grid search function includes cross validation with
five folds. The <a class="reference external" href="https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV">Scikit-Learn
documentation</a>
contains more information on how to set the different parameters.</p>
<p>If we take out the random noise, running the above codes results in <span class="math notranslate nohighlight">\(\lambda=0\)</span> yielding the best fit.</p>
</div>
<div class="section" id="randomized-grid-search">
<h2>Randomized Grid Search<a class="headerlink" href="#randomized-grid-search" title="Permalink to this headline"></a></h2>
<p>An alternative to the above manual grid set up, is to use a random
search where the parameters are tuned from a random distribution
(uniform below) for a fixed number of iterations. A model is
constructed and evaluated for each combination of chosen parameters.
We repeat the previous example but now with a random search. Note
that values of <span class="math notranslate nohighlight">\(\lambda\)</span> are now limited to be within <span class="math notranslate nohighlight">\(x\in
[0,1]\)</span>. This domain may not be the most relevant one for the specific
case under study.</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">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.linear_model</span> <span class="kn">import</span> <span class="n">Ridge</span>
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">GridSearchCV</span>
<span class="kn">from</span> <span class="nn">scipy.stats</span> <span class="kn">import</span> <span class="n">uniform</span> <span class="k">as</span> <span class="n">randuniform</span>
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">RandomizedSearchCV</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">2021</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="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="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="n">param_grid</span> <span class="o">=</span> <span class="p">{</span><span class="s1">&#39;alpha&#39;</span><span class="p">:</span> <span class="n">randuniform</span><span class="p">()}</span>
<span class="c1"># create and fit a ridge regression model, testing each alpha</span>
<span class="n">model</span> <span class="o">=</span> <span class="n">Ridge</span><span class="p">()</span>
<span class="n">gridsearch</span> <span class="o">=</span> <span class="n">RandomizedSearchCV</span><span class="p">(</span><span class="n">estimator</span><span class="o">=</span><span class="n">model</span><span class="p">,</span> <span class="n">param_distributions</span><span class="o">=</span><span class="n">param_grid</span><span class="p">,</span> <span class="n">n_iter</span><span class="o">=</span><span class="mi">100</span><span class="p">)</span>
<span class="n">gridsearch</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="nb">print</span><span class="p">(</span><span class="n">gridsearch</span><span class="p">)</span>
<span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">gridsearch</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"># summarize the results of the grid search</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;Best estimated lambda-value: </span><span class="si">{</span><span class="n">gridsearch</span><span class="o">.</span><span class="n">best_estimator_</span><span class="o">.</span><span class="n">alpha</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;MSE score: </span><span class="si">{</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="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;R2 score: </span><span class="si">{</span><span class="n">R2</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="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="section" id="wisconsin-cancer-data">
<h2>Wisconsin Cancer Data<a class="headerlink" href="#wisconsin-cancer-data" title="Permalink to this headline"></a></h2>
<p>We show here how we can use a simple regression case on the breast
cancer data using Logistic regression as our algorithm for
classification.</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">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</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.datasets</span> <span class="kn">import</span> <span class="n">load_breast_cancer</span>
<span class="kn">from</span> <span class="nn">sklearn.linear_model</span> <span class="kn">import</span> <span class="n">LogisticRegression</span>
<span class="c1"># Load the data</span>
<span class="n">cancer</span> <span class="o">=</span> <span class="n">load_breast_cancer</span><span class="p">()</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">cancer</span><span class="o">.</span><span class="n">data</span><span class="p">,</span><span class="n">cancer</span><span class="o">.</span><span class="n">target</span><span class="p">,</span><span class="n">random_state</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">X_train</span><span class="o">.</span><span class="n">shape</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">X_test</span><span class="o">.</span><span class="n">shape</span><span class="p">)</span>
<span class="c1"># Logistic Regression</span>
<span class="n">logreg</span> <span class="o">=</span> <span class="n">LogisticRegression</span><span class="p">(</span><span class="n">solver</span><span class="o">=</span><span class="s1">&#39;lbfgs&#39;</span><span class="p">)</span>
<span class="n">logreg</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="nb">print</span><span class="p">(</span><span class="s2">&quot;Test set accuracy with Logistic Regression: </span><span class="si">{:.2f}</span><span class="s2">&quot;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">logreg</span><span class="o">.</span><span class="n">score</span><span class="p">(</span><span class="n">X_test</span><span class="p">,</span><span class="n">y_test</span><span class="p">)))</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="section" id="using-the-correlation-matrix">
<h2>Using the correlation matrix<a class="headerlink" href="#using-the-correlation-matrix" title="Permalink to this headline"></a></h2>
<p>In addition to the above scores, we could also study the covariance (and the correlation matrix).
We use <strong>Pandas</strong> to compute the correlation matrix.</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">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</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.datasets</span> <span class="kn">import</span> <span class="n">load_breast_cancer</span>
<span class="kn">from</span> <span class="nn">sklearn.linear_model</span> <span class="kn">import</span> <span class="n">LogisticRegression</span>
<span class="n">cancer</span> <span class="o">=</span> <span class="n">load_breast_cancer</span><span class="p">()</span>
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
<span class="c1"># Making a data frame</span>
<span class="n">cancerpd</span> <span class="o">=</span> <span class="n">pd</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">(</span><span class="n">cancer</span><span class="o">.</span><span class="n">data</span><span class="p">,</span> <span class="n">columns</span><span class="o">=</span><span class="n">cancer</span><span class="o">.</span><span class="n">feature_names</span><span class="p">)</span>
<span class="n">fig</span><span class="p">,</span> <span class="n">axes</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">(</span><span class="mi">15</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="n">figsize</span><span class="o">=</span><span class="p">(</span><span class="mi">10</span><span class="p">,</span><span class="mi">20</span><span class="p">))</span>
<span class="n">malignant</span> <span class="o">=</span> <span class="n">cancer</span><span class="o">.</span><span class="n">data</span><span class="p">[</span><span class="n">cancer</span><span class="o">.</span><span class="n">target</span> <span class="o">==</span> <span class="mi">0</span><span class="p">]</span>
<span class="n">benign</span> <span class="o">=</span> <span class="n">cancer</span><span class="o">.</span><span class="n">data</span><span class="p">[</span><span class="n">cancer</span><span class="o">.</span><span class="n">target</span> <span class="o">==</span> <span class="mi">1</span><span class="p">]</span>
<span class="n">ax</span> <span class="o">=</span> <span class="n">axes</span><span class="o">.</span><span class="n">ravel</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="mi">30</span><span class="p">):</span>
<span class="n">_</span><span class="p">,</span> <span class="n">bins</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">histogram</span><span class="p">(</span><span class="n">cancer</span><span class="o">.</span><span class="n">data</span><span class="p">[:,</span><span class="n">i</span><span class="p">],</span> <span class="n">bins</span> <span class="o">=</span><span class="mi">50</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">hist</span><span class="p">(</span><span class="n">malignant</span><span class="p">[:,</span><span class="n">i</span><span class="p">],</span> <span class="n">bins</span> <span class="o">=</span> <span class="n">bins</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">0.5</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">hist</span><span class="p">(</span><span class="n">benign</span><span class="p">[:,</span><span class="n">i</span><span class="p">],</span> <span class="n">bins</span> <span class="o">=</span> <span class="n">bins</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">0.5</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">set_title</span><span class="p">(</span><span class="n">cancer</span><span class="o">.</span><span class="n">feature_names</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
<span class="n">ax</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">set_yticks</span><span class="p">(())</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s2">&quot;Feature magnitude&quot;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s2">&quot;Frequency&quot;</span><span class="p">)</span>
<span class="n">ax</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">legend</span><span class="p">([</span><span class="s2">&quot;Malignant&quot;</span><span class="p">,</span> <span class="s2">&quot;Benign&quot;</span><span class="p">],</span> <span class="n">loc</span> <span class="o">=</span><span class="s2">&quot;best&quot;</span><span class="p">)</span>
<span class="n">fig</span><span class="o">.</span><span class="n">tight_layout</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
<span class="kn">import</span> <span class="nn">seaborn</span> <span class="k">as</span> <span class="nn">sns</span>
<span class="n">correlation_matrix</span> <span class="o">=</span> <span class="n">cancerpd</span><span class="o">.</span><span class="n">corr</span><span class="p">()</span><span class="o">.</span><span class="n">round</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
<span class="c1"># use the heatmap function from seaborn to plot the correlation matrix</span>
<span class="c1"># annot = True to print the values inside the square</span>
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">(</span><span class="n">figsize</span><span class="o">=</span><span class="p">(</span><span class="mi">15</span><span class="p">,</span><span class="mi">8</span><span class="p">))</span>
<span class="n">sns</span><span class="o">.</span><span class="n">heatmap</span><span class="p">(</span><span class="n">data</span><span class="o">=</span><span class="n">correlation_matrix</span><span class="p">,</span> <span class="n">annot</span><span class="o">=</span><span class="kc">True</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="discussing-the-correlation-data">
<h2>Discussing the correlation data<a class="headerlink" href="#discussing-the-correlation-data" title="Permalink to this headline"></a></h2>
<p>In the above example we note two things. In the first plot we display
the overlap of benign and malignant tumors as functions of the various
features in the Wisconsing breast cancer data set. We see that for
some of the features we can distinguish clearly the benign and
malignant cases while for other features we cannot. This can point to
us which features may be of greater interest when we wish to classify
a benign or not benign tumour.</p>
<p>In the second figure we have computed the so-called correlation
matrix, which in our case with thirty features becomes a <span class="math notranslate nohighlight">\(30\times 30\)</span>
matrix.</p>
<p>We constructed this matrix using <strong>pandas</strong> via the statements</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="n">cancerpd</span> <span class="o">=</span> <span class="n">pd</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">(</span><span class="n">cancer</span><span class="o">.</span><span class="n">data</span><span class="p">,</span> <span class="n">columns</span><span class="o">=</span><span class="n">cancer</span><span class="o">.</span><span class="n">feature_names</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>
<p>and then</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="n">correlation_matrix</span> <span class="o">=</span> <span class="n">cancerpd</span><span class="o">.</span><span class="n">corr</span><span class="p">()</span><span class="o">.</span><span class="n">round</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>
<p>Diagonalizing this matrix we can in turn say something about which
features are of relevance and which are not. This leads us to
the classical Principal Component Analysis (PCA) theorem with
applications. This will be discussed later this semester (<a class="reference external" href="https://compphysics.github.io/MachineLearning/doc/pub/week43/html/week43-bs.html">week 43</a>).</p>
</div>
<div class="section" id="other-measures-in-classification-studies-cancer-data-again">
<h2>Other measures in classification studies: Cancer Data again<a class="headerlink" href="#other-measures-in-classification-studies-cancer-data-again" 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">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</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.datasets</span> <span class="kn">import</span> <span class="n">load_breast_cancer</span>
<span class="kn">from</span> <span class="nn">sklearn.linear_model</span> <span class="kn">import</span> <span class="n">LogisticRegression</span>
<span class="c1"># Load the data</span>
<span class="n">cancer</span> <span class="o">=</span> <span class="n">load_breast_cancer</span><span class="p">()</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">cancer</span><span class="o">.</span><span class="n">data</span><span class="p">,</span><span class="n">cancer</span><span class="o">.</span><span class="n">target</span><span class="p">,</span><span class="n">random_state</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">X_train</span><span class="o">.</span><span class="n">shape</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">X_test</span><span class="o">.</span><span class="n">shape</span><span class="p">)</span>
<span class="c1"># Logistic Regression</span>
<span class="n">logreg</span> <span class="o">=</span> <span class="n">LogisticRegression</span><span class="p">(</span><span class="n">solver</span><span class="o">=</span><span class="s1">&#39;lbfgs&#39;</span><span class="p">)</span>
<span class="n">logreg</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="kn">from</span> <span class="nn">sklearn.preprocessing</span> <span class="kn">import</span> <span class="n">LabelEncoder</span>
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">cross_validate</span>
<span class="c1">#Cross validation</span>
<span class="n">accuracy</span> <span class="o">=</span> <span class="n">cross_validate</span><span class="p">(</span><span class="n">logreg</span><span class="p">,</span><span class="n">X_test</span><span class="p">,</span><span class="n">y_test</span><span class="p">,</span><span class="n">cv</span><span class="o">=</span><span class="mi">10</span><span class="p">)[</span><span class="s1">&#39;test_score&#39;</span><span class="p">]</span>
<span class="nb">print</span><span class="p">(</span><span class="n">accuracy</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;Test set accuracy with Logistic Regression: </span><span class="si">{:.2f}</span><span class="s2">&quot;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">logreg</span><span class="o">.</span><span class="n">score</span><span class="p">(</span><span class="n">X_test</span><span class="p">,</span><span class="n">y_test</span><span class="p">)))</span>
<span class="kn">import</span> <span class="nn">scikitplot</span> <span class="k">as</span> <span class="nn">skplt</span>
<span class="n">y_pred</span> <span class="o">=</span> <span class="n">logreg</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">skplt</span><span class="o">.</span><span class="n">metrics</span><span class="o">.</span><span class="n">plot_confusion_matrix</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span> <span class="n">y_pred</span><span class="p">,</span> <span class="n">normalize</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
<span class="n">y_probas</span> <span class="o">=</span> <span class="n">logreg</span><span class="o">.</span><span class="n">predict_proba</span><span class="p">(</span><span class="n">X_test</span><span class="p">)</span>
<span class="n">skplt</span><span class="o">.</span><span class="n">metrics</span><span class="o">.</span><span class="n">plot_roc</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span> <span class="n">y_probas</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
<span class="n">skplt</span><span class="o">.</span><span class="n">metrics</span><span class="o">.</span><span class="n">plot_cumulative_gain</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span> <span class="n">y_probas</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="optimization-the-central-part-of-any-machine-learning-algortithm">
<h2>Optimization, the central part of any Machine Learning algortithm<a class="headerlink" href="#optimization-the-central-part-of-any-machine-learning-algortithm" title="Permalink to this headline"></a></h2>
<p><a class="reference external" href="https://www.uio.no/studier/emner/matnat/fys/FYS-STK3155/h20/forelesningsvideoer/OverarchingAimsWeek39.mp4?vrtx=view-as-webpage">Overview Video, why do we care about gradient methods?</a></p>
<p>Almost every problem in machine learning and data science starts with
a dataset <span class="math notranslate nohighlight">\(X\)</span>, a model <span class="math notranslate nohighlight">\(g(\beta)\)</span>, which is a function of the
parameters <span class="math notranslate nohighlight">\(\beta\)</span> and a cost function <span class="math notranslate nohighlight">\(C(X, g(\beta))\)</span> that allows
us to judge how well the model <span class="math notranslate nohighlight">\(g(\beta)\)</span> explains the observations
<span class="math notranslate nohighlight">\(X\)</span>. The model is fit by finding the values of <span class="math notranslate nohighlight">\(\beta\)</span> that minimize
the cost function. Ideally we would be able to solve for <span class="math notranslate nohighlight">\(\beta\)</span>
analytically, however this is not possible in general and we must use
some approximative/numerical method to compute the minimum.</p>
</div>
<div class="section" id="revisiting-our-logistic-regression-case">
<h2>Revisiting our Logistic Regression case<a class="headerlink" href="#revisiting-our-logistic-regression-case" title="Permalink to this headline"></a></h2>
<p>In our discussion on Logistic Regression we studied the
case of
two classes, with <span class="math notranslate nohighlight">\(y_i\)</span> either
<span class="math notranslate nohighlight">\(0\)</span> or <span class="math notranslate nohighlight">\(1\)</span>. Furthermore we assumed also that we have only two
parameters <span class="math notranslate nohighlight">\(\beta\)</span> in our fitting, that is we
defined probabilities</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{align*}
p(y_i=1|x_i,\boldsymbol{\beta}) &amp;= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\
p(y_i=0|x_i,\boldsymbol{\beta}) &amp;= 1 - p(y_i=1|x_i,\boldsymbol{\beta}),
\end{align*}
\end{split}\]</div>
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> are the weights we wish to extract from data, in our case <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span>.</p>
</div>
<div class="section" id="the-equations-to-solve">
<h2>The equations to solve<a class="headerlink" href="#the-equations-to-solve" title="Permalink to this headline"></a></h2>
<p>Our compact equations used a definition of a vector <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> with <span class="math notranslate nohighlight">\(n\)</span>
elements <span class="math notranslate nohighlight">\(y_i\)</span>, an <span class="math notranslate nohighlight">\(n\times p\)</span> matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> which contains the
<span class="math notranslate nohighlight">\(x_i\)</span> values and a vector <span class="math notranslate nohighlight">\(\boldsymbol{p}\)</span> of fitted probabilities
<span class="math notranslate nohighlight">\(p(y_i\vert x_i,\boldsymbol{\beta})\)</span>. We rewrote in a more compact form
the first derivative of the cost function as</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial \mathcal{C}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -\boldsymbol{X}^T\left(\boldsymbol{y}-\boldsymbol{p}\right).
\]</div>
<p>If we in addition define a diagonal matrix <span class="math notranslate nohighlight">\(\boldsymbol{W}\)</span> with elements
<span class="math notranslate nohighlight">\(p(y_i\vert x_i,\boldsymbol{\beta})(1-p(y_i\vert x_i,\boldsymbol{\beta})\)</span>, we can obtain a compact expression of the second derivative as</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial^2 \mathcal{C}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}\partial \boldsymbol{\beta}^T} = \boldsymbol{X}^T\boldsymbol{W}\boldsymbol{X}.
\]</div>
<p>This defines what is called the Hessian matrix.</p>
</div>
<div class="section" id="solving-using-newton-raphson-s-method">
<h2>Solving using Newton-Raphson’s method<a class="headerlink" href="#solving-using-newton-raphson-s-method" title="Permalink to this headline"></a></h2>
<p>If we can set up these equations, Newton-Raphson’s iterative method is normally the method of choice. It requires however that we can compute in an efficient way the matrices that define the first and second derivatives.</p>
<p>Our iterative scheme is then given by</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\beta}^{\mathrm{new}} = \boldsymbol{\beta}^{\mathrm{old}}-\left(\frac{\partial^2 \mathcal{C}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}\partial \boldsymbol{\beta}^T}\right)^{-1}_{\boldsymbol{\beta}^{\mathrm{old}}}\times \left(\frac{\partial \mathcal{C}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\right)_{\boldsymbol{\beta}^{\mathrm{old}}},
\]</div>
<p>or in matrix form as</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\beta}^{\mathrm{new}} = \boldsymbol{\beta}^{\mathrm{old}}-\left(\boldsymbol{X}^T\boldsymbol{W}\boldsymbol{X} \right)^{-1}\times \left(-\boldsymbol{X}^T(\boldsymbol{y}-\boldsymbol{p}) \right)_{\boldsymbol{\beta}^{\mathrm{old}}}.
\]</div>
<p>The right-hand side is computed with the old values of <span class="math notranslate nohighlight">\(\beta\)</span>.</p>
<p>If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement.</p>
</div>
<div class="section" id="brief-reminder-on-newton-raphson-s-method">
<h2>Brief reminder on Newton-Raphson’s method<a class="headerlink" href="#brief-reminder-on-newton-raphson-s-method" title="Permalink to this headline"></a></h2>
<p>Let us quickly remind ourselves how we derive the above method.</p>
<p>Perhaps the most celebrated of all one-dimensional root-finding
routines is Newton’s method, also called the Newton-Raphson
method. This method requires the evaluation of both the
function <span class="math notranslate nohighlight">\(f\)</span> and its derivative <span class="math notranslate nohighlight">\(f'\)</span> at arbitrary points.
If you can only calculate the derivative
numerically and/or your function is not of the smooth type, we
normally discourage the use of this method.</p>
</div>
<div class="section" id="the-equations">
<h2>The equations<a class="headerlink" href="#the-equations" title="Permalink to this headline"></a></h2>
<p>The Newton-Raphson formula consists geometrically of extending the
tangent line at a current point until it crosses zero, then setting
the next guess to the abscissa of that zero-crossing. The mathematics
behind this method is rather simple. Employing a Taylor expansion for
<span class="math notranslate nohighlight">\(x\)</span> sufficiently close to the solution <span class="math notranslate nohighlight">\(s\)</span>, we have</p>
<!-- Equation labels as ordinary links -->
<div id="eq:taylornr"></div>
<div class="math notranslate nohighlight">
\[
f(s)=0=f(x)+(s-x)f'(x)+\frac{(s-x)^2}{2}f''(x) +\dots.
\label{eq:taylornr} \tag{2}
\]</div>
<p>For small enough values of the function and for well-behaved
functions, the terms beyond linear are unimportant, hence we obtain</p>
<div class="math notranslate nohighlight">
\[
f(x)+(s-x)f'(x)\approx 0,
\]</div>
<p>yielding</p>
<div class="math notranslate nohighlight">
\[
s\approx x-\frac{f(x)}{f'(x)}.
\]</div>
<p>Having in mind an iterative procedure, it is natural to start iterating with</p>
<div class="math notranslate nohighlight">
\[
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.
\]</div>
</div>
<div class="section" id="simple-geometric-interpretation">
<h2>Simple geometric interpretation<a class="headerlink" href="#simple-geometric-interpretation" title="Permalink to this headline"></a></h2>
<p>The above is Newton-Raphson’s method. It has a simple geometric
interpretation, namely <span class="math notranslate nohighlight">\(x_{n+1}\)</span> is the point where the tangent from
<span class="math notranslate nohighlight">\((x_n,f(x_n))\)</span> crosses the <span class="math notranslate nohighlight">\(x\)</span>-axis. Close to the solution,
Newton-Raphson converges fast to the desired result. However, if we
are far from a root, where the higher-order terms in the series are
important, the Newton-Raphson formula can give grossly inaccurate
results. For instance, the initial guess for the root might be so far
from the true root as to let the search interval include a local
maximum or minimum of the function. If an iteration places a trial
guess near such a local extremum, so that the first derivative nearly
vanishes, then Newton-Raphson may fail totally</p>
</div>
<div class="section" id="extending-to-more-than-one-variable">
<h2>Extending to more than one variable<a class="headerlink" href="#extending-to-more-than-one-variable" title="Permalink to this headline"></a></h2>
<p>Newton’s method can be generalized to systems of several non-linear equations
and variables. Consider the case with two equations</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{array}{cc} f_1(x_1,x_2) &amp;=0\\
f_2(x_1,x_2) &amp;=0,\end{array}
\end{split}\]</div>
<p>which we Taylor expand to obtain</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{array}{cc} 0=f_1(x_1+h_1,x_2+h_2)=&amp;f_1(x_1,x_2)+h_1
\partial f_1/\partial x_1+h_2
\partial f_1/\partial x_2+\dots\\
0=f_2(x_1+h_1,x_2+h_2)=&amp;f_2(x_1,x_2)+h_1
\partial f_2/\partial x_1+h_2
\partial f_2/\partial x_2+\dots
\end{array}.
\end{split}\]</div>
<p>Defining the Jacobian matrix <span class="math notranslate nohighlight">\({\bf \boldsymbol{J}}\)</span> we have</p>
<div class="math notranslate nohighlight">
\[\begin{split}
{\bf \boldsymbol{J}}=\left( \begin{array}{cc}
\partial f_1/\partial x_1 &amp; \partial f_1/\partial x_2 \\
\partial f_2/\partial x_1 &amp;\partial f_2/\partial x_2
\end{array} \right),
\end{split}\]</div>
<p>we can rephrase Newton’s method as</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\left(\begin{array}{c} x_1^{n+1} \\ x_2^{n+1} \end{array} \right)=
\left(\begin{array}{c} x_1^{n} \\ x_2^{n} \end{array} \right)+
\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right),
\end{split}\]</div>
<p>where we have defined</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right)=
-{\bf \boldsymbol{J}}^{-1}
\left(\begin{array}{c} f_1(x_1^{n},x_2^{n}) \\ f_2(x_1^{n},x_2^{n}) \end{array} \right).
\end{split}\]</div>
<p>We need thus to compute the inverse of the Jacobian matrix and it
is to understand that difficulties may
arise in case <span class="math notranslate nohighlight">\({\bf \boldsymbol{J}}\)</span> is nearly singular.</p>
<p>It is rather straightforward to extend the above scheme to systems of
more than two non-linear equations. In our case, the Jacobian matrix is given by the Hessian that represents the second derivative of cost function.</p>
</div>
<div class="section" id="steepest-descent">
<h2>Steepest descent<a class="headerlink" href="#steepest-descent" title="Permalink to this headline"></a></h2>
<p>The basic idea of gradient descent is
that a function <span class="math notranslate nohighlight">\(F(\mathbf{x})\)</span>,
<span class="math notranslate nohighlight">\(\mathbf{x} \equiv (x_1,\cdots,x_n)\)</span>, decreases fastest if one goes from <span class="math notranslate nohighlight">\(\bf {x}\)</span> in the
direction of the negative gradient <span class="math notranslate nohighlight">\(-\nabla F(\mathbf{x})\)</span>.</p>
<p>It can be shown that if</p>
<div class="math notranslate nohighlight">
\[
\mathbf{x}_{k+1} = \mathbf{x}_k - \gamma_k \nabla F(\mathbf{x}_k),
\]</div>
<p>with <span class="math notranslate nohighlight">\(\gamma_k &gt; 0\)</span>.</p>
<p>For <span class="math notranslate nohighlight">\(\gamma_k\)</span> small enough, then <span class="math notranslate nohighlight">\(F(\mathbf{x}_{k+1}) \leq
F(\mathbf{x}_k)\)</span>. This means that for a sufficiently small <span class="math notranslate nohighlight">\(\gamma_k\)</span>
we are always moving towards smaller function values, i.e a minimum.</p>
</div>
<div class="section" id="more-on-steepest-descent">
<h2>More on Steepest descent<a class="headerlink" href="#more-on-steepest-descent" title="Permalink to this headline"></a></h2>
<p>The previous observation is the basis of the method of steepest
descent, which is also referred to as just gradient descent (GD). One
starts with an initial guess <span class="math notranslate nohighlight">\(\mathbf{x}_0\)</span> for a minimum of <span class="math notranslate nohighlight">\(F\)</span> and
computes new approximations according to</p>
<div class="math notranslate nohighlight">
\[
\mathbf{x}_{k+1} = \mathbf{x}_k - \gamma_k \nabla F(\mathbf{x}_k), \ \ k \geq 0.
\]</div>
<p>The parameter <span class="math notranslate nohighlight">\(\gamma_k\)</span> is often referred to as the step length or
the learning rate within the context of Machine Learning.</p>
</div>
<div class="section" id="the-ideal">
<h2>The ideal<a class="headerlink" href="#the-ideal" title="Permalink to this headline"></a></h2>
<p>Ideally the sequence <span class="math notranslate nohighlight">\(\{\mathbf{x}_k \}_{k=0}\)</span> converges to a global
minimum of the function <span class="math notranslate nohighlight">\(F\)</span>. In general we do not know if we are in a
global or local minimum. In the special case when <span class="math notranslate nohighlight">\(F\)</span> is a convex
function, all local minima are also global minima, so in this case
gradient descent can converge to the global solution. The advantage of
this scheme is that it is conceptually simple and straightforward to
implement. However the method in this form has some severe
limitations:</p>
<p>In machine learing we are often faced with non-convex high dimensional
cost functions with many local minima. Since GD is deterministic we
will get stuck in a local minimum, if the method converges, unless we
have a very good intial guess. This also implies that the scheme is
sensitive to the chosen initial condition.</p>
<p>Note that the gradient is a function of <span class="math notranslate nohighlight">\(\mathbf{x} =
(x_1,\cdots,x_n)\)</span> which makes it expensive to compute numerically.</p>
</div>
<div class="section" id="the-sensitiveness-of-the-gradient-descent">
<h2>The sensitiveness of the gradient descent<a class="headerlink" href="#the-sensitiveness-of-the-gradient-descent" title="Permalink to this headline"></a></h2>
<p>The gradient descent method
is sensitive to the choice of learning rate <span class="math notranslate nohighlight">\(\gamma_k\)</span>. This is due
to the fact that we are only guaranteed that <span class="math notranslate nohighlight">\(F(\mathbf{x}_{k+1}) \leq
F(\mathbf{x}_k)\)</span> for sufficiently small <span class="math notranslate nohighlight">\(\gamma_k\)</span>. The problem is to
determine an optimal learning rate. If the learning rate is chosen too
small the method will take a long time to converge and if it is too
large we can experience erratic behavior.</p>
<p>Many of these shortcomings can be alleviated by introducing
randomness. One such method is that of Stochastic Gradient Descent
(SGD), to be discussed next week.</p>
</div>
<div class="section" id="convex-functions">
<h2>Convex functions<a class="headerlink" href="#convex-functions" title="Permalink to this headline"></a></h2>
<p>Ideally we want our cost/loss function to be convex(concave).</p>
<p>First we give the definition of a convex set: A set <span class="math notranslate nohighlight">\(C\)</span> in
<span class="math notranslate nohighlight">\(\mathbb{R}^n\)</span> is said to be convex if, for all <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> in <span class="math notranslate nohighlight">\(C\)</span> and
all <span class="math notranslate nohighlight">\(t \in (0,1)\)</span> , the point <span class="math notranslate nohighlight">\((1 − t)x + ty\)</span> also belongs to
C. Geometrically this means that every point on the line segment
connecting <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> is in <span class="math notranslate nohighlight">\(C\)</span> as discussed below.</p>
<p>The convex subsets of <span class="math notranslate nohighlight">\(\mathbb{R}\)</span> are the intervals of
<span class="math notranslate nohighlight">\(\mathbb{R}\)</span>. Examples of convex sets of <span class="math notranslate nohighlight">\(\mathbb{R}^2\)</span> are the
regular polygons (triangles, rectangles, pentagons, etc…).</p>
</div>
<div class="section" id="convex-function">
<h2>Convex function<a class="headerlink" href="#convex-function" title="Permalink to this headline"></a></h2>
<p><strong>Convex function</strong>: Let <span class="math notranslate nohighlight">\(X \subset \mathbb{R}^n\)</span> be a convex set. Assume that the function <span class="math notranslate nohighlight">\(f: X \rightarrow \mathbb{R}\)</span> is continuous, then <span class="math notranslate nohighlight">\(f\)</span> is said to be convex if $<span class="math notranslate nohighlight">\(f(tx_1 + (1-t)x_2) \leq tf(x_1) + (1-t)f(x_2) \)</span><span class="math notranslate nohighlight">\( for all \)</span>x_1, x_2 \in X<span class="math notranslate nohighlight">\( and for all \)</span>t \in [0,1]<span class="math notranslate nohighlight">\(. If \)</span>\leq<span class="math notranslate nohighlight">\( is replaced with a strict inequaltiy in the definition, we demand \)</span>x_1 \neq x_2<span class="math notranslate nohighlight">\( and \)</span>t\in(0,1)<span class="math notranslate nohighlight">\( then \)</span>f<span class="math notranslate nohighlight">\( is said to be strictly convex. For a single variable function, convexity means that if you draw a straight line connecting \)</span>f(x_1)<span class="math notranslate nohighlight">\( and \)</span>f(x_2)<span class="math notranslate nohighlight">\(, the value of the function on the interval \)</span>[x_1,x_2]$ is always below the line as illustrated below.</p>
</div>
<div class="section" id="conditions-on-convex-functions">
<h2>Conditions on convex functions<a class="headerlink" href="#conditions-on-convex-functions" title="Permalink to this headline"></a></h2>
<p>In the following we state first and second-order conditions which
ensures convexity of a function <span class="math notranslate nohighlight">\(f\)</span>. We write <span class="math notranslate nohighlight">\(D_f\)</span> to denote the
domain of <span class="math notranslate nohighlight">\(f\)</span>, i.e the subset of <span class="math notranslate nohighlight">\(R^n\)</span> where <span class="math notranslate nohighlight">\(f\)</span> is defined. For more
details and proofs we refer to: [S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press](<a class="reference external" href="http://stanford.edu/boyd/cvxbook/">http://stanford.edu/boyd/cvxbook/</a>, 2004).</p>
<p><strong>First order condition.</strong></p>
<p>Suppose <span class="math notranslate nohighlight">\(f\)</span> is differentiable (i.e <span class="math notranslate nohighlight">\(\nabla f(x)\)</span> is well defined for
all <span class="math notranslate nohighlight">\(x\)</span> in the domain of <span class="math notranslate nohighlight">\(f\)</span>). Then <span class="math notranslate nohighlight">\(f\)</span> is convex if and only if <span class="math notranslate nohighlight">\(D_f\)</span>
is a convex set and $<span class="math notranslate nohighlight">\(f(y) \geq f(x) + \nabla f(x)^T (y-x) \)</span><span class="math notranslate nohighlight">\( holds
for all \)</span>x,y \in D_f<span class="math notranslate nohighlight">\(. This condition means that for a convex function
the first order Taylor expansion (right hand side above) at any point
a global under estimator of the function. To convince yourself you can
make a drawing of \)</span>f(x) = x^2+1<span class="math notranslate nohighlight">\( and draw the tangent line to \)</span>f(x)$ and
note that it is always below the graph.</p>
<p><strong>Second order condition.</strong></p>
<p>Assume that <span class="math notranslate nohighlight">\(f\)</span> is twice
differentiable, i.e the Hessian matrix exists at each point in
<span class="math notranslate nohighlight">\(D_f\)</span>. Then <span class="math notranslate nohighlight">\(f\)</span> is convex if and only if <span class="math notranslate nohighlight">\(D_f\)</span> is a convex set and its
Hessian is positive semi-definite for all <span class="math notranslate nohighlight">\(x\in D_f\)</span>. For a
single-variable function this reduces to <span class="math notranslate nohighlight">\(f''(x) \geq 0\)</span>. Geometrically this means that <span class="math notranslate nohighlight">\(f\)</span> has nonnegative curvature
everywhere.</p>
<p>This condition is particularly useful since it gives us an procedure for determining if the function under consideration is convex, apart from using the definition.</p>
</div>
<div class="section" id="more-on-convex-functions">
<h2>More on convex functions<a class="headerlink" href="#more-on-convex-functions" title="Permalink to this headline"></a></h2>
<p>The next result is of great importance to us and the reason why we are
going on about convex functions. In machine learning we frequently
have to minimize a loss/cost function in order to find the best
parameters for the model we are considering.</p>
<p>Ideally we want the
global minimum (for high-dimensional models it is hard to know
if we have local or global minimum). However, if the cost/loss function
is convex the following result provides invaluable information:</p>
<p><strong>Any minimum is global for convex functions.</strong></p>
<p>Consider the problem of finding <span class="math notranslate nohighlight">\(x \in \mathbb{R}^n\)</span> such that <span class="math notranslate nohighlight">\(f(x)\)</span>
is minimal, where <span class="math notranslate nohighlight">\(f\)</span> is convex and differentiable. Then, any point
<span class="math notranslate nohighlight">\(x^*\)</span> that satisfies <span class="math notranslate nohighlight">\(\nabla f(x^*) = 0\)</span> is a global minimum.</p>
<p>This result means that if we know that the cost/loss function is convex and we are able to find a minimum, we are guaranteed that it is a global minimum.</p>
</div>
<div class="section" id="some-simple-problems">
<h2>Some simple problems<a class="headerlink" href="#some-simple-problems" title="Permalink to this headline"></a></h2>
<ol class="simple">
<li><p>Show that <span class="math notranslate nohighlight">\(f(x)=x^2\)</span> is convex for <span class="math notranslate nohighlight">\(x \in \mathbb{R}\)</span> using the definition of convexity. Hint: If you re-write the definition, <span class="math notranslate nohighlight">\(f\)</span> is convex if the following holds for all <span class="math notranslate nohighlight">\(x,y \in D_f\)</span> and any <span class="math notranslate nohighlight">\(\lambda \in [0,1]\)</span> <span class="math notranslate nohighlight">\(\lambda f(x)+(1-\lambda)f(y)-f(\lambda x + (1-\lambda) y ) \geq 0\)</span>.</p></li>
<li><p>Using the second order condition show that the following functions are convex on the specified domain.</p></li>
</ol>
<ul class="simple">
<li><p><span class="math notranslate nohighlight">\(f(x) = e^x\)</span> is convex for <span class="math notranslate nohighlight">\(x \in \mathbb{R}\)</span>.</p></li>
<li><p><span class="math notranslate nohighlight">\(g(x) = -\ln(x)\)</span> is convex for <span class="math notranslate nohighlight">\(x \in (0,\infty)\)</span>.</p></li>
</ul>
<ol class="simple">
<li><p>Let <span class="math notranslate nohighlight">\(f(x) = x^2\)</span> and <span class="math notranslate nohighlight">\(g(x) = e^x\)</span>. Show that <span class="math notranslate nohighlight">\(f(g(x))\)</span> and <span class="math notranslate nohighlight">\(g(f(x))\)</span> is convex for <span class="math notranslate nohighlight">\(x \in \mathbb{R}\)</span>. Also show that if <span class="math notranslate nohighlight">\(f(x)\)</span> is any convex function than <span class="math notranslate nohighlight">\(h(x) = e^{f(x)}\)</span> is convex.</p></li>
<li><p>A norm is any function that satisfy the following properties</p></li>
</ol>
<ul class="simple">
<li><p><span class="math notranslate nohighlight">\(f(\alpha x) = |\alpha| f(x)\)</span> for all <span class="math notranslate nohighlight">\(\alpha \in \mathbb{R}\)</span>.</p></li>
<li><p><span class="math notranslate nohighlight">\(f(x+y) \leq f(x) + f(y)\)</span></p></li>
<li><p><span class="math notranslate nohighlight">\(f(x) \leq 0\)</span> for all <span class="math notranslate nohighlight">\(x \in \mathbb{R}^n\)</span> with equality if and only if <span class="math notranslate nohighlight">\(x = 0\)</span></p></li>
</ul>
<p>Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this).</p>
</div>
<div class="section" id="revisiting-our-first-homework">
<h2>Revisiting our first homework<a class="headerlink" href="#revisiting-our-first-homework" title="Permalink to this headline"></a></h2>
<p>We will use linear regression as a case study for the gradient descent
methods. Linear regression is a great test case for the gradient
descent methods discussed in the lectures since it has several
desirable properties such as:</p>
<ol class="simple">
<li><p>An analytical solution (recall homework set 1).</p></li>
<li><p>The gradient can be computed analytically.</p></li>
<li><p>The cost function is convex which guarantees that gradient descent converges for small enough learning rates</p></li>
</ol>
<p>We revisit an example similar to what we had in the first homework set. We had a function of the type</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="n">x</span> <span class="o">=</span> <span class="mi">2</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">m</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">y</span> <span class="o">=</span> <span class="mi">4</span><span class="o">+</span><span class="mi">3</span><span class="o">*</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">randn</span><span class="p">(</span><span class="n">m</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>
<p>with <span class="math notranslate nohighlight">\(x_i \in [0,1] \)</span> is chosen randomly using a uniform distribution. Additionally we have a stochastic noise chosen according to a normal distribution <span class="math notranslate nohighlight">\(\cal {N}(0,1)\)</span>.
The linear regression model is given by</p>
<div class="math notranslate nohighlight">
\[
h_\beta(x) = \boldsymbol{y} = \beta_0 + \beta_1 x,
\]</div>
<p>such that</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{y}_i = \beta_0 + \beta_1 x_i.
\]</div>
</div>
<div class="section" id="gradient-descent-example">
<h2>Gradient descent example<a class="headerlink" href="#gradient-descent-example" title="Permalink to this headline"></a></h2>
<p>Let <span class="math notranslate nohighlight">\(\mathbf{y} = (y_1,\cdots,y_n)^T\)</span>, <span class="math notranslate nohighlight">\(\mathbf{\boldsymbol{y}} = (\boldsymbol{y}_1,\cdots,\boldsymbol{y}_n)^T\)</span> and <span class="math notranslate nohighlight">\(\beta = (\beta_0, \beta_1)^T\)</span></p>
<p>It is convenient to write <span class="math notranslate nohighlight">\(\mathbf{\boldsymbol{y}} = X\beta\)</span> where <span class="math notranslate nohighlight">\(X \in \mathbb{R}^{100 \times 2} \)</span> is the design matrix given by (we keep the intercept here)</p>
<div class="math notranslate nohighlight">
\[\begin{split}
X \equiv \begin{bmatrix}
1 &amp; x_1 \\
\vdots &amp; \vdots \\
1 &amp; x_{100} &amp; \\
\end{bmatrix}.
\end{split}\]</div>
<p>The cost/loss/risk function is given by (</p>
<div class="math notranslate nohighlight">
\[
C(\beta) = \frac{1}{n}||X\beta-\mathbf{y}||_{2}^{2} = \frac{1}{n}\sum_{i=1}^{100}\left[ (\beta_0 + \beta_1 x_i)^2 - 2 y_i (\beta_0 + \beta_1 x_i) + y_i^2\right]
\]</div>
<p>and we want to find <span class="math notranslate nohighlight">\(\beta\)</span> such that <span class="math notranslate nohighlight">\(C(\beta)\)</span> is minimized.</p>
</div>
<div class="section" id="the-derivative-of-the-cost-loss-function">
<h2>The derivative of the cost/loss function<a class="headerlink" href="#the-derivative-of-the-cost-loss-function" title="Permalink to this headline"></a></h2>
<p>Computing <span class="math notranslate nohighlight">\(\partial C(\beta) / \partial \beta_0\)</span> and <span class="math notranslate nohighlight">\(\partial C(\beta) / \partial \beta_1\)</span> we can show that the gradient can be written as</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\nabla_{\beta} C(\beta) = \frac{2}{n}\begin{bmatrix} \sum_{i=1}^{100} \left(\beta_0+\beta_1x_i-y_i\right) \\
\sum_{i=1}^{100}\left( x_i (\beta_0+\beta_1x_i)-y_ix_i\right) \\
\end{bmatrix} = \frac{2}{n}X^T(X\beta - \mathbf{y}),
\end{split}\]</div>
<p>where <span class="math notranslate nohighlight">\(X\)</span> is the design matrix defined above.</p>
</div>
<div class="section" id="the-hessian-matrix">
<h2>The Hessian matrix<a class="headerlink" href="#the-hessian-matrix" title="Permalink to this headline"></a></h2>
<p>The Hessian matrix of <span class="math notranslate nohighlight">\(C(\beta)\)</span> is given by</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{H} \equiv \begin{bmatrix}
\frac{\partial^2 C(\beta)}{\partial \beta_0^2} &amp; \frac{\partial^2 C(\beta)}{\partial \beta_0 \partial \beta_1} \\
\frac{\partial^2 C(\beta)}{\partial \beta_0 \partial \beta_1} &amp; \frac{\partial^2 C(\beta)}{\partial \beta_1^2} &amp; \\
\end{bmatrix} = \frac{2}{n}X^T X.
\end{split}\]</div>
<p>This result implies that <span class="math notranslate nohighlight">\(C(\beta)\)</span> is a convex function since the matrix <span class="math notranslate nohighlight">\(X^T X\)</span> always is positive semi-definite.</p>
</div>
<div class="section" id="simple-program">
<h2>Simple program<a class="headerlink" href="#simple-program" title="Permalink to this headline"></a></h2>
<p>We can now write a program that minimizes <span class="math notranslate nohighlight">\(C(\beta)\)</span> using the gradient descent method with a constant learning rate <span class="math notranslate nohighlight">\(\gamma\)</span> according to</p>
<div class="math notranslate nohighlight">
\[
\beta_{k+1} = \beta_k - \gamma \nabla_\beta C(\beta_k), \ k=0,1,\cdots
\]</div>
<p>We can use the expression we computed for the gradient and let use a
<span class="math notranslate nohighlight">\(\beta_0\)</span> be chosen randomly and let <span class="math notranslate nohighlight">\(\gamma = 0.001\)</span>. Stop iterating
when <span class="math notranslate nohighlight">\(||\nabla_\beta C(\beta_k) || \leq \epsilon = 10^{-8}\)</span>. <strong>Note that the code below does not include the latter stop criterion</strong>.</p>
<p>And finally we can compare our solution for <span class="math notranslate nohighlight">\(\beta\)</span> with the analytic result given by
<span class="math notranslate nohighlight">\(\beta= (X^TX)^{-1} X^T \mathbf{y}\)</span>.</p>
</div>
<div class="section" id="id1">
<h2>Gradient Descent Example<a class="headerlink" href="#id1" title="Permalink to this headline"></a></h2>
<p>Here our simple example</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="c1"># Importing various packages</span>
<span class="kn">from</span> <span class="nn">random</span> <span class="kn">import</span> <span class="n">random</span><span class="p">,</span> <span class="n">seed</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">mpl_toolkits.mplot3d</span> <span class="kn">import</span> <span class="n">Axes3D</span>
<span class="kn">from</span> <span class="nn">matplotlib</span> <span class="kn">import</span> <span class="n">cm</span>
<span class="kn">from</span> <span class="nn">matplotlib.ticker</span> <span class="kn">import</span> <span class="n">LinearLocator</span><span class="p">,</span> <span class="n">FormatStrFormatter</span>
<span class="kn">import</span> <span class="nn">sys</span>
<span class="c1"># the number of datapoints</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="mi">2</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="mi">1</span><span class="p">)</span>
<span class="n">y</span> <span class="o">=</span> <span class="mi">4</span><span class="o">+</span><span class="mi">3</span><span class="o">*</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">randn</span><span class="p">(</span><span class="n">n</span><span class="p">,</span><span class="mi">1</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">c_</span><span class="p">[</span><span class="n">np</span><span class="o">.</span><span class="n">ones</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="mi">1</span><span class="p">)),</span> <span class="n">x</span><span class="p">]</span>
<span class="c1"># Hessian matrix</span>
<span class="n">H</span> <span class="o">=</span> <span class="p">(</span><span class="mf">2.0</span><span class="o">/</span><span class="n">n</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">X</span>
<span class="c1"># Get the eigenvalues</span>
<span class="n">EigValues</span><span class="p">,</span> <span class="n">EigVectors</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">eig</span><span class="p">(</span><span class="n">H</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;Eigenvalues of Hessian Matrix:</span><span class="si">{</span><span class="n">EigValues</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="n">beta_linreg</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">beta_linreg</span><span class="p">)</span>
<span class="n">beta</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">2</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">eta</span> <span class="o">=</span> <span class="mf">1.0</span><span class="o">/</span><span class="n">np</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">EigValues</span><span class="p">)</span>
<span class="n">Niterations</span> <span class="o">=</span> <span class="mi">1000</span>
<span class="k">for</span> <span class="nb">iter</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">Niterations</span><span class="p">):</span>
<span class="n">gradient</span> <span class="o">=</span> <span class="p">(</span><span class="mf">2.0</span><span class="o">/</span><span class="n">n</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="p">(</span><span class="n">X</span> <span class="o">@</span> <span class="n">beta</span><span class="o">-</span><span class="n">y</span><span class="p">)</span>
<span class="n">beta</span> <span class="o">-=</span> <span class="n">eta</span><span class="o">*</span><span class="n">gradient</span>
<span class="nb">print</span><span class="p">(</span><span class="n">beta</span><span class="p">)</span>
<span class="n">xnew</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="mi">0</span><span class="p">],[</span><span class="mi">2</span><span class="p">]])</span>
<span class="n">xbnew</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">c_</span><span class="p">[</span><span class="n">np</span><span class="o">.</span><span class="n">ones</span><span class="p">((</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">)),</span> <span class="n">xnew</span><span class="p">]</span>
<span class="n">ypredict</span> <span class="o">=</span> <span class="n">xbnew</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">beta</span><span class="p">)</span>
<span class="n">ypredict2</span> <span class="o">=</span> <span class="n">xbnew</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">beta_linreg</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">xnew</span><span class="p">,</span> <span class="n">ypredict</span><span class="p">,</span> <span class="s2">&quot;r-&quot;</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">xnew</span><span class="p">,</span> <span class="n">ypredict2</span><span class="p">,</span> <span class="s2">&quot;b-&quot;</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">y</span> <span class="p">,</span><span class="s1">&#39;ro&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">axis</span><span class="p">([</span><span class="mi">0</span><span class="p">,</span><span class="mf">2.0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span> <span class="mf">15.0</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="sa">r</span><span class="s1">&#39;$x$&#39;</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="sa">r</span><span class="s1">&#39;$y$&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="sa">r</span><span class="s1">&#39;Gradient descent example&#39;</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="and-a-corresponding-example-using-scikit-learn">
<h2>And a corresponding example using <strong>scikit-learn</strong><a class="headerlink" href="#and-a-corresponding-example-using-scikit-learn" 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="c1"># Importing various packages</span>
<span class="kn">from</span> <span class="nn">random</span> <span class="kn">import</span> <span class="n">random</span><span class="p">,</span> <span class="n">seed</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.linear_model</span> <span class="kn">import</span> <span class="n">SGDRegressor</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="mi">2</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="mi">1</span><span class="p">)</span>
<span class="n">y</span> <span class="o">=</span> <span class="mi">4</span><span class="o">+</span><span class="mi">3</span><span class="o">*</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">randn</span><span class="p">(</span><span class="n">n</span><span class="p">,</span><span class="mi">1</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">c_</span><span class="p">[</span><span class="n">np</span><span class="o">.</span><span class="n">ones</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="mi">1</span><span class="p">)),</span> <span class="n">x</span><span class="p">]</span>
<span class="n">beta_linreg</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="p">(</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="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">beta_linreg</span><span class="p">)</span>
<span class="n">sgdreg</span> <span class="o">=</span> <span class="n">SGDRegressor</span><span class="p">(</span><span class="n">max_iter</span> <span class="o">=</span> <span class="mi">50</span><span class="p">,</span> <span class="n">penalty</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">eta0</span><span class="o">=</span><span class="mf">0.1</span><span class="p">)</span>
<span class="n">sgdreg</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="o">.</span><span class="n">ravel</span><span class="p">())</span>
<span class="nb">print</span><span class="p">(</span><span class="n">sgdreg</span><span class="o">.</span><span class="n">intercept_</span><span class="p">,</span> <span class="n">sgdreg</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="section" id="gradient-descent-and-ridge">
<h2>Gradient descent and Ridge<a class="headerlink" href="#gradient-descent-and-ridge" title="Permalink to this headline"></a></h2>
<p>We have also discussed Ridge regression where the loss function contains a regularized term given by the <span class="math notranslate nohighlight">\(L_2\)</span> norm of <span class="math notranslate nohighlight">\(\beta\)</span>,</p>
<div class="math notranslate nohighlight">
\[
C_{\text{ridge}}(\beta) = \frac{1}{n}||X\beta -\mathbf{y}||^2 + \lambda ||\beta||^2, \ \lambda \geq 0.
\]</div>
<p>In order to minimize <span class="math notranslate nohighlight">\(C_{\text{ridge}}(\beta)\)</span> using GD we adjust the gradient as follows</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\nabla_\beta C_{\text{ridge}}(\beta) = \frac{2}{n}\begin{bmatrix} \sum_{i=1}^{100} \left(\beta_0+\beta_1x_i-y_i\right) \\
\sum_{i=1}^{100}\left( x_i (\beta_0+\beta_1x_i)-y_ix_i\right) \\
\end{bmatrix} + 2\lambda\begin{bmatrix} \beta_0 \\ \beta_1\end{bmatrix} = 2 (\frac{1}{n}X^T(X\beta - \mathbf{y})+\lambda \beta).
\end{split}\]</div>
<p>We can easily extend our program to minimize <span class="math notranslate nohighlight">\(C_{\text{ridge}}(\beta)\)</span> using gradient descent and compare with the analytical solution given by</p>
<div class="math notranslate nohighlight">
\[
\beta_{\text{ridge}} = \left(X^T X + n\lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y}.
\]</div>
</div>
<div class="section" id="the-hessian-matrix-for-ridge-regression">
<h2>The Hessian matrix for Ridge Regression<a class="headerlink" href="#the-hessian-matrix-for-ridge-regression" title="Permalink to this headline"></a></h2>
<p>The Hessian matrix of Ridge Regression for our simple example is given by</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\boldsymbol{H} \equiv \begin{bmatrix}
\frac{\partial^2 C(\beta)}{\partial \beta_0^2} &amp; \frac{\partial^2 C(\beta)}{\partial \beta_0 \partial \beta_1} \\
\frac{\partial^2 C(\beta)}{\partial \beta_0 \partial \beta_1} &amp; \frac{\partial^2 C(\beta)}{\partial \beta_1^2} &amp; \\
\end{bmatrix} = \frac{2}{n}X^T X+2\lambda\boldsymbol{I}.
\end{split}\]</div>
<p>This implies that the Hessian matrix is positive definite, hence the stationary point is a
minimum.
Note that the Ridge cost function is convex being a sum of two convex
functions. Therefore, the stationary point is a global
minimum of this function.</p>
</div>
<div class="section" id="program-example-for-gradient-descent-with-ridge-regression">
<h2>Program example for gradient descent with Ridge Regression<a class="headerlink" href="#program-example-for-gradient-descent-with-ridge-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">from</span> <span class="nn">random</span> <span class="kn">import</span> <span class="n">random</span><span class="p">,</span> <span class="n">seed</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">mpl_toolkits.mplot3d</span> <span class="kn">import</span> <span class="n">Axes3D</span>
<span class="kn">from</span> <span class="nn">matplotlib</span> <span class="kn">import</span> <span class="n">cm</span>
<span class="kn">from</span> <span class="nn">matplotlib.ticker</span> <span class="kn">import</span> <span class="n">LinearLocator</span><span class="p">,</span> <span class="n">FormatStrFormatter</span>
<span class="kn">import</span> <span class="nn">sys</span>
<span class="c1"># the number of datapoints</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="mi">2</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="mi">1</span><span class="p">)</span>
<span class="n">y</span> <span class="o">=</span> <span class="mi">4</span><span class="o">+</span><span class="mi">3</span><span class="o">*</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">randn</span><span class="p">(</span><span class="n">n</span><span class="p">,</span><span class="mi">1</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">c_</span><span class="p">[</span><span class="n">np</span><span class="o">.</span><span class="n">ones</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="mi">1</span><span class="p">)),</span> <span class="n">x</span><span class="p">]</span>
<span class="n">XT_X</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">X</span>
<span class="c1">#Ridge parameter lambda</span>
<span class="n">lmbda</span> <span class="o">=</span> <span class="mf">0.001</span>
<span class="n">Id</span> <span class="o">=</span> <span class="n">n</span><span class="o">*</span><span class="n">lmbda</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">XT_X</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span>
<span class="c1"># Hessian matrix</span>
<span class="n">H</span> <span class="o">=</span> <span class="p">(</span><span class="mf">2.0</span><span class="o">/</span><span class="n">n</span><span class="p">)</span><span class="o">*</span> <span class="n">XT_X</span><span class="o">+</span><span class="mi">2</span><span class="o">*</span><span class="n">lmbda</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">XT_X</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span>
<span class="c1"># Get the eigenvalues</span>
<span class="n">EigValues</span><span class="p">,</span> <span class="n">EigVectors</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">eig</span><span class="p">(</span><span class="n">H</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;Eigenvalues of Hessian Matrix:</span><span class="si">{</span><span class="n">EigValues</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="n">beta_linreg</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">XT_X</span><span class="o">+</span><span class="n">Id</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">beta_linreg</span><span class="p">)</span>
<span class="c1"># Start plain gradient descent</span>
<span class="n">beta</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">2</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">eta</span> <span class="o">=</span> <span class="mf">1.0</span><span class="o">/</span><span class="n">np</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">EigValues</span><span class="p">)</span>
<span class="n">Niterations</span> <span class="o">=</span> <span class="mi">100</span>
<span class="k">for</span> <span class="nb">iter</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">Niterations</span><span class="p">):</span>
<span class="n">gradients</span> <span class="o">=</span> <span class="mf">2.0</span><span class="o">/</span><span class="n">n</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="p">(</span><span class="n">X</span> <span class="o">@</span> <span class="p">(</span><span class="n">beta</span><span class="p">)</span><span class="o">-</span><span class="n">y</span><span class="p">)</span><span class="o">+</span><span class="mi">2</span><span class="o">*</span><span class="n">lmbda</span><span class="o">*</span><span class="n">beta</span>
<span class="n">beta</span> <span class="o">-=</span> <span class="n">eta</span><span class="o">*</span><span class="n">gradients</span>
<span class="nb">print</span><span class="p">(</span><span class="n">beta</span><span class="p">)</span>
<span class="n">ypredict</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta</span>
<span class="n">ypredict2</span> <span class="o">=</span> <span class="n">X</span> <span class="o">@</span> <span class="n">beta_linreg</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">ypredict</span><span class="p">,</span> <span class="s2">&quot;r-&quot;</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">ypredict2</span><span class="p">,</span> <span class="s2">&quot;b-&quot;</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">y</span> <span class="p">,</span><span class="s1">&#39;ro&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">axis</span><span class="p">([</span><span class="mi">0</span><span class="p">,</span><span class="mf">2.0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span> <span class="mf">15.0</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="sa">r</span><span class="s1">&#39;$x$&#39;</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="sa">r</span><span class="s1">&#39;$y$&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="sa">r</span><span class="s1">&#39;Gradient descent example for Ridge&#39;</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="using-gradient-descent-methods-limitations">
<h2>Using gradient descent methods, limitations<a class="headerlink" href="#using-gradient-descent-methods-limitations" title="Permalink to this headline"></a></h2>
<ul class="simple">
<li><p><strong>Gradient descent (GD) finds local minima of our function</strong>. Since the GD algorithm is deterministic, if it converges, it will converge to a local minimum of our cost/loss/risk function. Because in ML we are often dealing with extremely rugged landscapes with many local minima, this can lead to poor performance.</p></li>
<li><p><strong>GD is sensitive to initial conditions</strong>. One consequence of the local nature of GD is that initial conditions matter. Depending on where one starts, one will end up at a different local minima. Therefore, it is very important to think about how one initializes the training process. This is true for GD as well as more complicated variants of GD.</p></li>
<li><p><strong>Gradients are computationally expensive to calculate for large datasets</strong>. In many cases in statistics and ML, the cost/loss/risk function is a sum of terms, with one term for each data point. For example, in linear regression, <span class="math notranslate nohighlight">\(E \propto \sum_{i=1}^n (y_i - \mathbf{w}^T\cdot\mathbf{x}_i)^2\)</span>; for logistic regression, the square error is replaced by the cross entropy. To calculate the gradient we have to sum over <em>all</em> <span class="math notranslate nohighlight">\(n\)</span> data points. Doing this at every GD step becomes extremely computationally expensive. An ingenious solution to this, is to calculate the gradients using small subsets of the data called “mini batches”. This has the added benefit of introducing stochasticity into our algorithm.</p></li>
<li><p><strong>GD is very sensitive to choices of learning rates</strong>. GD is extremely sensitive to the choice of learning rates. If the learning rate is very small, the training process take an extremely long time. For larger learning rates, GD can diverge and give poor results. Furthermore, depending on what the local landscape looks like, we have to modify the learning rates to ensure convergence. Ideally, we would <em>adaptively</em> choose the learning rates to match the landscape.</p></li>
<li><p><strong>GD treats all directions in parameter space uniformly.</strong> Another major drawback of GD is that unlike Newton’s method, the learning rate for GD is the same in all directions in parameter space. For this reason, the maximum learning rate is set by the behavior of the steepest direction and this can significantly slow down training. Ideally, we would like to take large steps in flat directions and small steps in steep directions. Since we are exploring rugged landscapes where curvatures change, this requires us to keep track of not only the gradient but second derivatives. The ideal scenario would be to calculate the Hessian but this proves to be too computationally expensive.</p></li>
<li><p>GD can take exponential time to escape saddle points, even with random initialization. As we mentioned, GD is extremely sensitive to initial condition since it determines the particular local minimum GD would eventually reach. However, even with a good initialization scheme, through the introduction of randomness, GD can still take exponential time to escape saddle points.</p></li>
</ul>
</div>
<div class="section" id="challenge-yourself-the-coming-weekend">
<h2>Challenge yourself the coming weekend<a class="headerlink" href="#challenge-yourself-the-coming-weekend" title="Permalink to this headline"></a></h2>
<p>Write a code which implements gradient descent for a logistic regression example.</p>
</div>
<div class="section" id="lab-session-material-from-last-week-and-relevant-for-the-first-project">
<h2>Lab session: Material from last week and relevant for the first project<a class="headerlink" href="#lab-session-material-from-last-week-and-relevant-for-the-first-project" title="Permalink to this headline"></a></h2>
</div>
<div class="section" id="various-steps-in-cross-validation">
<h2>Various steps in cross-validation<a class="headerlink" href="#various-steps-in-cross-validation" title="Permalink to this headline"></a></h2>
<p>When the repetitive splitting of the data set is done randomly,
samples may accidently end up in a fast majority of the splits in
either training or test set. Such samples may have an unbalanced
influence on either model building or prediction evaluation. To avoid
this <span class="math notranslate nohighlight">\(k\)</span>-fold cross-validation structures the data splitting. The
samples are divided into <span class="math notranslate nohighlight">\(k\)</span> more or less equally sized exhaustive and
mutually exclusive subsets. In turn (at each split) one of these
subsets plays the role of the test set while the union of the
remaining subsets constitutes the training set. Such a splitting
warrants a balanced representation of each sample in both training and
test set over the splits. Still the division into the <span class="math notranslate nohighlight">\(k\)</span> subsets
involves a degree of randomness. This may be fully excluded when
choosing <span class="math notranslate nohighlight">\(k=n\)</span>. This particular case is referred to as leave-one-out
cross-validation (LOOCV).</p>
</div>
<div class="section" id="how-to-set-up-the-cross-validation-for-ridge-and-or-lasso">
<h2>How to set up the cross-validation for Ridge and/or Lasso<a class="headerlink" href="#how-to-set-up-the-cross-validation-for-ridge-and-or-lasso" title="Permalink to this headline"></a></h2>
<ul class="simple">
<li><p>Define a range of interest for the penalty parameter.</p></li>
<li><p>Divide the data set into training and test set comprising samples <span class="math notranslate nohighlight">\(\{1, \ldots, n\} \setminus i\)</span> and <span class="math notranslate nohighlight">\(\{ i \}\)</span>, respectively.</p></li>
<li><p>Fit the linear regression model by means of for example Ridge or Lasso regression for each <span class="math notranslate nohighlight">\(\lambda\)</span> in the grid using the training set, and the corresponding estimate of the error variance <span class="math notranslate nohighlight">\(\boldsymbol{\sigma}_{-i}^2(\lambda)\)</span>, as</p></li>
</ul>
<div class="math notranslate nohighlight">
\[
\begin{align*}
\boldsymbol{\beta}_{-i}(\lambda) &amp; = ( \boldsymbol{X}_{-i, \ast}^{T}
\boldsymbol{X}_{-i, \ast} + \lambda \boldsymbol{I}_{pp})^{-1}
\boldsymbol{X}_{-i, \ast}^{T} \boldsymbol{y}_{-i}
\end{align*}
\]</div>
<ul class="simple">
<li><p>Evaluate the prediction performance of these models on the test set by <span class="math notranslate nohighlight">\(C[y_i, \boldsymbol{X}_{i, \ast}; \boldsymbol{\beta}_{-i}(\lambda), \boldsymbol{\sigma}_{-i}^2(\lambda)]\)</span>. Or, by the prediction error <span class="math notranslate nohighlight">\(|y_i - \boldsymbol{X}_{i, \ast} \boldsymbol{\beta}_{-i}(\lambda)|\)</span>, the relative error, the error squared or the R2 score function.</p></li>
<li><p>Repeat the first three steps such that each sample plays the role of the test set once.</p></li>
<li><p>Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data.</p></li>
</ul>
</div>
<div class="section" id="cross-validation-in-brief">
<h2>Cross-validation in brief<a class="headerlink" href="#cross-validation-in-brief" title="Permalink to this headline"></a></h2>
<p>For the various values of <span class="math notranslate nohighlight">\(k\)</span></p>
<ol class="simple">
<li><p>shuffle the dataset randomly.</p></li>
<li><p>Split the dataset into <span class="math notranslate nohighlight">\(k\)</span> groups.</p></li>
<li><p>For each unique group:</p></li>
</ol>
<p>a. Decide which group to use as set for test data</p>
<p>b. Take the remaining groups as a training data set</p>
<p>c. Fit a model on the training set and evaluate it on the test set</p>
<p>d. Retain the evaluation score and discard the model</p>
<ol class="simple">
<li><p>Summarize the model using the sample of model evaluation scores</p></li>
</ol>
</div>
<div class="section" id="code-example-for-cross-validation-and-k-fold-cross-validation">
<h2>Code Example for Cross-validation and <span class="math notranslate nohighlight">\(k\)</span>-fold Cross-validation<a class="headerlink" href="#code-example-for-cross-validation-and-k-fold-cross-validation" title="Permalink to this headline"></a></h2>
<p>The code here uses Ridge regression with cross-validation (CV) resampling and <span class="math notranslate nohighlight">\(k\)</span>-fold CV in order to fit a specific polynomial.</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">KFold</span>
<span class="kn">from</span> <span class="nn">sklearn.linear_model</span> <span class="kn">import</span> <span class="n">Ridge</span>
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">cross_val_score</span>
<span class="kn">from</span> <span class="nn">sklearn.preprocessing</span> <span class="kn">import</span> <span class="n">PolynomialFeatures</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="c1"># Generate the data.</span>
<span class="n">nsamples</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">randn</span><span class="p">(</span><span class="n">nsamples</span><span class="p">)</span>
<span class="n">y</span> <span class="o">=</span> <span class="mi">3</span><span class="o">*</span><span class="n">x</span><span class="o">**</span><span class="mi">2</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">nsamples</span><span class="p">)</span>
<span class="c1">## Cross-validation on Ridge regression using KFold only</span>
<span class="c1"># Decide degree on polynomial to fit</span>
<span class="n">poly</span> <span class="o">=</span> <span class="n">PolynomialFeatures</span><span class="p">(</span><span class="n">degree</span> <span class="o">=</span> <span class="mi">6</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">500</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">5</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
<span class="c1"># Initialize a KFold instance</span>
<span class="n">k</span> <span class="o">=</span> <span class="mi">5</span>
<span class="n">kfold</span> <span class="o">=</span> <span class="n">KFold</span><span class="p">(</span><span class="n">n_splits</span> <span class="o">=</span> <span class="n">k</span><span class="p">)</span>
<span class="c1"># Perform the cross-validation to estimate MSE</span>
<span class="n">scores_KFold</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">k</span><span class="p">))</span>
<span class="n">i</span> <span class="o">=</span> <span class="mi">0</span>
<span class="k">for</span> <span class="n">lmb</span> <span class="ow">in</span> <span class="n">lambdas</span><span class="p">:</span>
<span class="n">ridge</span> <span class="o">=</span> <span class="n">Ridge</span><span class="p">(</span><span class="n">alpha</span> <span class="o">=</span> <span class="n">lmb</span><span class="p">)</span>
<span class="n">j</span> <span class="o">=</span> <span class="mi">0</span>
<span class="k">for</span> <span class="n">train_inds</span><span class="p">,</span> <span class="n">test_inds</span> <span class="ow">in</span> <span class="n">kfold</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="n">x</span><span class="p">):</span>
<span class="n">xtrain</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">train_inds</span><span class="p">]</span>
<span class="n">ytrain</span> <span class="o">=</span> <span class="n">y</span><span class="p">[</span><span class="n">train_inds</span><span class="p">]</span>
<span class="n">xtest</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">test_inds</span><span class="p">]</span>
<span class="n">ytest</span> <span class="o">=</span> <span class="n">y</span><span class="p">[</span><span class="n">test_inds</span><span class="p">]</span>
<span class="n">Xtrain</span> <span class="o">=</span> <span class="n">poly</span><span class="o">.</span><span class="n">fit_transform</span><span class="p">(</span><span class="n">xtrain</span><span class="p">[:,</span> <span class="n">np</span><span class="o">.</span><span class="n">newaxis</span><span class="p">])</span>
<span class="n">ridge</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">Xtrain</span><span class="p">,</span> <span class="n">ytrain</span><span class="p">[:,</span> <span class="n">np</span><span class="o">.</span><span class="n">newaxis</span><span class="p">])</span>
<span class="n">Xtest</span> <span class="o">=</span> <span class="n">poly</span><span class="o">.</span><span class="n">fit_transform</span><span class="p">(</span><span class="n">xtest</span><span class="p">[:,</span> <span class="n">np</span><span class="o">.</span><span class="n">newaxis</span><span class="p">])</span>
<span class="n">ypred</span> <span class="o">=</span> <span class="n">ridge</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">Xtest</span><span class="p">)</span>
<span class="n">scores_KFold</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</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">ypred</span> <span class="o">-</span> <span class="n">ytest</span><span class="p">[:,</span> <span class="n">np</span><span class="o">.</span><span class="n">newaxis</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">size</span><span class="p">(</span><span class="n">ypred</span><span class="p">)</span>
<span class="n">j</span> <span class="o">+=</span> <span class="mi">1</span>
<span class="n">i</span> <span class="o">+=</span> <span class="mi">1</span>
<span class="n">estimated_mse_KFold</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">scores_KFold</span><span class="p">,</span> <span class="n">axis</span> <span class="o">=</span> <span class="mi">1</span><span class="p">)</span>
<span class="c1">## Cross-validation using cross_val_score from sklearn along with KFold</span>
<span class="c1"># kfold is an instance initialized above as:</span>
<span class="c1"># kfold = KFold(n_splits = k)</span>
<span class="n">estimated_mse_sklearn</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">i</span> <span class="o">=</span> <span class="mi">0</span>
<span class="k">for</span> <span class="n">lmb</span> <span class="ow">in</span> <span class="n">lambdas</span><span class="p">:</span>
<span class="n">ridge</span> <span class="o">=</span> <span class="n">Ridge</span><span class="p">(</span><span class="n">alpha</span> <span class="o">=</span> <span class="n">lmb</span><span class="p">)</span>
<span class="n">X</span> <span class="o">=</span> <span class="n">poly</span><span class="o">.</span><span class="n">fit_transform</span><span class="p">(</span><span class="n">x</span><span class="p">[:,</span> <span class="n">np</span><span class="o">.</span><span class="n">newaxis</span><span class="p">])</span>
<span class="n">estimated_mse_folds</span> <span class="o">=</span> <span class="n">cross_val_score</span><span class="p">(</span><span class="n">ridge</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">np</span><span class="o">.</span><span class="n">newaxis</span><span class="p">],</span> <span class="n">scoring</span><span class="o">=</span><span class="s1">&#39;neg_mean_squared_error&#39;</span><span class="p">,</span> <span class="n">cv</span><span class="o">=</span><span class="n">kfold</span><span class="p">)</span>
<span class="c1"># cross_val_score return an array containing the estimated negative mse for every fold.</span>
<span class="c1"># we have to the the mean of every array in order to get an estimate of the mse of the model</span>
<span class="n">estimated_mse_sklearn</span><span class="p">[</span><span class="n">i</span><span class="p">]</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="o">-</span><span class="n">estimated_mse_folds</span><span class="p">)</span>
<span class="n">i</span> <span class="o">+=</span> <span class="mi">1</span>
<span class="c1">## Plot and compare the slightly different ways to perform cross-validation</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">estimated_mse_sklearn</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">&#39;cross_val_score&#39;</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">estimated_mse_KFold</span><span class="p">,</span> <span class="s1">&#39;r--&#39;</span><span class="p">,</span> <span class="n">label</span> <span class="o">=</span> <span class="s1">&#39;KFold&#39;</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">&#39;log10(lambda)&#39;</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">&#39;mse&#39;</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="exercisesweek38.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 38</p>
</div>
</a>
<a class='right-next' id="next-link" href="exercisesweek39.html" title="next page">
<div class="prev-next-info">
<p class="prev-next-subtitle">next</p>
<p class="prev-next-title">Exercises week 39</p>
</div>
<i class="fas fa-angle-right"></i>
</a>
</div>
</div>
</div>
<footer class="footer">
<p>
By Morten Hjorth-Jensen<br/>
&copy; Copyright 2021.<br/>
</p>
</footer>
</main>
</div>
</div>
<script src="_static/js/index.be7d3bbb2ef33a8344ce.js"></script>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink