Statistics Machine Learning Python Draft
Statistics Machine Learning Python Draft
Statistics Machine Learning Python Draft
Python
Release 0.3 beta
1 Introduction 1
1.1 Python ecosystem for data-science . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Introduction to Machine Learning . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Data analysis methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Python language 9
2.1 Import libraries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Basic operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Data types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 Execution control statements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.5 Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.6 List comprehensions, iterators, etc. . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.7 Regular expression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.8 System programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.9 Scripts and argument parsing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.10 Networking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.11 Modules and packages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.12 Object Oriented Programming (OOP) . . . . . . . . . . . . . . . . . . . . . . . . 32
2.13 Style guide for Python programming . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.14 Documenting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.15 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3 Scientific Python 37
3.1 Numpy: arrays and matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2 Pandas: data manipulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.3 Matplotlib: data visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4 Statistics 75
4.1 Univariate statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.2 Lab 1: Brain volumes study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
4.3 Multivariate statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
4.4 Time Series in python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
i
5.6 Resampling Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219
5.7 Ensemble learning: bagging, boosting and stacking . . . . . . . . . . . . . . . . . 232
5.8 Gradient descent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 247
ii
CHAPTER
ONE
INTRODUCTION
• Interpreted
• Garbage collector (do not prevent from memory leak)
• Dynamically-typed language (Java is statically typed)
1.1.2 Anaconda
Anaconda is a python distribution that ships most of python tools and libraries
Installation
1. Download anaconda (Python 3.x) http://continuum.io/downloads
2. Install it, on Linux
bash Anaconda3-2.4.1-Linux-x86_64.sh
export PATH="${HOME}/anaconda3/bin:$PATH"
Install additional packages. Those commands install qt back-end (Fix a temporary issue to run
spyder)
1
Statistics and Machine Learning in Python, Release 0.3 beta
conda list
Environments
• A conda environment is a directory that contains a specific collection of conda packages
that you have installed.
• Control packages environment for a specific purpose: collaborating with someone else,
delivering an application to your client,
• Switch between environments
List of all environments
:: conda info –envs
1. Create new environment
2. Activate
3. Install new package
Miniconda
Anaconda without the collection of (>700) packages. With Miniconda you download only the
packages you want with the conda command: conda install PACKAGENAME
1. Download anaconda (Python 3.x) https://conda.io/miniconda.html
2. Install it, on Linux
bash Miniconda3-latest-Linux-x86_64.sh
export PATH=${HOME}/miniconda3/bin:$PATH
2 Chapter 1. Introduction
Statistics and Machine Learning in Python, Release 0.3 beta
1.1.3 Commands
python: python interpreter. On the dos/unix command line execute wholes file:
python file.py
Interactive mode:
python
ipython
For neuroimaging:
1.1.4 Libraries
scipy.org: https://www.scipy.org/docs.html
Numpy: Basic numerical operation. Matrix operation plus some basic solvers.:
import numpy as np
X = np.array([[1, 2], [3, 4]])
#v = np.array([1, 2]).reshape((2, 1))
v = np.array([1, 2])
np.dot(X, v) # no broadcasting
X * v # broadcasting
np.dot(v, X)
X - X.mean(axis=0)
import scipy
import scipy.linalg
scipy.linalg.svd(X, full_matrices=False)
Matplotlib: visualization:
import numpy as np
import matplotlib.pyplot as plt
#%matplotlib qt
x = np.linspace(0, 10, 50)
sinus = np.sin(x)
plt.plot(x, sinus)
plt.show()
4 Chapter 1. Introduction
Statistics and Machine Learning in Python, Release 0.3 beta
• Linear model.
• Non parametric statistics.
• Linear algebra: matrix operations, inversion, eigenvalues.
6 Chapter 1. Introduction
Statistics and Machine Learning in Python, Release 0.3 beta
6. (i) Sanity check (basic descriptive statistics); (ii) data cleaning (impute missing data,
recoding); Final Quality Control (QC) perform descriptive statistics and think ! (re-
move possible confounding variable, etc.).
7. Explore data (visualization, PCA) and perform basic univariate statistics for association
between the target an input variables.
8. Perform more complex multivariate-machine learning.
9. Model validation using a left-out-sample strategy (cross-validation, etc.).
10. Apply on new data.
8 Chapter 1. Introduction
CHAPTER
TWO
PYTHON LANGUAGE
# import a function
from math import sqrt
sqrt(25) # no longer have to reference the module
# define an alias
import numpy as np
# Numbers
10 + 4 # add (returns 14)
10 - 4 # subtract (returns 6)
10 * 4 # multiply (returns 40)
10 ** 4 # exponent (returns 10000)
10 / 4 # divide (returns 2 because both types are 'int')
10 / float(4) # divide (returns 2.5)
5 % 4 # modulo (returns 1) - also known as the remainder
9
Statistics and Machine Learning in Python, Release 0.3 beta
# Boolean operations
# comparisons (these return True)
5 > 3
5 >= 3
5 != 3
5 == 5
Out:
True
Out:
True
2.3.1 Lists
Different objects categorized along a certain ordered sequence, lists are ordered, iterable, mu-
table (adding or removing objects changes the list size), can contain multiple data types.
# create an empty list (two ways)
empty_list = []
empty_list = list()
# create a list
simpsons = ['homer', 'marge', 'bart']
# examine a list
simpsons[0] # print element 0 ('homer')
len(simpsons) # returns the length (3)
# sort a list in place (modifies but does not return the list)
simpsons.sort()
simpsons.sort(reverse=True) # sort in reverse
simpsons.sort(key=len) # sort by a key
# return a sorted list (but does not modify the original list)
sorted(simpsons)
(continues on next page)
# examine objects
id(num) == id(same_num) # returns True
id(num) == id(new_num) # returns False
num is same_num # returns True
num is new_num # returns False
num == same_num # returns True
num == new_num # returns True (their contents are equivalent)
# conatenate +, replicate *
[1, 2, 3] + [4, 5, 6]
["a"] * 2 + ["b"] * 3
Out:
['a', 'a', 'b', 'b', 'b']
2.3.2 Tuples
Like lists, but their size cannot change: ordered, iterable, immutable, can contain multiple data
types
# create a tuple
digits = (0, 1, 'two') # create a tuple directly
digits = tuple([0, 1, 'two']) # create a tuple from a list
zero = (0,) # trailing comma is required to indicate it's a tuple
# examine a tuple
digits[2] # returns 'two'
len(digits) # returns 3
digits.count(0) # counts the number of instances of that value (1)
digits.index(1) # returns the index of the first instance of that value (1)
# concatenate tuples
digits = digits + (3, 4)
# create a single tuple with elements repeated (also works with lists)
(3, 4) * 2 # returns (3, 4, 3, 4)
2.3.3 Strings
# create a string
s = str(42) # convert another data type into a string
s = 'I like you'
# examine a string
s[0] # returns 'I'
len(s) # returns 10
# concatenate strings
s3 = 'The meaning of life is'
s4 = '42'
s3 + ' ' + s4 # returns 'The meaning of life is 42'
s3 + ' ' + str(42) # same thing
# string formatting
# more examples: http://mkaz.com/2012/10/10/python-string-format/
'pi is {:.2f}'.format(3.14159) # returns 'pi is 3.14'
Out:
'pi is 3.14'
Out:
first line
second line
Out:
sequece of bytes are not strings, should be decoded before some operations
print(s.decode('utf-8').split())
Out:
2.3.5 Dictionaries
Dictionaries are structures which can contain multiple data types, and is ordered with key-value
pairs: for each (unique) key, the dictionary outputs one value. Keys can be strings, numbers, or
tuples, while the corresponding values can be any Python object. Dictionaries are: unordered,
iterable, mutable
# examine a dictionary
family['dad'] # returns 'homer'
len(family) # returns 3
family.keys() # returns list: ['dad', 'mom', 'size']
family.values() # returns list: ['homer', 'marge', 6]
family.items() # returns list of tuples:
# [('dad', 'homer'), ('mom', 'marge'), ('size', 6)]
'mom' in family # returns True
'marge' in family # returns False (only checks keys)
Out:
Error 'grandma'
2.3.6 Sets
Like dictionaries, but with unique keys only (no corresponding values). They are: unordered, it-
erable, mutable, can contain multiple data types made up of unique elements (strings, numbers,
or tuples)
# create a set
languages = {'python', 'r', 'java'} # create a set directly
snakes = set(['cobra', 'viper', 'python']) # create a set from a list
# examine a set
len(languages) # returns 3
'python' in languages # returns True
# set operations
languages & snakes # returns intersection: {'python'}
languages | snakes # returns union: {'cobra', 'r', 'java', 'viper', 'python'}
languages - snakes # returns set difference: {'r', 'java'}
snakes - languages # returns set difference: {'cobra', 'viper'}
try:
languages.remove('c') # try to remove a non-existing element (throws an error)
except KeyError as e:
print("Error", e)
Out:
Error 'c'
[0, 1, 2, 9]
x = 3
# if statement
if x > 0:
print('positive')
# if/else statement
if x > 0:
print('positive')
else:
print('zero or negative')
# if/elif/else statement
if x > 0:
print('positive')
elif x == 0:
print('zero')
else:
print('negative')
Out:
positive
positive
positive
positive
2.4.2 Loops
Loops are a set of instructions which repeat until termination conditions are met. This can
include iterating through all values in an object, go through a range of values, etc
# range returns a list of integers
range(0, 3) # returns [0, 1, 2]: includes first value but excludes second value
range(3) # same thing: starting at zero is the default
range(0, 5, 2) # returns [0, 2, 4]: third argument specifies the 'stride'
# for loop
fruits = ['apple', 'banana', 'cherry']
for i in range(len(fruits)):
print(fruits[i].upper())
# use range when iterating over a large sequence to avoid actually creating the integer␣
˓→list in memory
v = 0
for i in range(10 ** 6):
v += 1
Out:
APPLE
BANANA
CHERRY
APPLE
BANANA
CHERRY
quote = """
our incomes are like our shoes; if too small they gall and pinch us
but if too large they cause us to stumble and to trip
"""
# use enumerate if you need to access the index value within the loop
for index, fruit in enumerate(fruits):
print(index, fruit)
# for/else loop
for fruit in fruits:
if fruit == 'banana':
print("Found the banana!")
break # exit the loop and skip the 'else' block
else:
# this block executes ONLY if the for loop completes without hitting 'break'
print("Can't find the banana")
# while loop
count = 0
while count < 5:
print("This will print 5 times")
count += 1 # equivalent to 'count = count + 1'
Out:
dad homer
mom marge
size 6
0 apple
1 banana
2 cherry
Can't find the banana
Found the banana!
This will print 5 times
This will print 5 times
This will print 5 times
This will print 5 times
This will print 5 times
key = 'c'
try:
dct[key]
except:
print("Key %s is missing. Add it with empty value" % key)
dct['c'] = []
print(dct)
Out:
2.5 Functions
Functions are sets of instructions launched when called upon, they can have multiple input
values and a return value
2.5. Functions 19
Statistics and Machine Learning in Python, Release 0.3 beta
#
def add(a, b):
return a + b
add(2, 3)
add("deux", "trois")
# default arguments
def power_this(x, power=2):
return x ** power
power_this(2) # 4
power_this(2, 3) # 8
# return values can be assigned into multiple variables using tuple unpacking
min_num, max_num = min_max(nums) # min_num = 1, max_num = 3
Out:
this is text
3
3
Process which affects whole lists without iterating through loops. For more: http://
python-3-patterns-idioms-test.readthedocs.io/en/latest/Comprehensions.html
# for loop to create a list of cubes
nums = [1, 2, 3, 4, 5]
cubes = []
for num in nums:
cubes.append(num**3)
# set comprehension
fruits = ['apple', 'banana', 'cherry']
unique_lengths = {len(fruit) for fruit in fruits} # {5, 6}
# dictionary comprehension
fruit_lengths = {fruit:len(fruit) for fruit in fruits} # {'apple': 5, 'banana
˓→': 6, 'cherry': 6} (continues on next page)
Out:
import re
Out:
Method/Attribute Purpose
match(string) Determine if the RE matches at the beginning of the string.
search(string) Scan through a string, looking for any location where this RE matches.
findall(string) Find all substrings where the RE matches, and returns them as a list.
finditer(string) Find all substrings where the RE matches, and returns them as an itera-
tor.
regex.sub("SUB-", "toto")
Out:
'toto'
Out:
'helloworld'
import os
Out:
/home/ed203246/git/pystatsml/python_lang
Temporary directory
import tempfile
tmpdir = tempfile.gettempdir()
Join paths
mytmpdir = os.path.join(tmpdir, "foobar")
# list containing the names of the entries in the directory given by path.
os.listdir(tmpdir)
Out:
['tracker-extract-files.16094', 'pymp-b2pv57wx', 'pymp-ewr8p5l2', 'systemd-private-
˓→cd5e03034b9b4cc4abcad7be9bb39d90-fwupd.service-waEboi', 'systemd-private-
˓→cd5e03034b9b4cc4abcad7be9bb39d90-systemd-timesyncd.service-XpeGLh', 'dropbox-antifreeze-
˓→MFb6ch', 'snap.chromium', 'snap.libreoffice', 'config-err-OKsGNU', '.font-unix', 'plop2
˓→', '.X11-unix', 'spyder-ed203246', 'net-export', 'VMwareDnD', '.org.chromium.Chromium.
Create a directory
if not os.path.exists(mytmpdir):
os.mkdir(mytmpdir)
# Write
lines = ["Dans python tout est bon", "Enfin, presque"]
# Read
## read one line at a time (entire file does not have to fit into memory)
f = open(filename, "r")
f.readline() # one string per line (including newlines)
f.readline() # next line
f.close()
## read one line at a time (entire file does not have to fit into memory)
f = open(filename, 'r')
f.readline() # one string per line (including newlines)
f.readline() # next line
f.close()
## use list comprehension to duplicate readlines without reading entire file at once
f = open(filename, 'r')
[line for line in f]
f.close()
Out:
/tmp/foobar/myfile.txt
Walk
import os
WD = os.path.join(tmpdir, "foobar")
Out:
import tempfile
import glob
tmpdir = tempfile.gettempdir()
Out:
['/tmp/plop2/myfile.txt', '/tmp/foobar/myfile.txt']
['myfile', 'myfile']
import shutil
shutil.copy(src, dst)
try:
shutil.copytree(src, dst)
shutil.rmtree(dst)
shutil.move(src, dst)
except (FileExistsError, FileNotFoundError) as e:
pass
Out:
• For more advanced use cases, the underlying Popen interface can be used directly.
• Run the command described by args.
• Wait for command to complete
• return a CompletedProcess instance.
• Does not capture stdout or stderr by default. To do so, pass PIPE for the stdout and/or
stderr arguments.
import subprocess
# Capture output
out = subprocess.run(["ls", "-a", "/"], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
# out.stdout is a sequence of bytes that should be decoded into a utf-8 string
print(out.stdout.decode('utf-8').split("\n")[:5])
Out:
0
['.', '..', 'bin', 'boot', 'cdrom']
Process
A process is a name given to a program instance that has been loaded into memory
and managed by the operating system.
Process = address space + execution context (thread of control)
Process address space (segments):
• Code.
• Data (static/global).
• Heap (dynamic memory allocation).
• Stack.
Execution context:
• Data registers.
• Stack pointer (SP).
• Program counter (PC).
• Working Registers.
OS Scheduling of processes: context switching (ie. save/load Execution context)
Pros/cons
• Context switching expensive.
• (potentially) complex data sharing (not necessary true).
• Cooperating processes - no need for memory protection (separate address
spaces).
• Relevant for parrallel computation with memory allocation.
Threads
• Threads share the same address space (Data registers): access to code, heap
and (global) data.
• Separate execution stack, PC and Working Registers.
Pros/cons
• Faster context switching only SP, PC and Working Registers.
• Can exploit fine-grain concurrency
• Simple data sharing through the shared address space.
• Precautions have to be taken or two threads will write to the same memory at
the same time. This is what the global interpreter lock (GIL) is for.
• Relevant for GUI, I/O (Network, disk) concurrent operation
In Python
• The threading module uses threads.
import time
import threading
startime = time.time()
# Will execute both in parallel
thread1.start()
thread2.start()
# Joins threads back to the parent process
thread1.join()
thread2.join()
print("Threading ellapsed time ", time.time() - startime)
print(out_list[:10])
Out:
Multiprocessing
import multiprocessing
startime = time.time()
p1.start()
p2.start()
p1.join()
p2.join()
print("Multiprocessing ellapsed time ", time.time() - startime)
Out:
import multiprocessing
import time
startime = time.time()
p1.start()
p2.start()
p1.join()
p2.join()
print(out_list[:10])
Out:
import os
import os.path
import argparse
import re
import pandas as pd
if __name__ == "__main__":
# parse command line options
output = "word_count.csv"
parser = argparse.ArgumentParser()
parser.add_argument('-i', '--input',
help='list of input files.',
nargs='+', type=str)
parser.add_argument('-o', '--output',
help='output csv file (default %s)' % output,
(continues on next page)
if options.input is None :
parser.print_help()
raise SystemExit("Error: input files are missing")
else:
filenames = [f for f in options.input if os.path.isfile(f)]
# Match words
regex = re.compile("[a-zA-Z]+")
count = dict()
for filename in filenames:
fd = open(filename, "r")
for line in fd:
for word in regex.findall(line.lower()):
if not word in count:
count[word] = 1
else:
count[word] += 1
fd = open(options.output, "w")
# Pandas
df = pd.DataFrame([[k, count[k]] for k in count], columns=["word", "count"])
df.to_csv(options.output, index=False)
2.10 Networking
# TODO
2.10.1 FTP
Out:
2.10.2 HTTP
# TODO
2.10.3 Sockets
# TODO
2.10.4 xmlrpc
# TODO
A module is a Python file. A package is a directory which MUST contain a special file called
__init__.py
To import, extend variable PYTHONPATH:
export PYTHONPATH=path_to_parent_python_module:${PYTHONPATH}
Or
import sys
sys.path.append("path_to_parent_python_module")
The __init__.py file can be empty. But you can set which modules the package exports as the
API, while keeping other modules internal, by overriding the __all__ variable, like so:
parentmodule/__init__.py file:
import parentmodule.submodule1
import parentmodule.function1
Sources
• http://python-textbok.readthedocs.org/en/latest/Object_Oriented_Programming.html
Principles
• Encapsulate data (attributes) and code (methods) into objects.
• Class = template or blueprint that can be used to create objects.
• An object is a specific instance of a class.
• Inheritance: OOP allows classes to inherit commonly used state and behaviour from other
classes. Reduce code duplication
• Polymorphism: (usually obtained through polymorphism) calling code is agnostic as to
whether an object belongs to a parent class or one of its descendants (abstraction, modu-
larity). The same method called on 2 objects of 2 different classes will behave differently.
import math
class Shape2D:
def area(self):
raise NotImplementedError()
# Inheritance + Encapsulation
class Square(Shape2D):
def __init__(self, width):
self.width = width
def area(self):
return self.width ** 2
class Disk(Shape2D):
def __init__(self, radius):
self.radius = radius
def area(self):
return math.pi * self.radius ** 2
# Polymorphism
print([s.area() for s in shapes])
s = Shape2D()
try:
(continues on next page)
Out:
[4, 28.274333882308138]
NotImplementedError
See PEP 8
• Spaces (four) are the preferred indentation method.
• Two blank lines for top level function or classes definition.
• One blank line to indicate logical sections.
• Never use: from lib import *
• Bad: Capitalized_Words_With_Underscores
• Function and Variable Names: lower_case_with_underscores
• Class Names: CapitalizedWords (aka: CamelCase)
2.14 Documenting
Parameters
----------
a : float
First operand.
b : float, optional
Second operand. The default is 2.
Returns
-------
Sum of operands.
print(help(my_function))
Out:
my_function(a, b=2)
This function ...
Parameters
----------
a : float
First operand.
b : float, optional
Second operand. The default is 2.
Returns
-------
Sum of operands.
Example
-------
>>> my_function(3)
5
None
"""
Created on Thu Nov 14 12:08:41 CET 2019
@author: [email protected]
Some description
"""
2.15 Exercises
Create a function that acts as a simple calulator If the operation is not specified, default to
addition If the operation is misspecified, return an prompt message Ex: calc(4,5,"multiply")
returns 20 Ex: calc(3,5) returns 8 Ex: calc(1, 2, "something") returns error message
Given a list of numbers, return a list where all adjacent duplicate elements have been reduced
to a single element. Ex: [1, 2, 2, 3, 2] returns [1, 2, 3, 2]. You may create a new list or
modify the passed in list.
Remove all duplicate values (adjacent or not) Ex: [1, 2, 2, 3, 2] returns [1, 2, 3]
2.15. Exercises 35
Statistics and Machine Learning in Python, Release 0.3 beta
THREE
SCIENTIFIC PYTHON
NumPy is an extension to the Python programming language, adding support for large, multi-
dimensional (numerical) arrays and matrices, along with a large library of high-level mathe-
matical functions to operate on these arrays.
Sources:
• Kevin Markham: https://github.com/justmarkham
import numpy as np
Create ndarrays from lists. note: every element must be the same type (will be converted if
possible)
Out:
np.zeros(10)
np.zeros((3, 6))
np.ones(10)
np.linspace(0, 1, 5) # 0 to 1 (inclusive) with 5 points
np.logspace(0, 3, 4) # 10^0 to 10^3 (inclusive) with 4 points
Out:
37
Statistics and Machine Learning in Python, Release 0.3 beta
int_array = np.arange(5)
float_array = int_array.astype(float)
arr1.dtype # float64
arr2.dtype # int32
arr2.ndim # 2
arr2.shape # (2, 4) - axis 0 is rows, axis 1 is columns
arr2.size # 8 - total number of elements
len(arr2) # 2 - size of first dimension (aka axis)
Out:
3.1.3 Reshaping
Out:
(2, 5)
[[0. 1.]
[2. 3.]
[4. 5.]
[6. 7.]
[8. 9.]]
Add an axis
a = np.array([0, 1])
a_col = a[:, np.newaxis]
print(a_col)
#or
a_col = a[:, None]
Out:
[[0]
[1]]
Transpose
print(a_col.T)
Out:
[[0 1]]
arr_flt = arr.flatten()
arr_flt[0] = 33
print(arr_flt)
print(arr)
Out:
[33. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[[0. 1. 2. 3. 4.]
[5. 6. 7. 8. 9.]]
arr_flt = arr.ravel()
arr_flt[0] = 33
print(arr_flt)
print(arr)
Out:
[33. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[[33. 1. 2. 3. 4.]
[ 5. 6. 7. 8. 9.]]
Numpy internals: By default Numpy use C convention, ie, Row-major language: The matrix is
stored by rows. In C, the last index changes most rapidly as one moves through the array as
stored in memory.
For 2D arrays, sequential move in the memory will:
• iterate over rows (axis 0)
– iterate over columns (axis 1)
For 3D arrays, sequential move in the memory will:
• iterate over plans (axis 0)
– iterate over rows (axis 1)
Out:
[ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23]
x = x.reshape(2, 3, 4)
print(x)
Out:
[[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
[[12 13 14 15]
[16 17 18 19]
[20 21 22 23]]]
print(x[0, :, :])
Out:
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
print(x[:, 0, :])
Out:
[[ 0 1 2 3]
[12 13 14 15]]
print(x[:, :, 0])
Out:
[[ 0 4 8]
[12 16 20]]
Ravel
print(x.ravel())
Out:
[ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23]
a = np.array([0, 1])
b = np.array([2, 3])
ab = np.stack((a, b)).T
print(ab)
# or
np.hstack((a[:, None], b[:, None]))
Out:
[[0 2]
[1 3]]
array([[0, 2],
[1, 3]])
3.1.6 Selection
Single item
Out:
3.0
Slicing
Syntax: start:stop:step with start (default 0) stop (default last) step (default 1)
Out:
[[1. 2. 3.]
[6. 7. 8.]]
arr2[0, 0] = 33
print(arr2)
print(arr)
Out:
[[33. 2. 3.]
[ 6. 7. 8.]]
[[ 0. 33. 2. 3. 4.]
[ 5. 6. 7. 8. 9.]]
print(arr[0, ::-1])
# The rule of thumb here can be: in the context of lvalue indexing (i.e. the indices are␣
˓→placed in the left hand side value of an assignment), no view or copy of the array is␣
˓→created (because there is no need to). However, with regular values, the above rules␣
˓→for creating views does apply.
Out:
[ 4. 3. 2. 33. 0.]
Out:
[[33. 2. 3.]
[ 6. 7. 8.]]
[[44. 2. 3.]
[ 6. 7. 8.]]
[[ 0. 33. 2. 3. 4.]
[ 5. 6. 7. 8. 9.]]
print(arr2)
arr2[0] = 44
print(arr2)
print(arr)
Out:
[33. 6. 7. 8. 9.]
[44. 6. 7. 8. 9.]
[[ 0. 33. 2. 3. 4.]
[ 5. 6. 7. 8. 9.]]
However, In the context of lvalue indexing (left hand side value of an assignment) Fancy autho-
rizes the modification of the original array
arr[arr > 5] = 0
print(arr)
Out:
[[0. 0. 2. 3. 4.]
[5. 0. 0. 0. 0.]]
Out:
nums = np.arange(5)
nums * 10 # multiply each element by 10
nums = np.sqrt(nums) # square root of each element
np.ceil(nums) # also floor, rint (round to nearest int)
np.isnan(nums) # checks for NaN
nums + np.arange(5) # add element-wise
np.maximum(nums, np.array([1, -2, 3, -4, 5])) # compare element-wise
# random numbers
np.random.seed(12234) # Set the seed
np.random.rand(2, 3) # 2 x 3 matrix in [0, 1]
np.random.randn(10) # random normals (mean 0, sd 1)
np.random.randint(0, 2, 10) # 10 randomly picked 0 or 1
Out:
array([0, 0, 0, 1, 1, 0, 1, 1, 1, 1])
3.1.8 Broadcasting
Rules
Starting with the trailing axis and working backward, Numpy compares arrays dimensions.
• If two dimensions are equal then continues
• If one of the operand has dimension 1 stretches it to match the largest one
• When one of the shapes runs out of dimensions (because it has less dimensions than
the other shape), Numpy will use 1 in the comparison process until the other shape’s
dimensions run out as well.
a = np.array([[ 0, 0, 0],
[10, 10, 10],
[20, 20, 20],
[30, 30, 30]])
b = np.array([0, 1, 2])
print(a + b)
Out:
[[ 0 1 2]
[10 11 12]
[20 21 22]
[30 31 32]]
Examples
Shapes of operands A, B and result:
A (2d array): 5 x 4
B (1d array): 1
Result (2d array): 5 x 4
A (2d array): 5 x 4
B (1d array): 4
Result (2d array): 5 x 4
A (3d array): 15 x 3 x 5
B (3d array): 15 x 1 x 5
Result (3d array): 15 x 3 x 5
A (3d array): 15 x 3 x 5
B (2d array): 3 x 5
Result (3d array): 15 x 3 x 5
A (3d array): 15 x 3 x 5
B (2d array): 3 x 1
Result (3d array): 15 x 3 x 5
3.1.9 Exercises
• For each column find the row index of the minimum value.
• Write a function standardize(X) that return an array whose columns are centered and
scaled (by std-dev).
Total running time of the script: ( 0 minutes 0.011 seconds)
It is often said that 80% of data analysis is spent on the cleaning and small, but important,
aspect of data manipulation and cleaning with Pandas.
Sources:
• Kevin Markham: https://github.com/justmarkham
• Pandas doc: http://pandas.pydata.org/pandas-docs/stable/index.html
Data structures
• Series is a one-dimensional labeled array capable of holding any data type (inte-
gers, strings, floating point numbers, Python objects, etc.). The axis labels are col-
lectively referred to as the index. The basic method to create a Series is to call
pd.Series([1,3,5,np.nan,6,8])
• DataFrame is a 2-dimensional labeled data structure with columns of potentially different
types. You can think of it like a spreadsheet or SQL table, or a dict of Series objects. It
stems from the R data.frame() object.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
print(user3)
Out:
Concatenate DataFrame
user1.append(user2)
users = pd.concat([user1, user2, user3])
print(users)
Out:
Join DataFrame
Out:
name height
0 alice 165
1 john 180
2 eric 175
3 julie 171
print(merge_inter)
Out:
Out:
Reshaping by pivoting
Out:
Out:
3.2.3 Summarizing
Out:
<class 'pandas.core.frame.DataFrame'>
Int64Index: 6 entries, 0 to 5
Data columns (total 5 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 name 6 non-null object
1 age 6 non-null int64
2 gender 6 non-null object
3 job 6 non-null object
4 height 4 non-null float64
dtypes: float64(1), int64(1), object(3)
memory usage: 288.0+ bytes
df = users.copy()
df.iloc[0] # first row
df.iloc[0, 0] # first item of first row
df.iloc[0, 0] = 55
for i in range(users.shape[0]):
row = df.iloc[i]
row.age *= 100 # setting a copy, and not the original frame data.
Out:
/home/ed203246/anaconda3/lib/python3.7/site-packages/pandas/core/generic.py:5168:␣
˓→SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame
df = users.copy()
df.loc[0] # first row
df.loc[0, "age"] # first item of first row
df.loc[0, "age"] = 55
for i in range(df.shape[0]):
df.loc[i, "age"] *= 10
print(df) # df is modified
Out:
Out:
3.2.7 Sorting
df = users.copy()
print(df)
Out:
print(df.describe())
Out:
age height
count 6.000000 4.000000
mean 33.666667 172.750000
std 14.895189 6.344289
min 19.000000 165.000000
25% 23.000000 169.500000
50% 29.500000 173.000000
75% 41.250000 176.250000
max 58.000000 180.000000
print(df.describe(include='all'))
print(df.describe(include=['object'])) # limit to one (or more) types
Out:
print(df.groupby("job").mean())
print(df.groupby("job")["age"].mean())
print(df.groupby("job").describe(include='all'))
Out:
age height
job
engineer 33.000000 NaN
manager 58.000000 NaN
scientist 44.000000 171.000000
student 22.333333 173.333333
job
engineer 33.000000
manager 58.000000
scientist 44.000000
student 22.333333
(continues on next page)
engineer 1 1 peter 1 NaN NaN NaN NaN NaN NaN NaN 1.0 NaN NaN ..
˓→. NaN NaN NaN 0.0 NaN NaN NaN NaN NaN NaN NaN NaN ␣
˓→NaN NaN
manager 1 1 paul 1 NaN NaN NaN NaN NaN NaN NaN 1.0 NaN NaN ..
˓→. NaN NaN NaN 0.0 NaN NaN NaN NaN NaN NaN NaN NaN ␣
˓→NaN NaN
scientist 1 1 julie 1 NaN NaN NaN NaN NaN NaN NaN 1.0 NaN NaN ..
˓→. NaN NaN NaN 1.0 NaN NaN NaN 171.000000 NaN 171.0 171.0 171.0 ␣
˓→171.0 171.0
student 3 3 eric 1 NaN NaN NaN NaN NaN NaN NaN 3.0 NaN NaN ..
˓→. NaN NaN NaN 3.0 NaN NaN NaN 173.333333 7.637626 165.0 170.0 175.0 ␣
˓→177.5 180.0
[4 rows x 44 columns]
Groupby in a loop
Out:
df = users.append(df.iloc[0], ignore_index=True)
Out:
0 False
1 False
2 False
3 False
4 False
5 False
6 True
dtype: bool
Missing data
Out:
name 0
age 0
gender 0
job 0
height 2
dtype: int64
print(df)
Out:
name age gender job height
0 alice 19 F student 165.00
1 john 26 M student 180.00
(continues on next page)
df = users.copy()
print(df.columns)
df.columns = ['age', 'genre', 'travail', 'nom', 'taille']
df['travail'].str.contains("etu|inge")
Out:
0 NaN
1 NaN
2 NaN
3 NaN
4 NaN
5 NaN
Name: travail, dtype: object
Assume random variable follows the normal distribution Exclude data outside 3 standard-
deviations: - Probability that a sample lies within 1 sd: 68.27% - Probability that a sample
lies within 3 sd: 99.73% (68.27 + 2 * 15.73)
size_outlr_mean = size.copy()
size_outlr_mean[((size - size.mean()).abs() > 3 * size.std())] = size.mean()
print(size_outlr_mean.mean())
Out:
251.72690773066762
Median absolute deviation (MAD), based on the median, is a robust non-parametric statistics.
https://en.wikipedia.org/wiki/Median_absolute_deviation
Out:
179.66768954321032 184.46486745682998
csv
url = 'https://raw.github.com/neurospin/pystatsml/master/datasets/salary_table.csv'
salary = pd.read_csv(url)
Excel
pd.read_excel(xls_filename, sheet_name='users')
# Multiple sheets
with pd.ExcelWriter(xls_filename) as writer:
users.to_excel(writer, sheet_name='users', index=False)
df.to_excel(writer, sheet_name='salary', index=False)
pd.read_excel(xls_filename, sheet_name='users')
pd.read_excel(xls_filename, sheet_name='salary')
SQL (SQLite)
import pandas as pd
import sqlite3
Connect
conn = sqlite3.connect(db_filename)
url = 'https://raw.github.com/neurospin/pystatsml/master/datasets/salary_table.csv'
salary = pd.read_csv(url)
Push modifications
cur = conn.cursor()
values = (100, 14000, 5, 'Bachelor', 'N')
cur.execute("insert into salary values (?, ?, ?, ?, ?)", values)
conn.commit()
Out:
3.2.13 Exercises
Data Frame
4. For each species compute the mean of numerical columns and store it in a stats table
like:
Missing data
df = users.copy()
df.loc[[0, 2], "age"] = None
df.loc[[1, 3], "gender"] = None
1. Write a function fillmissing_with_mean(df) that fill all missing value of numerical col-
umn with the mean of the current columns.
2. Save the original users and “imputed” frame in a single excel file “users.xlsx” with 2 sheets:
original, imputed.
Total running time of the script: ( 0 minutes 1.137 seconds)
import numpy as np
import matplotlib.pyplot as plt
plt.plot(x, sinus)
plt.show()
# Rapid multiplot
cosinus = np.cos(x)
plt.plot(x, sinus, "-b", x, sinus, "ob", x, cosinus, "-r", x, cosinus, "or")
plt.xlabel('this is x!')
plt.ylabel('this is y!')
plt.title('My First Plot')
plt.show()
# Step by step
plt.plot(x, sinus, label='sinus', color='blue', linestyle='--', linewidth=2)
plt.plot(x, cosinus, label='cosinus', color='red', linestyle='-', linewidth=2)
plt.legend()
plt.show()
Load dataset
import pandas as pd
try:
salary = pd.read_csv("../datasets/salary_table.csv")
except:
url = 'https://raw.github.com/neurospin/pystatsml/master/datasets/salary_table.csv'
salary = pd.read_csv(url)
df = salary
<matplotlib.collections.PathCollection at 0x7f2d82e0be10>
## Figure size
plt.figure(figsize=(6,5))
## Set labels
plt.xlabel('Experience')
plt.ylabel('Salary')
plt.legend(loc=4) # lower right
plt.show()
# Prefer vectorial format (SVG: Scalable Vector Graphics) can be edited with
# Inkscape, Adobe Illustrator, Blender, etc.
plt.plot(x, sinus)
plt.savefig("sinus.svg")
plt.close()
# Or pdf
plt.plot(x, sinus)
(continues on next page)
3.3.4 Seaborn
Boxplot
Box plots are non-parametric: they display variation in samples of a statistical population with-
out making any assumptions of the underlying statistical distribution.
<AxesSubplot:xlabel='education', ylabel='salary'>
<AxesSubplot:xlabel='management', ylabel='salary'>
Fig. 2: title
i = 0
for edu, d in salary.groupby(['education']):
sns.distplot(d.salary[d.management == "Y"], color="b", bins=10, label="Manager",␣
˓→ax=axes[i])
sns.distplot(d.salary[d.management == "N"], color="r", bins=10, label="Employee",␣
˓→ax=axes[i])
axes[i].set_title(edu)
axes[i].set_ylabel('Density')
i += 1
ax = plt.legend()
/home/ed203246/anaconda3/lib/python3.7/site-packages/seaborn/distributions.
˓→py:2551: FutureWarning: distplot is a deprecated function and will be␣
˓→removed in a future version. Please adapt your code to use either displot (a␣
warnings.warn(msg, FutureWarning)
/home/ed203246/anaconda3/lib/python3.7/site-packages/seaborn/distributions.
˓→py:2551: FutureWarning: distplot is a deprecated function and will be␣
˓→removed in a future version. Please adapt your code to use either displot (a␣
warnings.warn(msg, FutureWarning)
/home/ed203246/anaconda3/lib/python3.7/site-packages/seaborn/distributions.
˓→py:2551: FutureWarning: distplot is a deprecated function and will be␣
˓→removed in a future version. Please adapt your code to use either displot (a␣
˓→removed in a future version. Please adapt your code to use either displot (a␣
warnings.warn(msg, FutureWarning)
/home/ed203246/anaconda3/lib/python3.7/site-packages/seaborn/distributions.
˓→py:2551: FutureWarning: distplot is a deprecated function and will be␣
˓→removed in a future version. Please adapt your code to use either displot (a␣
warnings.warn(msg, FutureWarning)
/home/ed203246/anaconda3/lib/python3.7/site-packages/seaborn/distributions.
˓→py:2551: FutureWarning: distplot is a deprecated function and will be␣
˓→removed in a future version. Please adapt your code to use either displot (a␣
warnings.warn(msg, FutureWarning)
ax = sns.violinplot(x="salary", data=salary)
Tune bandwidth
Tips dataset One waiter recorded information about each tip he received over a period of a few
months working in one restaurant. He collected several variables:
ax = sns.violinplot(x=tips["total_bill"])
Group by day
g = sns.PairGrid(salary, hue="management")
g.map_diag(plt.hist)
g.map_offdiag(plt.scatter)
ax = g.add_legend()
ax = sns.pointplot(x="timepoint", y="signal",
hue="region", style="event",
data=fmri)
# version 0.9
# sns.lineplot(x="timepoint", y="signal",
# hue="region", style="event",
# data=fmri)
FOUR
STATISTICS
Mean
The estimator 𝑥
¯ on a sample of size 𝑛: 𝑥 = 𝑥1 , ..., 𝑥𝑛 is given by
1 ∑︁
𝑥
¯= 𝑥𝑖
𝑛
𝑖
75
Statistics and Machine Learning in Python, Release 0.3 beta
Variance
Note here the subtracted 1 degree of freedom (df) in the divisor. In standard statistical practice,
𝑑𝑓 = 1 provides an unbiased estimator of the variance of a hypothetical infinite population.
With 𝑑𝑓 = 0 it instead provides a maximum likelihood estimate of the variance for normally
distributed variables.
Standard deviation
√︀
𝑆𝑡𝑑(𝑋) = 𝑉 𝑎𝑟(𝑋)
√︀
The estimator is simply 𝜎𝑥 = 𝜎𝑥2 .
Covariance
Cov(𝑋, 𝑋) = Var(𝑋)
Cov(𝑋, 𝑌 ) = Cov(𝑌, 𝑋)
Cov(𝑐𝑋, 𝑌 ) = 𝑐 Cov(𝑋, 𝑌 )
Cov(𝑋 + 𝑐, 𝑌 ) = Cov(𝑋, 𝑌 )
Correlation
𝐶𝑜𝑣(𝑋, 𝑌 )
𝐶𝑜𝑟(𝑋, 𝑌 ) =
𝑆𝑡𝑑(𝑋)𝑆𝑡𝑑(𝑌 )
The estimator is
𝜎𝑥𝑦
𝜌𝑥𝑦 = .
𝜎𝑥 𝜎𝑦
76 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
The standard error (SE) is the standard deviation (of the sampling distribution) of a statistic:
𝑆𝑡𝑑(𝑋)
𝑆𝐸(𝑋) = √ .
𝑛
Exercises
• Generate 2 random samples: 𝑥 ∼ 𝑁 (1.78, 0.1) and 𝑦 ∼ 𝑁 (1.66, 0.1), both of size 10.
• Compute 𝑥¯, 𝜎𝑥 , 𝜎𝑥𝑦 (xbar, xvar, xycov) using only the np.sum() operation. Explore
the np. module to find out which numpy functions performs the same computations and
compare them (using assert) with your previous results.
Normal distribution
The normal distribution, noted 𝒩 (𝜇, 𝜎) with parameters: 𝜇 mean (location) and 𝜎 > 0 std-dev.
Estimators: 𝑥
¯ and 𝜎𝑥 .
The normal distribution, noted 𝒩 , is useful because of the central limit theorem (CLT) which
states that: given certain conditions, the arithmetic mean of a sufficiently large number of iter-
ates of independent random variables, each with a well-defined expected value and well-defined
variance, will be approximately normally distributed, regardless of the underlying distribution.
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
%matplotlib inline
mu = 0 # mean
variance = 2 #variance
sigma = np.sqrt(variance) #standard deviation",
x = np.linspace(mu-3*variance,mu+3*variance, 100)
plt.plot(x, norm.pdf(x, mu, sigma))
[<matplotlib.lines.Line2D at 0x7f3c95a58b50>]
The chi-square or 𝜒2𝑛 distribution with 𝑛 degrees of freedom (df) is the distribution of a sum of
the squares of 𝑛 independent standard normal random variables 𝒩 (0, 1). Let 𝑋 ∼ 𝒩 (𝜇, 𝜎 2 ),
then, 𝑍 = (𝑋 − 𝜇)/𝜎 ∼ 𝒩 (0, 1), then:
• The squared standard 𝑍 2 ∼ 𝜒21 (one df).
∑︀𝑛
• The distribution of sum of squares of 𝑛 normal random variables: 𝑖 𝑍𝑖2 ∼ 𝜒2𝑛
The sum of two 𝜒2 RV with 𝑝 and 𝑞 df is a 𝜒2 RV with 𝑝 + 𝑞 df. This is useful when sum-
ming/subtracting sum of squares.
The 𝜒2 -distribution is used to model errors measured as sum of squares or the distribution of
the sample variance.
The 𝐹 -distribution, 𝐹𝑛,𝑝 , with 𝑛 and 𝑝 degrees of freedom is the ratio of two independent 𝜒2
variables. Let 𝑋 ∼ 𝜒2𝑛 and 𝑌 ∼ 𝜒2𝑝 then:
𝑋/𝑛
𝐹𝑛,𝑝 =
𝑌 /𝑝
The 𝐹 -distribution plays a central role in hypothesis testing answering the question: Are two
variances equals?, is the ratio or two errors significantly large ?.
import numpy as np
from scipy.stats import f
import matplotlib.pyplot as plt
%matplotlib inline
78 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
Let 𝑀 ∼ 𝒩 (0, 1) and 𝑉 ∼ 𝜒2𝑛 . The 𝑡-distribution, 𝑇𝑛 , with 𝑛 degrees of freedom is the ratio:
𝑀
𝑇𝑛 = √︀
𝑉 /𝑛
The distribution of the difference between an estimated parameter and its true (or assumed)
value divided by the standard deviation of the estimated parameter (standard error) follow a
𝑡-distribution. Is this parameters different from a given value?
Examples
• Test a proportion: Biased coin ? 200 heads have been found over 300 flips, is it coins
biased ?
• Test the association between two variables.
– Exemple height and sex: In a sample of 25 individuals (15 females, 10 males), is
female height is different from male height ?
– Exemple age and arterial hypertension: In a sample of 25 individuals is age height
correlated with arterial hypertension ?
Steps
1. Model the data
2. Fit: estimate the model parameters (frequency, mean, correlation, regression coeficient)
3. Compute a test statistic from model the parameters.
4. Formulate the null hypothesis: What would be the (distribution of the) test statistic if the
observations are the result of pure chance.
5. Compute the probability (𝑝-value) to obtain a larger value for the test statistic by chance
(under the null hypothesis).
Biased coin ? 2 heads have been found over 3 flips, is it coins biased ?
1. Model the data: number of heads follow a Binomial disctribution.
2. Compute model parameters: N=3, P = the frequency of number of heads over the number
of flip: 2/3.
3. Compute a test statistic, same as frequency.
4. Under the null hypothesis the distribution of the number of tail is:
80 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
1 2 3 count #heads
0
H 1
H 1
H 1
H H 2
H H 2
H H 2
H H H 3
8 possibles configurations, probabilities of differents values for 𝑝 are: 𝑥 measure the number of
success.
• 𝑃 (𝑥 = 0) = 1/8
• 𝑃 (𝑥 = 1) = 3/8
• 𝑃 (𝑥 = 2) = 3/8
• 𝑃 (𝑥 = 3) = 1/8
Text(0.5, 0, 'Distribution of the number of head over 3 flip under the null hypothesis')
3. Compute the probability (𝑝-value) to observe a value larger or equal that 2 under the null
hypothesis ? This probability is the 𝑝-value:
Biased coin ? 60 heads have been found over 100 flips, is it coins biased ?
1. Model the data: number of heads follow a Binomial disctribution.
2. Compute model parameters: N=100, P=60/100.
3. Compute a test statistic, same as frequency.
4. Compute a test statistic: 60/100.
5. Under the null hypothesis the distribution of the number of tail (𝑘) follow the binomial
distribution of parameters N=100, P=0.5:
(︂ )︂
100
𝑃 𝑟(𝑋 = 𝑘|𝐻0 ) = 𝑃 𝑟(𝑋 = 𝑘|𝑛 = 100, 𝑝 = 0.5) = 0.5𝑘 (1 − 0.5)(100−𝑘) .
𝑘
100 (︂ )︂
∑︁ 100
𝑃 (𝑋 = 𝑘 ≥ 60|𝐻0 ) = 0.5𝑘 (1 − 0.5)(100−𝑘)
𝑘
𝑘=60
60 (︂ )︂
∑︁ 100
=1− 0.5𝑘 (1 − 0.5)(100−𝑘) , the cumulative distribution function.
𝑘
𝑘=1
import scipy.stats
import matplotlib.pyplot as plt
_ = plt.legend()
0.01760010010885238
82 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
The one-sample 𝑡-test is used to determine whether a sample comes from a population with a
specific mean. For example you want to test if the average height of a population is 1.75 𝑚.
Assumptions
In testing the null hypothesis that the population mean is equal to a specified value 𝜇0 = 1.75,
one uses the statistic:
difference of means √
𝑡= 𝑛 (4.8)
std-dev √of noise
𝑡 = effect size 𝑛 (4.9)
¯ − 𝜇0 √
𝑥
𝑡= 𝑛 (4.10)
𝑠𝑥
4 Compute the probability of the test statistic under the null hypotheis. This require to have
the distribution of the t statistic under 𝐻0 .
Example
Given the following samples, we will test whether its true mean is 1.75.
Warning, when computing the std or the variance, set ddof=1. The default value, ddof=0, leads
to the biased estimator of the variance.
import numpy as np
x = [1.83, 1.83, 1.73, 1.82, 1.83, 1.73, 1.99, 1.85, 1.68, 1.87]
print(xbar)
84 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
1.816
2.3968766311585883
The :math:`p`-value is the probability to observe a value 𝑡 more extreme than the observed
one 𝑡𝑜𝑏𝑠 under the null hypothesis 𝐻0 : 𝑃 (𝑡 > 𝑡𝑜𝑏𝑠 |𝐻0 )
import numpy as np
import scipy.stats as stats
n = 50
x = np.random.normal(size=n)
y = 2 * x + np.random.normal(size=n)
86 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
0.9421665358030606 1.9724847076189737e-24
The two-sample 𝑡-test (Snedecor and Cochran, 1989) is used to determine if two population
means are equal. There are several variations on this test. If data are paired (e.g. 2 measures,
before and after treatment for each individual) use the one-sample 𝑡-test of the difference. The
variances of the two samples may be assumed to be equal (a.k.a. homoscedasticity) or unequal
(a.k.a. heteroscedasticity).
Assumptions
Assume that the two random variables are normally distributed: 𝑦1 ∼ 𝒩 (𝜇1 , 𝜎1 ), 𝑦2 ∼
𝒩 (𝜇2 , 𝜎2 ).
3. 𝑡-test
difference of means
𝑡= (4.11)
standard dev of error
difference of means
= (4.12)
its standard error
𝑦¯1 − 𝑦¯2 √
= √︀∑︀ 𝑛−2 (4.13)
𝜀2
𝑦¯1 − 𝑦¯2
= (4.14)
𝑠𝑦¯1 −𝑦¯2
𝑠2𝑦1 𝑠2𝑦
𝑠2𝑦¯1 −𝑦¯2 = 𝑠2𝑦¯1 + 𝑠2𝑦¯2 = + 2 (4.15)
𝑛1 𝑛2
thus (4.16)
√︃
𝑠2𝑦1 𝑠2𝑦
𝑠𝑦¯1 −𝑦¯2 = + 2 (4.17)
𝑛1 𝑛2
To compute the 𝑝-value one needs the degrees of freedom associated with this variance estimate.
It is approximated using the Welch–Satterthwaite equation:
(︂ )︂2
𝑠2𝑦1 𝑠2𝑦2
𝑛1 + 𝑛2
𝜈≈ 𝑠4𝑦1 𝑠4𝑦2
.
𝑛21 (𝑛1 −1)
+ 𝑛22 (𝑛2 −1)
88 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
If we assume equal variance (ie, 𝑠2𝑦1 = 𝑠2𝑦1 = 𝑠2 ), where 𝑠2 is an estimator of the common
variance of the two samples:
then
√︃ √︂
𝑠2 𝑠2 1 1
𝑠𝑦¯1 −𝑦¯2 = + =𝑠 +
𝑛1 𝑛2 𝑛1 𝑛2
Therefore, the 𝑡 statistic, that is used to test whether the means are different is:
𝑦¯ − 𝑦¯2
𝑡= √︁1 ,
𝑠 · 𝑛11 + 1
𝑛2
𝑦¯1 − 𝑦¯2 √
𝑡= √ · 𝑛 (4.20)
𝑠 2
√
≈ effect size · 𝑛 (4.21)
difference of means √
≈ · 𝑛 (4.22)
standard deviation of the noise
Example
Given the following two samples, test whether their means are equal using the standard t-test,
assuming equal variance.
height = np.array([ 1.83, 1.83, 1.73, 1.82, 1.83, 1.73, 1.99, 1.85, 1.68, 1.87,
1.66, 1.71, 1.73, 1.64, 1.70, 1.60, 1.79, 1.73, 1.62, 1.77])
Ttest_indResult(statistic=3.5511519888466885, pvalue=0.00228208937112721)
Analysis of variance (ANOVA) provides a statistical test of whether or not the means of several
(k) groups are equal, and therefore generalizes the 𝑡-test to more than two groups. ANOVAs
are useful for comparing (testing) three or more means (groups or variables) for statistical
significance. It is conceptually similar to multiple two-sample 𝑡-tests, but is less conservative.
Here we will consider one-way ANOVA with one independent variable, ie one-way anova.
Wikipedia:
• Test if any group is on average superior, or inferior, to the others versus the null hypothesis
that all four strategies yield the same mean response
• Detect any of several possible differences.
• The advantage of the ANOVA 𝐹 -test is that we do not need to pre-specify which strategies
are to be compared, and we do not need to adjust for making multiple comparisons.
• The disadvantage of the ANOVA 𝐹 -test is that if we reject the null hypothesis, we do not
know which strategies can be said to be significantly different from the others.
Assumptions
1. The samples are randomly selected in an independent manner from the k populations.
2. All k populations have distributions that are approximately normal. Check by plotting
groups distribution.
3. The k population variances are equal. Check by plotting groups distribution.
Is there a difference in Petal Width in species from iris dataset. Let 𝑦1 , 𝑦2 and 𝑦3 be Petal Width
in three species.
Here we assume (see assumptions) that the three populations were sampled from three random
variables that are normally distributed. I.e., 𝑌1 ∼ 𝑁 (𝜇1 , 𝜎1 ), 𝑌2 ∼ 𝑁 (𝜇2 , 𝜎2 ) and 𝑌3 ∼ 𝑁 (𝜇3 , 𝜎3 ).
3. 𝐹 -test
Explained variance
𝐹 = (4.23)
Unexplained variance
Between-group variability 𝑠2
= = 2𝐵 . (4.24)
Within-group variability 𝑠𝑊
90 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
where 𝑦¯𝑖· denotes the sample mean in the 𝑖th group, 𝑛𝑖 is the number of observations in the 𝑖th
group, 𝑦¯ denotes the overall mean of the data, and 𝐾 denotes the number of groups.
The “unexplained variance”, or “within-group variability” is
∑︁
𝑠2𝑊 = (𝑦𝑖𝑗 − 𝑦¯𝑖· )2 /(𝑁 − 𝐾),
𝑖𝑗
where 𝑦𝑖𝑗 is the 𝑗th observation in the 𝑖th out of 𝐾 groups and 𝑁 is the overall sample size.
This 𝐹 -statistic follows the 𝐹 -distribution with 𝐾 − 1 and 𝑁 − 𝐾 degrees of freedom under the
null hypothesis. The statistic will be large if the between-group variability is large relative to
the within-group variability, which is unlikely to happen if the population means of the groups
all have the same value.
Note that when there are only two groups for the one-way ANOVA F-test, 𝐹 = 𝑡2 where 𝑡 is the
Student’s 𝑡 statistic.
# Group means
means = iris.groupby("Species").mean().reset_index()
print(means)
# Plot groups
ax = sns.violinplot(x="Species", y="Sepal_Length", data=iris)
ax = sns.swarmplot(x="Species", y="Sepal_Length", data=iris,
color="white")
ax = sns.swarmplot(x="Species", y="Sepal_Length", color="black", data=means, size=10)
# ANOVA
lm = ols('Sepal_Length ~ Species', data=iris).fit()
sm.stats.anova_lm(lm, typ=2) # Type 2 ANOVA DataFrame
Computes the chi-square, 𝜒2 , statistic and 𝑝-value for the hypothesis test of independence of
frequencies in the observed contingency table (cross-table). The observed frequencies are tested
against an expected contingency table obtained by computing expected frequencies based on
the marginal sums under the assumption of independence.
Example: 20 participants: 10 exposed to some chemical product and 10 non exposed (exposed
= 1 or 0). Among the 20 participants 10 had cancer 10 not (cancer = 1 or 0). 𝜒2 tests the
association between those two variables.
import numpy as np
import pandas as pd
import scipy.stats as stats
# Dataset:
# 15 samples:
# 10 first exposed
exposed = np.array([1] * 10 + [0] * 10)
# 8 first with cancer, 10 without, the last two with.
cancer = np.array([1] * 8 + [0] * 10 + [1] * 2)
92 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
Observed table:
---------------
cancer 0 1
exposed
0 8 2
1 2 8
Statistics:
-----------
Chi2 = 5.000000, pval = 0.025347
Expected table:
---------------
[[5. 5.]
[5. 5.]]
cancer_marg = crosstab.sum(axis=1)
cancer_freq = cancer_marg / cancer_marg.sum()
print('Expected frequencies:')
print(np.outer(exposed_freq, cancer_freq))
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
x = np.array([44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 46, 47, 48, 60.1])
y = np.array([2.6, 3.1, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 4, 4.1, 4.5, 3.8])
plt.plot(x, y, "bo")
# Non-Parametric Spearman
cor, pval = stats.spearmanr(x, y)
print("Non-Parametric Spearman cor test, cor: %.4f, pval: %.4f" % (cor, pval))
94 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
Source: https://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test
The Wilcoxon signed-rank test is a non-parametric statistical hypothesis test used when com-
paring two related samples, matched samples, or repeated measurements on a single sample
to assess whether their population mean ranks differ (i.e. it is a paired difference test). It is
equivalent to one-sample test of the difference of paired samples.
It can be used as an alternative to the paired Student’s 𝑡-test, 𝑡-test for matched pairs, or the 𝑡-
test for dependent samples when the population cannot be assumed to be normally distributed.
When to use it? Observe the data distribution: - presence of outliers - the distribution of the
residuals is not Gaussian
It has a lower sensitivity compared to 𝑡-test. May be problematic to use when the sample size is
small.
Null hypothesis 𝐻0 : difference between the pairs follows a symmetric distribution around zero.
# create an outlier
bv1[0] -= 10
# Paired t-test
print(stats.ttest_rel(bv0, bv1))
# Wilcoxon
print(stats.wilcoxon(bv0, bv1))
Ttest_relResult(statistic=0.7508746146807028, pvalue=0.4619272512287318)
WilcoxonResult(statistic=28.0, pvalue=0.002712249755859375)
# create an outlier
bv1[0] -= 10
# Two-samples t-test
print(stats.ttest_ind(bv0, bv1))
# Wilcoxon
print(stats.mannwhitneyu(bv0, bv1))
Ttest_indResult(statistic=0.5170316045185565, pvalue=0.6081306040701799)
MannwhitneyuResult(statistic=45.0, pvalue=1.4624324663684467e-05)
96 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
Assumptions
Using the dataset “salary”, explore the association between the dependant variable (e.g. Salary)
and the independent variable (e.g.: Experience is quantitative), considering only non-managers.
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
url = 'https://raw.github.com/neurospin/pystatsml/master/datasets/salary_table.csv'
salary = pd.read_csv(url)
salary = salary[salary.management == 'N']
Model the data on some hypothesis e.g.: salary is a linear function of the experience.
salary𝑖 = 𝛽0 + 𝛽 experience𝑖 + 𝜖𝑖 ,
more generally
𝑦𝑖 = 𝛽0 + 𝛽 𝑥𝑖 + 𝜖𝑖
This can be rewritten in the matrix form using the design matrix made of values of independant
variable and the intercept:
⎡ ⎤ ⎡ ⎤ ⎡ ⎤
𝑦1 1 𝑥1 𝜖1
⎢𝑦2 ⎥ ⎢1 𝑥2 ⎥ [︂ ]︂ ⎢𝜖2 ⎥
⎢𝑦3 ⎥ = ⎢1 𝑥3 ⎥ 𝛽0 + ⎢𝜖3 ⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎣𝑦4 ⎦ ⎣1 𝑥4 ⎦ 𝛽1
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎣𝜖4 ⎦
𝑦5 1 𝑥5 𝜖5
Recall from calculus that an extreme point can be found by computing where the derivative is
zero, i.e. to find the intercept, we perform the steps:
𝜕𝑆𝑆𝐸 ∑︁
= (𝑦𝑖 − 𝛽 𝑥𝑖 − 𝛽0 ) = 0
𝜕𝛽0
𝑖
∑︁ ∑︁
𝑦𝑖 = 𝛽 𝑥𝑖 + 𝑛 𝛽0
𝑖 𝑖
𝑛 𝑦¯ = 𝑛 𝛽 𝑥
¯ + 𝑛 𝛽0
𝛽0 = 𝑦¯ − 𝛽 𝑥
¯
To find the regression coefficient, we perform the steps:
𝜕𝑆𝑆𝐸 ∑︁
= 𝑥𝑖 (𝑦𝑖 − 𝛽 𝑥𝑖 − 𝛽0 ) = 0
𝜕𝛽
𝑖
Plug in 𝛽0 :
∑︁
𝑥𝑖 (𝑦𝑖 − 𝛽 𝑥𝑖 − 𝑦¯ + 𝛽 𝑥
¯) = 0
𝑖
∑︁ ∑︁ ∑︁
𝑥𝑖 𝑦𝑖 − 𝑦¯ 𝑥𝑖 = 𝛽 (𝑥𝑖 − 𝑥
¯)
𝑖 𝑖 𝑖
print("Using seaborn")
import seaborn as sns
sns.regplot(x="experience", y="salary", data=salary)
98 Chapter 4. Statistics
Statistics and Machine Learning in Python, Release 0.3 beta
Using seaborn
<AxesSubplot:xlabel='experience', ylabel='salary'>
R-squared
The goodness of fit of a statistical model describes how well it fits a set of observations. Mea-
sures of goodness of fit typically summarize the discrepancy between observed values and the
values expected under the model in question. We will consider the explained variance also
known as the coefficient of determination, denoted 𝑅2 pronounced R-squared.
The total sum of squares, 𝑆𝑆tot is the sum of the sum of squares explained by the regression,
𝑆𝑆reg , plus the sum of squares of residuals unexplained by the regression, 𝑆𝑆res , also called the
SSE, i.e. such that
Fig. 4: title
The mean of 𝑦 is
1 ∑︁
𝑦¯ = 𝑦𝑖 .
𝑛
𝑖
The total sum of squares is the total squared sum of deviations from the mean of 𝑦, i.e.
∑︁
𝑆𝑆tot = (𝑦𝑖 − 𝑦¯)2
𝑖
The regression sum of squares, also called the explained sum of squares:
∑︁
𝑆𝑆reg = 𝑦𝑖 − 𝑦¯)2 ,
(ˆ
𝑖
where 𝑦ˆ𝑖 = 𝛽𝑥𝑖 + 𝛽0 is the estimated value of salary 𝑦ˆ𝑖 given a value of experience 𝑥𝑖 .
The sum of squares of the residuals, also called the residual sum of squares (RSS) is:
∑︁
𝑆𝑆res = (𝑦𝑖 − 𝑦ˆ𝑖 )2 .
𝑖
𝑅2 is the explained sum of squares of errors. It is the variance explain by the regression divided
by the total variance, i.e.
explained SS 𝑆𝑆reg 𝑆𝑆𝑟𝑒𝑠
𝑅2 = = =1− .
total SS 𝑆𝑆𝑡𝑜𝑡 𝑆𝑆𝑡𝑜𝑡
Test
ˆ 2 = 𝑆𝑆res /(𝑛 − 2) be an estimator of the variance of 𝜖. The 2 in the denominator stems
Let 𝜎
from the 2 estimated parameters: intercept and coefficient.
𝑆𝑆res
• Unexplained variance: ^2
𝜎
∼ 𝜒2𝑛−2
𝑆𝑆
• Explained variance: 𝜎^ 2reg ∼ 𝜒21 . The single degree of freedom comes from the difference
between 𝑆𝑆
^2
𝜎
tot
(∼ 𝜒2𝑛−1 ) and 𝑆𝑆
^2
𝜎
res
(∼ 𝜒2𝑛−2 ), i.e. (𝑛 − 1) − (𝑛 − 2) degree of freedom.
The Fisher statistics of the ratio of two variances:
Explained variance 𝑆𝑆reg /1
𝐹 = = ∼ 𝐹 (1, 𝑛 − 2)
Unexplained variance 𝑆𝑆res /(𝑛 − 2)
Using the 𝐹 -distribution, compute the probability of observing a value greater than 𝐹 under
𝐻0 , i.e.: 𝑃 (𝑥 > 𝐹 |𝐻0 ), i.e. the survival function (1 − Cumulative Distribution Function) at 𝑥 of
the given 𝐹 -distribution.
y_mu = np.mean(y)
ss_tot = np.sum((y - y_mu) ** 2)
ss_res = np.sum(res ** 2)
r-squared: 0.932
MSE 459923.068
MAE 602.265
Multiple regression
Theory
or, simplified
𝑃 −1
𝛽𝑗 𝑥𝑗𝑖 + 𝜀𝑖 .
∑︁
𝑦𝑖 = 𝛽0 +
𝑗=1
Extending each sample with an intercept, 𝑥𝑖 := [1, 𝑥𝑖 ] ∈ 𝑅𝑃 +1 allows us to use a more general
notation based on linear algebra and write it as a simple dot product:
𝑦𝑖 = x𝑇𝑖 𝛽 + 𝜀𝑖 ,
where 𝛽 ∈ 𝑅𝑃 +1 is a vector of weights that define the 𝑃 + 1 parameters of the model. From
now we have 𝑃 regressors + the intercept.
Using the matrix notation:
⎡ ⎤ ⎡ ⎤ ⎡ ⎤
𝑦1 1 𝑥11 ... 𝑥1𝑃 ⎡ 𝛽 ⎤ 𝜖1
0
⎢𝑦2 ⎥ ⎢1 𝑥21 ... 𝑥2𝑃 ⎥ ⎢ 𝛽 ⎥ ⎢𝜖2 ⎥
⎥ ⎢
⎢ ⎥ ⎢
⎢𝑦3 ⎥ = ⎢1 ⎢ 1⎥ ⎢ ⎥
𝑥31 ... 𝑥3𝑃 ⎥
⎥⎢ .. ⎥ + ⎢𝜖3 ⎥
𝑥4𝑃 ⎦ .
⎢ ⎥ ⎢ ⎦ ⎣ ⎥
⎣𝑦4 ⎦ ⎣1 𝑥41 ...
⎣
𝜖4 ⎦
𝑦5 1 𝑥5 ... 𝑥5 𝛽 𝑃 𝜖5
Let 𝑋 = [𝑥𝑇0 , ..., 𝑥𝑇𝑁 ] be the (𝑁 × 𝑃 + 1) design matrix of 𝑁 samples of 𝑃 input features with
one column of one and let be 𝑦 = [𝑦1 , ..., 𝑦𝑁 ] be a vector of the 𝑁 targets.
𝑦 = 𝑋𝛽 + 𝜀
Using the matrix notation, the mean squared error (MSE) loss can be rewritten:
1
𝑀 𝑆𝐸(𝛽) = ||𝑦 − 𝑋𝛽||22 .
𝑁
The 𝛽 that minimises the MSE can be found by:
(︂ )︂
1
∇𝛽 ||𝑦 − 𝑋𝛽||22 =0 (4.25)
𝑁
1
∇𝛽 (𝑦 − 𝑋𝛽)𝑇 (𝑦 − 𝑋𝛽) = 0 (4.26)
𝑁
1
∇𝛽 (𝑦 𝑇 𝑦 − 2𝛽 𝑇 𝑋 𝑇 𝑦 + 𝛽 𝑇 𝑋 𝑇 𝑋𝛽) = 0 (4.27)
𝑁
−2𝑋 𝑇 𝑦 + 2𝑋 𝑇 𝑋𝛽 = 0 (4.28)
𝑇 𝑇
𝑋 𝑋𝛽 = 𝑋 𝑦 (4.29)
𝛽 = (𝑋 𝑇 𝑋)−1 𝑋 𝑇 𝑦, (4.30)
import numpy as np
from scipy import linalg
np.random.seed(seed=42) # make the example reproducible
# Dataset
N, P = 50, 4
X = np.random.normal(size= N * P).reshape((N, P))
## Our model needs an intercept so we add a column of 1s:
X[:, 0] = 1
print(X[:5, :])
Sources: http://statsmodels.sourceforge.net/devel/examples/
Multiple regression
import statsmodels.api as sm
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly␣
˓→specified.
Use R language syntax for data.frame. For an additive model: 𝑦𝑖 = 𝛽 0 + 𝑥1𝑖 𝛽 1 + 𝑥2𝑖 𝛽 2 + 𝜖𝑖 ≡ y ~
x1 + x2.
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly␣
˓→specified.
Analysis of covariance (ANCOVA) is a linear model that blends ANOVA and linear regression.
ANCOVA evaluates whether population means of a dependent variable (DV) are equal across
levels of a categorical independent variable (IV) often called a treatment, while statistically
controlling for the effects of other quantitative or continuous variables that are not of primary
interest, known as covariates (CV).
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
try:
df = pd.read_csv("../datasets/salary_table.csv")
except:
url = 'https://raw.github.com/neurospin/pystatsml/master/datasets/salary_table.csv'
df = pd.read_csv(url)
Normality assumption of the residuals can be rejected (p-value < 0.05). There is an efect of the
“management” factor, take it into account.
One-way AN(C)OVA
sum_sq df F PR(>F)
management 5.755739e+08 1.0 183.593466 4.054116e-17
experience 3.334992e+08 1.0 106.377768 3.349662e-13
Residual 1.348070e+08 43.0 NaN NaN
Jarque-Bera normality test p-value 0.004
Distribution of residuals is still not normal but closer to normality. Both management and
experience are significantly associated with salary.
Two-way AN(C)OVA
sum_sq df F PR(>F)
education 9.152624e+07 2.0 43.351589 7.672450e-11
management 5.075724e+08 1.0 480.825394 2.901444e-24
experience 3.380979e+08 1.0 320.281524 5.546313e-21
Residual 4.328072e+07 41.0 NaN NaN
Jarque-Bera normality test p-value 0.506
Normality assumtion cannot be rejected. Assume it. Education, management and experience
are significantly associated with salary.
oneway is nested within twoway. Comparing two nested models tells us if the additional predic-
tors (i.e. education) of the full model significantly decrease the residuals. Such comparison can
be done using an 𝐹 -test on residuals:
Factor coding
See http://statsmodels.sourceforge.net/devel/contrasts.html
By default Pandas use “dummy coding”. Explore:
print(twoway.model.data.param_names)
print(twoway.model.data.exog[:10, :])
import numpy as np
np.random.seed(seed=42) # make example reproducible
# Dataset
n_samples, n_features = 100, 1000
n_info = int(n_features/10) # number of features with information
n1, n2 = int(n_samples/2), n_samples - int(n_samples/2)
snr = .5
Y = np.random.randn(n_samples, n_features)
grp = np.array(["g1"] * n1 + ["g2"] * n2)
#
import scipy.stats as stats
import matplotlib.pyplot as plt
tvals, pvals = np.full(n_features, np.NAN), np.full(n_features, np.NAN)
for j in range(n_features):
tvals[j], pvals[j] = stats.ttest_ind(Y[grp=="g1", j], Y[grp=="g2", j],
equal_var=True)
axis[2].hist([pvals[n_info:], pvals[:n_info]],
stacked=True, bins=100, label=["Negatives", "Positives"])
axis[2].set_xlabel("p-value histogram")
axis[2].set_ylabel("density")
axis[2].legend()
plt.tight_layout()
Note that under the null hypothesis the distribution of the p-values is uniform.
Statistical measures:
• True Positive (TP) equivalent to a hit. The test correctly concludes the presence of an
effect.
• True Negative (TN). The test correctly concludes the absence of an effect.
• False Positive (FP) equivalent to a false alarm, Type I error. The test improperly con-
cludes the presence of an effect. Thresholding at 𝑝-value < 0.05 leads to 47 FP.
• False Negative (FN) equivalent to a miss, Type II error. The test improperly concludes the
absence of an effect.
The Bonferroni correction is based on the idea that if an experimenter is testing 𝑃 hypothe-
ses, then one way of maintaining the familywise error rate (FWER) is to test each individual
hypothesis at a statistical significance level of 1/𝑃 times the desired maximum overall level.
So, if the desired significance level for the whole family of tests is 𝛼 (usually 0.05), then the
Bonferroni correction would test each individual hypothesis at a significance level of 𝛼/𝑃 . For
example, if a trial is testing 𝑃 = 8 hypotheses with a desired 𝛼 = 0.05, then the Bonferroni
correction would test each individual hypothesis at 𝛼 = 0.05/8 = 0.00625.
FDR-controlling procedures are designed to control the expected proportion of rejected null
hypotheses that were incorrect rejections (“false discoveries”). FDR-controlling procedures pro-
vide less stringent control of Type I errors compared to the familywise error rate (FWER) con-
trolling procedures (such as the Bonferroni correction), which control the probability of at least
one Type I error. Thus, FDR-controlling procedures have greater power, at the cost of increased
rates of Type I errors.
4.1.12 Exercises
Load the dataset: birthwt Risk Factors Associated with Low Infant Birth Weight at https://raw.
github.com/neurospin/pystatsml/master/datasets/birthwt.csv
1. Test the association of mother’s (bwt) age and birth weight using the correlation test and
linear regeression.
2. Test the association of mother’s weight (lwt) and birth weight using the correlation test
and linear regeression.
3. Produce two scatter plot of: (i) age by birth weight; (ii) mother’s weight by birth weight.
Conclusion ?
Considering the salary and the experience of the salary table. https://raw.github.com/
neurospin/pystatsml/master/datasets/salary_table.csv
Compute:
• Estimate the model paramters 𝛽, 𝛽0 using scipy stats.linregress(x,y)
• Compute the predicted values 𝑦ˆ
Compute:
• 𝑦¯: y_mu
• 𝑆𝑆tot : ss_tot
• 𝑆𝑆reg : ss_reg
• 𝑆𝑆res : ss_res
• Check partition of variance formula based on sum of squares by using assert np.
allclose(val1, val2, atol=1e-05)
• Compute 𝑅2 and compare it with the r_value above
• Compute the 𝐹 score
• Compute the 𝑝-value:
• Plot the 𝐹 (1, 𝑛) distribution for 100 𝑓 values within [10, 25]. Draw 𝑃 (𝐹 (1, 𝑛) > 𝐹 ),
i.e. color the surface defined by the 𝑥 values larger than 𝐹 below the 𝐹 (1, 𝑛).
• 𝑃 (𝐹 (1, 𝑛) > 𝐹 ) is the 𝑝-value, compute it.
Multiple regression
import numpy as np
from scipy import linalg
np.random.seed(seed=42) # make the example reproducible
# Dataset
N, P = 50, 4
X = np.random.normal(size= N * P).reshape((N, P))
## Our model needs an intercept so we add a column of 1s:
X[:, 0] = 1
print(X[:5, :])
Given the following two sample, test whether their means are equals.
𝑦 =𝑔+𝜀
Where the noise 𝜀 ∼ 𝑁 (1, 1) and 𝑔 ∈ {0, 1} is a group indicator variable with 50 ones and 50
zeros.
• Write a function tstat(y, g) that compute the two samples t-test of y splited in two
groups defined by g.
• Sample the t-statistic distribution under the null hypothesis using random permutations.
• Assess the p-value.
Write a function univar_stat(df, target, variables) that computes the parametric statistics
and 𝑝-values between the target variable (provided as as string) and all variables (provided
as a list of string) of the pandas DataFrame df. The target is a quantitative variable but vari-
ables may be quantitative or qualitative. The function returns a DataFrame with four columns:
variable, test, value, p_value.
Apply it to the salary dataset available at https://raw.github.com/neurospin/pystatsml/master/
datasets/salary_table.csv, with target being S: salaries for IT staff in a corporation.
Multiple comparisons
This exercise has 2 goals: apply you knowledge of statistics using vectorized numpy operations.
Given the dataset provided for multiple comparisons, compute the two-sample 𝑡-test (assuming
equal variance) for each (column) feature of the Y array given the two groups defined by grp
variable. You should return two vectors of size n_features: one for the 𝑡-values and one for the
𝑝-values.
ANOVA
# dataset
mu_k = np.array([1, 2, 3]) # means of 3 samples
sd_k = np.array([1, 1, 1]) # sd of 3 samples
n_k = np.array([10, 20, 30]) # sizes of 3 samples
grp = [0, 1, 2] # group labels
n = np.sum(n_k)
label = np.hstack([[k] * n_k[k] for k in [0, 1, 2]])
(continues on next page)
y = np.zeros(n)
for k in grp:
y[label == k] = np.random.normal(mu_k[k], sd_k[k], n_k[k])
The study provides the brain volumes of grey matter (gm), white matter (wm) and cerebrospinal
fluid) (csf) of 808 anatomical MRI scans.
import os
import os.path
import pandas as pd
import tempfile
import urllib.request
WD = os.path.join(tempfile.gettempdir(), "brainvol")
os.makedirs(WD, exist_ok=True)
#os.chdir(WD)
Fetch data
• Demographic data demo.csv (columns: participant_id, site, group, age, sex) and tissue
volume data: group is Control or Patient. site is the recruiting site.
• Gray matter volume gm.csv (columns: participant_id, session, gm_vol)
• White matter volume wm.csv (columns: participant_id, session, wm_vol)
• Cerebrospinal Fluid csf.csv (columns: participant_id, session, csf_vol)
base_url = 'https://raw.github.com/neurospin/pystatsml/master/datasets/brain_volumes/%s'
data = dict()
for file in ["demo.csv", "gm.csv", "wm.csv", "csf.csv"]:
(continues on next page)
Out:
brain_vol = brain_vol.dropna()
assert brain_vol.shape == (766, 9)
import os
import pandas as pd
import seaborn as sns
import statsmodels.formula.api as smfrmla
import statsmodels.api as sm
(continues on next page)
Descriptive statistics Most of participants have several MRI sessions (column session) Select
on rows from session one “ses-01”
desc_glob_num = brain_vol1.describe()
print(desc_glob_num)
Out:
Out:
Out:
Out:
gm_vol
count meanstd min 25% 50% 75% max
group
Control 86.00 0.72 0.09 0.48 0.66 0.71 0.78 1.03
Patient 157.00 0.70 0.08 0.53 0.65 0.70 0.76 0.90
4.2.3 Statistics
Objectives:
1. Site effect of gray matter atrophy
2. Test the association between the age and gray matter atrophy in the control and patient
population independently.
3. Test for differences of atrophy between the patients and the controls
4. Test for interaction between age and clinical status, ie: is the brain atrophy process in
patient population faster than in the control population.
5. The effect of the medication in the patient population.
import statsmodels.api as sm
import statsmodels.formula.api as smfrmla
import scipy.stats
import seaborn as sns
Out:
/home/ed203246/anaconda3/lib/python3.7/site-packages/seaborn/_decorators.py:43:␣
˓→FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12,␣
˓→the only valid positional argument will be `data`, and passing other arguments without␣
˓→an explicit keyword will result in an error or misinterpretation.
FutureWarning
(continues on next page)
<AxesSubplot:xlabel='site', ylabel='gm_f'>
Out:
print(sm.stats.anova_lm(anova, typ=2))
Out:
2. Test the association between the age and gray matter atrophy in the control and patient
population independently.
Plot
Out:
/home/ed203246/anaconda3/lib/python3.7/site-packages/seaborn/_decorators.py:43:␣
˓→FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12,␣
˓→the only valid positional argument will be `data`, and passing other arguments without␣
Out:
Out:
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly␣
˓→specified.
Age explains 10.57% of the grey matter fraction variance
--- In patient population ---
OLS Regression Results
==============================================================================
(continues on next page)
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly␣
˓→specified.
Age explains 27.96% of the grey matter fraction variance
Before testing for differences of atrophy between the patients ans the controls Preliminary tests
for age x group effect (patients would be older or younger than Controls)
Plot
Out:
/home/ed203246/anaconda3/lib/python3.7/site-packages/seaborn/_decorators.py:43:␣
˓→FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12,␣
˓→the only valid positional argument will be `data`, and passing other arguments without␣
<AxesSubplot:xlabel='group', ylabel='age'>
Out:
Ttest_indResult(statistic=-1.2155557697674162, pvalue=0.225343592508479)
Out:
OLS Regression Results
==============================================================================
Dep. Variable: age R-squared: 0.006
(continues on next page)
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly␣
˓→specified.
No significant difference in age between patients and controls
Preliminary tests for sex x group (more/less males in patients than in Controls)
Out:
3. Test for differences of atrophy between the patients and the controls
Out:
sum_sq df F PR(>F)
group 0.00 1.00 0.01 0.92
Residual 0.46 241.00 nan nan
No significant difference in age between patients and controls
print(sm.stats.anova_lm(smfrmla.ols(
"gm_f ~ group + age + site", data=brain_vol1).fit(), typ=2))
print("No significant difference in age between patients and controls")
Out:
sum_sq df F PR(>F)
group 0.00 1.00 1.82 0.18
site 0.11 5.00 19.79 0.00
age 0.09 1.00 86.86 0.00
Residual 0.25 235.00 nan nan
No significant difference in age between patients and controls
4. Test for interaction between age and clinical status, ie: is the brain atrophy process in
patient population faster than in the control population.
print("%.3f%% of grey matter loss per year (almost %.1f%% per decade)" %\
(ancova.params.age * 100, ancova.params.age * 100 * 10))
Out:
sum_sq df F PR(>F)
site 0.11 5.00 20.28 0.00
age 0.10 1.00 89.37 0.00
group:age 0.00 1.00 3.28 0.07
Residual 0.25 235.00 nan nan
= Parameters =
Intercept 0.52
site[T.S3] 0.01
site[T.S4] 0.03
site[T.S5] 0.01
site[T.S7] 0.06
site[T.S8] 0.02
age -0.00
group[T.Patient]:age -0.00
dtype: float64
-0.148% of grey matter loss per year (almost -1.5% per decade)
grey matter loss in patients is accelerated by -0.232% per decade
Multivariate statistics includes all statistical techniques for analyzing samples made of two or
more variables. The data set (a 𝑁 × 𝑃 matrix X) is a collection of 𝑁 independent samples
column vectors [x1 , . . . , x𝑖 , . . . , x𝑁 ] of length 𝑃
⎡ 𝑇 ⎤ ⎡ ⎤ ⎡ ⎤
−x1 − 𝑥11 · · · 𝑥1𝑗 · · · 𝑥1𝑃 𝑥11 . . . 𝑥1𝑃
⎢ .. ⎥ ⎢ .. .. .. ⎥ ⎢ .. .. ⎥
⎢ . ⎥ ⎢ . . . ⎥⎥ ⎢ . . ⎥
⎢
⎢ 𝑇 ⎥ ⎢ ⎥
X = ⎢−x𝑖 −⎥ = ⎢ 𝑥𝑖1 · · · 𝑥𝑖𝑗 · · · 𝑥𝑖𝑃 ⎥ = ⎢
⎢ ⎥ ⎢ ⎥ ⎢ X ⎥ .
⎢ .. ⎥ ⎢ .. .. .. ⎥ ⎢ .. .. ⎥
⎥
⎣ . ⎦ ⎣ . . . ⎦ ⎣ . . ⎦
𝑇
−x𝑃 − 𝑥𝑁 1 · · · 𝑥𝑁 𝑗 · · · 𝑥𝑁 𝑃 𝑥𝑁 1 . . . 𝑥𝑁 𝑃 𝑁 ×𝑃
Source: Wikipedia
Algebraic definition
The dot product, denoted ’‘·” of two 𝑃 -dimensional vectors a = [𝑎1 , 𝑎2 , ..., 𝑎𝑃 ] and a =
[𝑏1 , 𝑏2 , ..., 𝑏𝑃 ] is defined as
⎡ ⎤
𝑏1
⎢ .. ⎥
]︀ ⎢ . ⎥
⎢ ⎥
∑︁
𝑇 𝑇
[︀
a·b=a b= 𝑎𝑖 𝑏𝑖 = 𝑎1 . . . a . . . 𝑎𝑃 ⎢ ⎢ b ⎥.
⎥
.
𝑖
⎣ .. ⎦
⎢ ⎥
𝑏𝑃
The Euclidean norm of a vector can be computed using the dot product, as
√
‖a‖2 = a · a.
In Euclidean space, a Euclidean vector is a geometrical object that possesses both a magnitude
and a direction. A vector can be pictured as an arrow. Its magnitude is its length, and its
direction is the direction that the arrow points. The magnitude of a vector a is denoted by ‖a‖2 .
The dot product of two Euclidean vectors a and b is defined by
a · b = 0.
At the other extreme, if they are codirectional, then the angle between them is 0° and
a · b = ‖a‖2 ‖b‖2
a · a = ‖a‖22 .
The scalar projection (or scalar component) of a Euclidean vector a in the direction of a Eu-
clidean vector b is given by
𝑎𝑏 = ‖a‖2 cos 𝜃,
Fig. 5: Projection.
import numpy as np
np.random.seed(42)
a = np.random.randn(10)
b = np.random.randn(10)
np.dot(a, b)
-4.085788532659924
𝑥𝑖𝑃 𝑥¯𝑃
• The covariance matrix ΣXX is a symmetric positive semi-definite matrix whose element
in the 𝑗, 𝑘 position is the covariance between the 𝑗 𝑡ℎ and 𝑘 𝑡ℎ elements of a random vector
i.e. the 𝑗 𝑡ℎ and 𝑘 𝑡ℎ columns of X.
• The covariance matrix generalizes the notion of covariance to multiple dimensions.
• The covariance matrix describe the shape of the sample distribution around the mean
assuming an elliptical distribution:
where
𝑁
1 1 ∑︁
𝑠𝑗𝑘 = 𝑠𝑘𝑗 = xj 𝑇 x k = 𝑥𝑖𝑗 𝑥𝑖𝑘
𝑁 −1 𝑁 −1
𝑖=1
np.random.seed(42)
colors = sns.color_palette()
# Generate dataset
for i in range(len(mean)):
X[i] = np.random.multivariate_normal(mean[i], Cov[i], n_samples)
# Plot
for i in range(len(mean)):
# Points
plt.scatter(X[i][:, 0], X[i][:, 1], color=colors[i], label="class %i" % i)
# Means
plt.scatter(mean[i][0], mean[i][1], marker="o", s=200, facecolors='w',
edgecolors=colors[i], linewidth=2)
# Ellipses representing the covariance matrices
pystatsml.plot_utils.plot_cov_ellipse(Cov[i], pos=mean[i], facecolor='none',
linewidth=2, edgecolor=colors[i])
plt.axis('equal')
_ = plt.legend(loc='upper left')
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
url = 'https://python-graph-gallery.com/wp-content/uploads/mtcars.csv'
df = pd.read_csv(url)
f, ax = plt.subplots(figsize=(5.5, 4.5))
cmap = sns.color_palette("RdBu_r", 11)
# Draw the heatmap with the mask and correct aspect ratio
_ = sns.heatmap(corr, mask=None, cmap=cmap, vmax=1, center=0,
square=True, linewidths=.5, cbar_kws={"shrink": .5})
lab=0
reordered = np.concatenate(clusters)
R = corr.loc[reordered, reordered]
f, ax = plt.subplots(figsize=(5.5, 4.5))
# Draw the heatmap with the mask and correct aspect ratio
_ = sns.heatmap(R, mask=None, cmap=cmap, vmax=1, center=0,
square=True, linewidths=.5, cbar_kws={"shrink": .5})
[['mpg', 'cyl', 'disp', 'hp', 'wt', 'qsec', 'vs', 'carb'], ['am', 'gear'], ['drat']]
In statistics, precision is the reciprocal of the variance, and the precision matrix is the matrix
inverse of the covariance matrix.
It is related to partial correlations that measures the degree of association between two vari-
ables, while controlling the effect of other variables.
import numpy as np
print(Pcor.round(2))
# Precision matrix:
[[ 6.79 -3.21 -3.21 0. 0. 0. ]
[-3.21 6.79 -3.21 0. 0. 0. ]
[-3.21 -3.21 6.79 0. 0. 0. ]
[ 0. -0. -0. 5.26 -4.74 -0. ]
[ 0. 0. 0. -4.74 5.26 0. ]
[ 0. 0. 0. 0. 0. 1. ]]
# Partial correlations:
[[ nan 0.47 0.47 -0. -0. -0. ]
[ nan nan 0.47 -0. -0. -0. ]
[ nan nan nan -0. -0. -0. ]
[ nan nan nan nan 0.9 0. ]
[ nan nan nan nan nan -0. ]
[ nan nan nan nan nan nan]]
• The Mahalanobis distance is a measure of the distance between two points x and 𝜇 where
the dispersion (i.e. the covariance structure) of the samples is taken into account.
• The dispersion is considered through covariance matrix.
This is formally expressed as
√︁
𝐷𝑀 (x, 𝜇) = (x − 𝜇)𝑇 Σ−1 (x − 𝜇).
Intuitions
• Distances along the principal directions of dispersion are contracted since they correspond
to likely dispersion of points.
• Distances othogonal to the principal directions of dispersion are dilated since they corre-
spond to unlikely dispersion of points.
For example
√
𝐷𝑀 (1) = 1𝑇 Σ−1 1.
ones = np.ones(Cov.shape[0])
d_euc = np.sqrt(np.dot(ones, ones))
d_mah = np.sqrt(np.dot(np.dot(ones, Prec), ones))
The first dot product that distances along the principal directions of dispersion are contracted:
print(np.dot(ones, Prec))
import numpy as np
import scipy
import matplotlib.pyplot as plt
import seaborn as sns
import pystatsml.plot_utils
%matplotlib inline
np.random.seed(40)
colors = sns.color_palette()
Covi = scipy.linalg.inv(Cov)
dm_m_x1 = scipy.spatial.distance.mahalanobis(mean, x1, Covi)
dm_m_x2 = scipy.spatial.distance.mahalanobis(mean, x2, Covi)
# Plot distances
vm_x1 = (x1 - mean) / d2_m_x1
vm_x2 = (x2 - mean) / d2_m_x2
jitter = .1
plt.plot([mean[0] - jitter, d2_m_x1 * vm_x1[0] - jitter],
[mean[1], d2_m_x1 * vm_x1[1]], color='k')
plt.plot([mean[0] - jitter, d2_m_x2 * vm_x2[0] - jitter],
[mean[1], d2_m_x2 * vm_x2[1]], color='k')
plt.legend(loc='lower right')
plt.text(-6.1, 3,
'Euclidian: d(m, x1) = %.1f<d(m, x2) = %.1f' % (d2_m_x1, d2_m_x2), color='k')
plt.text(-6.1, 3.5,
'Mahalanobis: d(m, x1) = %.1f>d(m, x2) = %.1f' % (dm_m_x1, dm_m_x2), color='r')
plt.axis('equal')
print('Euclidian d(m, x1) = %.2f < d(m, x2) = %.2f' % (d2_m_x1, d2_m_x2))
print('Mahalanobis d(m, x1) = %.2f > d(m, x2) = %.2f' % (dm_m_x1, dm_m_x2))
If the covariance matrix is the identity matrix, the Mahalanobis distance reduces to the Eu-
clidean distance. If the covariance matrix is diagonal, then the resulting distance measure is
called a normalized Euclidean distance.
More generally, the Mahalanobis distance is a measure of the distance between a point x and a
distribution 𝒩 (x|𝜇, Σ). It is a multi-dimensional generalization of the idea of measuring how
many standard deviations away x is from the mean. This distance is zero if x is at the mean,
and grows as x moves away from the mean: along each principal component axis, it measures
the number of standard deviations from x to the mean of the distribution.
The distribution, or probability density function (PDF) (sometimes just density), of a continuous
random variable is a function that describes the relative likelihood for this random variable to
take on a given value.
The multivariate normal distribution, or multivariate Gaussian distribution, of a 𝑃 -dimensional
random vector x = [𝑥1 , 𝑥2 , . . . , 𝑥𝑃 ]𝑇 is
1 1
𝒩 (x|𝜇, Σ) = exp{− (x − 𝜇)𝑇 Σ−1 (x − 𝜇)}.
(2𝜋)𝑃/2 |Σ|1/2 2
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D
# x, y grid
x, y = np.mgrid[-3:3:.1, -3:3:.1]
X = np.stack((x.ravel(), y.ravel())).T
norm = multivariate_normal_pdf(X, mean, sigma).reshape(x.shape)
# Do it with scipy
norm_scpy = multivariate_normal(mu, sigma).pdf(np.stack((x, y), axis=2))
assert np.allclose(norm, norm_scpy)
# Plot
fig = plt.figure(figsize=(10, 7))
ax = fig.gca(projection='3d')
surf = ax.plot_surface(x, y, norm, rstride=3,
cstride=3, cmap=plt.cm.coolwarm,
linewidth=1, antialiased=False
)
ax.set_zlim(0, 0.2)
ax.zaxis.set_major_locator(plt.LinearLocator(10))
ax.zaxis.set_major_formatter(plt.FormatStrFormatter('%.02f'))
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('p(x)')
4.3.8 Exercises
Two libraries:
• Pandas: https://pandas.pydata.org/pandas-docs/stable/timeseries.html
• scipy http://www.statsmodels.org/devel/tsa.html
4.4.1 Stationarity
A TS is said to be stationary if its statistical properties such as mean, variance remain constant
over time.
• constant mean
• constant variance
• an autocovariance that does not depend on time.
what is making a TS non-stationary. There are 2 major reasons behind non-stationaruty of a
TS:
1. Trend – varying mean over time. For eg, in this case we saw that on average, the number
of passengers was growing over time.
2. Seasonality – variations at specific time-frames. eg people might have a tendency to buy
cars in a particular month because of pay increment or festivals.
import pandas as pd
import numpy as np
# String as index
prices = {'apple': 4.99,
'banana': 1.99,
'orange': 3.99}
ser = pd.Series(prices)
print(ser)
0 1
1 3
dtype: int64
apple 4.99
banana 1.99
orange 3.99
dtype: float64
a 1
b 2
dtype: int64
2
source: https://www.datacamp.com/community/tutorials/time-series-analysis-tutorial
Get Google Trends data of keywords such as ‘diet’ and ‘gym’ and see how they vary over time
while learning about trends and seasonality in time series data.
In the Facebook Live code along session on the 4th of January, we checked out Google trends
data of keywords ‘diet’, ‘gym’ and ‘finance’ to see how they vary over time. We asked ourselves
if there could be more searches for these terms in January when we’re all trying to turn over a
new leaf?
In this tutorial, you’ll go through the code that we put together during the session step by step.
You’re not going to do much mathematics but you are going to do the following:
• Read data
• Recode data
• Exploratory Data Analysis
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
try:
url = "https://raw.githubusercontent.com/datacamp/datacamp_facebook_live_ny_
˓→resolution/master/datasets/multiTimeline.csv"
df = pd.read_csv(url, skiprows=2)
except:
df = pd.read_csv("../datasets/multiTimeline.csv", skiprows=2)
print(df.head())
# Rename columns
df.columns = ['month', 'diet', 'gym', 'finance']
# Describe
print(df.describe())
Next, you’ll turn the ‘month’ column into a DateTime data type and make it the index of the
DataFrame.
Note that you do this because you saw in the result of the .info() method that the ‘Month’
column was actually an of data type object. Now, that generic data type encapsulates everything
from strings to integers, etc. That’s not exactly what you want when you want to be looking
at time series data. That’s why you’ll use .to_datetime() to convert the ‘month’ column in your
DataFrame to a DateTime.
Be careful! Make sure to include the inplace argument when you’re setting the index of the
DataFrame df so that you actually alter the original index and set it to the ‘month’ column.
df.month = pd.to_datetime(df.month)
df.set_index('month', inplace=True)
print(df.head())
You can use a built-in pandas visualization method .plot() to plot your data as 3 line plots on a
single figure (one for each column, namely, ‘diet’, ‘gym’, and ‘finance’).
df.plot()
plt.xlabel('Year');
Text(0.5, 0, 'Year')
Note that this data is relative. As you can read on Google trends:
Numbers represent search interest relative to the highest point on the chart for the given region
and time. A value of 100 is the peak popularity for the term. A value of 50 means that the term
is half as popular. Likewise a score of 0 means the term was less than 1% as popular as the
peak.
Rolling average, for each time point, take the average of the points on either side of it. Note
that the number of points is specified by a window size.
Remove Seasonality with pandas Series.
See: http://pandas.pydata.org/pandas-docs/stable/timeseries.html A: ‘year end frequency’
year frequency
diet = df['diet']
diet_resamp_yr = diet.resample('A').mean()
diet_roll_yr = diet.rolling(12).mean()
<matplotlib.legend.Legend at 0x7f2a045af0d0>
x = np.asarray(df[['diet']])
win = 12
(continues on next page)
[<matplotlib.lines.Line2D at 0x7f2a0021ee90>]
gym = df['gym']
Text(0.5, 0, 'Year')
Detrending
Text(0.5, 0, 'Year')
df.diff().plot()
plt.xlabel('Year')
Text(0.5, 0, 'Year')
df.plot()
plt.xlabel('Year');
print(df.corr())
sns.heatmap(df.corr(), cmap="coolwarm")
<AxesSubplot:>
‘diet’ and ‘gym’ are negatively correlated! Remember that you have a seasonal and a trend
component. From the correlation coefficient, ‘diet’ and ‘gym’ are negatively correlated:
• trends components are negatively correlated.
• seasonal components would positively correlated and their
The actual correlation coefficient is actually capturing both of those.
df.diff().plot()
plt.xlabel('Year');
print(df.diff().corr())
sns.heatmap(df.diff().corr(), cmap="coolwarm")
<AxesSubplot:>
x = gym
plt.subplot(411)
plt.plot(x, label='Original')
plt.legend(loc='best')
plt.subplot(412)
plt.plot(trend, label='Trend')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(seasonal,label='Seasonality')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(residual, label='Residuals')
plt.legend(loc='best')
plt.tight_layout()
4.4.10 Autocorrelation
A time series is periodic if it repeats itself at equally spaced intervals, say, every 12 months.
Autocorrelation Function (ACF): It is a measure of the correlation between the TS with a lagged
version of itself. For instance at lag 5, ACF would compare series at time instant t1. . . t2 with
series at instant t1-5. . . t2-5 (t1-5 and t2 being end points).
Plot
x = df["diet"].astype(float)
autocorrelation_plot(x)
<AxesSubplot:xlabel='Lag', ylabel='Autocorrelation'>
/home/ed203246/anaconda3/lib/python3.7/site-packages/statsmodels/tsa/stattools.py:666:␣
˓→FutureWarning: fft=True will become the default after the release of the 0.12 release␣
˓→of statsmodels. To suppress this warning, explicitly set fft=False.
FutureWarning,
ACF peaks every 12 months: Time series is correlated with itself shifted by 12 months.
4.4.11 Time Series Forecasting with Python using Autoregressive Moving Average
(ARMA) models
Source:
• https://www.packtpub.com/mapt/book/big_data_and_business_intelligence/
9781783553358/7/ch07lvl1sec77/arma-models
• http://en.wikipedia.org/wiki/Autoregressive%E2%80%93moving-average_model
• ARIMA: https://www.analyticsvidhya.com/blog/2016/02/
time-series-forecasting-codes-python/
ARMA models are often used to forecast a time series. These models combine autoregressive
and moving average models. In moving average models, we assume that a variable is the sum
of the mean of the time series and a linear combination of noise components.
The autoregressive and moving average models can have different orders. In general, we can
define an ARMA model with p autoregressive terms and q moving average terms as follows:
𝑝
∑︁ 𝑞
∑︁
𝑥𝑡 = 𝑎𝑖 𝑥𝑡−𝑖 + 𝑏𝑖 𝜀𝑡−𝑖 + 𝜀𝑡
𝑖 𝑖
Choosing p and q
Plot the partial autocorrelation functions for an estimate of p, and likewise using the autocorre-
lation functions for an estimate of q.
Partial Autocorrelation Function (PACF): This measures the correlation between the TS with a
lagged version of itself but after eliminating the variations already explained by the intervening
comparisons. Eg at lag 5, it will check the correlation but remove the effects already explained
by lags 1 to 4.
x = df["gym"].astype(float)
#Plot ACF:
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(x_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(x_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function (q=1)')
#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(x_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(x_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function (p=1)')
plt.tight_layout()
/home/ed203246/anaconda3/lib/python3.7/site-packages/statsmodels/tsa/stattools.py:666:␣
˓→FutureWarning: fft=True will become the default after the release of the 0.12 release␣
˓→of statsmodels. To suppress this warning, explicitly set fft=False.
FutureWarning,
In this plot, the two dotted lines on either sides of 0 are the confidence interevals. These can be
used to determine the p and q values as:
• p: The lag value where the PACF chart crosses the upper confidence interval for the first
time, in this case p=1.
• q: The lag value where the ACF chart crosses the upper confidence interval for the first
time, in this case q=1.
1. Define the model by calling ARMA() and passing in the p and q parameters.
2. The model is prepared on the training data by calling the fit() function.
3. Predictions can be made by calling the predict() function and specifying the index of the
time or times to be predicted.
print(model.summary())
plt.plot(x)
plt.plot(model.predict(), color='red')
plt.title('RSS: %.4f'% sum((model.fittedvalues-x)**2))
/home/ed203246/anaconda3/lib/python3.7/site-packages/statsmodels/tsa/arima_model.py:472:␣
˓→FutureWarning:
To silence this warning and continue using ARMA and ARIMA until they are
removed, use:
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
FutureWarning)
warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)
/home/ed203246/anaconda3/lib/python3.7/site-packages/statsmodels/tsa/base/tsa_model.
˓→py:527: ValueWarning: No frequency information was provided, so inferred frequency MS␣
˓→will be used.
% freq, ValueWarning)
FIVE
MACHINE LEARNING
5.1.1 Introduction
In machine learning and statistics, dimensionality reduction or dimension reduction is the pro-
cess of reducing the number of features under consideration, and can be divided into feature
selection (not addressed here) and feature extraction.
Feature extraction starts from an initial set of measured data and builds derived values (fea-
tures) intended to be informative and non-redundant, facilitating the subsequent learning and
generalization steps, and in some cases leading to better human interpretations. Feature extrac-
tion is related to dimensionality reduction.
The input matrix X, of dimension 𝑁 × 𝑃 , is
⎡ ⎤
𝑥11 . . . 𝑥1𝑃
⎢ ⎥
⎢ .. .. ⎥
⎢ ⎥
⎢ . X . ⎥
⎢ ⎥
⎣ ⎦
𝑥𝑁 1 . . . 𝑥 𝑁 𝑃
where the rows represent the samples and columns represent the variables.
The goal is to learn a transformation that extracts a few relevant features. This is generally
done by exploiting the covariance ΣXX between the input features.
Decompose the data matrix X𝑁 ×𝑃 into a product of a mixing matrix U𝑁 ×𝐾 and a dictionary
matrix V𝑃 ×𝐾 .
X = UV𝑇 ,
X ≈ X̂ = UV𝑇 ,
159
Statistics and Machine Learning in Python, Release 0.3 beta
X = UDV𝑇 ,
where
⎡ ⎤ ⎡ ⎤
𝑥11 𝑥1𝑃 𝑢11⎡ 𝑢1𝐾
⎤⎡ ⎤
⎢ ⎥ ⎢ ⎥ 𝑑1 0 𝑣11 𝑣1𝑃
⎢ ⎥ ⎢ ⎥
⎢
⎢ X ⎥=⎢
⎥ ⎢ U ⎥
⎥
⎣ D ⎦ ⎣ V𝑇 ⎦.
⎣ ⎦ ⎣ ⎦ 0 𝑑𝐾 𝑣𝐾1 𝑣𝐾𝑃
𝑥𝑁 1 𝑥𝑁 𝑃 𝑢𝑁 1 𝑢𝑁 𝐾
U: right-singular
• V = [v1 , · · · , v𝐾 ] is a 𝑃 × 𝐾 orthogonal matrix.
• It is a dictionary of patterns to be combined (according to the mixing coefficients) to
reconstruct the original samples.
• V perfoms the initial rotations (projection) along the 𝐾 = min(𝑁, 𝑃 ) principal compo-
nent directions, also called loadings.
• Each v𝑗 performs the linear combination of the variables that has maximum sample vari-
ance, subject to being uncorrelated with the previous v𝑗−1 .
D: singular values
• D is a 𝐾 × 𝐾 diagonal matrix made of the singular values of X with 𝑑1 ≥ 𝑑2 ≥ · · · ≥
𝑑𝐾 ≥ 0.
• D scale the projection along the coordinate axes by 𝑑1 , 𝑑2 , · · · , 𝑑𝐾 .
V transforms correlated variables (X) into a set of uncorrelated ones (UD) that better expose
the various relationships among the original data items.
X = UDV𝑇 , (5.1)
𝑇
XV = UDV V, (5.2)
XV = UDI, (5.3)
XV = UD (5.4)
At the same time, SVD is a method for identifying and ordering the dimensions along which
data points exhibit the most variation.
import numpy as np
import scipy
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
np.random.seed(42)
# dataset
n_samples = 100
experience = np.random.normal(size=n_samples)
salary = 1500 + experience + np.random.normal(size=n_samples, scale=.5)
X = np.column_stack([experience, salary])
plt.figure(figsize=(9, 3))
plt.subplot(131)
plt.scatter(U[:, 0], U[:, 1], s=50)
plt.axis('equal')
plt.title("U: Rotated and scaled data")
plt.subplot(132)
# Project data
PC = np.dot(X, Vh.T)
plt.scatter(PC[:, 0], PC[:, 1], s=50)
plt.axis('equal')
plt.title("XV: Rotated data")
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.subplot(133)
plt.scatter(X[:, 0], X[:, 1], s=50)
for i in range(Vh.shape[0]):
plt.arrow(x=0, y=0, dx=Vh[i, 0], dy=Vh[i, 1], head_width=0.2,
head_length=0.2, linewidth=2, fc='r', ec='r')
plt.text(Vh[i, 0], Vh[i, 1],'v%i' % (i+1), color="r", fontsize=15,
horizontalalignment='right', verticalalignment='top')
plt.axis('equal')
plt.ylim(-4, 4)
plt.tight_layout()
Sources:
• C. M. Bishop Pattern Recognition and Machine Learning, Springer, 2006
• Everything you did and didn’t know about PCA
• Principal Component Analysis in 3 Simple Steps
Principles
• Principal components analysis is the main method used for linear dimension reduction.
• The idea of principal component analysis is to find the 𝐾 principal components di-
rections (called the loadings) V𝐾×𝑃 that capture the variation in the data as much as
possible.
• It converts a set of 𝑁 𝑃 -dimensional observations N𝑁 ×𝑃 of possibly correlated variables
into a set of 𝑁 𝐾-dimensional samples C𝑁 ×𝐾 , where the 𝐾 < 𝑃 . The new variables are
linearly uncorrelated. The columns of C𝑁 ×𝐾 are called the principal components.
• The dimension reduction is obtained by using only 𝐾 < 𝑃 components that exploit corre-
lation (covariance) among the original variables.
• PCA is mathematically defined as an orthogonal linear transformation V𝐾×𝑃 that trans-
forms the data to a new coordinate system such that the greatest variance by some projec-
tion of the data comes to lie on the first coordinate (called the first principal component),
the second greatest variance on the second coordinate, and so on.
C𝑁 ×𝐾 = X𝑁 ×𝑃 V𝑃 ×𝐾
• PCA can be thought of as fitting a 𝑃 -dimensional ellipsoid to the data, where each axis of
the ellipsoid represents a principal component. If some axis of the ellipse is small, then the
variance along that axis is also small, and by omitting that axis and its corresponding prin-
cipal component from our representation of the dataset, we lose only a commensurately
small amount of information.
• Finding the 𝐾 largest axes of the ellipse will permit to project the data onto a space having
dimensionality 𝐾 < 𝑃 while maximizing the variance of the projected data.
Dataset preprocessing
Centering
Consider a data matrix, X , with column-wise zero empirical mean (the sample mean of each
column has been shifted to zero), ie. X is replaced by X − 1x̄𝑇 .
Standardizing
Optionally, standardize the columns, i.e., scale them by their standard-deviation. Without stan-
dardization, a variable with a high variance will capture most of the effect of the PCA. The
principal direction will be aligned with this variable. Standardization will, however, raise noise
variables to the save level as informative variables.
The covariance matrix of centered standardized data is the correlation matrix.
To begin with, consider the projection onto a one-dimensional space (𝐾 = 1). We can define
the direction of this space using a 𝑃 -dimensional vector v, which for convenience (and without
loss of generality) we shall choose to be a unit vector so that ‖v‖2 = 1 (note that we are only
interested in the direction defined by v, not in the magnitude of v itself). PCA consists of two
mains steps:
Projection in the directions that capture the greatest variance
Each 𝑃 -dimensional data point x𝑖 is then projected onto v, where the coordinate (in the co-
ordinate system of v) is a scalar value, namely x𝑇𝑖 v. I.e., we want to find the vector v that
maximizes these coordinates along v, which we will see corresponds to maximizing the vari-
ance of the projected data. This is equivalently expressed as
1 ∑︁ (︀ 𝑇 )︀2
v = arg max x𝑖 v .
‖v‖=1 𝑁
𝑖
where SXX is a biased estiamte of the covariance matrix of the data, i.e.
1 𝑇
SXX = X X.
𝑁
We now maximize the projected variance v𝑇 SXX v with respect to v. Clearly, this has to be a
constrained maximization to prevent ‖v2 ‖ → ∞. The appropriate constraint comes from the
normalization condition ‖v‖2 ≡ ‖v‖22 = v𝑇 v = 1. To enforce this constraint, we introduce a
Lagrange multiplier that we shall denote by 𝜆, and then make an unconstrained maximization
of
By setting the gradient with respect to v equal to zero, we see that this quantity has a stationary
point when
SXX v = 𝜆v.
v𝑇 SXX v = 𝜆,
and so the variance will be at a maximum when v is equal to the eigenvector corresponding to
the largest eigenvalue, 𝜆. This eigenvector is known as the first principal component.
We can define additional principal components in an incremental fashion by choosing each new
direction to be that which maximizes the projected variance amongst all possible directions that
are orthogonal to those already considered. If we consider the general case of a 𝐾-dimensional
projection space, the optimal linear projection for which the variance of the projected data is
maximized is now defined by the 𝐾 eigenvectors, v1 , . . . , vK , of the data covariance matrix
SXX that corresponds to the 𝐾 largest eigenvalues, 𝜆1 ≥ 𝜆2 ≥ · · · ≥ 𝜆𝐾 .
Back to SVD
X𝑇 X = (UDV𝑇 )𝑇 (UDV𝑇 )
= VD𝑇 U𝑇 UDV𝑇
= VD2 V𝑇
V𝑇 X𝑇 XV = D2
1 1
V𝑇 X𝑇 XV = D2
𝑁 −1 𝑁 −1
1
V𝑇 SXX V = D2
𝑁 −1
.
Considering only the 𝑘 𝑡ℎ right-singular vectors v𝑘 associated to the singular value 𝑑𝑘
1
vk 𝑇 SXX vk = 𝑑2 ,
𝑁 −1 𝑘
It turns out that if you have done the singular value decomposition then you already have
the Eigenvalue decomposition for X𝑇 X. Where - The eigenvectors of SXX are equivalent to
the right singular vectors, V, of X. - The eigenvalues, 𝜆𝑘 , of SXX , i.e. the variances of the
components, are equal to 𝑁 1−1 times the squared singular values, 𝑑𝑘 .
Moreover computing PCA with SVD do not require to form the matrix X𝑇 X, so computing the
SVD is now the standard way to calculate a principal components analysis from a data matrix,
unless only a handful of components are required.
PCA outputs
The SVD or the eigendecomposition of the data covariance matrix provides three main quanti-
ties:
1. Principal component directions or loadings are the eigenvectors of X𝑇 X. The V𝐾×𝑃
or the right-singular vectors of an SVD of X are called principal component directions of
X. They are generally computed using the SVD of X.
C𝑁 ×𝐾 = UDV𝑇𝑁 ×𝑃 V𝑃 ×𝐾 (5.5)
C𝑁 ×𝐾 = UD𝑇𝑁 ×𝐾 I𝐾×𝐾 (5.6)
C𝑁 ×𝐾 = UD𝑇𝑁 ×𝐾 (5.7)
(5.8)
Thus c𝑗 = Xv𝑗 = u𝑗 𝑑𝑗 , for 𝑗 = 1, . . . 𝐾. Hence u𝑗 is simply the projection of the row vectors of
X, i.e., the input predictor vectors, on the direction v𝑗 , scaled by 𝑑𝑗 .
⎡ ⎤
𝑥1,1 𝑣1,1 + . . . + 𝑥1,𝑃 𝑣1,𝑃
⎢ 𝑥2,1 𝑣1,1 + . . . + 𝑥2,𝑃 𝑣1,𝑃 ⎥
c1 = ⎢ ..
⎢ ⎥
⎥
⎣ . ⎦
𝑥𝑁,1 𝑣1,1 + . . . + 𝑥𝑁,𝑃 𝑣1,𝑃
1
𝑣𝑎𝑟(c𝑘 ) = (Xv𝑘 )2 (5.9)
𝑁 −1
1
= (u𝑘 𝑑𝑘 )2 (5.10)
𝑁 −1
1
= 𝑑2 (5.11)
𝑁 −1 𝑘
We must choose 𝐾 * ∈ [1, . . . , 𝐾], the number of required components. This can be done by
calculating the explained variance ratio of the 𝐾 * first components and by choosing 𝐾 * such
that the cumulative explained variance ratio is greater than some given threshold (e.g., ≈
90%). This is expressed as
∑︀𝐾 *
𝑗 𝑣𝑎𝑟(c𝑘 )
cumulative explained variance(c𝑘 ) = ∑︀𝐾 .
𝑗 𝑣𝑎𝑟(c𝑘 )
PCs
Plot the samples projeted on first the principal components as e.g. PC1 against PC2.
PC directions
Exploring the loadings associated with a component provides the contribution of each original
variable in the component.
Remark: The loadings (PC directions) are the coefficients of multiple regression of PC on origi-
nal variables:
c = Xv (5.12)
X𝑇 c = X𝑇 Xv (5.13)
𝑇 −1 𝑇
(X X) X c=v (5.14)
Another way to evaluate the contribution of the original variables in each PC can be obtained
by computing the correlation between the PCs and the original variables, i.e. columns of X,
denoted x𝑗 , for 𝑗 = 1, . . . , 𝑃 . For the 𝑘 𝑡ℎ PC, compute and plot the correlations with all original
variables
𝑐𝑜𝑟(c𝑘 , x𝑗 ), 𝑗 = 1 . . . 𝐾, 𝑗 = 1 . . . 𝐾.
import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
np.random.seed(42)
# dataset
n_samples = 100
experience = np.random.normal(size=n_samples)
salary = 1500 + experience + np.random.normal(size=n_samples, scale=.5)
X = np.column_stack([experience, salary])
PC = pca.transform(X)
plt.subplot(121)
plt.scatter(X[:, 0], X[:, 1])
plt.xlabel("x1"); plt.ylabel("x2")
plt.subplot(122)
plt.scatter(PC[:, 0], PC[:, 1])
plt.xlabel("PC1 (var=%.2f)" % pca.explained_variance_ratio_[0])
plt.ylabel("PC2 (var=%.2f)" % pca.explained_variance_ratio_[1])
plt.axis('equal')
plt.tight_layout()
[0.93646607 0.06353393]
Resources:
• http://www.stat.pitt.edu/sungkyu/course/2221Fall13/lec8_mds_combined.pdf
• https://en.wikipedia.org/wiki/Multidimensional_scaling
• Hastie, Tibshirani and Friedman (2009). The Elements of Statistical Learning: Data Mining,
Inference, and Prediction. New York: Springer, Second Edition.
The purpose of MDS is to find a low-dimensional projection of the data in which the pairwise
distances between data points is preserved, as closely as possible (in a least-squares sense).
• Let D be the (𝑁 × 𝑁 ) pairwise distance matrix where 𝑑𝑖𝑗 is a distance between points 𝑖
and 𝑗.
• The MDS concept can be extended to a wide variety of data types specified in terms of a
similarity matrix.
Given the dissimilarity (distance) matrix D𝑁 ×𝑁 = [𝑑𝑖𝑗 ], MDS attempts to find 𝐾-dimensional
projections of the 𝑁 points x1 , . . . , x𝑁 ∈ R𝐾 , concatenated in an X𝑁 ×𝐾 matrix, so that 𝑑𝑖𝑗 ≈
‖x𝑖 − x𝑗 ‖ are as close as possible. This can be obtained by the minimization of a loss function
called the stress function
∑︁
stress(X) = (𝑑𝑖𝑗 − ‖x𝑖 − x𝑗 ‖)2 .
𝑖̸=𝑗
The Sammon mapping performs better at preserving small distances compared to the least-
squares scaling.
Example
The eurodist datset provides the road distances (in kilometers) between 21 cities in Europe.
Given this matrix of pairwise (non-Euclidean) distances D = [𝑑𝑖𝑗 ], MDS can be used to recover
the coordinates of the cities in some Euclidean referential whose orientation is arbitrary.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
print(df.iloc[:5, :5])
city = df["city"]
D = np.array(df.iloc[:, 1:]) # Distance matrix
for i in range(len(city)):
plt.text(Xr[i, 0], Xr[i, 1], city[i])
plt.axis('equal')
(-1894.0919178069155,
2914.3554370871234,
-1712.9733697197494,
2145.437068788015)
We must choose 𝐾 * ∈ {1, . . . , 𝐾} the number of required components. Plotting the values of
the stress function, obtained using 𝑘 ≤ 𝑁 − 1 components. In general, start with 1, . . . 𝐾 ≤ 4.
Choose 𝐾 * where you can clearly distinguish an elbow in the stress curve.
Thus, in the plot below, we choose to retain information accounted for by the first two compo-
nents, since this is where the elbow is in the stress curve.
print(stress)
plt.plot(k_range, stress)
plt.xlabel("k")
plt.ylabel("stress")
Sources:
• Scikit-learn documentation
• Wikipedia
Nonlinear dimensionality reduction or manifold learning cover unsupervised methods that
attempt to identify low-dimensional manifolds within the original 𝑃 -dimensional space that
represent high data density. Then those methods provide a mapping from the high-dimensional
space to the low-dimensional embedding.
Isomap
ax = fig.add_subplot(121, projection='3d')
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=color, cmap=plt.cm.Spectral)
ax.view_init(4, -72)
plt.title('2D "S shape" manifold in 3D')
Y = manifold.Isomap(n_neighbors=10, n_components=2).fit_transform(X)
ax = fig.add_subplot(122)
plt.scatter(Y[:, 0], Y[:, 1], c=color, cmap=plt.cm.Spectral)
plt.title("Isomap")
plt.xlabel("First component")
plt.ylabel("Second component")
plt.axis('tight')
5.1.6 Exercises
PCA
MDS
5.2 Clustering
Wikipedia: Cluster analysis or clustering is the task of grouping a set of objects in such a way
that objects in the same group (called a cluster) are more similar (in some sense or another)
to each other than to those in other groups (clusters). Clustering is one of the main task of
exploratory data mining, and a common technique for statistical data analysis, used in many
fields, including machine learning, pattern recognition, image analysis, information retrieval,
and bioinformatics.
Sources: http://scikit-learn.org/stable/modules/clustering.html
which represents the sum of the squares of the Euclidean distances of each data point to its
assigned vector 𝜇𝑘 . Our goal is to find values for the {𝑟𝑖𝑘 } and the {𝜇𝑘 } so as to minimize the
function 𝐽. We can do this through an iterative procedure in which each iteration involves two
successive steps corresponding to successive optimizations with respect to the 𝑟𝑖𝑘 and the 𝜇𝑘
. First we choose some initial values for the 𝜇𝑘 . Then in the first phase we minimize 𝐽 with
respect to the 𝑟𝑖𝑘 , keeping the 𝜇𝑘 fixed. In the second phase we minimize 𝐽 with respect to
the 𝜇𝑘 , keeping 𝑟𝑖𝑘 fixed. This two-stage optimization process is then repeated until conver-
gence. We shall see that these two stages of updating 𝑟𝑖𝑘 and 𝜇𝑘 correspond respectively to the
expectation (E) and maximization (M) steps of the expectation-maximisation (EM) algorithm,
and to emphasize this we shall use the terms E step and M step in the context of the 𝐾-means
algorithm.
Consider first the determination of the 𝑟𝑖𝑘 . Because 𝐽 in is a linear function of 𝑟𝑖𝑘 , this opti-
mization can be performed easily to give a closed form solution. The terms involving different
𝑖 are independent and so we can optimize for each 𝑖 separately by choosing 𝑟𝑖𝑘 to be 1 for
whichever value of 𝑘 gives the minimum value of ||𝑥𝑖 − 𝜇𝑘 ||2 . In other words, we simply assign
the 𝑖th data point to the closest cluster centre. More formally, this can be expressed as
{︃
1, if 𝑘 = arg min𝑗 ||𝑥𝑖 − 𝜇𝑗 ||2 .
𝑟𝑖𝑘 = (5.15)
0, otherwise.
Now consider the optimization of the 𝜇𝑘 with the 𝑟𝑖𝑘 held fixed. The objective function 𝐽 is a
quadratic function of 𝜇𝑘 , and it can be minimized by setting its derivative with respect to 𝜇𝑘 to
zero giving
∑︁
2 𝑟𝑖𝑘 (𝑥𝑖 − 𝜇𝑘 ) = 0
𝑖
The denominator in this expression is equal to the number of points assigned to cluster 𝑘, and so
this result has a simple interpretation, namely set 𝜇𝑘 equal to the mean of all of the data points
𝑥𝑖 assigned to cluster 𝑘. For this reason, the procedure is known as the 𝐾-means algorithm.
The two phases of re-assigning data points to clusters and re-computing the cluster means are
repeated in turn until there is no further change in the assignments (or until some maximum
number of iterations is exceeded). Because each phase reduces the value of the objective func-
tion 𝐽, convergence of the algorithm is assured. However, it may converge to a local rather than
global minimum of 𝐽.
iris = datasets.load_iris()
X = iris.data[:, :2] # use only 'sepal length and sepal width'
y_iris = iris.target
km2 = cluster.KMeans(n_clusters=2).fit(X)
km3 = cluster.KMeans(n_clusters=3).fit(X)
km4 = cluster.KMeans(n_clusters=4).fit(X)
plt.figure(figsize=(9, 3))
plt.subplot(131)
plt.scatter(X[:, 0], X[:, 1], c=km2.labels_)
plt.title("K=2, J=%.2f" % km2.inertia_)
plt.subplot(132)
plt.scatter(X[:, 0], X[:, 1], c=km3.labels_)
plt.title("K=3, J=%.2f" % km3.inertia_)
plt.subplot(133)
plt.scatter(X[:, 0], X[:, 1], c=km4.labels_)#.astype(np.float))
plt.title("K=4, J=%.2f" % km4.inertia_)
Exercises
1. Analyse clusters
• Analyse the plot above visually. What would a good value of 𝐾 be?
• If you instead consider the inertia, the value of 𝐽, what would a good value of 𝐾 be?
• Explain why there is such difference.
• For 𝐾 = 2 why did 𝐾-means clustering not find the two “natural” clusters? See
the assumptions of 𝐾-means: http://scikit-learn.org/stable/auto_examples/cluster/plot_
kmeans_assumptions.html#example-cluster-plot-kmeans-assumptions-py
Write a function kmeans(X, K) that return an integer vector of the samples’ labels.
The Gaussian mixture model (GMM) is a simple linear superposition of Gaussian components
over the data, aimed at providing a rich class of density models. We turn to a formulation of
Gaussian mixtures in terms of discrete latent variables: the 𝐾 hidden classes to be discovered.
Differences compared to 𝐾-means:
• Whereas the 𝐾-means algorithm performs a hard assignment of data points to clusters, in
which each data point is associated uniquely with one cluster, the GMM algorithm makes
a soft assignment based on posterior probabilities.
• Whereas the classic 𝐾-means is only based on Euclidean distances, classic GMM use a
Mahalanobis distances that can deal with non-spherical distributions. It should be noted
that Mahalanobis could be plugged within an improved version of 𝐾-Means clustering.
The Mahalanobis distance is unitless and scale-invariant, and takes into account the cor-
relations of the data set.
where:
• The 𝑝(𝑘) are ∑︀
the mixing coefficients also know as the class probability of class 𝑘, and they
sum to one: 𝐾 𝑘=1 𝑝(𝑘) = 1.
To compute the classes parameters: 𝑝(𝑘), 𝜇𝑘 , Σ𝑘 we sum over all samples, by weighting each
sample 𝑖 by its responsibility or contribution
∑︀ to class 𝑘: 𝑝(𝑘 | 𝑥𝑖 ) such that for each point its
contribution to all classes sum to one 𝑘 𝑝(𝑘 | 𝑥𝑖 ) = 1. This contribution is the conditional
probability of class 𝑘 given 𝑥: 𝑝(𝑘 | 𝑥) (sometimes called the posterior). It can be computed
using Bayes’ rule:
𝑝(𝑥 | 𝑘)𝑝(𝑘)
𝑝(𝑘 | 𝑥) = (5.16)
𝑝(𝑥)
𝒩 (𝑥 | 𝜇𝑘 , Σ𝑘 )𝑝(𝑘)
= ∑︀𝐾 (5.17)
𝑘=1 𝒩 (𝑥 | 𝜇𝑘 , Σ𝑘 )𝑝(𝑘)
Since the class parameters, 𝑝(𝑘), 𝜇𝑘 and Σ𝑘 , depend on the responsibilities 𝑝(𝑘 | 𝑥) and the
responsibilities depend on class parameters, we need a two-step iterative algorithm: the
expectation-maximization (EM) algorithm. We discuss this algorithm next.
### The expectation-maximization (EM) algorithm for Gaussian mixtures
Given a Gaussian mixture model, the goal is to maximize the likelihood function with respect
to the parameters (comprised of the means and covariances of the components and the mixing
coefficients).
Initialize the means 𝜇𝑘 , covariances Σ𝑘 and mixing coefficients 𝑝(𝑘)
1. E step. For each sample 𝑖, evaluate the responsibilities for each class 𝑘 using the current
parameter values
𝒩 (𝑥𝑖 | 𝜇𝑘 , Σ𝑘 )𝑝(𝑘)
𝑝(𝑘 | 𝑥𝑖 ) = ∑︀𝐾
𝑘=1 𝒩 (𝑥𝑖 | 𝜇𝑘 , Σ𝑘 )𝑝(𝑘)
2. M step. For each class, re-estimate the parameters using the current responsibilities
𝑁
1 ∑︁
𝜇new
𝑘 = 𝑝(𝑘 | 𝑥𝑖 )𝑥𝑖 (5.18)
𝑁𝑘
𝑖=1
𝑁
1 ∑︁
Σnew
𝑘 = 𝑝(𝑘 | 𝑥𝑖 )(𝑥𝑖 − 𝜇new new 𝑇
𝑘 )(𝑥𝑖 − 𝜇𝑘 ) (5.19)
𝑁𝑘
𝑖=1
new 𝑁𝑘
𝑝 (𝑘) = (5.20)
𝑁
3. Evaluate the log-likelihood
𝑁
{︃ 𝐾 }︃
∑︁ ∑︁
ln 𝒩 (𝑥|𝜇𝑘 , Σ𝑘 )𝑝(𝑘) ,
𝑖=1 𝑘=1
and check for convergence of either the parameters or the log-likelihood. If the convergence
criterion is not satisfied return to step 1.
import numpy as np
from sklearn import datasets
import matplotlib.pyplot as plt
import seaborn as sns # nice color
import sklearn
from sklearn.mixture import GaussianMixture
import pystatsml.plot_utils
colors = sns.color_palette()
iris = datasets.load_iris()
X = iris.data[:, :2] # 'sepal length (cm)''sepal width (cm)'
y_iris = iris.target
plt.figure(figsize=(9, 3))
plt.subplot(131)
plt.scatter(X[:, 0], X[:, 1], c=[colors[lab] for lab in gmm2.predict(X)])#, color=colors)
for i in range(gmm2.covariances_.shape[0]):
pystatsml.plot_utils.plot_cov_ellipse(cov=gmm2.covariances_[i, :], pos=gmm2.means_[i,␣
˓→:],
facecolor='none', linewidth=2, edgecolor=colors[i])
plt.scatter(gmm2.means_[i, 0], gmm2.means_[i, 1], edgecolor=colors[i],
marker="o", s=100, facecolor="w", linewidth=2)
plt.title("K=2")
plt.subplot(132)
plt.scatter(X[:, 0], X[:, 1], c=[colors[lab] for lab in gmm3.predict(X)])
for i in range(gmm3.covariances_.shape[0]):
pystatsml.plot_utils.plot_cov_ellipse(cov=gmm3.covariances_[i, :], pos=gmm3.means_[i,␣
˓→:],
facecolor='none', linewidth=2, edgecolor=colors[i])
plt.scatter(gmm3.means_[i, 0], gmm3.means_[i, 1], edgecolor=colors[i],
marker="o", s=100, facecolor="w", linewidth=2)
(continues on next page)
plt.subplot(133)
plt.scatter(X[:, 0], X[:, 1], c=[colors[lab] for lab in gmm4.predict(X)]) # .astype(np.
˓→float))
for i in range(gmm4.covariances_.shape[0]):
pystatsml.plot_utils.plot_cov_ellipse(cov=gmm4.covariances_[i, :], pos=gmm4.means_[i,␣
˓→:],
In statistics, the Bayesian information criterion (BIC) is a criterion for model selection among
a finite set of models; the model with the lowest BIC is preferred. It is based, in part, on the
likelihood function and it is closely related to the Akaike information criterion (AIC).
X = iris.data
y_iris = iris.target
bic = list()
#print(X)
ks = np.arange(1, 10)
for k in ks:
gmm = GaussianMixture(n_components=k, covariance_type='full')
gmm.fit(X)
bic.append(gmm.bic(X))
k_chosen = ks[np.argmin(bic)]
plt.plot(ks, bic)
plt.xlabel("k")
(continues on next page)
Choose k= 2
Hierarchical clustering is an approach to clustering that build hierarchies of clusters in two main
approaches:
• Agglomerative: A bottom-up strategy, where each observation starts in their own cluster,
and pairs of clusters are merged upwards in the hierarchy.
• Divisive: A top-down strategy, where all observations start out in the same cluster, and
then the clusters are split recursively downwards in the hierarchy.
In order to decide which clusters to merge or to split, a measure of dissimilarity between clusters
is introduced. More specific, this comprise a distance measure and a linkage criterion. The
distance measure is just what it sounds like, and the linkage criterion is essentially a function of
the distances between points, for instance the minimum distance between points in two clusters,
the maximum distance between points in two clusters, the average distance between points in
two clusters, etc. One particular linkage criterion, the Ward criterion, will be discussed next.
Ward clustering
Ward clustering belongs to the family of agglomerative hierarchical clustering algorithms. This
means that they are based on a “bottoms up” approach: each sample starts in its own cluster,
and pairs of clusters are merged as one moves up the hierarchy.
In Ward clustering, the criterion for choosing the pair of clusters to merge at each step is the
minimum variance criterion. Ward’s minimum variance criterion minimizes the total within-
cluster variance by each merge. To implement this method, at each step: find the pair of
clusters that leads to minimum increase in total within-cluster variance after merging. This
increase is a weighted squared distance between cluster centers.
The main advantage of agglomerative hierarchical clustering over 𝐾-means clustering is that
you can benefit from known neighborhood information, for example, neighboring pixels in an
image.
iris = datasets.load_iris()
X = iris.data[:, :2] # 'sepal length (cm)''sepal width (cm)'
y_iris = iris.target
plt.figure(figsize=(9, 3))
plt.subplot(131)
plt.scatter(X[:, 0], X[:, 1], c=ward2.labels_)
plt.title("K=2")
plt.subplot(132)
plt.scatter(X[:, 0], X[:, 1], c=ward3.labels_)
plt.title("K=3")
plt.subplot(133)
plt.scatter(X[:, 0], X[:, 1], c=ward4.labels_) # .astype(np.float))
plt.title("K=4")
5.2.5 Exercises
Perform clustering of the iris dataset based on all variables using Gaussian mixture models. Use
PCA to visualize clusters.
Linear regression models the output, or target variable 𝑦 ∈ R as a linear combination of the 𝑃 -
dimensional input x ∈ R𝑃 . Let X be the 𝑁 × 𝑃 matrix with each row an input vector (with a 1
in the first position), and similarly let y be the 𝑁 -dimensional vector of outputs in the training
set, the linear model will predict the y given x using the parameter vector, or weight vector
w ∈ R𝑃 according to
y = Xw + 𝜀,
where 𝜀 ∈ R𝑁 are the residuals, or the errors of the prediction. The w is found by minimizing
an objective function, which is the loss function, 𝐿(w), i.e. the error measured on the data.
This error is the sum of squared errors (SSE) loss.
Minimizing the SSE is the Ordinary Least Square OLS regression as objective function. which is
a simple ordinary least squares (OLS) minimization whose analytic solution is:
𝐿(w, X, y) ∑︁
𝜕 =2 x𝑖 (x𝑖 · w − 𝑦𝑖 )
𝜕w
𝑖
Scikit learn offer many models for supervised learning, and they all follow the same application
programming interface (API), namely:
model = Estimator()
model.fit(X, y)
predictions = model.predict(X)
%matplotlib inline
#import warnings
#warnings.filterwarnings(action='once')
lr = lm.LinearRegression().fit(X, y)
y_pred = lr.predict(X)
print("R-squared =", metrics.r2_score(y, y_pred))
# Plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
XX = np.column_stack([xx1.ravel(), xx2.ravel()])
yy = lr.predict(XX)
ax.plot_surface(xx1, xx2, yy.reshape(xx1.shape), color='None')
ax.set_xlabel('TV')
ax.set_ylabel('Radio')
_ = ax.set_zlabel('Sales')
R-squared = 0.8971942610828956
Coefficients = [0.04575482 0.18799423]
5.3.3 Overfitting
In statistics and machine learning, overfitting occurs when a statistical model describes random
errors or noise instead of the underlying relationships. Overfitting generally occurs when a
model is excessively complex, such as having too many parameters relative to the number
of observations. A model that has been overfit will generally have poor predictive performance,
as it can exaggerate minor fluctuations in the data.
A learning algorithm is trained using some set of training samples. If the learning algorithm has
the capacity to overfit the training samples the performance on the training sample set will
improve while the performance on unseen test sample set will decline.
The overfitting phenomenon has three main explanations: - excessively complex models, - mul-
ticollinearity, and - high dimensionality.
Model complexity
Complex learners with too many parameters relative to the number of observations may overfit
the training dataset.
Multicollinearity
Predictors are highly correlated, meaning that one can be linearly predicted from the others.
In this situation the coefficient estimates of the multiple regression may change erratically in
response to small changes in the model or the data. Multicollinearity does not reduce the
predictive power or reliability of the model as a whole, at least not within the sample data
set; it only affects computations regarding individual predictors. That is, a multiple regression
model with correlated predictors can indicate how well the entire bundle of predictors predicts
the outcome variable, but it may not give valid results about any individual predictor, or about
which predictors are redundant with respect to others. In case of perfect multicollinearity the
predictor matrix is singular and therefore cannot be inverted. Under these circumstances, for a
general linear model y = Xw + 𝜀, the ordinary least-squares estimator, w𝑂𝐿𝑆 = (X𝑇 X)−1 X𝑇 y,
does not exist.
An example where correlated predictor may produce an unstable model follows:
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
X = np.column_stack([bv, tax])
beta_star = np.array([.1, 0]) # true solution
'''
Since tax and bv are correlated, there is an infinite number of linear combinations
leading to the same prediction.
'''
High dimensionality
High dimensions means a large number of input features. Linear predictor associate one pa-
rameter to each input feature, so a high-dimensional situation (𝑃 , number of features, is large)
with a relatively small number of samples 𝑁 (so-called large 𝑃 small 𝑁 situation) generally
lead to an overfit of the training data. Thus it is generally a bad idea to add many input features
into the learner. This phenomenon is called the curse of dimensionality.
One of the most important criteria to use when choosing a learning algorithm is based on the
relative size of 𝑃 and 𝑁 .
• Remenber that the “covariance” matrix X𝑇 X used in the linear model is a 𝑃 × 𝑃 matrix of
rank min(𝑁, 𝑃 ). Thus if 𝑃 > 𝑁 the equation system is overparameterized and admit an
infinity of solutions that might be specific to the learning dataset. See also ill-conditioned
or singular matrices.
• The sampling density of 𝑁 samples in an 𝑃 -dimensional space is proportional to 𝑁 1/𝑃 .
Thus a high-dimensional space becomes very sparse, leading to poor estimations of sam-
ples densities.
• Another consequence of the sparse sampling in high dimensions is that all sample points
are close to an edge of the sample. Consider 𝑁 data points uniformly distributed in a
𝑃 -dimensional unit ball centered at the origin. Suppose we consider a nearest-neighbor
estimate at the origin. The median distance from the origin to the closest data point is
given by the expression
)︂1/𝑃
1𝑁
(︂
𝑑(𝑃, 𝑁 ) = 1 − .
2
A more complicated expression exists for the mean distance to the closest point. For N = 500,
P = 10 , 𝑑(𝑃, 𝑁 ) ≈ 0.52, more than halfway to the boundary. Hence most data points are
closer to the boundary of the sample space than to any other data point. The reason that
this presents a problem is that prediction is much more difficult near the edges of the training
sample. One must extrapolate from neighboring sample points rather than interpolate between
them. (Source: T Hastie, R Tibshirani, J Friedman. The Elements of Statistical Learning: Data
Mining, Inference, and Prediction. Second Edition, 2009.)
• Structural risk minimization provides a theoretical background of this phenomenon. (See
VC dimension.)
• See also bias–variance trade-off.
def fit_on_increasing_size(model):
n_samples = 100
n_features_ = np.arange(10, 800, 20)
r2_train, r2_test, snr = [], [], []
for n_features in n_features_:
# Sample the dataset (* 2 nb of samples)
n_features_info = int(n_features/10)
np.random.seed(42) # Make reproducible
X = np.random.randn(n_samples * 2, n_features)
beta = np.zeros(n_features)
beta[:n_features_info] = 1
Xbeta = np.dot(X, beta)
eps = np.random.randn(n_samples * 2)
y = Xbeta + eps
# Split the dataset into train and test sample
Xtrain, Xtest = X[:n_samples, :], X[n_samples:, :]
ytrain, ytest = y[:n_samples], y[n_samples:]
# fit/predict
lr = model.fit(Xtrain, ytrain)
y_pred_train = lr.predict(Xtrain)
y_pred_test = lr.predict(Xtest)
snr.append(Xbeta.std() / eps.std())
r2_train.append(metrics.r2_score(ytrain, y_pred_train))
r2_test.append(metrics.r2_score(ytest, y_pred_test))
return n_features_, np.array(r2_train), np.array(r2_test), np.array(snr)
argmax = n_features[np.argmax(r2_test)]
# plot
fig, axis = plt.subplots(1, 2, figsize=(9, 3))
Exercises
Regarding linear models, overfitting generally leads to excessively complex solutions (coeffi-
cient vectors), accounting for noise or spurious correlations within predictors. Regularization
aims to alleviate this phenomenon by constraining (biasing or reducing) the capacity of the
learning algorithm in order to promote simple solutions. Regularization penalizes “large” solu-
tions forcing the coefficients to be small, i.e. to shrink them toward zeros.
The objective function 𝐽(w) to minimize with respect to w is composed of a loss function 𝐿(w)
for goodness-of-fit and a penalty term Ω(w) (regularization to avoid overfitting). This is a
trade-off where the respective contribution of the loss and the penalty terms is controlled by
the regularization parameter 𝜆.
Therefore the loss function 𝐿(w) (generally the SSE) is combined with a penalty function
Ω(w) leading to the general form:
The respective contribution of the loss and the penalty is controlled by the regularization
parameter 𝜆.
Ridge regression impose a ℓ2 penalty on the coefficients, i.e. it penalizes with the Euclidean
norm of the coefficients while minimizing SSE. The objective function becomes:
𝑁
∑︁
Ridge(w) = (𝑦𝑖 − x𝑇𝑖 w)2 + 𝜆‖w‖22 (5.25)
𝑖
= ‖y − xw‖22 + 𝜆‖w‖22 . (5.26)
The w that minimises 𝐹𝑅𝑖𝑑𝑔𝑒 (w) can be found by the following derivation:
∇w Ridge(w) = 0 (5.27)
∇w (y − Xw) (y − Xw) + 𝜆w𝑇 w = 0
𝑇
(︀ )︀
(5.28)
∇w (y𝑇 y − 2w𝑇 X𝑇 y + w𝑇 X𝑇 Xw + 𝜆w𝑇 w) = 0
(︀ )︀
(5.29)
𝑇 𝑇
−2X y + 2X Xw + 2𝜆w = 0 (5.30)
−X𝑇 y + (X𝑇 X + 𝜆I)w = 0 (5.31)
𝑇 𝑇
(X X + 𝜆I)w = x y (5.32)
w = (X𝑇 X + 𝜆I)−1 x𝑇 y (5.33)
• The solution adds a positive constant to the diagonal of X𝑇 X before inversion. This makes
the problem nonsingular, even if X𝑇 X is not of full rank, and was the main motivation
behind ridge regression.
• Increasing 𝜆 shrinks the w coefficients toward 0.
• This approach penalizes the objective function by the Euclidian (:math:`ell_2`) norm
of the coefficients such that solutions with large coefficients become unattractive.
The gradient of the loss:
𝐿(w, X, y) ∑︁
𝜕 = 2( x𝑖 (x𝑖 · w − 𝑦𝑖 ) + 𝜆w)
𝜕w
𝑖
# lambda is alpha!
mod = lm.Ridge(alpha=10)
argmax = n_features[np.argmax(r2_test)]
# plot
fig, axis = plt.subplots(1, 2, figsize=(9, 3))
Exercice
Lasso regression penalizes the coefficients by the ℓ1 norm. This constraint will reduce (bias)
the capacity of the learning algorithm. To add such a penalty forces the coefficients to be small,
i.e. it shrinks them toward zero. The objective function to minimize becomes:
This penalty forces some coefficients to be exactly zero, providing a feature selection property.
# lambda is alpha !
mod = lm.Lasso(alpha=.1)
argmax = n_features[np.argmax(r2_test)]
# plot
fig, axis = plt.subplots(1, 2, figsize=(9, 3))
Occam’s razor
Occam’s razor (also written as Ockham’s razor, and lex parsimoniae in Latin, which means
law of parsimony) is a problem solving principle attributed to William of Ockham (1287-1347),
who was an English Franciscan friar and scholastic philosopher and theologian. The principle
can be interpreted as stating that among competing hypotheses, the one with the fewest
assumptions should be selected.
Principle of parsimony
The penalty based on the ℓ1 norm promotes sparsity (scattered, or not dense): it forces many
coefficients to be exactly zero. This also makes the coefficient vector scattered.
The figure bellow illustrates the OLS loss under a constraint acting on the ℓ1 norm of the coef-
ficient vector. I.e., it illustrates the following optimization problem:
minimize ‖y − Xw‖22
w
subject to ‖w‖1 ≤ 1.
Optimization issues
Section to be completed
• No more closed-form solution.
• Convex but not differentiable.
• Requires specific optimization algorithms, such as the fast iterative shrinkage-thresholding
algorithm (FISTA): Amir Beck and Marc Teboulle, A Fast Iterative Shrinkage-Thresholding
Algorithm for Linear Inverse Problems SIAM J. Imaging Sci., 2009.
The ridge penalty shrinks the coefficients toward zero. The figure illustrates: the OLS solution
on the left. The ℓ1 and ℓ2 penalties in the middle pane. The penalized OLS in the right pane.
The right pane shows how the penalties shrink the coefficients toward zero. The black points
are the minimum found in each case, and the white points represents the true solution used to
generate the data.
The Elastic-net estimator combines the ℓ1 and ℓ2 penalties, and results in the problem to
Rational
• If there are groups of highly correlated variables, Lasso tends to arbitrarily select only
one from each group. These models are difficult to interpret because covariates that are
strongly associated with the outcome are not included in the predictive model. Conversely,
the elastic net encourages a grouping effect, where strongly correlated predictors tend to
be in or out of the model together.
• Studies on real world data and simulation studies show that the elastic net often outper-
forms the lasso, while enjoying a similar sparsity of representation.
argmax = n_features[np.argmax(r2_test)]
# plot
fig, axis = plt.subplots(1, 2, figsize=(9, 3))
This geometric method does not make any probabilistic assumptions, instead it relies on dis-
tances. It looks for the linear projection of the data points onto a vector, 𝑤, that maximizes
the between/within variance ratio, denoted 𝐹 (𝑤). Under a few assumptions, it will provide the
same results as linear discriminant analysis (LDA), explained below.
Suppose two classes of observations, 𝐶0 and 𝐶1 , have means 𝜇0 and 𝜇1 and the same total
within-class scatter (“covariance”) matrix,
∑︁ ∑︁
𝑆𝑊 = (𝑥𝑖 − 𝜇0 )(𝑥𝑖 − 𝜇0 )𝑇 + (𝑥𝑗 − 𝜇1 )(𝑥𝑗 − 𝜇1 )𝑇 (5.36)
𝑖∈𝐶0 𝑗∈𝐶1
𝑇
= 𝑋𝑐 𝑋𝑐 , (5.37)
where 𝑋0 and 𝑋1 are the (𝑁0 × 𝑃 ) and (𝑁1 × 𝑃 ) matrices of samples of classes 𝐶0 and 𝐶1 .
Let 𝑆𝐵 being the scatter “between-class” matrix, given by
𝑆𝐵 = (𝜇1 − 𝜇0 )(𝜇1 − 𝜇0 )𝑇 .
2
𝜎between
𝐹Fisher (𝑤) = 2 (5.38)
𝜎within
(𝑤𝑇 𝜇1 − 𝑤𝑇 𝜇0 )2
= (5.39)
𝑤𝑇 𝑋𝑐𝑇 𝑋𝑐 𝑤
(𝑤𝑇 (𝜇1 − 𝜇0 ))2
= (5.40)
𝑤𝑇 𝑋𝑐𝑇 𝑋𝑐 𝑤
𝑤𝑇 (𝜇1 − 𝜇0 )(𝜇1 − 𝜇0 )𝑇 𝑤
= (5.41)
𝑤𝑇 𝑋𝑐𝑇 𝑋𝑐 𝑤
𝑤 𝑇 𝑆𝐵 𝑤
= 𝑇 . (5.42)
𝑤 𝑆𝑊 𝑤
In the two-class case, the maximum separation occurs by a projection on the (𝜇1 − 𝜇0 ) using
the Mahalanobis metric 𝑆𝑊 −1 , so that
𝑤 ∝ 𝑆𝑊 −1 (𝜇1 − 𝜇0 ).
Demonstration
∇𝑤 𝐹Fisher (𝑤) = 0
(︂ 𝑇 )︂
𝑤 𝑆𝐵 𝑤
∇𝑤 =0
𝑤 𝑇 𝑆𝑊 𝑤
(𝑤𝑇 𝑆𝑊 𝑤)(2𝑆𝐵 𝑤) − (𝑤𝑇 𝑆𝐵 𝑤)(2𝑆𝑊 𝑤) = 0
(𝑤𝑇 𝑆𝑊 𝑤)(𝑆𝐵 𝑤) = (𝑤𝑇 𝑆𝐵 𝑤)(𝑆𝑊 𝑤)
𝑤 𝑇 𝑆𝐵 𝑤
𝑆𝐵 𝑤 = (𝑆𝑊 𝑤)
𝑤 𝑇 𝑆𝑊 𝑤
𝑆𝐵 𝑤 = 𝜆(𝑆𝑊 𝑤)
𝑆𝑊 −1 𝑆𝐵 𝑤 = 𝜆𝑤.
Since we do not care about the magnitude of 𝑤, only its direction, we replaced the scalar factor
(𝑤𝑇 𝑆𝐵 𝑤)/(𝑤𝑇 𝑆𝑊 𝑤) by 𝜆.
In the multiple-class case, the solutions 𝑤 are determined by the eigenvectors of 𝑆𝑊 −1 𝑆𝐵 that
correspond to the 𝐾 − 1 largest eigenvalues.
However, in the two-class case (in which 𝑆𝐵 = (𝜇1 − 𝜇0 )(𝜇1 − 𝜇0 )𝑇 ) it is easy to show that
𝑤 = 𝑆𝑊 −1 (𝜇1 − 𝜇0 ) is the unique eigenvector of 𝑆𝑊 −1 𝑆𝐵 :
𝑆𝑊 −1 (𝜇1 − 𝜇0 )(𝜇1 − 𝜇0 )𝑇 𝑤 = 𝜆𝑤
𝑆𝑊 −1 (𝜇1 − 𝜇0 )(𝜇1 − 𝜇0 )𝑇 𝑆𝑊 −1 (𝜇1 − 𝜇0 ) = 𝜆𝑆𝑊 −1 (𝜇1 − 𝜇0 ),
𝑤 ∝ 𝑆𝑊 −1 (𝜇1 − 𝜇0 ).
warnings.filterwarnings('ignore')
%matplotlib inline
Exercise
import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
# Dataset
n_samples, n_features = 100, 2
mean0, mean1 = np.array([0, 0]), np.array([0, 2])
Cov = np.array([[1, .8],[.8, 1]])
np.random.seed(42)
X0 = np.random.multivariate_normal(mean0, Cov, n_samples)
X1 = np.random.multivariate_normal(mean1, Cov, n_samples)
X = np.vstack([X0, X1])
y = np.array([0] * X0.shape[0] + [1] * X1.shape[0])
errors = y_pred_lda != y
print("Nb errors=%i, error rate=%.2f" % (errors.sum(), errors.sum() / len(y_pred_lda)))
Logistic regression is called a generalized linear models. ie.: it is a linear model with a link
function that maps the output of linear multiple regression to the posterior probability of class
1 𝑝(1|𝑥) using the logistic sigmoid function:
1
𝑝(1|𝑤, 𝑥𝑖 ) =
1 + exp(−𝑤 · 𝑥𝑖 )
x = np.linspace(-6, 6, 100)
plt.subplot(121)
(continues on next page)
x = np.linspace(-3, 3, 100)
plt.subplot(122)
plt.plot(x, logistic_loss(x), label='Logistic loss')
plt.plot(x, np.maximum(0, 1 - x), label='Hinge loss')
plt.legend()
plt.title('Losses')
plt.grid(True)
The Loss function for sample 𝑖 is the negative log of the probability:
{︃
− log(𝑝(1|𝑤, 𝑥𝑖 )) if 𝑦𝑖 = 1
𝐿(𝑤, 𝑥𝑖 , 𝑦𝑖 ) =
− log(1 − 𝑝(1|𝑤, 𝑥𝑖 ) if 𝑦𝑖 = 0
For the whole dataset 𝑋, 𝑦 = {𝑥𝑖 , 𝑦𝑖 } the loss function to minimize 𝐿(𝑤, 𝑋, 𝑦) is the negative
negative log likelihood (nll) that can be simplied using a 0/1 coding of the label in the case of
binary classification:
𝐿(𝑤, 𝑋, 𝑦) ∑︁
𝜕 = 𝑥𝑖 (𝑦𝑖 − 𝑝(1|𝑤, 𝑥𝑖 ))
𝜕𝑤
𝑖
Logistic regression is a discriminative model since it focuses only on the posterior probability
of each class 𝑝(𝐶𝑘 |𝑥). It only requires to estimate the 𝑃 weight of the 𝑤 vector. Thus it should
be favoured over LDA with many input features. In small dimension and balanced situations it
would provide similar predictions than LDA.
However imbalanced group sizes cannot be explicitly controlled. It can be managed using a
reweighting of the input samples.
logreg.fit(X, y)
y_pred_logreg = logreg.predict(X)
errors = y_pred_logreg != y
print("Nb errors=%i, error rate=%.2f" % (errors.sum(), errors.sum() / len(y_pred_logreg)))
print(logreg.coef_)
Exercise
Explore the Logistic Regression parameters and proposes a solution in cases of highly im-
balanced training dataset 𝑁1 ≫ 𝑁0 when we know that in reality both classes have the same
probability 𝑝(𝐶1 ) = 𝑝(𝐶0 ).
5.4.4 Overfitting
When the matrix 𝑆𝑊 is not full rank or 𝑃 ≫ 𝑁 , the The Fisher most discriminant projection
estimate of the is not unique. This can be solved using a biased version of 𝑆𝑊 :
𝑆𝑊 𝑅𝑖𝑑𝑔𝑒 = 𝑆𝑊 + 𝜆𝐼
where 𝐼 is the 𝑃 × 𝑃 identity matrix. This leads to the regularized (ridge) estimator of the
Fisher’s linear discriminant analysis:
Increasing 𝜆 will:
• Shrinks the coefficients toward zero.
• The covariance will converge toward the diagonal matrix, reducing the contribution of
the pairwise covariances.
The objective function to be minimized is now the combination of the logistic loss (negative
log likelyhood) − log ℒ(𝑤) with a penalty of the L2 norm of the weights vector. In the two-class
case, using the 0/1 coding we obtain:
# Dataset
# Build a classification task using 3 informative features
from sklearn import datasets
X, y = datasets.make_classification(n_samples=100,
n_features=20,
n_informative=3,
n_redundant=0,
n_repeated=0,
(continues on next page)
import numpy as np
import matplotlib.pyplot as plt
lr.fit(X, y)
y_pred_lr = lr.predict(X)
errors = y_pred_lr != y
print("Nb errors=%i, error rate=%.2f" % (errors.sum(), errors.sum() / len(y)))
print(lr.coef_)
(1, 20)
Diff 0.0
Nb errors=26, error rate=0.26
[[-0.12899737 0.7579822 -0.01228473 -0.11412421 0.25491221 0.4329847
0.14564739 0.16763962 0.85071394 0.02116803 -0.1611039 -0.0146019
-0.03399884 0.43127728 -0.05831644 -0.0812323 0.15877844 0.29387389
0.54659524 0.03376169]]
The objective function to be minimized is now the combination of the logistic loss − log ℒ(𝑤)
with a penalty of the L1 norm of the weights vector. In the two-class case, using the 0/1 coding
we obtain:
lrl1.fit(X, y)
y_pred_lrl1 = lrl1.predict(X)
errors = y_pred_lrl1 != y
(continues on next page)
Support Vector Machine seek for separating hyperplane with maximum margin to enforce ro-
bustness against noise. Like logistic regression it is a discriminative method that only focuses
of predictions.
Here we present the non separable case of Maximum Margin Classifiers with ±1 coding (ie.:
𝑦𝑖 {−1, +1}). In the next figure the legend aply to samples of “dot” class.
Here we introduced the slack variables: 𝜉𝑖 , with 𝜉𝑖 = 0 for points that are on or inside the
correct margin boundary and 𝜉𝑖 = |𝑦𝑖 − (𝑤 𝑐𝑑𝑜𝑡 · 𝑥𝑖 )| for other points. Thus:
1. If 𝑦𝑖 (𝑤 · 𝑥𝑖 ) ≥ 1 then the point lies outside the margin but on the correct side of the
decision boundary. In this case 𝜉𝑖 = 0. The constraint is thus not active for this point. It
does not contribute to the prediction.
2. If 1 > 𝑦𝑖 (𝑤 · 𝑥𝑖 ) ≥ 0 then the point lies inside the margin and on the correct side of the
decision boundary. In this case 0 < 𝜉𝑖 ≤ 1. The constraint is active for this point. It does
contribute to the prediction as a support vector.
3. If 0 < 𝑦𝑖 (𝑤 · 𝑥𝑖 )) then the point is on the wrong side of the decision boundary (missclassi-
fication). In this case 0 < 𝜉𝑖 > 1. The constraint is active for this point. It does contribute
to the prediction as a support vector.
This loss is called the hinge loss, defined as:
max(0, 1 − 𝑦𝑖 (𝑤 · 𝑥𝑖 ))
So linear SVM is closed to Ridge logistic regression, using the hinge loss instead of the logistic
loss. Both will provide very similar predictions.
svmlin = svm.LinearSVC()
# Remark: by default LinearSVC uses squared_hinge as loss
svmlin.fit(X, y)
y_pred_svmlin = svmlin.predict(X)
errors = y_pred_svmlin != y
print("Nb errors=%i, error rate=%.2f" % (errors.sum(), errors.sum() / len(y_pred_svmlin)))
print(svmlin.coef_)
Linear SVM for classification (also called SVM-C or SVC) with l1-regularization
∑︀𝑁
min 𝐹Lasso linear SVM (𝑤) = ||𝑤||1 + 𝐶 𝑖 𝜉𝑖
with ∀𝑖 𝑦𝑖 (𝑤 · 𝑥𝑖 ) ≥ 1 − 𝜉𝑖
svmlinl1.fit(X, y)
y_pred_svmlinl1 = svmlinl1.predict(X)
errors = y_pred_svmlinl1 != y
print("Nb errors=%i, error rate=%.2f" % (errors.sum(), errors.sum() / len(y_pred_
˓→svmlinl1)))
print(svmlinl1.coef_)
## Exercise
Compare predictions of Logistic regression (LR) and their SVM counterparts, ie.: L2 LR vs L2
SVM and L1 LR vs L1 SVM
• Compute the correlation between pairs of weights vectors.
• Compare the predictions of two classifiers using their decision function:
– Give the equation of the decision function for a linear classifier, assuming that their
is no intercept.
– Compute the correlation decision function.
– Plot the pairwise decision function of the classifiers.
• Conclude on the differences between Linear SVM and logistic regression.
The objective function to be minimized is now the combination of the logistic loss log 𝐿(𝑤) or
the hinge loss with combination of L1 and L2 penalties. In the two-class case, using the 0/1
coding we obtain:
X, y = datasets.make_classification(n_samples=100,
n_features=20,
n_informative=3,
n_redundant=0,
n_repeated=0,
n_classes=2,
random_state=0,
shuffle=False)
SGDClassifier(class_weight='balanced', penalty='elasticnet')
Exercise
source: https://en.wikipedia.org/wiki/Sensitivity_and_specificity
Imagine a study evaluating a new test that screens people for a disease. Each person taking the
test either has or does not have the disease. The test outcome can be positive (classifying the
person as having the disease) or negative (classifying the person as not having the disease). The
test results for each subject may or may not match the subject’s actual status. In that setting:
• True positive (TP): Sick people correctly identified as sick
• False positive (FP): Healthy people incorrectly identified as sick
• True negative (TN): Healthy people correctly identified as healthy
• False negative (FN): Sick people incorrectly identified as healthy
• Accuracy (ACC):
ACC = (TP + TN) / (TP + FP + FN + TN)
• Sensitivity (SEN) or recall of the positive class or true positive rate (TPR) or hit rate:
SEN = TP / P = TP / (TP+FN)
• Specificity (SPC) or recall of the negative class or true negative rate:
SPC = TN / N = TN / (TN+FP)
• Precision or positive predictive value (PPV):
PPV = TP / (TP + FP)
• Balanced accuracy (bACC):is a useful performance measure is the balanced accuracy
which avoids inflated performance estimates on imbalanced datasets (Brodersen, et
al. (2010). “The balanced accuracy and its posterior distribution”). It is defined as the
arithmetic mean of sensitivity and specificity, or the average accuracy obtained on either
class:
bACC = 1/2 * (SEN + SPC)
• F1 Score (or F-score) which is a weighted average of precision and recall are usefull to
deal with imballaced datasets
The four outcomes can be formulated in a 2×2 contingency table or confusion matrix https:
//en.wikipedia.org/wiki/Sensitivity_and_specificity
For more precision see: http://scikit-learn.org/stable/modules/model_evaluation.html
metrics.accuracy_score(y_true, y_pred)
# Balanced accuracy
b_acc = recalls.mean()
P-value associated to classification rate. Compared the number of correct classifications (=ac-
curacy ×𝑁 ) to the null hypothesis of Binomial distribution of parameters 𝑝 (typically 50% of
chance level) and 𝑁 (Number of observations).
Is 65% of accuracy a significant prediction rate among 70 observations?
Since this is an exact, two-sided test of the null hypothesis, the p-value can be divided by 2
since we test that the accuracy is superior to the chance level.
import scipy.stats
acc, N = 0.65, 70
pval = scipy.stats.binom_test(x=int(acc * N), n=N, p=0.5) / 2
print(pval)
0.01123144774625465
Some classifier may have found a good discriminative projection 𝑤. However if the threshold
to decide the final predicted class is poorly adjusted, the performances will highlight an high
specificity and a low sensitivity or the contrary.
In this case it is recommended to use the AUC of a ROC analysis which basically provide a mea-
sure of overlap of the two classes when points are projected on the discriminative axis. For more
detail on ROC and AUC see:https://en.wikipedia.org/wiki/Receiver_operating_characteristic.
print("Predictions:", y_pred)
metrics.accuracy_score(y_true, y_pred)
Predictions: [0 0 0 0 0 0 0 0]
Recalls: [1. 0.]
AUC: 1.0
Learning with discriminative (logistic regression, SVM) methods is generally based on minimiz-
ing the misclassification of training samples, which may be unsuitable for imbalanced datasets
where the recognition might be biased in favor of the most numerous class. This problem
can be addressed with a generative approach, which typically requires more parameters to be
determined leading to reduced performances in high dimension.
Dealing with imbalanced class may be addressed by three main ways (see Japkowicz and
Stephen (2002) for a review), resampling, reweighting and one class learning.
In sampling strategies, either the minority class is oversampled or majority class is undersam-
pled or some combination of the two is deployed. Undersampling (Zhang and Mani, 2003) the
majority class would lead to a poor usage of the left-out samples. Sometime one cannot afford
such strategy since we are also facing a small sample size problem even for the majority class.
Informed oversampling, which goes beyond a trivial duplication of minority class samples, re-
quires the estimation of class conditional distributions in order to generate synthetic samples.
Here generative models are required. An alternative, proposed in (Chawla et al., 2002) generate
samples along the line segments joining any/all of the k minority class nearest neighbors. Such
procedure blindly generalizes the minority area without regard to the majority class, which may
be particularly problematic with high-dimensional and potentially skewed class distribution.
import numpy as np
from sklearn import linear_model
from sklearn import datasets
from sklearn import metrics
import matplotlib.pyplot as plt
# dataset
X, y = datasets.make_classification(n_samples=500,
n_features=5,
n_informative=2,
n_redundant=0,
n_repeated=0,
n_classes=2,
random_state=1,
shuffle=False)
5.4.13 Exercise
Write a class FisherLinearDiscriminant that implements the Fisher’s linear discriminant anal-
ysis. This class must be compliant with the scikit-learn API by providing two methods: - fit(X,
y) which fits the model and returns the object itself; - predict(X) which returns a vector of the
predicted values. Apply the object on the dataset presented for the LDA.
SVM are based kernel methods require only a user-specified kernel function 𝐾(𝑥𝑖 , 𝑥𝑗 ), i.e., a
similarity function over pairs of data points (𝑥𝑖 , 𝑥𝑗 ) into kernel (dual) space on which learning
algorithms operate linearly, i.e. every operation on points is a linear combination of 𝐾(𝑥𝑖 , 𝑥𝑗 ).
Outline of the SVM algorithm:
1. Map points 𝑥 into kernel space using a kernel function: 𝑥 → 𝐾(𝑥, .).
2. Learning algorithms operates linearly by dot product into high-kernel space 𝐾(., 𝑥𝑖 ) ·
𝐾(., 𝑥𝑗 ).
• Using the kernel trick (Mercer’s Theorem) replaces dot product in high dimensional
space by a simpler operation such that 𝐾(., 𝑥𝑖 ) · 𝐾(., 𝑥𝑗 ) = 𝐾(𝑥𝑖 , 𝑥𝑗 ). Thus we only
need to compute a similarity measure for each pairs of point and store in a 𝑁 × 𝑁
Gram matrix.
• Finally, The learning process consist of estimating the 𝛼𝑖 of the decision function that
maximises the hinge loss (of 𝑓 (𝑥)) plus some penalty when applied on all training
points.
(︃ 𝑁 )︃
∑︁
𝑓 (𝑥) = sign 𝛼𝑖 𝑦𝑖 𝐾(𝑥𝑖 , 𝑥) .
𝑖
One of the most commonly used kernel is the Radial Basis Function (RBF) Kernel. For a pair of
points 𝑥𝑖 , 𝑥𝑗 the RBF kernel is defined as:
‖𝑥𝑖 − 𝑥𝑗 ‖2
(︂ )︂
𝐾(𝑥𝑖 , 𝑥𝑗 ) = exp − (5.50)
2𝜎 2
= exp −𝛾 ‖𝑥𝑖 − 𝑥𝑗 ‖2
(︀ )︀
(5.51)
Where 𝜎 (or 𝛾) defines the kernel width parameter. Basically, we consider a Gaussian function
centered on each training sample 𝑥𝑖 . it has a ready interpretation as a similarity measure as it
decreases with squared Euclidean distance between the two feature vectors.
Non linear SVM also exists for regression problems.
%matplotlib inline
import warnings
warnings.filterwarnings(action='once')
import numpy as np
from sklearn.svm import SVC
from sklearn import datasets
import matplotlib.pyplot as plt
# dataset
X, y = datasets.make_classification(n_samples=10, n_features=2,n_redundant=0,
n_classes=2,
(continues on next page)
clf.decision_function(X)
# Usefull internals:
# Array of support vectors
clf.support_vectors_
/home/ed203246/anaconda3/lib/python3.7/site-packages/ipykernel/ipkernel.py:287:␣
˓→DeprecationWarning: should_run_async will not call transform_cell automatically␣
˓→in the future. Please pass the result to transformed_cell argument and any␣
and should_run_async(code)
/home/ed203246/anaconda3/lib/python3.7/importlib/_bootstrap.py:219:␣
˓→RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility.␣
#Errors: 0
True
A tree can be “learned” by splitting the training dataset into subsets based on an features value
test.
Each internal node represents a “test” on an feature resulting on the split of the current sample.
At each step the algorithm selects the feature and a cutoff value that maximises a given metric.
Different metrics exist for regression tree (target is continuous) or classification tree (the target
is qualitative).
This process is repeated on each derived subset in a recursive manner called recursive partition-
ing. The recursion is completed when the subset at a node has all the same value of the target
variable, or when splitting no longer adds value to the predictions. This general principle is
implemented by many recursive partitioning tree algorithms.
Decision trees are simple to understand and interpret however they tend to overfit the data.
However decision trees tend to overfit the training set. Leo Breiman propose random forest to
deal with this issue.
A single decision tree is usually overfits the data it is learning from because it learn from only
one pathway of decisions. Predictions from a single decision tree usually don’t make accurate
predictions on new data.
A random forest is a meta estimator that fits a number of decision tree learners on various
sub-samples of the dataset and use averaging to improve the predictive accuracy and control
over-fitting.
Random forest models reduce the risk of overfitting by introducing randomness by:
• building multiple trees (n_estimators)
• drawing observations with replacement (i.e., a bootstrapped sample)
• splitting nodes on the best split among a random subset of the features selected at every
node
#Errors: 0
/home/ed203246/anaconda3/lib/python3.7/site-packages/ipykernel/ipkernel.py:287:␣
˓→DeprecationWarning: should_run_async will not call transform_cell automatically␣
˓→in the future. Please pass the result to transformed_cell argument and any␣
and should_run_async(code)
/home/ed203246/anaconda3/lib/python3.7/importlib/_bootstrap.py:219:␣
˓→RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility.␣
Extra Trees is like Random Forest, in that it builds multiple trees and splits nodes using random
subsets of features, but with two key differences: it does not bootstrap observations (meaning
it samples without replacement), and nodes are split on random splits, not best splits. So, in
summary, ExtraTrees: builds multiple trees with bootstrap = False by default, which means it
samples without replacement nodes are split based on random splits among a random subset of
the features selected at every node In Extra Trees, randomness doesn’t come from bootstrapping
of data, but rather comes from the random splits of all observations. ExtraTrees is named for
(Extremely Randomized Trees).
warnings.filterwarnings('ignore')
%matplotlib inline
Machine learning algorithms overfit taining data. Predictive performances MUST be evaluated
on independant hold-out dataset.
1. Training dataset: Dataset used to fit the model (set the model parameters like weights).
The training error can be easily calculated by applying the statistical learning method to
the observations used in its training. But because of overfitting, the training error rate
can dramatically underestimate the error that would be obtained on new samples.
2. Validation dataset: Dataset used to provide an unbiased evaluation of a model fit on the
training dataset while tuning model hyperparameters. The validation error is the aver-
age error that results from a learning method to predict the response on a new (validation)
samples that is, on samples that were not used in training the method.
3. Test Dataset: Dataset used to provide an unbiased evaluation of a final model fit on the
training dataset. It is only used once a model is completely trained(using the train and
validation sets).
What is the Difference Between Test and Validation Datasets? by Jason Brownlee
Thus the original dataset is generally split in a training, validation and a test data sets. Large
training+validation set (80%) small test set (20%) might provide a poor estimation of the pre-
dictive performances (same argument stands for train vs validation samples). On the contrary,
large test set and small training set might produce a poorly estimated learner. This is why, on
situation where we cannot afford such split, it recommended to use cross-validation scheme to
estimate the predictive power of a learning algorithm.
Cross-Validation scheme randomly divides the set of observations into 𝐾 groups, or folds, of
approximately equal size. The first fold is treated as a validation set, and the method 𝑓 () is
fitted on the remaining union of 𝐾 − 1 folds: (𝑓 (𝑋 −𝐾 , 𝑦 −𝐾 )).
The measure of performance (the score function 𝒮), either a error measure or an correct predic-
tion measure is an average of a loss error or correct prediction measure, noted ℒ, between a true
target value and the predicted target value. The score function is evaluated of the on the obser-
vations in the held-out fold. For each sample 𝑖 we consider the model estimated 𝑓 (𝑋 −𝑘(𝑖) , 𝑦 −𝑘(𝑖)
on the data set without the group 𝑘 that contains 𝑖 noted −𝑘(𝑖). This procedure is repeated 𝐾
times; each time, a different group of observations is treated as a test set. Then we compare the
predicted value (𝑓−𝑘(𝑖) (𝑥𝑖 ) = 𝑦ˆ𝑖 ) with true value 𝑦𝑖 using a Error or Loss function ℒ(𝑦, 𝑦ˆ).
For 10-fold we can either average over 10 values (Macro measure) or concatenate the 10 ex-
periments and compute the micro measures.
Two strategies micro vs macro estimates:
Micro measure: average(individual scores): compute a score 𝒮 for each sample and average
over all samples. It is simillar to average score(concatenation): an averaged score computed
over all concatenated samples.
𝑁
1 ∑︁ (︁ )︁
𝒮(𝑓 ) = ℒ 𝑦𝑖 , 𝑓 (𝑥−𝑘(𝑖) , 𝑦 −𝑘(𝑖) ) .
𝑁
𝑖
Macro measure mean(CV scores) (the most commonly used method): compute a score 𝒮 on
each each fold 𝑘 and average accross folds:
𝐾
1 ∑︁
𝒮(𝑓 ) = 𝒮𝑘 (𝑓 ).
𝐾
𝑘
𝐾
1 ∑︁ 1 ∑︁ (︁ )︁
𝒮(𝑓 ) = ℒ 𝑦𝑖 , 𝑓 (𝑥−𝑘(𝑖) , 𝑦 −𝑘(𝑖) ) .
𝐾 𝑁𝑘
𝑘 𝑖∈𝑘
These two measures (an average of average vs. a global average) are generaly similar. They
may differ slightly is folds are of different sizes.
This validation scheme is known as the K-Fold CV. Typical choices of 𝐾 are 5 or 10, [Kohavi
1995]. The extreme case where 𝐾 = 𝑁 is known as leave-one-out cross-validation, LOO-CV.
CV for regression
Usually the error function ℒ() is the r-squared score. However other function could be used.
%matplotlib inline
import warnings
warnings.filterwarnings(action='once')
import numpy as np
from sklearn import datasets
import sklearn.linear_model as lm
import sklearn.metrics as metrics
from sklearn.model_selection import KFold
X, y = datasets.make_regression(n_samples=100, n_features=100,
n_informative=10, random_state=42)
estimator = lm.Ridge(alpha=10)
cv = KFold(n_splits=5, random_state=42)
r2_train, r2_test = list(), list()
Train r2:0.99
Test r2:0.73
# provide a cv
cv = KFold(n_splits=5, random_state=42)
scores = cross_val_score(estimator=estimator, X=X, y=y, cv=cv)
print("Test r2:%.2f" % scores.mean())
Test r2:0.73
Test r2:0.73
CV for classification
With classification problems it is essential to sample folds where each set contains approxi-
mately the same percentage of samples of each target class as the complete set. This is called
stratification. In this case, we will use StratifiedKFold with is a variation of k-fold which
returns stratified folds.
Usually the error function 𝐿() are, at least, the sensitivity and the specificity. However other
function could be used.
import numpy as np
from sklearn import datasets
import sklearn.linear_model as lm
import sklearn.metrics as metrics
from sklearn.model_selection import StratifiedKFold
X, y = datasets.make_classification(n_samples=100, n_features=100,
n_informative=10, random_state=42)
cv = StratifiedKFold(n_splits=5)
# Or vector of test predictions (for both macro and micro measures, not for training␣
˓→samples)
y_test_pred = np.zeros(len(y))
== Macro measures ==
Train SPC:1.00; SEN:1.00
Test SPC:0.78; SEN:0.82
Test ACC:0.80, ballanced ACC:0.80 Folds: [0.9, 0.7, 0.95, 0.7, 0.75]
Test ACC:0.80 Folds: [0.9, 0.7, 0.95, 0.7, 0.75]
== Micro measures ==
Test SPC:0.78; SEN:0.82
Test ACC:0.80
Test ACC:0.80
Note that with Scikit-learn user-friendly function we average the scores’ average obtained on
individual folds which may provide slightly different results that the overall average presented
earlier.
Dataset
import numpy as np
from sklearn import datasets
import sklearn.linear_model as lm
import sklearn.metrics as metrics
from sklearn.model_selection import StratifiedKFold
X, y = datasets.make_classification(n_samples=20, n_features=5, n_informative=2, random_
˓→state=42)
cv = StratifiedKFold(n_splits=5)
print(coefs_cv.mean(axis=0))
print("Std Err of the coef")
print(coefs_cv.std(axis=0) / np.sqrt(coefs_cv.shape[0]))
parallel = Parallel(n_jobs=5)
cv_ret = parallel(
delayed(_split_fit_predict)(
clone(estimator), X, y, train, test)
for train, test in cv.split(X, y))
• the test set (test) is used for assessment of the generalization error of the final chosen
model.
Model selection of the best hyper parameters over a grid of possible values
For each possible values of hyper parameters 𝛼𝑘 :
1. Fit the learner on training set: 𝑓 (𝑋𝑡𝑟𝑎𝑖𝑛 , 𝑦𝑡𝑟𝑎𝑖𝑛 , 𝛼𝑘 )
2. Evaluate the model on the validation set and keep the parameter(s) that minimises the
error measure
𝛼* = arg min 𝐿(𝑓 (𝑋𝑡𝑟𝑎𝑖𝑛 ), 𝑦𝑣𝑎𝑙 , 𝛼𝑘 )
3. Refit the learner on all training + validation data using the best hyper parameters: 𝑓 * ≡
𝑓 (𝑋𝑡𝑟𝑎𝑖𝑛∪𝑣𝑎𝑙 , 𝑦𝑡𝑟𝑎𝑖𝑛∪𝑣𝑎𝑙 , 𝛼* )
4. ** Model assessment ** of 𝑓 * on the test set: 𝐿(𝑓 * (𝑋𝑡𝑒𝑠𝑡 ), 𝑦𝑡𝑒𝑠𝑡 )
Most of time, we cannot afford such three-way split. Thus, again we will use CV, but in this case
we need two nested CVs.
One outer CV loop, for model assessment. This CV performs 𝐾 splits of the dataset into
training plus validation (𝑋−𝐾 , 𝑦−𝐾 ) set and a test set 𝑋𝐾 , 𝑦𝐾
One inner CV loop, for model selection. For each run of the outer loop, the inner loop loop
performs 𝐿 splits of dataset (𝑋−𝐾 , 𝑦−𝐾 ) into training set: (𝑋−𝐾,−𝐿 , 𝑦−𝐾,−𝐿 ) and a validation
set: (𝑋−𝐾,𝐿 , 𝑦−𝐾,𝐿 ).
Note that the inner CV loop combined with the learner form a new learner with an automatic
model (parameter) selection procedure. This new learner can be easily constructed using Scikit-
learn. The learned is wrapped inside a GridSearchCV class.
Then the new learned can be plugged into the classical outer CV loop.
import numpy as np
from sklearn import datasets
import sklearn.linear_model as lm
from sklearn.model_selection import GridSearchCV
import sklearn.metrics as metrics
from sklearn.model_selection import KFold
# Dataset
noise_sd = 10
X, y, coef = datasets.make_regression(n_samples=50, n_features=100, noise=noise_sd,
n_informative=2, random_state=42, coef=True)
# Use this to tune the noise parameter such that snr < 5
print("SNR:", np.std(np.dot(X, coef)) / noise_sd)
(continues on next page)
# Warp
model = GridSearchCV(lm.ElasticNet(max_iter=10000), param_grid, cv=5)
SNR: 2.6358469446381614
Sklearn will automatically select a grid of parameters, most of time use the defaults values.
n_jobs is the number of CPUs to use during the cross validation. If -1, use all the CPUs.
1) Biased usage: fit on all data, ommit outer CV loop
model.fit(X, y)
print("Train r2:%.2f" % metrics.r2_score(y, model.predict(X)))
print(model.best_params_)
Train r2:0.96
{'alpha': 1.0, 'l1_ratio': 0.9}
cv = KFold(n_splits=5, random_state=42)
r2_train, r2_test = list(), list()
alphas = list()
r2_test.append(metrics.r2_score(y_test, model.predict(X_test)))
r2_train.append(metrics.r2_score(y_train, model.predict(X_train)))
alphas.append(model.best_params_)
Train r2:1.00
Test r2:0.55
Selected alphas: [{'alpha': 0.001, 'l1_ratio': 0.9}, {'alpha': 0.001, 'l1_ratio': 0.9}, {
˓→'alpha': 0.001, 'l1_ratio': 0.9}, {'alpha': 0.01, 'l1_ratio': 0.9}, {'alpha': 0.001,
˓→'l1_ratio': 0.9}]
Test r2:0.55
# Dataset
X, y, coef = datasets.make_regression(n_samples=50, n_features=100, noise=10,
n_informative=2, random_state=42, coef=True)
X, y = datasets.make_classification(n_samples=100, n_features=100,
n_informative=10, random_state=42)
A permutation test is a type of non-parametric randomization test in which the null distribution
of a test statistic is estimated by randomly permuting the observations.
Permutation tests are highly attractive because they make no assumptions other than that the
observations are independent and identically distributed under the null hypothesis.
1. Compute a observed statistic 𝑡𝑜𝑏𝑠 on the data.
2. Use randomization to compute the distribution of 𝑡 under the null hypothesis: Perform 𝑁
random permutation of the data. For each sample of permuted data, 𝑖 the data compute
the statistic 𝑡𝑖 . This procedure provides the distribution of 𝑡 under the null hypothesis 𝐻0 :
𝑃 (𝑡|𝐻0 )
3. Compute the p-value = 𝑃 (𝑡 > 𝑡𝑜𝑏𝑠 |𝐻0 ) |{𝑡𝑖 > 𝑡𝑜𝑏𝑠 }|, where 𝑡𝑖 ’s include 𝑡𝑜𝑏𝑠 .
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
#%matplotlib qt
np.random.seed(42)
x = np.random.normal(loc=10, scale=1, size=100)
y = x + np.random.normal(loc=-3, scale=3, size=100) # snr = 1/2
# Plot
# Re-weight to obtain distribution
(continues on next page)
Exercise
Given the logistic regression presented above and its validation given a 5 folds CV.
1. Compute the p-value associated with the prediction accuracy using a permutation test.
2. Compute the p-value associated with the prediction accuracy using a parametric test.
5.6.6 Bootstrapping
import numpy as np
from sklearn import datasets
import sklearn.linear_model as lm
import sklearn.metrics as metrics
import pandas as pd
# Regression dataset
n_features = 5
n_features_info = 2
n_samples = 100
X = np.random.randn(n_samples, n_features)
beta = np.zeros(n_features)
beta[:n_features_info] = 1
Xbeta = np.dot(X, beta)
eps = np.random.randn(n_samples)
y = Xbeta + eps
# Bootstrap loop
nboot = 100 # !! Should be at least 1000
scores_names = ["r2"]
scores_boot = np.zeros((nboot, len(scores_names)))
coefs_boot = np.zeros((nboot, X.shape[1]))
orig_all = np.arange(X.shape[0])
for boot_i in range(nboot):
boot_tr = np.random.choice(orig_all, size=len(orig_all), replace=True)
boot_te = np.setdiff1d(orig_all, boot_tr, assume_unique=False)
Xtr, ytr = X[boot_tr, :], y[boot_tr]
Xte, yte = X[boot_te, :], y[boot_te]
model.fit(Xtr, ytr)
y_pred = model.predict(Xte).ravel()
scores_boot[boot_i, :] = metrics.r2_score(yte, y_pred)
coefs_boot[boot_i, :] = model.coef_
coefs_boot = pd.DataFrame(coefs_boot)
coefs_stat = coefs_boot.describe(percentiles=[.975, .5, .025])
print("Coefficients distribution")
print(coefs_stat)
These methods are Ensemble learning techniques. These models are machine learning
paradigms where multiple models (often called “weak learners”) are trained to solve the same
problem and combined to get better results. The main hypothesis is that when weak models
are correctly combined we can obtain more accurate and/or robust models.
a high bias (low degree of freedom models, for example) or because they have too much
variance to be robust (high degree of freedom models, for example). Then, the idea of ensem-
ble methods is to combining several of them together in order to create a strong learner (or
ensemble model) that achieves better performances.
Usually, ensemble models are used in order to :
• decrease the variance for bagging (Bootstrap Aggregating) technique
• reduce bias for the boosting technique
• improving the predictive force for stacking technique.
To understand these techniques, first, we will explore what is boostrapping and its different
hypothesis.
5.7.2 Bootstrapping
distribution so that sampling from the dataset is a good approximation of sampling from the
real distribution (representativity).
• Second, the size N of the dataset should be large enough compared to the size B of the
bootstrap samples so that samples are not too much correlated (independence).
Bootstrap samples are often used, for example, to evaluate variance or confidence intervals
of a statistical estimators. By definition, a statistical estimator is a function of some observa-
tions and, so, a random variable with variance coming from these observations. In order to
estimate the variance of such an estimator, we need to evaluate it on several independent sam-
ples drawn from the distribution of interest. In most of the cases, considering truly independent
samples would require too much data compared to the amount really available. We can then use
bootstrapping to generate several bootstrap samples that can be considered as being “almost-
representative” and “almost-independent” (almost i.i.d. samples). These bootstrap samples will
allow us to approximate the variance of the estimator, by evaluating its value for each of them.
Bootstrapping is often used to evaluate variance or confidence interval of some statistical esti-
mators.
5.7.3 Bagging
In parallel methods we fit the different considered learners independently from each others
and, so, it is possible to train them concurrently. The most famous such approach is “bagging”
(standing for “bootstrap aggregating”) that aims at producing an ensemble model that is more
robust than the individual models composing it.
When training a model, no matter if we are dealing with a classification or a regression problem,
we obtain a function that takes an input, returns an output and that is defined with respect to
the training dataset.
The idea of bagging is then simple: we want to fit several independent models and “average”
their predictions in order to obtain a model with a lower variance. However, we can’t, in
practice, fit fully independent models because it would require too much data. So, we rely on
the good “approximate properties” of bootstrap samples (representativity and independence)
to fit models that are almost independent.
First, we create multiple bootstrap samples so that each new bootstrap sample will act as
another (almost) independent dataset drawn from true distribution. Then, we can fit a weak
learner for each of these samples and finally aggregate them such that we kind of “aver-
age” their outputs and, so, obtain an ensemble model with less variance that its components.
Roughly speaking, as the bootstrap samples are approximatively independent and identically
distributed (i.i.d.), so are the learned base models. Then, “averaging” weak learners outputs
do not change the expected answer but reduce its variance.
So, assuming that we have L bootstrap samples (approximations of L independent datasets) of
size B denoted
and then aggregate them into some kind of averaging process in order to get an ensemble model
with a lower variance. For example, we can define our strong model such that
There are several possible ways to aggregate the multiple models fitted in parallel. - For a
regression problem, the outputs of individual models can literally be averaged to obtain the
output of the ensemble model. - For classification problem the class outputted by each model
can be seen as a vote and the class that receives the majority of the votes is returned by
the ensemble model (this is called hard-voting). Still for a classification problem, we can
also consider the probabilities of each classes returned by all the models, average these
probabilities and keep the class with the highest average probability (this is called soft-
voting). –> Averages or votes can either be simple or weighted if any relevant weights can be
used.
Finally, we can mention that one of the big advantages of bagging is that it can be parallelised.
As the different models are fitted independently from each others, intensive parallelisation tech-
niques can be used if required.
Bagging consists in fitting several base models on different bootstrap samples and build an
ensemble model that “average” the results of these weak learners.
Question : - Can you name an algorithms based on Bagging technique , Hint : leaf
###### Examples
Here, we are trying some example of stacking
• Bagged Decision Trees for Classification
import pandas
from sklearn import model_selection
from sklearn.ensemble import BaggingClassifier
from sklearn.tree import DecisionTreeClassifier
names = ['preg', 'plas', 'pres', 'skin', 'test', 'mass', 'pedi', 'age', 'class']
dataframe = pandas.read_csv("https://raw.githubusercontent.com/jbrownlee/Datasets/master/
˓→pima-indians-diabetes.data.csv",names=names)
array = dataframe.values
x = array[:,0:8]
y = array[:,8]
(continues on next page)
import pandas
from sklearn import model_selection
from sklearn.ensemble import RandomForestClassifier
names = ['preg', 'plas', 'pres', 'skin', 'test', 'mass', 'pedi', 'age', 'class']
dataframe = pandas.read_csv("https://raw.githubusercontent.com/jbrownlee/Datasets/master/
˓→pima-indians-diabetes.data.csv",names=names)
array = dataframe.values
x = array[:,0:8]
y = array[:,8]
Both of these algorithms will print, Accuracy: 0.77 (+/- 0.07). They are equivalent.
5.7.4 Boosting
In sequential methods the different combined weak models are no longer fitted indepen-
dently from each others. The idea is to fit models iteratively such that the training of model at
a given step depends on the models fitted at the previous steps. “Boosting” is the most famous
of these approaches and it produces an ensemble model that is in general less biased than the
weak learners that compose it.
Boosting methods work in the same spirit as bagging methods: we build a family of models
that are aggregated to obtain a strong learner that performs better.
However, unlike bagging that mainly aims at reducing variance, boosting is a technique
that consists in fitting sequentially multiple weak learners in a very adaptative way: each
model in the sequence is fitted giving more importance to observations in the dataset that
were badly handled by the previous models in the sequence. Intuitively, each new model
focus its efforts on the most difficult observations to fit up to now, so that we obtain, at the
end of the process, a strong learner with lower bias (even if we can notice that boosting can
also have the effect of reducing variance).
–> Boosting, like bagging, can be used for regression as well as for classification problems.
Being mainly focused at reducing bias, the base models that are often considered for boosting
are* *models with low variance but high bias. For example, if we want to usetreesas our base
models, we will choosemost of the time shallow decision trees with only a few depths.**
Another important reason that motivates the use of low variance but high bias models as weak
learners for boosting is that these models are in general less computationally expensive to fit
(few degrees of freedom when parametrised). Indeed, as computations to fit the different mod-
els can’t be done in parallel (unlike bagging), it could become too expensive to fit sequentially
several complex models.
Once the weak learners have been chosen, we still need to define how they will be sequentially
fitted and how they will be aggregated. We will discuss these questions in the two follow-
ing subsections, describing more especially two important boosting algorithms: adaboost and
gradient boosting.
In a nutshell, these two meta-algorithms differ on how they create and aggregate the weak
learners during the sequential process. Adaptive boosting updates the weights attached to
each of the training dataset observations whereas gradient boosting updates the value of
these observations. This main difference comes from the way both methods try to solve the
optimisation problem of finding the best model that can be written as a weighted sum of weak
learners.
Boosting consists in, iteratively, fitting a weak learner, aggregate it to the ensemble model and
“update” the training dataset to better take into account the strengths and weakness of the
current ensemble model when fitting the next base model.
1/ Adaptative boosting
In adaptative boosting (often called “adaboost”), we try to define our ensemble model as a
weighted sum of L weak learners
Finding the best ensemble model with this form is a difficult optimisation problem. Then,
instead of trying to solve it in one single shot (finding all the coefficients and weak learners
that give the best overall additive model), we make use of an iterative optimisation process
that is much more tractable, even if it can lead to a sub-optimal solution. More especially,
we add the weak learners one by one, looking at each iteration for the best possible pair
(coefficient, weak learner) to add to the current ensemble model. In other words, we define
recurrently the (s_l)’s such that
where c_l and w_l are chosen such that s_l is the model that fit the best the training data and,
so, that is the best possible improvement over s_(l-1). We can then denote
where E(.) is the fitting error of the given model and e(.,.) is the loss/error function. Thus,
instead of optimising “globally” over all the L models in the sum, we approximate the optimum
by optimising “locally” building and adding the weak learners to the strong model one by one.
More especially, when considering a binary classification, we can show that the adaboost algo-
rithm can be re-written into a process that proceeds as follow. First, it updates the observations
weights in the dataset and train a new weak learner with a special focus given to the obser-
vations misclassified by the current ensemble model. Second, it adds the weak learner to the
weighted sum according to an update coefficient that expresse the performances of this weak
model: the better a weak learner performs, the more it contributes to the strong learner.
So, assume that we are facing a binary classification problem, with N observations in our
dataset and we want to use adaboost algorithm with a given family of weak models. At the
very beginning of the algorithm (first model of the sequence), all the observations have
the same weights 1/N. Then, we repeat L times (for the L learners in the sequence) the
following steps:
fit the best possible weak model with the current observations weights
compute the value of the update coefficient that is some kind of scalar evaluation metric of the
weak learner that indicates how much this weak learner should be taken into account into the
ensemble model
update the strong learner by adding the new weak learner multiplied by its update coefficient
compute new observations weights that expresse which observations we would like to focus
on at the next iteration (weights of observations wrongly predicted by the aggregated model
increase and weights of the correctly predicted observations decrease)
Repeating these steps, we have then build sequentially our L models and aggregate them into
a simple linear combination weighted by coefficients expressing the performance of each
learner.
Notice that there exists variants of the initial adaboost algorithm such that LogitBoost (classifi-
cation) or L2Boost (regression) that mainly differ by their choice of loss function.
Adaboost updates weights of the observations at each iteration. Weights of well classified obser-
vations decrease relatively to weights of misclassified observations. Models that perform better
have higher weights in the final ensemble model.
2/ Gradient boosting
In gradient boosting, the ensemble model we try to build is also a weighted sum of weak
learners
Just as we mentioned for adaboost, finding the optimal model under this form is too difficult
and an iterative approach is required. The main difference with adaptative boosting is in
the definition of the sequential optimisation process. Indeed, gradient boosting casts the
problem into a gradient descent one: at each iteration we fit a weak learner to the opposite
of the gradient of the current fitting error with respect to the current ensemble model.
Let’s try to clarify this last point. First, theoretical gradient descent process over the ensemble
model can be written
where E(.) is the fitting error of the given model, c_l is a coefficient corresponding to the step
size and
This entity is the opposite of the gradient of the fitting error with respect to the ensemble
model at step l-1. This opposite of the gradient is a function that can, in practice, only be
evaluated for observations in the training dataset (for which we know inputs and outputs):
these evaluations are called pseudo-residuals attached to each observations. Moreover, even if
we know for the observations the values of these pseudo-residuals, we don’t want to add to
our ensemble model any kind of function: we only want to add a new instance of weak model.
So, the natural thing to do is to fit a weak learner to the pseudo-residuals computed for each
observation. Finally, the coefficient c_l is computed following a one dimensional optimisation
process (line-search to obtain the best step size c_l).
So, assume that we want to use gradient boosting technique with a given family of weak models.
At the very beginning of the algorithm (first model of the sequence), the pseudo-residuals are
set equal to the observation values. Then, we repeat L times (for the L models of the sequence)
the following steps:
fit the best possible weak model to pseudo-residuals (approximate the opposite of the gradient
with respect to the current strong learner)
compute the value of the optimal step size that defines by how much we update the ensemble
model in the direction of the new weak learner
update the ensemble model by adding the new weak learner multiplied by the step size (make
a step of gradient descent)
compute new pseudo-residuals that indicate, for each observation, in which direction we would
like to update next the ensemble model predictions
Repeating these steps, we have then build sequentially our L models and aggregate them fol-
lowing a gradient descent approach. Notice that, while adaptative boosting tries to solve at
each iteration exactly the “local” optimisation problem (find the best weak learner and its
coefficient to add to the strong model), gradient boosting uses instead a gradient descent
approach and can more easily be adapted to large number of loss functions. Thus, gradi-
ent boosting can be considered as a generalization of adaboost to arbitrary differentiable
loss functions.
Note There is an algorithm which gained huge popularity after a Kaggle’s competitions. It is
XGBoost (Extreme Gradient Boosting). This is a gradient boosting algorithm which has more
flexibility (varying number of terminal nodes and left weights) parameters to avoid sub-
learners correlations. Having these important qualities, XGBOOST is one of the most used
algorithm in data science. LIGHTGBM is a recent implementation of this algorithm. It was
published by Microsoft and it gives us the same scores (if parameters are equivalents) but it
runs quicker than a classic XGBOOST.
Gradient boosting updates values of the observations at each iteration. Weak learners are
trained to fit the pseudo-residuals that indicate in which direction to correct the current en-
semble model predictions to lower the error.
Examples
Here, we are trying an example of Boosting and compare it to a Bagging. Both of algorithms
take the same weak learners to build the macro-model
• Adaboost Classifier
breast_cancer = load_breast_cancer()
x = pd.DataFrame(breast_cancer.data, columns=breast_cancer.feature_names)
y = pd.Categorical.from_codes(breast_cancer.target, breast_cancer.target_names)
# Transforming string Target to an int
encoder = LabelEncoder()
binary_encoded_y = pd.Series(encoder.fit_transform(y))
breast_cancer = load_breast_cancer()
x = pd.DataFrame(breast_cancer.data, columns=breast_cancer.feature_names)
y = pd.Categorical.from_codes(breast_cancer.target, breast_cancer.target_names)
# Transforming string Target to an int
encoder = LabelEncoder()
binary_encoded_y = pd.Series(encoder.fit_transform(y))
Comparaison
Stacking mainly differ from bagging and boosting on two points : - First stacking often con-
siders heterogeneous weak learners (different learning algorithms are combined) whereas
bagging and boosting consider mainly homogeneous weak learners. - Second, stacking learns
to combine the base models using a meta-model whereas bagging and boosting combine weak
learners following deterministic algorithms.
As we already mentioned, the idea of stacking is to learn several different weak learners and
combine them by training a meta-model to output predictions based on the multiple predic-
tions returned by these weak models. So, we need to define two things in order to build our
stacking model: the L learners we want to fit and the meta-model that combines them.
For example, for a classification problem, we can choose as weak learners a KNN classifier, a
logistic regression and a SVM, and decide to learn a neural network as meta-model. Then, the
neural network will take as inputs the outputs of our three weak learners and will learn to
return final predictions based on it.
So, assume that we want to fit a stacking ensemble composed of L weak learners. Then we have
to follow the steps thereafter:
• split the training data in two folds
• choose L weak learners and fit them to data of the first fold
• for each of the L weak learners, make predictions for observations in the second fold
• fit the meta-model on the second fold, using predictions made by the weak learners
as inputs
In the previous steps, we split the dataset in two folds because predictions on data that have
been used for the training of the weak learners are not relevant for the training of the meta-
model.
Stacking consists in training a meta-model to produce outputs based on the outputs returned
by some lower layer weak learners.
A possible extension of stacking is multi-level stacking. It consists in doing stacking with
multiple layers. As an example,
Multi-level stacking considers several layers of stacking: some meta-models are trained on out-
puts returned by lower layer meta-models and so on. Here we have represented a 3-layers
stacking model.
Examples
Here, we are trying an example of Stacking and compare it to a Bagging & a Boosting. We note
that, many other applications (datasets) would show more difference between these techniques.
from sklearn.ensemble import AdaBoostClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.datasets import load_breast_cancer
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import accuracy_score
from sklearn.metrics import f1_score
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
breast_cancer = load_breast_cancer()
x = pd.DataFrame(breast_cancer.data, columns=breast_cancer.feature_names)
y = pd.Categorical.from_codes(breast_cancer.target, breast_cancer.target_names)
boosting_clf_ada_boost= AdaBoostClassifier(
DecisionTreeClassifier(max_depth=1),
n_estimators=3
)
bagging_clf_rf = RandomForestClassifier(n_estimators=200, max_depth=1,random_state=2020)
clf_logistic_reg = LogisticRegression(solver='liblinear',random_state=2020)
'''
This is a test class for stacking !
Please fill Free to change it to fit your needs
We suppose that at least the First N-1 Classifiers have
a predict_proba function.
'''
def __init__(self,classifiers):
if(len(classifiers) < 2):
raise numberOfClassifierException("You must fit your classifier with 2␣
˓→classifiers at least");
else:
self._classifiers = classifiers
def fit(self,data_x,data_y):
stacked_data_x = data_x.copy()
for classfier in self._classifiers[:-1]:
classfier.fit(data_x,data_y)
stacked_data_x = np.column_stack((stacked_data_x,classfier.predict_proba(data_
˓→ x)))
last_classifier = self._classifiers[-1]
last_classifier.fit(stacked_data_x,data_y)
def predict(self,data_x):
stacked_data_x = data_x.copy()
for classfier in self._classifiers[:-1]:
prob_predictions = classfier.predict_proba(data_x)
(continues on next page)
last_classifier = self._classifiers[-1]
return last_classifier.predict(stacked_data_x)
bagging_clf_rf.fit(train_x, train_y)
boosting_clf_ada_boost.fit(train_x, train_y)
classifers_list = [clf_rf,clf_ada_boost,clf_logistic_reg]
clf_stacking = Stacking(classifers_list)
clf_stacking.fit(train_x,train_y)
predictions_bagging = bagging_clf_rf.predict(test_x)
predictions_boosting = boosting_clf_ada_boost.predict(test_x)
predictions_stacking = clf_stacking.predict(test_x)
Comparaison
5.8.1 Introduction
Consider the 3-dimensional graph below in the context of a cost function. Our goal is to move
from the mountain in the top right corner (high cost) to the dark blue sea in the bottom
left (low cost). The arrows represent the direction of steepest descent (negative gradient)
from any given point–the direction that decreases the cost function as quickly as possible
Learning rate
The size of these steps is called the learning rate. With a high learning rate we can cover
more ground each step, but we risk overshooting the lowest point since the slope of the hill is
constantly changing. With a very low learning rate, we can confidently move in the direction
of the negative gradient since we are recalculating it so frequently. A low learning rate is
more precise, but calculating the gradient is time-consuming, so it will take us a very long
time to get to the bottom.
impacts of learning rate choice.
Cost function
A Loss Function (Error function) tells us “how good” our model is at making predictions for
a given set of parameters. The cost function has its own curve and its own gradients. The
slope of this curve tells us how to update our parameters to make the model more accurate.
To solve for the gradient, we iterate through our data points using our
:math:beta_1 and :math:beta_0 values and compute the
partial derivatives. This new gradient tells us the slope of our cost function at our cur-
rent position (current parameter values) and the direction we should move to update our
parameters. The size of our update is controlled by the learning rate.
Pseudocode of this algorithm
m : 1
b : 1
m_deriv : 0
b_deriv : 0
data_length : length(X)
loop i : 1 --> number_iterations:
loop i : 1 -> data_length :
m_deriv : m_deriv -X[i] * ((m*X[i] + b) - Y[i])
b_deriv : b_deriv - ((m*X[i] + b) - Y[i])
m : m - (m_deriv / data_length) * learning_rate
b : b - (b_deriv / data_length) * learning_rate
return m, b
There are three variants of gradient descent, which differ in how much data we use to
compute the gradient of the objective function. Depending on the amount of data, we make
a trade-off between the accuracy of the parameter update and the time it takes to perform
an update.
Batch gradient descent, known also as Vanilla gradient descent, computes the gradient of the
cost function with respect to the parameters 𝜃 for the entire training dataset :
𝜃 = 𝜃 − 𝜂 · ∇𝜃 𝐽(𝜃)
As we need to calculate the gradients for the whole dataset to perform just one update, batch
gradient descent can be very slow and is intractable for datasets that don’t fit in memory.
Batch gradient descent also doesn’t allow us to update our model online.
Stochastic gradient descent (SGD) in contrast performs a parameter update for each training
example 𝑥(𝑖) and label 𝑦 (𝑖)
• Choose an initial vector of parameters 𝑤 and learning rate 𝜂.
• Repeat until an approximate minimum is obtained:
– Randomly shuffle examples in the training set.
– For 𝑖 ∈ 1, . . . , 𝑛
(𝑖) (𝑖)
* 𝜃 = 𝜃 − 𝜂 · ∇𝜃 𝐽(𝜃; 𝑥 ; 𝑦 )
Batch gradient descent performs redundant computations for large datasets, as it recom-
putes gradients for similar examples before each parameter update. SGD does away with
this redundancy by performing one update at a time. It is therefore usually much faster
and can also be used to learn online. SGD performs frequent updates with a high variance
that cause the objective function to fluctuate heavily as in the image below.
SGD fluctuation.
While batch gradient descent converges to the minimum of the basin the parameters are
placed in, SGD’s fluctuation, on the one hand, enables it to jump to new and potentially
better local minima. On the other hand, this ultimately complicates convergence to the
exact minimum, as SGD will keep overshooting. However, it has been shown that when we
slowly decrease the learning rate, SGD shows the same convergence behaviour as batch
gradient descent, almost certainly converging to a local or the global minimum for non-
convex and convex optimization respectively.
Mini-batch gradient descent finally takes the best of both worlds and performs an update for
every mini-batch of n training examples:
This way, it :
• reduces the variance of the parameter updates, which can lead to more stable con-
vergence.
• can make use of highly optimized matrix optimizations common to state-of-the-art deep
learning libraries that make computing the gradient very efficient. Common mini-batch
sizes range between 50 and 256, but can vary for different applications.
Mini-batch gradient descent is typically the algorithm of choice when training a neural
network.
Vanilla mini-batch gradient descent, however, does not guarantee good convergence, but offers
a few challenges that need to be addressed:
• Choosing a proper learning rate can be difficult. A learning rate that is too small leads
to painfully slow convergence, while a learning rate that is too large can hinder conver-
gence and cause the loss function to fluctuate around the minimum or even to diverge.
• Learning rate schedules try to adjust the learning rate during training by e.g. an-
nealing, i.e. reducing the learning rate according to a pre-defined schedule or when the
change in objective between epochs falls below a threshold. These schedules and thresh-
olds, however, have to be defined in advance and are thus unable to adapt to a dataset’s
characteristics.
• Additionally, the same learning rate applies to all parameter updates. If our data is sparse
and our features have very different frequencies, we might not want to update all of
them to the same extent, but perform a larger update for rarely occurring features.
• Another key challenge of minimizing highly non-convex error functions common for
neural networks is avoiding getting trapped in their numerous suboptimal local min-
ima. These saddle points (local minimas) are usually surrounded by a plateau of the
same error, which makes it notoriously hard for SGD to escape, as the gradient is close
to zero in all dimensions.
In the following, we will outline some algorithms that are widely used by the deep learning
community to deal with the aforementioned challenges.
Momentum
SGD has trouble navigating ravines (areas where the surface curves much more steeply in
one dimension than in another), which are common around local optima. In these scenarios,
SGD oscillates across the slopes of the ravine while only making hesitant progress along
the bottom towards the local optimum as in the image below.
𝑣𝑡 = 𝜌𝑣𝑡−1 + ∇𝜃 𝐽(𝜃)
(5.52)
𝜃 = 𝜃 − 𝑣𝑡
vx = 0
while True:
dx = gradient(J, x)
(continues on next page)
Note: The momentum term :math:`rho` is usually set to 0.9 or a similar value.
Essentially, when using momentum, we push a ball down a hill. The ball accumulates mo-
mentum as it rolls downhill, becoming faster and faster on the way (until it reaches its
terminal velocity if there is air resistance, i.e. :math:`rho` <1 ).
The same thing happens to our parameter updates: The momentum term increases for dimen-
sions whose gradients point in the same directions and reduces updates for dimensions
whose gradients change directions. As a result, we gain faster convergence and reduced
oscillation.
• Added element-wise scaling of the gradient based on the historical sum of squares in each
dimension.
• “Per-parameter learning rates” or “adaptive learning rates”
grad_squared = 0
while True:
dx = gradient(J, x)
grad_squared += dx * dx
x -= learning_rate * dx / (np.sqrt(grad_squared) + 1e-7)
grad_squared = 0
while True:
dx = gradient(J, x)
grad_squared += decay_rate * grad_squared + (1 - decay_rate) * dx * dx
x -= learning_rate * dx / (np.sqrt(grad_squared) + 1e-7)
However, a ball that rolls down a hill, blindly following the slope, is highly unsatisfactory. We’d
like to have a smarter ball, a ball that has a notion of where it is going so that it knows to
slow down before the hill slopes up again. Nesterov accelerated gradient (NAG) is a way to
give our momentum term this kind of prescience. We know that we will use our momentum
term 𝛾𝑣𝑡−1 to move the parameters 𝜃.
Computing 𝜃 − 𝛾𝑣𝑡−1 thus gives us an approximation of the next position of the parameters
(the gradient is missing for the full update), a rough idea where our parameters are going to
be. We can now effectively look ahead by calculating the gradient not w.r.t. to our current
parameters 𝜃 but w.r.t. the approximate future position of our parameters:
Again, we set the momentum term 𝛾 to a value of around 0.9. While Momentum first com-
putes the current gradient and then takes a big jump in the direction of the updated
accumulated gradient , NAG first makes a big jump in the direction of the previous ac-
cumulated gradient, measures the gradient and then makes a correction, which results
in the complete NAG update. This anticipatory update prevents us from going too fast and
results in increased responsiveness, which has significantly increased the performance of
RNNs on a number of tasks
Adam
Adaptive Moment Estimation (Adam) is a method that computes adaptive learning rates for
each parameter. In addition to storing an exponentially decaying average of past squared
gradients :math:`v_t`, Adam also keeps an exponentially decaying average of past gradi-
ents :math:`m_t`, similar to momentum. Whereas momentum can be seen as a ball running
down a slope, Adam behaves like a heavy ball with friction, which thus prefers flat minima
in the error surface. We compute the decaying averages of past and past squared gradients 𝑚𝑡
and 𝑣𝑡 respectively as follows:
𝑚𝑡 and 𝑣𝑡 are estimates of the first moment (the mean) and the second moment (the uncentered
variance) of the gradients respectively, hence the name of the method. Adam (almost)
first_moment = 0
second_moment = 0
while True:
dx = gradient(J, x)
# Momentum:
first_moment = beta1 * first_moment + (1 - beta1) * dx
# AdaGrad/RMSProp
second_moment = beta2 * second_moment + (1 - beta2) * dx * dx
x -= learning_rate * first_moment / (np.sqrt(second_moment) + 1e-7)
As 𝑚𝑡 and 𝑣𝑡 are initialized as vectors of 0’s, the authors of Adam observe that they are biased
towards zero, especially during the initial time steps, and especially when the decay rates are
small (i.e. 𝛽1 and 𝛽2 are close to 1). They counteract these biases by computing bias-corrected
first and second moment estimates:
𝑚𝑡
𝑚
ˆ𝑡 = (5.55)
1 − 𝛽1𝑡
𝑣𝑡
𝑣ˆ𝑡 = (5.56)
1 − 𝛽2𝑡
They then use these to update the parameters (Adam update rule):
𝜂
𝜃𝑡+1 = 𝜃𝑡 − √ 𝑚
ˆ𝑡
𝑣ˆ𝑡 + 𝜖
• 𝑚
ˆ 𝑡 Accumulate gradient: velocity.
• 𝑣ˆ𝑡 Element-wise scaling of the gradient based on the historical sum of squares in each
dimension.
• Choose Adam as default optimizer
• Default values of 0.9 for 𝛽1 , 0.999 for 𝛽2 , and 10−7 for 𝜖.
• learning rate in a range between 1𝑒 − 3 and 5𝑒 − 4
SIX
DEEP LEARNING
6.1 Backpropagation
%matplotlib inline
Y = max(XW(1) , 0)W(2)
A fully-connected ReLU network with one hidden layer and no biases, trained to predict y from
x using Euclidean error.
Chaine rule
𝑥 → 𝑧 (1) = 𝑥𝑇 𝑤(1) → ℎ(1) = max(𝑧 (1) , 0) → 𝑧 (2) = ℎ(1)𝑇 𝑤(2) → 𝐿(𝑧 (2) , 𝑦) = (𝑧 (2) − 𝑦)2
𝑤(1) ↗ 𝑤(2) ↗
𝜕𝑧 (1) 𝜕ℎ(1) 1 if 𝑧 (1) >0 𝜕𝑧
(2) 𝜕𝐿
=𝑥 = {else 0 = ℎ(1) = 2(𝑧 (2) − 𝑦)
𝜕𝑤(1) 𝜕𝑧 (1) 𝜕𝑤(2) 𝜕𝑧 (2)
𝜕𝑧 (1) 𝜕𝑧 (2)
= 𝑤(1) = 𝑤(2)
𝜕𝑥 𝜕ℎ(1)
Backward: compute gradient of the loss given each parameters vectors applying chaine rule
from the loss downstream to the parameters:
For 𝑤(2) :
257
Statistics and Machine Learning in Python, Release 0.3 beta
𝜕𝐿 𝜕𝐿 𝜕𝑧 (2)
= (6.1)
𝜕𝑤(2) 𝜕𝑧 (2) 𝜕𝑤(2)
=2(𝑧 (2) − 𝑦)ℎ(1) (6.2)
For 𝑤(1) :
Given a function 𝑧 = 𝑥 with 𝑧 the output, 𝑥 the input and 𝑤 the coeficients.
• Scalar to Scalar: 𝑥 ∈ R, 𝑧 ∈ R, 𝑤 ∈ R
Regular derivative:
𝜕𝑧
=𝑥∈R
𝜕𝑤
If 𝑤 changes by a small amount, how much will 𝑧 change?
• Vector to Scalar: 𝑥 ∈ R𝑁 , 𝑧 ∈ R, 𝑤 ∈ R𝑁
𝜕𝑧
Derivative is Gradient of partial derivative: 𝜕𝑤 ∈ R𝑁
⎡ 𝜕𝑧 ⎤
𝜕𝑤1
⎢ .. ⎥
⎢ . ⎥
𝜕𝑧 ⎢ 𝜕𝑧 ⎥
= ∇𝑤 𝑧 = ⎢
⎢ 𝜕𝑤𝑖
⎥ (6.5)
𝜕𝑤 ⎢ ..
⎥
⎥
⎣ . ⎦
𝜕𝑧
𝜕𝑤𝑁
For each element 𝑤𝑖 of 𝑤, if it changes by a small amount then how much will y change?
• Vector to Vector: 𝑤 ∈ R𝑁 , 𝑧 ∈ R𝑀
Derivative is Jacobian of partial derivative:
TO COMPLETE
𝜕𝑧
𝜕𝑤 ∈ R𝑁 ×𝑀
Backpropagation summary
Backpropagation algorithm in a graph: 1. Forward pass, for each node compute local partial
derivatives of ouput given inputs 2. Backward pass: apply chain rule from the end to each
parameters - Update parameter with gradient descent using the current upstream gradient and
the current local gradient - Compute upstream gradient for the backward nodes
Think locally and remember that at each node: - For the loss the gradient is the error - At each
step, the upstream gradient is obtained by multiplying the upstream gradient (an error) with
the current parameters (vector of matrix). - At each step, the current local gradient equal the
input, therfore the current update is the current upstream gradient time the input.
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import sklearn.model_selection
iris = sns.load_dataset("iris")
#g = sns.pairplot(iris, hue="species")
df = iris[iris.species != "setosa"]
g = sns.pairplot(df, hue="species")
df['species_n'] = iris.species.map({'versicolor':1, 'virginica':2})
# Scale
from sklearn.preprocessing import StandardScaler
scalerx, scalery = StandardScaler(), StandardScaler()
X_iris = scalerx.fit_transform(X_iris)
Y_iris = StandardScaler().fit_transform(Y_iris)
/home/edouard/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:5:␣
˓→SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
"""
This implementation uses numpy to manually compute the forward pass, loss, and backward
pass.
W1 = np.random.randn(D_in, H)
W2 = np.random.randn(H, D_out)
learning_rate = lr
for t in range(nite):
# Forward pass: compute predicted y
z1 = X.dot(W1)
h1 = np.maximum(z1, 0)
Y_pred = h1.dot(W2)
# Update weights
W1 -= learning_rate * grad_w1
W2 -= learning_rate * grad_w2
losses_tr.append(loss)
losses_val.append(loss_val)
if t % 10 == 0:
print(t, loss, loss_val)
0 15126.224825529907 2910.260853330454
10 71.5381374591153 104.97056197642135
20 50.756938353833334 80.02800827986354
30 46.546510744624236 72.85211241738614
40 44.41413064447564 69.31127324764276
[<matplotlib.lines.Line2D at 0x7f960cf5e9b0>,
<matplotlib.lines.Line2D at 0x7f960cf5eb00>]
source
Numpy is a great framework, but it cannot utilize GPUs to accelerate its numerical compu-
tations. For modern deep neural networks, GPUs often provide speedups of 50x or greater,
so unfortunately numpy won’t be enough for modern deep learning. Here we introduce the
most fundamental PyTorch concept: the Tensor. A PyTorch Tensor is conceptually identical to a
numpy array: a Tensor is an n-dimensional array, and PyTorch provides many functions for op-
erating on these Tensors. Behind the scenes, Tensors can keep track of a computational graph
and gradients, but they’re also useful as a generic tool for scientific computing. Also unlike
numpy, PyTorch Tensors can utilize GPUs to accelerate their numeric computations. To run a
PyTorch Tensor on GPU, you simply need to cast it to a new datatype. Here we use PyTorch
Tensors to fit a two-layer network to random data. Like the numpy example above we need to
manually implement the forward and backward passes through the network:
import torch
dtype = torch.float
device = torch.device("cpu")
# device = torch.device("cuda:0") # Uncomment this to run on GPU
learning_rate = lr
for t in range(nite):
# Forward pass: compute predicted y
z1 = X.mm(W1)
h1 = z1.clamp(min=0)
y_pred = h1.mm(W2)
losses_tr.append(loss)
losses_val.append(loss_val)
if t % 10 == 0:
print(t, loss, loss_val)
0 8086.1591796875 5429.57275390625
10 225.77589416503906 331.83734130859375
20 86.46501159667969 117.72447204589844
30 52.375606536865234 73.84156036376953
40 43.16458511352539 64.0667495727539
[<matplotlib.lines.Line2D at 0x7f960033c470>,
<matplotlib.lines.Line2D at 0x7f960033c5c0>]
source
A fully-connected ReLU network with one hidden layer and no biases, trained to predict y
from x by minimizing squared Euclidean distance. This implementation computes the forward
pass using operations on PyTorch Tensors, and uses PyTorch autograd to compute gradients.
A PyTorch Tensor represents a node in a computational graph. If x is a Tensor that has x.
requires_grad=True then x.grad is another Tensor holding the gradient of x with respect to
some scalar value.
import torch
dtype = torch.float
device = torch.device("cpu")
# device = torch.device("cuda:0") # Uncomment this to run on GPU
learning_rate = lr
for t in range(nite):
# Forward pass: compute predicted y using operations on Tensors; these
# are exactly the same operations we used to compute the forward pass using
# Tensors, but we do not need to keep references to intermediate values since
# we are not implementing the backward pass by hand.
y_pred = X.mm(W1).clamp(min=0).mm(W2)
# Use autograd to compute the backward pass. This call will compute the
# gradient of loss with respect to all Tensors with requires_grad=True.
# After this call w1.grad and w2.grad will be Tensors holding the gradient
# of the loss with respect to w1 and w2 respectively.
loss.backward()
y_pred = X_val.mm(W1).clamp(min=0).mm(W2)
if t % 10 == 0:
print(t, loss.item(), loss_val.item())
losses_tr.append(loss.item())
losses_val.append(loss_val.item())
0 8307.1806640625 2357.994873046875
10 111.97289276123047 250.04209899902344
20 65.83244323730469 201.63694763183594
30 53.70908737182617 183.17051696777344
40 48.719329833984375 173.3616943359375
[<matplotlib.lines.Line2D at 0x7f95ff2ad978>,
<matplotlib.lines.Line2D at 0x7f95ff2adac8>]
source
This implementation uses the nn package from PyTorch to build the network. PyTorch autograd
makes it easy to define computational graphs and take gradients, but raw autograd can be a bit
too low-level for defining complex neural networks; this is where the nn package can help. The
nn package defines a set of Modules, which you can think of as a neural network layer that has
produces output from input and may have some trainable weights.
import torch
X = torch.from_numpy(X)
Y = torch.from_numpy(Y)
X_val = torch.from_numpy(X_val)
Y_val = torch.from_numpy(Y_val)
learning_rate = lr
for t in range(nite):
# Forward pass: compute predicted y by passing x to the model. Module objects
# override the __call__ operator so you can call them like functions. When
# doing so you pass a Tensor of input data to the Module and it produces
# a Tensor of output data.
y_pred = model(X)
# Compute and print loss. We pass Tensors containing the predicted and true
# values of y, and the loss function returns a Tensor containing the
# loss.
loss = loss_fn(y_pred, Y)
# Backward pass: compute gradient of the loss with respect to all the learnable
# parameters of the model. Internally, the parameters of each Module are stored
# in Tensors with requires_grad=True, so this call will compute gradients for
# all learnable parameters in the model.
loss.backward()
if t % 10 == 0:
print(t, loss.item(), loss_val.item())
losses_tr.append(loss.item())
losses_val.append(loss_val.item())
lr=1e-4, nite=50)
0 82.32025146484375 91.3389892578125
10 50.322200775146484 63.563087463378906
20 40.825225830078125 57.13555145263672
30 37.53572082519531 55.74506378173828
40 36.191200256347656 55.499732971191406
[<matplotlib.lines.Line2D at 0x7f95ff296668>,
<matplotlib.lines.Line2D at 0x7f95ff2967b8>]
This implementation uses the nn package from PyTorch to build the network. Rather than man-
ually updating the weights of the model as we have been doing, we use the optim package to
define an Optimizer that will update the weights for us. The optim package defines many op-
timization algorithms that are commonly used for deep learning, including SGD+momentum,
RMSProp, Adam, etc.
import torch
X = torch.from_numpy(X)
Y = torch.from_numpy(Y)
X_val = torch.from_numpy(X_val)
Y_val = torch.from_numpy(Y_val)
# Use the optim package to define an Optimizer that will update the weights of
(continues on next page)
# Before the backward pass, use the optimizer object to zero all of the
# gradients for the variables it will update (which are the learnable
# weights of the model). This is because by default, gradients are
# accumulated in buffers( i.e, not overwritten) whenever .backward()
# is called. Checkout docs of torch.autograd.backward for more details.
optimizer.zero_grad()
with torch.no_grad():
y_pred = model(X_val)
loss_val = loss_fn(y_pred, Y_val)
if t % 10 == 0:
print(t, loss.item(), loss_val.item())
losses_tr.append(loss.item())
losses_val.append(loss_val.item())
0 92.271240234375 83.96189880371094
10 64.25907135009766 59.872535705566406
20 47.6252555847168 50.228126525878906
30 40.33802032470703 50.60377502441406
40 38.19448471069336 54.03163528442383
[<matplotlib.lines.Line2D at 0x7f95ff200080>,
<matplotlib.lines.Line2D at 0x7f95ff2001d0>]
%matplotlib inline
import os
import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torch.optim import lr_scheduler
import torchvision
from torchvision import transforms
from torchvision import datasets
from torchvision import models
#
from pathlib import Path
import matplotlib.pyplot as plt
# Device configuration
device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')
print(device)
cuda:0
Hyperparameters
train_loader = torch.utils.data.DataLoader(
datasets.MNIST('data', train=True, download=True,
transform=transforms.Compose([
transforms.ToTensor(),
transforms.Normalize((0.1307,), (0.3081,)) # Mean and Std of␣
˓→the MNIST dataset
])),
batch_size=batch_size_train, shuffle=True)
val_loader = torch.utils.data.DataLoader(
datasets.MNIST('data', train=False, transform=transforms.Compose([
transforms.ToTensor(),
transforms.Normalize((0.1307,), (0.3081,)) # Mean and Std of the MNIST dataset
])),
batch_size=batch_size_test, shuffle=True)
return train_loader, val_loader
So one test data batch is a tensor of shape: . This means we have 1000 examples of 28x28 pixels
in grayscale (i.e. no rgb channels, hence the one). We can plot some of them using matplotlib.
𝑓 (𝑥) = 𝜎(𝑥𝑇 𝑤)
𝑓 (𝑥) = softmax(𝑥𝑇 𝑊 + 𝑏)
X_train = train_loader.dataset.data.numpy()
#print(X_train.shape)
X_train = X_train.reshape((X_train.shape[0], -1))
(continues on next page)
X_test = val_loader.dataset.data.numpy()
X_test = X_test.reshape((X_test.shape[0], -1))
y_test = val_loader.dataset.targets.numpy()
print(X_train.shape, y_train.shape)
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)
coef = clf.coef_.copy()
plt.figure(figsize=(10, 5))
scale = np.abs(coef).max()
for i in range(10):
l1_plot = plt.subplot(2, 5, i + 1)
l1_plot.imshow(coef[i].reshape(28, 28), interpolation='nearest',
cmap=plt.cm.RdBu, vmin=-scale, vmax=scale)
l1_plot.set_xticks(())
l1_plot.set_yticks(())
l1_plot.set_xlabel('Class %i' % i)
plt.suptitle('Classification vector for...')
plt.show()
mlp.fit(X_train, y_train)
print("Training set score: %f" % mlp.score(X_train, y_train))
print("Test set score: %f" % mlp.score(X_test, y_test))
plt.show()
/home/ed203246/anaconda3/lib/python3.7/site-packages/sklearn/neural_network/multilayer_
˓→perceptron.py:562: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (5)␣
˓→reached and the optimization hasn't converged yet.
% self.max_iter, ConvergenceWarning)
class TwoLayerMLP(nn.Module):
# %load train_val_model.py
# %load train_val_model.py
import numpy as np
import torch
import time
import copy
best_model_wts = copy.deepcopy(model.state_dict())
best_acc = 0.0
running_loss = 0.0
running_corrects = 0
# forward
# track history if only in train
with torch.set_grad_enabled(phase == 'train'):
outputs = model(inputs)
_, preds = torch.max(outputs, 1)
loss = criterion(outputs, labels)
# statistics
running_loss += loss.item() * inputs.size(0)
running_corrects += torch.sum(preds == labels.data)
#nsamples = dataloaders[phase].dataset.data.shape[0]
epoch_loss = running_loss / nsamples
epoch_acc = running_corrects.double() / nsamples
losses[phase].append(epoch_loss)
accuracies[phase].append(epoch_acc)
if log_interval is not None and epoch % log_interval == 0:
print('{} Loss: {:.4f} Acc: {:.2f}%'.format(
phase, epoch_loss, 100 * epoch_acc))
print(next(model.parameters()).is_cuda)
torch.save(model.state_dict(), 'models/mod-%s.pth' % model.__class__.__name__)
True
torch.Size([50, 784])
torch.Size([50])
torch.Size([10, 50])
torch.Size([10])
Total number of parameters = 39760
Epoch 0/0
----------
train Loss: 0.4472 Acc: 87.65%
val Loss: 0.3115 Acc: 91.25%
Use the model to make new predictions. Consider the device, ie, load data on device
example_data.to(device) from prediction, then move back to cpu example_data.cpu().
with torch.no_grad():
output = model(example_data).cpu()
example_data = example_data.cpu()
(continues on next page)
# print(output.is_cuda)
# Softmax predictions
preds = output.argmax(dim=1)
Continue training from checkpoints: reload the model and run 10 more epochs
Epoch 0/9
----------
train Loss: 0.3088 Acc: 91.12%
val Loss: 0.2877 Acc: 91.92%
Epoch 2/9
----------
train Loss: 0.2847 Acc: 91.97%
val Loss: 0.2797 Acc: 92.05%
Epoch 4/9
----------
train Loss: 0.2743 Acc: 92.30%
val Loss: 0.2797 Acc: 92.11%
Epoch 6/9
----------
train Loss: 0.2692 Acc: 92.46%
(continues on next page)
Epoch 8/9
----------
train Loss: 0.2643 Acc: 92.66%
val Loss: 0.2684 Acc: 92.44%
• Define a MultiLayerMLP([D_in, 512, 256, 128, 64, D_out]) class that take the size of
the layers as parameters of the constructor.
• Add some non-linearity with relu acivation function
class MLP(nn.Module):
Epoch 0/9
----------
train Loss: 1.2111 Acc: 61.88%
val Loss: 0.3407 Acc: 89.73%
Epoch 2/9
----------
train Loss: 0.1774 Acc: 94.74%
val Loss: 0.1510 Acc: 95.47%
Epoch 4/9
----------
train Loss: 0.0984 Acc: 97.16%
val Loss: 0.1070 Acc: 96.76%
Epoch 6/9
----------
train Loss: 0.0636 Acc: 98.14%
val Loss: 0.0967 Acc: 96.98%
Epoch 8/9
----------
train Loss: 0.0431 Acc: 98.75%
val Loss: 0.0822 Acc: 97.55%
Reduce the size of the training dataset by considering only 10 minibatche for size16.
train_size = 10 * 16
# Stratified sub-sampling
targets = train_loader.dataset.targets.numpy()
nclasses = len(set(targets))
Train size= 160 Train label count= {0: 16, 1: 16, 2: 16, 3: 16, 4: 16, 5: 16, 6: 16, 7:␣
˓→16, 8: 16, 9: 16}
Batch sizes= [16, 16, 16, 16, 16, 16, 16, 16, 16, 16]
Datasets shape {'train': torch.Size([60000, 28, 28]), 'val': torch.Size([10000, 28, 28])}
N input features 784 N output 10
Epoch 0/99
----------
train Loss: 2.3066 Acc: 9.38%
val Loss: 2.3058 Acc: 10.34%
Epoch 20/99
----------
train Loss: 2.1213 Acc: 58.13%
val Loss: 2.1397 Acc: 51.34%
Epoch 40/99
----------
train Loss: 0.4651 Acc: 88.75%
val Loss: 0.8372 Acc: 73.63%
Epoch 60/99
----------
train Loss: 0.0539 Acc: 100.00%
val Loss: 0.8384 Acc: 75.46%
Epoch 80/99
----------
train Loss: 0.0142 Acc: 100.00%
val Loss: 0.9417 Acc: 75.55%
Epoch 0/99
----------
train Loss: 2.2523 Acc: 20.62%
val Loss: 2.0853 Acc: 45.51%
Epoch 20/99
----------
train Loss: 0.0010 Acc: 100.00%
val Loss: 1.0113 Acc: 78.08%
Epoch 40/99
----------
train Loss: 0.0002 Acc: 100.00%
val Loss: 1.1456 Acc: 78.12%
Epoch 60/99
----------
train Loss: 0.0001 Acc: 100.00%
val Loss: 1.2630 Acc: 77.98%
Epoch 80/99
----------
train Loss: 0.0000 Acc: 100.00%
val Loss: 1.3446 Acc: 77.87%
(continues on next page)
The CIFAR-10 dataset consists of 60000 32x32 colour images in 10 classes, with 6000 images
per class. There are 50000 training images and 10000 test images.
The dataset is divided into five training batches and one test batch, each with 10000 images.
The test batch contains exactly 1000 randomly-selected images from each class. The training
batches contain the remaining images in random order, but some training batches may contain
more images from one class than another. Between them, the training batches contain exactly
5000 images from each class.
Here are the classes in the dataset, as well as 10 random images from each: - airplane
- automobile
- bird
- cat
- deer
- dog
- frog
- horse
- ship
- truck
import numpy as np
import torch
import torch.nn as nn
import torchvision
import torchvision.transforms as transforms
# Device configuration
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
# Hyper-parameters
num_epochs = 5
learning_rate = 0.001
# CIFAR-10 dataset
train_dataset = torchvision.datasets.CIFAR10(root='data/',
train=True,
transform=transform,
download=True)
val_dataset = torchvision.datasets.CIFAR10(root='data/',
train=False,
transform=transforms.ToTensor())
# Data loader
train_loader = torch.utils.data.DataLoader(dataset=train_dataset,
batch_size=100,
shuffle=True)
val_loader = torch.utils.data.DataLoader(dataset=val_dataset,
batch_size=100,
shuffle=False)
Epoch 0/49
----------
train Loss: 2.0171 Acc: 24.24%
val Loss: 1.8761 Acc: 30.84%
Epoch 10/49
----------
train Loss: 1.5596 Acc: 43.70%
val Loss: 1.5853 Acc: 43.07%
Epoch 20/49
----------
train Loss: 1.4558 Acc: 47.59%
val Loss: 1.4210 Acc: 48.88%
Epoch 30/49
----------
train Loss: 1.3904 Acc: 49.79%
val Loss: 1.3890 Acc: 50.16%
Epoch 40/49
----------
train Loss: 1.3497 Acc: 51.24%
val Loss: 1.3625 Acc: 51.41%
6.3.1 Outline
2. Architecures
3. Train and test functions
4. CNN models
5. MNIST
6. CIFAR-10
Sources:
Deep learning - cs231n.stanford.edu
CNN - Stanford cs231n
Pytorch - WWW tutorials - github tutorials - github examples
MNIST and pytorch: - MNIST nextjournal.com/gkoehler/pytorch-mnist - MNIST
github/pytorch/examples - MNIST kaggle
6.3.2 Architectures
Sources:
• cv-tricks.com
• [zhenye-na.github.io(]https://zhenye-na.github.io/2018/12/01/cnn-deep-leearning-ai-
week2.html)
LeNet
Fig. 1: LeNet
AlexNet
Fig. 2: AlexNet
• Deeper, bigger,
• Featured Convolutional Layers stacked on top of each other (previously it was common to
only have a single CONV layer always immediately followed by a POOL layer).
• ReLu(Rectified Linear Unit) for the non-linear part, instead of a Tanh or Sigmoid.
The advantage of the ReLu over sigmoid is that it trains much faster than the latter because the
derivative of sigmoid becomes very small in the saturating region and therefore the updates to
the weights almost vanish. This is called vanishing gradient problem.
• Dropout: reduces the over-fitting by using a Dropout layer after every FC layer. Dropout
layer has a probability,(p), associated with it and is applied at every neuron of the response
map separately. It randomly switches off the activation with the probability p.
Fig. 4: Dropout
fitting. Another view of DropOut being helpful is that since neurons are randomly chosen, they
tend to avoid developing co-adaptations among themselves thereby enabling them to develop
meaningful features, independent of others.
• Data augmentation is carried out to reduce over-fitting. This Data augmentation includes
mirroring and cropping the images to increase the variation in the training data-set.
GoogLeNet. (Szegedy et al. from Google 2014) was a Convolutional Network . Its main contri-
bution was the development of an
• Inception Module that dramatically reduced the number of parameters in the network
(4M, compared to AlexNet with 60M).
• There are also several followup versions to the GoogLeNet, most recently Inception-v4.
VGGNet. (Karen Simonyan and Andrew Zisserman 2014)
• 16 CONV/FC layers and, appealingly, features an extremely homogeneous architecture.
• Only performs 3x3 convolutions and 2x2 pooling from the beginning to the end. Replace
large kernel-sized filters(11 and 5 in the first and second convolutional layer, respectively)
with multiple 3X3 kernel-sized filters one after another.
With a given receptive field(the effective area size of input image on which output depends),
multiple stacked smaller size kernel is better than the one with a larger size kernel because
multiple non-linear layers increases the depth of the network which enables it to learn more
complex features, and that too at a lower cost. For example, three 3X3 filters on top of each
other with stride 1 ha a receptive size of 7, but the number of parameters involved is 3*(9^2)
in comparison to 49^2 parameters of kernels with a size of 7.
• Lot more memory and parameters (140M)
ResNet. (Kaiming He et al. 2015)
Resnet block variants (Source):
• Skip connections
Fig. 6: VGGNet
Fig. 9: ResNet 18
• Batch normalization.
• State of the art CNN models and are the default choice (as of May 10, 2016). In particular,
also see more
• Recent developments that tweak the original architecture from Kaiming He et al. Identity
Mappings in Deep Residual Networks (published March 2016).
Models in pytorch
%matplotlib inline
import os
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from torch.optim import lr_scheduler
import torchvision
import torchvision.transforms as transforms
from torchvision import models
#
from pathlib import Path
import matplotlib.pyplot as plt
# Device configuration
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
# %load train_val_model.py
import numpy as np
import torch
import time
import copy
best_model_wts = copy.deepcopy(model.state_dict())
best_acc = 0.0
running_loss = 0.0
running_corrects = 0
# forward
# track history if only in train
with torch.set_grad_enabled(phase == 'train'):
outputs = model(inputs)
_, preds = torch.max(outputs, 1)
loss = criterion(outputs, labels)
# statistics
running_loss += loss.item() * inputs.size(0)
running_corrects += torch.sum(preds == labels.data)
#nsamples = dataloaders[phase].dataset.data.shape[0]
epoch_loss = running_loss / nsamples
epoch_acc = running_corrects.double() / nsamples
losses[phase].append(epoch_loss)
accuracies[phase].append(epoch_acc)
if log_interval is not None and epoch % log_interval == 0:
print('{} Loss: {:.4f} Acc: {:.2f}%'.format(
phase, epoch_loss, 100 * epoch_acc))
LeNet-5
class LeNet5(nn.Module):
"""
layers: (nb channels in input layer,
nb channels in 1rst conv,
nb channels in 2nd conv,
nb neurons for 1rst FC: TO BE TUNED,
nb neurons for 2nd FC,
nb neurons for 3rd FC,
nb neurons output FC TO BE TUNED)
"""
def __init__(self, layers = (1, 6, 16, 1024, 120, 84, 10), debug=False):
super(LeNet5, self).__init__()
self.layers = layers
self.debug = debug
self.conv1 = nn.Conv2d(layers[0], layers[1], 5, padding=2)
self.conv2 = nn.Conv2d(layers[1], layers[2], 5)
self.fc1 = nn.Linear(layers[3], layers[4])
self.fc2 = nn.Linear(layers[4], layers[5])
self.fc3 = nn.Linear(layers[5], layers[6])
class MiniVGGNet(torch.nn.Module):
def __init__(self, layers=(1, 16, 32, 1024, 120, 84, 10), debug=False):
super(MiniVGGNet, self).__init__()
self.layers = layers
self.debug = debug
# Conv block 1
(continues on next page)
# Conv block 2
self.conv21 = nn.Conv2d(in_channels=layers[1], out_channels=layers[2], kernel_
˓→size=3,
stride=1, padding=0, bias=True)
self.conv22 = nn.Conv2d(in_channels=layers[2], out_channels=layers[2], kernel_
˓→size=3,
stride=1, padding=1, bias=True)
x = F.relu(self.conv21(x))
x = F.relu(self.conv22(x))
x = F.max_pool2d(x, 2)
if self.debug:
print("### DEBUG: Shape of last convnet=", x.shape[1:], ". FC size=", np.
˓→prod(x.shape[1:]))
x = x.view(-1, self.layers[3])
x = F.relu(self.fc1(x))
x = F.relu(self.fc2(x))
x = self.fc3(x)
ResNet-like Model:
# ---------------------------------------------------------------------------- #
# An implementation of https://arxiv.org/pdf/1512.03385.pdf #
# See section 4.2 for the model architecture on CIFAR-10 #
# Some part of the code was referenced from below #
# https://github.com/pytorch/vision/blob/master/torchvision/models/resnet.py #
# ---------------------------------------------------------------------------- #
import torch.nn as nn
# 3x3 convolution
def conv3x3(in_channels, out_channels, stride=1):
return nn.Conv2d(in_channels, out_channels, kernel_size=3,
(continues on next page)
# Residual block
class ResidualBlock(nn.Module):
def __init__(self, in_channels, out_channels, stride=1, downsample=None):
super(ResidualBlock, self).__init__()
self.conv1 = conv3x3(in_channels, out_channels, stride)
self.bn1 = nn.BatchNorm2d(out_channels)
self.relu = nn.ReLU(inplace=True)
self.conv2 = conv3x3(out_channels, out_channels)
self.bn2 = nn.BatchNorm2d(out_channels)
self.downsample = downsample
# ResNet
class ResNet(nn.Module):
def __init__(self, block, layers, num_classes=10):
super(ResNet, self).__init__()
self.in_channels = 16
self.conv = conv3x3(3, 16)
self.bn = nn.BatchNorm2d(16)
self.relu = nn.ReLU(inplace=True)
self.layer1 = self.make_layer(block, 16, layers[0])
self.layer2 = self.make_layer(block, 32, layers[1], 2)
self.layer3 = self.make_layer(block, 64, layers[2], 2)
self.avg_pool = nn.AvgPool2d(8)
self.fc = nn.Linear(64, num_classes)
ResNet9
• DAWNBench on cifar10
• ResNet9: train to 94% CIFAR10 accuracy in 100 seconds
train_loader = torch.utils.data.DataLoader(
datasets.MNIST('data', train=True, download=True,
transform=transforms.Compose([
transforms.ToTensor(),
transforms.Normalize((0.1307,), (0.3081,))
])),
batch_size=batch_size_train, shuffle=True)
test_loader = torch.utils.data.DataLoader(
datasets.MNIST('data', train=False, transform=transforms.Compose([
transforms.ToTensor(),
transforms.Normalize((0.1307,), (0.3081,))
])),
batch_size=batch_size_test, shuffle=True)
return train_loader, test_loader
LeNet
Dry run in debug mode to get the shape of the last convnet layer.
model = LeNet5((1, 6, 16, 1, 120, 84, 10), debug=True)
batch_idx, (data_example, target_example) = next(enumerate(train_loader))
print(model)
_ = model(data_example)
LeNet5(
(conv1): Conv2d(1, 6, kernel_size=(5, 5), stride=(1, 1), padding=(2, 2))
(conv2): Conv2d(6, 16, kernel_size=(5, 5), stride=(1, 1))
(fc1): Linear(in_features=1, out_features=120, bias=True)
(fc2): Linear(in_features=120, out_features=84, bias=True)
(fc3): Linear(in_features=84, out_features=10, bias=True)
)
### DEBUG: Shape of last convnet= torch.Size([16, 5, 5]) . FC size= 400
torch.Size([6, 1, 5, 5])
torch.Size([6])
torch.Size([16, 6, 5, 5])
torch.Size([16])
torch.Size([120, 400])
torch.Size([120])
torch.Size([84, 120])
torch.Size([84])
torch.Size([10, 84])
(continues on next page)
Epoch 2/4
----------
train Loss: 0.0829 Acc: 97.46%
val Loss: 0.0561 Acc: 98.29%
Epoch 4/4
----------
train Loss: 0.0545 Acc: 98.32%
val Loss: 0.0547 Acc: 98.23%
MiniVGGNet
print(model)
_ = model(data_example)
MiniVGGNet(
(conv11): Conv2d(1, 16, kernel_size=(3, 3), stride=(1, 1))
(conv12): Conv2d(16, 16, kernel_size=(3, 3), stride=(1, 1))
(conv21): Conv2d(16, 32, kernel_size=(3, 3), stride=(1, 1))
(conv22): Conv2d(32, 32, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(continues on next page)
torch.Size([16, 1, 3, 3])
torch.Size([16])
torch.Size([16, 16, 3, 3])
torch.Size([16])
torch.Size([32, 16, 3, 3])
torch.Size([32])
torch.Size([32, 32, 3, 3])
torch.Size([32])
torch.Size([120, 800])
torch.Size([120])
torch.Size([84, 120])
torch.Size([84])
torch.Size([10, 84])
torch.Size([10])
Total number of parameters = 123502
Epoch 0/4
----------
train Loss: 1.4950 Acc: 45.62%
val Loss: 0.2489 Acc: 91.97%
Epoch 2/4
----------
train Loss: 0.0838 Acc: 97.30%
val Loss: 0.0653 Acc: 97.81%
Epoch 4/4
----------
train Loss: 0.0484 Acc: 98.51%
val Loss: 0.0443 Acc: 98.60%
Reduce the size of the training dataset by considering only 10 minibatche for size16.
train_loader, val_loader = load_mnist(16, 1000)
train_size = 10 * 16
# Stratified sub-sampling
targets = train_loader.dataset.targets.numpy()
nclasses = len(set(targets))
Train size= 160 Train label count= {0: 16, 1: 16, 2: 16, 3: 16, 4: 16, 5: 16, 6: 16, 7:␣
˓→16, 8: 16, 9: 16}
Batch sizes= [16, 16, 16, 16, 16, 16, 16, 16, 16, 16]
Datasets shape {'train': torch.Size([60000, 28, 28]), 'val': torch.Size([10000, 28, 28])}
N input features 784 N output 10
LeNet5
Epoch 0/99
----------
train Loss: 2.3009 Acc: 12.50%
val Loss: 2.2938 Acc: 10.90%
Epoch 20/99
----------
train Loss: 0.3881 Acc: 90.00%
val Loss: 0.6946 Acc: 78.96%
Epoch 40/99
----------
train Loss: 0.0278 Acc: 100.00%
val Loss: 0.5514 Acc: 86.55%
Epoch 60/99
----------
train Loss: 0.0068 Acc: 100.00%
val Loss: 0.6252 Acc: 86.77%
Epoch 80/99
----------
train Loss: 0.0033 Acc: 100.00%
val Loss: 0.6711 Acc: 86.83%
MiniVGGNet
Epoch 0/99
----------
train Loss: 2.3051 Acc: 10.00%
val Loss: 2.3061 Acc: 9.82%
Epoch 20/99
----------
train Loss: 2.2842 Acc: 10.00%
val Loss: 2.2860 Acc: 9.82%
Epoch 40/99
----------
train Loss: 0.6609 Acc: 75.63%
val Loss: 0.7725 Acc: 72.44%
Epoch 60/99
----------
train Loss: 0.2788 Acc: 93.75%
val Loss: 0.7235 Acc: 80.34%
Epoch 80/99
----------
train Loss: 0.0023 Acc: 100.00%
val Loss: 0.7130 Acc: 86.56%
(continues on next page)
import numpy as np
import torch
import torch.nn as nn
import torchvision
import torchvision.transforms as transforms
# Device configuration
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
# Hyper-parameters
num_epochs = 5
learning_rate = 0.001
# CIFAR-10 dataset
train_dataset = torchvision.datasets.CIFAR10(root='data/',
train=True,
transform=transform,
download=True)
val_dataset = torchvision.datasets.CIFAR10(root='data/',
train=False,
transform=transforms.ToTensor())
# Data loader
train_loader = torch.utils.data.DataLoader(dataset=train_dataset,
batch_size=100,
shuffle=True)
val_loader = torch.utils.data.DataLoader(dataset=val_dataset,
batch_size=100,
shuffle=False)
LeNet
LeNet5(
(conv1): Conv2d(3, 6, kernel_size=(5, 5), stride=(1, 1), padding=(2, 2))
(conv2): Conv2d(6, 16, kernel_size=(5, 5), stride=(1, 1))
(fc1): Linear(in_features=1, out_features=120, bias=True)
(fc2): Linear(in_features=120, out_features=84, bias=True)
(fc3): Linear(in_features=84, out_features=10, bias=True)
)
(continues on next page)
torch.Size([6, 3, 5, 5])
torch.Size([6])
torch.Size([16, 6, 5, 5])
torch.Size([16])
torch.Size([120, 576])
torch.Size([120])
torch.Size([84, 120])
torch.Size([84])
torch.Size([10, 84])
torch.Size([10])
Total number of parameters = 83126
Epoch 0/24
----------
train Loss: 2.3047 Acc: 10.00%
val Loss: 2.3043 Acc: 10.00%
Epoch 5/24
----------
train Loss: 2.3019 Acc: 10.01%
val Loss: 2.3015 Acc: 10.36%
Epoch 10/24
----------
train Loss: 2.2989 Acc: 12.97%
val Loss: 2.2979 Acc: 11.93%
Epoch 15/24
----------
train Loss: 2.2854 Acc: 10.34%
val Loss: 2.2808 Acc: 10.26%
Epoch 20/24
----------
train Loss: 2.1966 Acc: 15.87%
(continues on next page)
Epoch 0/24
----------
train Loss: 2.1439 Acc: 19.81%
val Loss: 1.9415 Acc: 29.59%
Epoch 5/24
----------
train Loss: 1.3457 Acc: 51.62%
val Loss: 1.2294 Acc: 55.74%
Epoch 10/24
----------
train Loss: 1.1607 Acc: 58.39%
val Loss: 1.1031 Acc: 60.68%
Epoch 15/24
----------
train Loss: 1.0710 Acc: 62.08%
(continues on next page)
Epoch 20/24
----------
train Loss: 1.0078 Acc: 64.25%
val Loss: 0.9505 Acc: 66.62%
Epoch 0/24
----------
train Loss: 1.8857 Acc: 29.70%
val Loss: 1.6223 Acc: 40.22%
Epoch 5/24
----------
train Loss: 1.3564 Acc: 50.88%
val Loss: 1.2271 Acc: 55.97%
Epoch 10/24
----------
train Loss: 1.2169 Acc: 56.35%
(continues on next page)
Epoch 15/24
----------
train Loss: 1.1296 Acc: 59.67%
val Loss: 1.0458 Acc: 63.05%
Epoch 20/24
----------
train Loss: 1.0830 Acc: 61.16%
val Loss: 1.0047 Acc: 64.49%
MiniVGGNet
MiniVGGNet(
(conv11): Conv2d(3, 16, kernel_size=(3, 3), stride=(1, 1))
(conv12): Conv2d(16, 16, kernel_size=(3, 3), stride=(1, 1))
(conv21): Conv2d(16, 32, kernel_size=(3, 3), stride=(1, 1))
(conv22): Conv2d(32, 32, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
(fc1): Linear(in_features=1, out_features=120, bias=True)
(fc2): Linear(in_features=120, out_features=84, bias=True)
(fc3): Linear(in_features=84, out_features=10, bias=True)
)
### DEBUG: Shape of last convnet= torch.Size([32, 6, 6]) . FC size= 1152
Epoch 0/24
----------
train Loss: 2.3027 Acc: 10.12%
val Loss: 2.3004 Acc: 10.31%
Epoch 5/24
----------
train Loss: 1.4790 Acc: 45.88%
val Loss: 1.3726 Acc: 50.32%
Epoch 10/24
----------
train Loss: 1.1115 Acc: 60.74%
val Loss: 1.0193 Acc: 64.00%
Epoch 15/24
----------
train Loss: 0.8937 Acc: 68.41%
val Loss: 0.8297 Acc: 71.18%
Epoch 20/24
----------
train Loss: 0.7848 Acc: 72.14%
val Loss: 0.7136 Acc: 75.42%
Adam
Epoch 0/24
----------
train Loss: 1.8366 Acc: 31.34%
val Loss: 1.5805 Acc: 41.62%
Epoch 5/24
----------
train Loss: 1.1755 Acc: 57.79%
val Loss: 1.1027 Acc: 60.83%
Epoch 10/24
----------
train Loss: 0.9741 Acc: 65.53%
val Loss: 0.8994 Acc: 68.29%
Epoch 15/24
----------
train Loss: 0.8611 Acc: 69.74%
val Loss: 0.8465 Acc: 70.90%
Epoch 20/24
----------
train Loss: 0.7916 Acc: 71.90%
val Loss: 0.7513 Acc: 74.03%
(continues on next page)
ResNet
Epoch 0/24
----------
train Loss: 1.4402 Acc: 46.84%
val Loss: 1.7289 Acc: 38.18%
Epoch 5/24
----------
train Loss: 0.6337 Acc: 77.88%
val Loss: 0.8672 Acc: 71.34%
Epoch 10/24
----------
train Loss: 0.4851 Acc: 83.11%
val Loss: 0.5754 Acc: 80.47%
Epoch 15/24
----------
(continues on next page)
Epoch 20/24
----------
train Loss: 0.3470 Acc: 87.99%
val Loss: 0.4696 Acc: 84.20%
Training complete in 7m 5s
Best val Acc: 85.60%
Sources:
• cs231n @ Stanford
• Sasank Chilamkurthy
Quote cs231n @ Stanford:
In practice, very few people train an entire Convolutional Network from scratch (with random
initialization), because it is relatively rare to have a dataset of sufficient size. Instead, it is
common to pretrain a ConvNet on a very large dataset (e.g. ImageNet, which contains 1.2
million images with 1000 categories), and then use the ConvNet either as an initialization or a
fixed feature extractor for the task of interest.
These two major transfer learning scenarios look as follows:
• ConvNet as fixed feature extractor:
– Take a ConvNet pretrained on ImageNet,
– Remove the last fully-connected layer (this layer’s outputs are the 1000 class scores
for a different task like ImageNet)
– Treat the rest of the ConvNet as a fixed feature extractor for the new dataset.
In practice:
– Freeze the weights for all of the network except that of the final fully connected layer.
This last fully connected layer is replaced with a new one with random weights and
only this layer is trained.
• Finetuning the convnet:
fine-tune the weights of the pretrained network by continuing the backpropagation. It is possi-
ble to fine-tune all the layers of the ConvNet
Instead of random initializaion, we initialize the network with a pretrained network, like the
one that is trained on imagenet 1000 dataset. Rest of the training looks as usual.
%matplotlib inline
import os
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from torch.optim import lr_scheduler
import torchvision
import torchvision.transforms as transforms
from torchvision import models
#
from pathlib import Path
import matplotlib.pyplot as plt
# Device configuration
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
since = time.time()
best_model_wts = copy.deepcopy(model.state_dict())
best_acc = 0.0
running_loss = 0.0
running_corrects = 0
# forward
# track history if only in train
with torch.set_grad_enabled(phase == 'train'):
outputs = model(inputs)
_, preds = torch.max(outputs, 1)
loss = criterion(outputs, labels)
# statistics
running_loss += loss.item() * inputs.size(0)
running_corrects += torch.sum(preds == labels.data)
#nsamples = dataloaders[phase].dataset.data.shape[0]
epoch_loss = running_loss / nsamples
(continues on next page)
losses[phase].append(epoch_loss)
accuracies[phase].append(epoch_acc)
if log_interval is not None and epoch % log_interval == 0:
print('{} Loss: {:.4f} Acc: {:.4f}'.format(
phase, epoch_loss, epoch_acc))
# CIFAR-10 dataset
train_dataset = torchvision.datasets.CIFAR10(root='data/',
train=True,
transform=transform,
download=True)
test_dataset = torchvision.datasets.CIFAR10(root='data/',
train=False,
transform=transforms.ToTensor())
# Data loader
(continues on next page)
val_loader = torch.utils.data.DataLoader(dataset=test_dataset,
batch_size=100,
shuffle=False)
model_ft = model_ft.to(device)
criterion = nn.CrossEntropyLoss()
epochs = np.arange(len(losses['train']))
_ = plt.plot(epochs, losses['train'], '-b', epochs, losses['val'], '--r')
Epoch 0/24
----------
train Loss: 1.2478 Acc: 0.5603
val Loss: 0.9084 Acc: 0.6866
(continues on next page)
Epoch 5/24
----------
train Loss: 0.5801 Acc: 0.7974
val Loss: 0.5918 Acc: 0.7951
Epoch 10/24
----------
train Loss: 0.4765 Acc: 0.8325
val Loss: 0.5257 Acc: 0.8178
Epoch 15/24
----------
train Loss: 0.4555 Acc: 0.8390
val Loss: 0.5205 Acc: 0.8201
Epoch 20/24
----------
train Loss: 0.4557 Acc: 0.8395
val Loss: 0.5183 Acc: 0.8212
Adam optimizer
model_ft = models.resnet18(pretrained=True)
num_ftrs = model_ft.fc.in_features
# Here the size of each output sample is set to 10.
model_ft.fc = nn.Linear(num_ftrs, D_out)
model_ft = model_ft.to(device)
criterion = nn.CrossEntropyLoss()
(continues on next page)
epochs = np.arange(len(losses['train']))
_ = plt.plot(epochs, losses['train'], '-b', epochs, losses['val'], '--r')
Epoch 0/24
----------
train Loss: 1.0491 Acc: 0.6407
val Loss: 0.8981 Acc: 0.6881
Epoch 5/24
----------
train Loss: 0.5495 Acc: 0.8135
val Loss: 0.9076 Acc: 0.7147
Epoch 10/24
----------
train Loss: 0.3352 Acc: 0.8834
val Loss: 0.4148 Acc: 0.8613
Epoch 15/24
----------
train Loss: 0.2819 Acc: 0.9017
val Loss: 0.4019 Acc: 0.8646
Epoch 20/24
----------
train Loss: 0.2719 Acc: 0.9050
val Loss: 0.4025 Acc: 0.8675
Freeze all the network except the final layer: requires_grad == False to freeze the parameters
so that the gradients are not computed in backward().
model_conv = torchvision.models.resnet18(pretrained=True)
for param in model_conv.parameters():
param.requires_grad = False
model_conv = model_conv.to(device)
criterion = nn.CrossEntropyLoss()
epochs = np.arange(len(losses['train']))
_ = plt.plot(epochs, losses['train'], '-b', epochs, losses['val'], '--r')
Epoch 0/24
----------
train Loss: 1.9107 Acc: 0.3265
val Loss: 1.7982 Acc: 0.3798
Epoch 10/24
----------
train Loss: 1.6411 Acc: 0.4278
val Loss: 1.6737 Acc: 0.4269
Epoch 15/24
----------
train Loss: 1.6315 Acc: 0.4299
val Loss: 1.6724 Acc: 0.4218
Epoch 20/24
----------
train Loss: 1.6400 Acc: 0.4274
val Loss: 1.6755 Acc: 0.4250
Adam optimizer
model_conv = torchvision.models.resnet18(pretrained=True)
for param in model_conv.parameters():
param.requires_grad = False
model_conv = model_conv.to(device)
(continues on next page)
criterion = nn.CrossEntropyLoss()
epochs = np.arange(len(losses['train']))
_ = plt.plot(epochs, losses['train'], '-b', epochs, losses['val'], '--r')
---------------------------------------------------------------------------
<ipython-input-16-dde92868b554> in <module>
19
20 model, losses, accuracies = train_val_model(model_conv, criterion, optimizer_conv,
---> 21 exp_lr_scheduler, dataloaders, num_epochs=25, log_interval=5)
22
23 epochs = np.arange(len(losses['train']))
SEVEN
• genindex
• modindex
• search
329