Learning Algorithms¶

Inspired by Sam Lau, who co-authored the Learning Data Science book.

In this lesson, we'll introduce machine learning from the ground up. By the end of this lesson, students will be able to:

  • Describe the difference between traditional algorithms and machine learning algorithms.
  • Identify the components of a machine learning model and dataset features and labels.
  • Apply sklearn to train a decision tree for classification and regression tasks.

A while back, we discussed data visualization using the Puget Sound Clean Air Agency's EPA-grade air quality sensors (AQS). However, these sensors are typically expensive, costing anywhere between \$15,000 and \$40,000 each, making it hard to deploy a large number of these sensors. Furthermore, EPA-grade AQS measurements also undergo calibration and accuracy checks that lead to delays of one or two hours, leading to data that is very accurate but not necessarily timely.

In contrast, "PurpleAir makes sensors that empower Community Scientists who collect hyper-local air quality data and share it with the public." In this lesson, we'll learn how we can use more accurate but less timely AQS measurements to calibrate the less accurate but more timely PurpleAir sensor (PAS) measurements so that we can provide the best information to the general public. The concepts in this lesson are actually used in the real world when you visit the EPA AirNow Fire and Smoke Map: the PAS data in this map are calibrated using the approach we will learn today.

In [1]:
import pandas as pd
import seaborn as sns

sns.set_theme()

Our dataset includes over 12,000 matched observations where we've paired each AQS measurement with a nearby PAS measurement, along with 3 other variables that experts have identified as potentially impacting PAS measurement quality. The dataset includes 5 columns:

  • The very accurate EPA-grade air quality sensor (AQS) measurement of the PM2.5 value.
  • The temperature in degrees celsius.
  • The relative humidity as a percentage between 0% and 100%.
  • The dew point, where a higher dew point means more moisture in the air.
  • The less accurate but more timely PurpleAir sensor (PAS) measurement of the PM2.5 value.

How can we use the PAS measurement together with the temperature, humidity, and dew point to predict the AQS measurement?

In [2]:
sensor_data = pd.read_csv("sensor_data.csv")
sensor_data
Out[2]:
AQS temp humidity dew PAS
0 6.7 18.027263 38.564815 3.629662 8.616954
1 3.8 16.115280 49.404315 5.442318 3.493916
2 4.0 19.897634 29.972222 1.734051 3.799601
3 4.7 21.378334 32.474513 4.165624 4.369691
4 3.2 18.443822 43.898226 5.867611 3.191071
... ... ... ... ... ...
12092 5.5 -12.101337 54.188889 -19.555834 2.386120
12093 16.8 4.159967 56.256030 -3.870659 32.444987
12094 15.6 1.707895 65.779221 -4.083768 25.297018
12095 14.0 -14.380144 48.206481 -23.015378 8.213208
12096 5.8 5.081813 52.200000 -4.016401 9.436011

12097 rows × 5 columns

Let's use data visualization to understand this dataset by scatter-plotting PAS to AQS. (Why do we place PAS on the x-axis and AQS on the y-axis?)

In [3]:
sns.relplot(sensor_data, x="PAS", y="AQS")
Out[3]:
<seaborn.axisgrid.FacetGrid at 0x7c584597a080>
No description has been provided for this image

Guessing game¶

Let's start with a simpler task: How can we use the PAS measurement alone to predict the AQS measurement? To do this, let's choose a line that best describes the trend, or the regression line that is produced by calling lmplot.

In [5]:
sns.lmplot(sensor_data, x="PAS", y="AQS")
Out[5]:
<seaborn.axisgrid.FacetGrid at 0x7c5902633070>
No description has been provided for this image
In [9]:
def plot_line(slope, intercept=0):
    grid = sns.relplot(sensor_data, x="PAS", y="AQS")
    grid.facet_axis(0, 0).plot([0, 260], [intercept, slope * 260 + intercept], c="orange")
    grid.set(title=f"Slope = {slope:.2f}, intercept = {intercept:.2f}")
    return grid


plot_line(0.45)
Out[9]:
<seaborn.axisgrid.FacetGrid at 0x7c583d5077f0>
No description has been provided for this image
In [10]:
sensor_data["PAS"]
Out[10]:
0         8.616954
1         3.493916
2         3.799601
3         4.369691
4         3.191071
           ...    
12092     2.386120
12093    32.444987
12094    25.297018
12095     8.213208
12096     9.436011
Name: PAS, Length: 12097, dtype: float64
In [11]:
sensor_data["AQS"]
Out[11]:
0         6.7
1         3.8
2         4.0
3         4.7
4         3.2
         ... 
12092     5.5
12093    16.8
12094    15.6
12095    14.0
12096     5.8
Name: AQS, Length: 12097, dtype: float64
In [12]:
sensor_data["PAS"] * 0.45
Out[12]:
0         3.877629
1         1.572262
2         1.709820
3         1.966361
4         1.435982
           ...    
12092     1.073754
12093    14.600244
12094    11.383658
12095     3.695944
12096     4.246205
Name: PAS, Length: 12097, dtype: float64
In [16]:
# "Error" represents the difference between our guess and the actual label
error_point45 = sensor_data["PAS"] * 0.45 - sensor_data["AQS"]
error_point45
Out[16]:
0        -2.822371
1        -2.227738
2        -2.290180
3        -2.733639
4        -1.764018
           ...    
12092    -4.426246
12093    -2.199756
12094    -4.216342
12095   -10.304056
12096    -1.553795
Length: 12097, dtype: float64
In [18]:
# Metric 1: "Mean Error" --- it can cancel out positives with negatives
error_point45.mean()
Out[18]:
-2.559337341646414
In [19]:
# Metric 2: "Mean Absolute Error"
error_point45.abs().mean()
Out[19]:
2.8898755399452463
In [20]:
# Metric 3: "Mean Squared Error" --- tend to penalize big errors moreso than small errors
(error_point45 ** 2).mean()
# Supervised machine learning: using a metric to tell the quality of your learning
Out[20]:
13.182456516471547
In [21]:
(error_slope2 ** 2).mean()
Out[21]:
693.5889006554958

What differentiates machine learning from just human guessing is the use of an algorithm to find the best line, which requires a metric for the quality of a trend.

We can visualize all of our guesses so far by plotting them against their mean squared errors on what's called a loss surface.

In [23]:
def plot_loss(slopes):
    from sklearn.metrics import mean_squared_error as mse
    losses = [mse(sensor_data["PAS"] * s, sensor_data["AQS"]) for s in slopes]
    grid = sns.relplot(x=slopes, y=losses)
    grid.set(title="Loss surface", xlabel="Slope", ylabel="MSE", xlim=[0, 1], ylim=[0, None])
    return grid


plot_loss([1, 0.5, 0.25, 0.45, 0.6, 0.75])
Out[23]:
<seaborn.axisgrid.FacetGrid at 0x7c58538c4940>
No description has been provided for this image

Gradient descent¶

So how do we write a machine learning algorithm that can optimize this metric and find the minimum mean squared error in the loss surface so that it selects the best possible line? Machine learning scientists can apply concepts from linear algebra to solve this system by selecting a random initial theta (slope) value and then rolling down the hill toward the minimum value at the bottom of the bowl. We can express this using numpy, a numeric computation module for Python that is a building block for pandas and sklearn (as we'll see later).

$$ \nabla_{\!\theta}\; \text{MSE}(\boldsymbol{\theta}, \mathbf{X}, \mathbf{y}) = -\frac{2}{n}(\mathbf{X}^\top \mathbf{y} - \mathbf{X}^\top \mathbf{X} \boldsymbol{\theta}) $$

In [29]:
import numpy as np


def grad_mse(theta, X, y):
    return np.array(- 2 / len(X) * (X.T @ y - X.T @ X * theta))


thetas = [np.random.random()]
print("Random initial theta value:", thetas[-1])

Random initial theta value: 0.964483396122684

We can then take a small step in the opposite direction of the gradient to roll down the hill until we converge on a good guess. To make this a machine learning algorithm, we simply put the update step in a loop until the theta values no longer noticeably change.

In [35]:
plot_line(thetas[-1])
plot_loss(thetas)

# Take a small step in the opposite direction of the gradient to roll downhill
thetas.append(thetas[-1] - 0.002 * grad_mse(thetas[-1], sensor_data["PAS"], sensor_data["AQS"]))
No description has been provided for this image
No description has been provided for this image

Linear regression models¶

What we've just described is the gradient descent algorithm for fitting a linear regression model. A linear regression model is a machine learning model that is used to predict a numeric value (like AQS measurements) using a linear combination of coefficients and features (columns from the training dataset). scikit-learn provides an easy way to do define a linear regression model, fit our training dataset X to the target values y, and examine the coefficients to look inside the model.

In [37]:
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
X = sensor_data[["PAS"]]
y = sensor_data["AQS"]
reg = LinearRegression().fit(X, y)

print("Model:", " + ".join([f"{reg.intercept_:.2f}"] + [f"{coef:.2f}({X.columns[i]})" for i, coef in enumerate(reg.coef_)]))
print("Error:", mean_squared_error(y, reg.predict(X)))
plot_line(reg.coef_[0], reg.intercept_)

Model: 3.18 + 0.40(PAS)
Error: 6.343125879204806
Out[37]:
<seaborn.axisgrid.FacetGrid at 0x7c583d640850>
No description has been provided for this image

This procedure is more or less how lmplot works!

In [38]:
sns.lmplot(sensor_data, x="PAS", y="AQS")
Out[38]:
<seaborn.axisgrid.FacetGrid at 0x7c583d6422f0>
No description has been provided for this image

But the advantage of designing our own model is that we can combine other variables to reduce the mean squared error loss. The final model that the EPA uses only takes into account the sensor measurement and the relative humidity, but not any other variables. Later, we'll learn why they made this decision.

In [39]:
X = sensor_data[["PAS", "humidity"]]
y = sensor_data["AQS"]
reg = LinearRegression().fit(X, y)

print("Model:", " + ".join([f"{reg.intercept_:.2f}"] + [f"{coef:.2f}({X.columns[i]})" for i, coef in enumerate(reg.coef_)]))
print("Error:", mean_squared_error(y, reg.predict(X)))
plot_line(reg.coef_[0], reg.intercept_)

Model: 6.25 + 0.43(PAS) + -0.07(humidity)
Error: 5.144725853340175
Out[39]:
<seaborn.axisgrid.FacetGrid at 0x7c583d6dc4f0>
No description has been provided for this image

Classification versus regression¶

Everything we've seen so far fall under the category of regression, where we aim to predict a numeric target value (one column) from one or more features (one or more other columns). The other main category of problems is classification, which is just like regression except we aim to predict a categorical target value. For example, we might want to answer the question: How can we predict whether a house belongs in NY or SF from its beds, baths, price, year of construction, square footage, price per square foot, and elevation?

In [40]:
homes = pd.read_csv("homes.csv")
homes
Out[40]:
beds bath price year_built sqft price_per_sqft elevation city
0 2.0 1.0 999000 1960 1000 999 10 NY
1 2.0 2.0 2750000 2006 1418 1939 0 NY
2 2.0 2.0 1350000 1900 2150 628 9 NY
3 1.0 1.0 629000 1903 500 1258 9 NY
4 0.0 1.0 439000 1930 500 878 10 NY
... ... ... ... ... ... ... ... ...
487 5.0 2.5 1800000 1890 3073 586 76 SF
488 2.0 1.0 695000 1923 1045 665 106 SF
489 3.0 2.0 1650000 1922 1483 1113 106 SF
490 1.0 1.0 649000 1983 850 764 163 SF
491 3.0 2.0 995000 1956 1305 762 216 SF

492 rows × 8 columns

In [47]:
# How does it work?
for bed_count in homes["beds"].unique():
    # try dividing the dataset by that bed count
Out[47]:
array([ 2. ,  1. ,  0. ,  3. ,  4. ,  5. , 10. ,  7. ,  8. ,  0.5,  6. ])

Let's learn about decision trees, a machine learning algorithm that can be used for classification (and also, as it turns out, regression too). Decision trees learn a nested if-then-else logical hierarchy to fit a training dataset. In the following visualization, the color and opacity of each box represents whether that subset of homes is more likely to be in NY or SF. The notes inside each node indicate information about the values:

  1. The first line shows the condition. If the condition is true, go left; if not, go right.
  2. The second line shows the percentage of samples represented by that node.
  3. The third line shows the proportion of homes within that sample that belong in ["NY", "SF"].
  4. The fourth line shows the majority class in that category, corresponding to the bigger number on line 3.
In [45]:
from sklearn.tree import DecisionTreeClassifier, plot_tree
X = homes.drop("city", axis=1)
y = homes["city"]
clf = DecisionTreeClassifier(max_depth=2).fit(X, y)

import matplotlib.pyplot as plt
plt.figure(dpi=300)
plot_tree(
    clf,
    feature_names=X.columns,
    class_names=["NY", "SF"],
    label="root",
    filled=True,
    impurity=False,
    proportion=True,
    rounded=False
) and None # Hide return value of plot_tree
No description has been provided for this image

We can also use this dataset for regression too. Write a one-line expression to train a DecisionTreeRegressor model to predict the price of a home in this dataset from all other variables.

In [55]:
from sklearn.tree import DecisionTreeRegressor

reg = DecisionTreeRegressor().fit(homes.drop(["price", "city"], axis=1), homes["price"])
reg
Out[55]:
DecisionTreeRegressor()
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
DecisionTreeRegressor()
In [56]:
plt.figure(dpi=300)
plot_tree(
    reg,
    feature_names=X.columns,
    label="root",
    filled=True,
    impurity=False,
    proportion=True,
    rounded=False,
    max_depth=3
) and None # Hide return value of plot_tree
No description has been provided for this image

Consider each of the following tasks and answer whether they would be best modeled as classification or regression.

Predict whether an email is spam or not spam.

Classification, since the target value is the category "spam" or "not spam".

Predict the number of views a video will receive based on subscriber count.

Regression, since the target value is a number.

Predict the next word to appear in a sentence.

Classification, since the target value is to choose from the dictionary of all possible next words.

If you're curious about how ChatGPT works, Jay Mody has a good introduction to GPT in 60 Lines of NumPy.

In [54]:
# How to incorporate the city value? It's a string so we have to convert it to columns with numbers
pd.concat([homes.drop(["price", "city"], axis=1), pd.get_dummies(homes["city"])], axis=1)
Out[54]:
beds bath year_built sqft price_per_sqft elevation NY SF
0 2.0 1.0 1960 1000 999 10 True False
1 2.0 2.0 2006 1418 1939 0 True False
2 2.0 2.0 1900 2150 628 9 True False
3 1.0 1.0 1903 500 1258 9 True False
4 0.0 1.0 1930 500 878 10 True False
... ... ... ... ... ... ... ... ...
487 5.0 2.5 1890 3073 586 76 False True
488 2.0 1.0 1923 1045 665 106 False True
489 3.0 2.0 1922 1483 1113 106 False True
490 1.0 1.0 1983 850 764 163 False True
491 3.0 2.0 1956 1305 762 216 False True

492 rows × 8 columns