从零开始的数据科学第二版（四）
原文：
zh.annas-archive.org/md5/48ab308fc34189a6d7d26b91b72a6df9译者：飞龙
第十三章：贝叶斯朴素分类
为心灵保持天真，为思想保持成熟。
阿纳托尔·法朗士
如果人们不能进行社交，那么社交网络就没什么用了。因此，DataSciencester 拥有一个受欢迎的功能，允许会员发送消息给其他会员。虽然大多数会员是负责任的公民，只发送受欢迎的“最近好吗？”消息，但一些不法分子坚持不懈地向其他成员发送关于致富计划、无需处方的药物和盈利数据科学证书项目的垃圾邮件。您的用户已经开始抱怨，因此消息副总裁要求您使用数据科学找出一种过滤这些垃圾邮件的方法。
一个非常愚蠢的垃圾邮件过滤器
假设一个“宇宙”，由从所有可能消息中随机选择的消息组成。设S为事件“消息是垃圾邮件”，B为事件“消息包含词bitcoin”。贝叶斯定理告诉我们，包含词bitcoin的消息是垃圾邮件的条件概率是：
P ( S | B ) = [ P ( B | S ) P ( S ) ] / [ P ( B | S ) P ( S ) + P ( B | ¬ S ) P ( ¬ S ) ]
分子是消息是垃圾邮件且包含bitcoin的概率，而分母只是消息包含bitcoin的概率。因此，您可以将这个计算视为简单地表示为垃圾邮件的bitcoin消息的比例。
如果我们有大量的已知是垃圾邮件的消息和大量的已知不是垃圾邮件的消息，那么我们可以很容易地估计P(B|S)和P(B|¬S)。如果我们进一步假设任何消息都有相等的可能性是垃圾邮件或不是垃圾邮件（所以P(S) = P(¬S) = 0.5），那么：
P ( S | B ) = P ( B | S ) / [ P ( B | S ) + P ( B | ¬ S ) ]
例如，如果垃圾邮件中有 50%的消息包含bitcoin这个词，而非垃圾邮件中只有 1%的消息包含，那么包含bitcoin的任意邮件是垃圾邮件的概率是：
0 . 5 / ( 0 . 5 + 0 . 01 ) = 98 %
一个更复杂的垃圾邮件过滤器
现在想象我们有一个包含许多词w[1] ... w[n]的词汇表。为了将其转化为概率理论的领域，我们将X[i]写成事件“消息包含词w[i]”。还想象一下（通过某种未指定的过程），我们已经得出了一个估计值P(X[i]|S)，表示垃圾邮件消息包含第i个词的概率，以及类似的估计P(X[i]|¬S)，表示非垃圾邮件消息包含第i个词的概率。
贝叶斯方法的关键在于做出（大胆的）假设，即每个单词的存在（或不存在）在消息是垃圾邮件与否的条件下是独立的。直观地说，这个假设意味着知道某个垃圾邮件消息是否包含词bitcoin并不会告诉你同一消息是否包含词rolex。用数学术语来说，这意味着：
P ( X 1 = x 1 , . . . , X n = x n | S ) = P ( X 1 = x 1 | S ) × ⋯ × P ( X n = x n | S )
这是一个极端的假设。（这个技术名称中带有天真也是有原因的。）假设我们的词汇表仅仅包括比特币和劳力士这两个词，而一半的垃圾邮件是关于“赚取比特币”，另一半是关于“正品劳力士”。在这种情况下，朴素贝叶斯估计一封垃圾邮件包含比特币和劳力士的概率是：
P ( X 1 = 1 , X 2 = 1 | S ) = P ( X 1 = 1 | S ) P ( X 2 = 1 | S ) = . 5 × . 5 = . 25
因为我们假设了比特币和劳力士实际上从不会同时出现的知识。尽管这种假设的不现实性，这种模型通常表现良好，并且在实际的垃圾邮件过滤器中历史悠久地使用。
我们用于“仅比特币”垃圾邮件过滤器的相同贝叶斯定理推理告诉我们可以使用以下方程来计算一封邮件是垃圾邮件的概率：
P ( S | X = x ) = P ( X = x | S ) / [ P ( X = x | S ) + P ( X = x | ¬ S ) ]
朴素贝叶斯假设允许我们通过简单地将每个词汇单词的概率估计相乘来计算右侧的每个概率。
实际操作中，通常要避免将大量概率相乘，以防止下溢问题，即计算机无法处理太接近 0 的浮点数。从代数中记得log ( a b ) = log a + log b 和exp( log x ) = x ，我们通常将p 1 ⋯ p n 计算为等效的（但更友好于浮点数的）形式：
exp( log ( p 1 ) + ⋯ + log ( p n ) )
唯一剩下的挑战是估计P ( Xi | S ) 和P ( Xi | ¬ S ) ，即垃圾邮件（或非垃圾邮件）包含单词w i 的概率。如果我们有大量标记为垃圾邮件和非垃圾邮件的“训练”邮件，一个明显的第一次尝试是简单地估计P ( X i | S ) 为仅仅是包含单词w i 的垃圾邮件的比例。
不过，这会造成一个很大的问题。想象一下，在我们的训练集中，词汇表中的单词data只出现在非垃圾邮件中。然后我们会估计P ( data | S ) = 0。结果是，我们的朴素贝叶斯分类器将始终为包含单词data的任何消息分配垃圾邮件概率 0，即使是像“免费比特币和正品劳力士手表数据”这样的消息也是如此。为了避免这个问题，我们通常使用某种平滑技术。
特别是，我们将选择一个伪计数——k——并估计在垃圾邮件中看到第i个单词的概率为：
P ( X i | S ) = ( k + number of spams containing w i ) / ( 2 k + number of spams )
我们对P ( X i | ¬ S )也是类似的。也就是说，当计算第i个单词的垃圾邮件概率时，我们假设我们还看到了包含该单词的k个额外的非垃圾邮件和k个额外的不包含该单词的非垃圾邮件。
例如，如果data出现在 0/98 封垃圾邮件中，如果k为 1，则我们将P(data|S)估计为 1/100 = 0.01，这使得我们的分类器仍然可以为包含单词data的消息分配一些非零的垃圾邮件概率。
实施
现在我们有了构建分类器所需的所有组件。首先，让我们创建一个简单的函数，将消息标记为不同的单词。我们首先将每条消息转换为小写，然后使用re.findall提取由字母、数字和撇号组成的“单词”。最后，我们将使用set获取不同的单词：
from typing import Set
import re
def tokenize(text: str) -> Set[str]:
text = text.lower() # Convert to lowercase,
all_words = re.findall("[a-z0-9']+", text) # extract the words, and
return set(all_words) # remove duplicates.
assert tokenize("Data Science is science") == {"data", "science", "is"}
我们还将为我们的训练数据定义一个类型：
from typing import NamedTuple
class Message(NamedTuple):
text: str
is_spam: bool
由于我们的分类器需要跟踪来自训练数据的标记、计数和标签，我们将其构建为一个类。按照惯例，我们将非垃圾邮件称为ham邮件。
构造函数只需要一个参数，即计算概率时使用的伪计数。它还初始化一个空的标记集合、计数器来跟踪在垃圾邮件和非垃圾邮件中看到每个标记的频率，以及它被训练的垃圾邮件和非垃圾邮件的数量：
from typing import List, Tuple, Dict, Iterable
import math
from collections import defaultdict
class NaiveBayesClassifier:
def __init__(self, k: float = 0.5) -> None:
self.k = k # smoothing factor
self.tokens: Set[str] = set()
self.token_spam_counts: Dict[str, int] = defaultdict(int)
self.token_ham_counts: Dict[str, int] = defaultdict(int)
self.spam_messages = self.ham_messages = 0
接下来，我们将为其提供一个方法，让它训练一堆消息。首先，我们增加spam_messages和ham_messages计数。然后我们对每个消息文本进行标记化，对于每个标记，根据消息类型递增token_spam_counts或token_ham_counts：
def train(self, messages: Iterable[Message]) -> None:
for message in messages:
# Increment message counts
if message.is_spam:
self.spam_messages += 1
else:
self.ham_messages += 1
# Increment word counts
for token in tokenize(message.text):
self.tokens.add(token)
if message.is_spam:
self.token_spam_counts[token] += 1
else:
self.token_ham_counts[token] += 1
最终，我们将想要预测P(垃圾邮件 | 标记)。正如我们之前看到的，要应用贝叶斯定理，我们需要知道每个词汇表中的标记对于垃圾邮件和非垃圾邮件的P(标记 | 垃圾邮件)和P(标记 | 非垃圾邮件)。因此，我们将创建一个“私有”辅助函数来计算这些：
def _probabilities(self, token: str) -> Tuple[float, float]:
"""returns P(token | spam) and P(token | ham)"""
spam = self.token_spam_counts[token]
ham = self.token_ham_counts[token]
p_token_spam = (spam + self.k) / (self.spam_messages + 2 * self.k)
p_token_ham = (ham + self.k) / (self.ham_messages + 2 * self.k)
return p_token_spam, p_token_ham
最后，我们准备编写我们的predict方法。如前所述，我们不会将许多小概率相乘，而是将对数概率相加：
def predict(self, text: str) -> float:
text_tokens = tokenize(text)
log_prob_if_spam = log_prob_if_ham = 0.0
# Iterate through each word in our vocabulary
for token in self.tokens:
prob_if_spam, prob_if_ham = self._probabilities(token)
# If *token* appears in the message,
# add the log probability of seeing it
if token in text_tokens:
log_prob_if_spam += math.log(prob_if_spam)
log_prob_if_ham += math.log(prob_if_ham)
# Otherwise add the log probability of _not_ seeing it,
# which is log(1 - probability of seeing it)
else:
log_prob_if_spam += math.log(1.0 - prob_if_spam)
log_prob_if_ham += math.log(1.0 - prob_if_ham)
prob_if_spam = math.exp(log_prob_if_spam)
prob_if_ham = math.exp(log_prob_if_ham)
return prob_if_spam / (prob_if_spam + prob_if_ham)
现在我们有了一个分类器。
测试我们的模型
让我们通过为其编写一些单元测试来确保我们的模型有效。
messages = [Message("spam rules", is_spam=True),
Message("ham rules", is_spam=False),
Message("hello ham", is_spam=False)]
model = NaiveBayesClassifier(k=0.5)
model.train(messages)
首先，让我们检查它是否正确计数：
assert model.tokens == {"spam", "ham", "rules", "hello"}
assert model.spam_messages == 1
assert model.ham_messages == 2
assert model.token_spam_counts == {"spam": 1, "rules": 1}
assert model.token_ham_counts == {"ham": 2, "rules": 1, "hello": 1}
现在让我们做出预测。我们还将（费力地）手动执行我们的朴素贝叶斯逻辑，并确保获得相同的结果：
text = "hello spam"
probs_if_spam = [
(1 + 0.5) / (1 + 2 * 0.5), # "spam" (present)
1 - (0 + 0.5) / (1 + 2 * 0.5), # "ham" (not present)
1 - (1 + 0.5) / (1 + 2 * 0.5), # "rules" (not present)
(0 + 0.5) / (1 + 2 * 0.5) # "hello" (present)
]
probs_if_ham = [
(0 + 0.5) / (2 + 2 * 0.5), # "spam" (present)
1 - (2 + 0.5) / (2 + 2 * 0.5), # "ham" (not present)
1 - (1 + 0.5) / (2 + 2 * 0.5), # "rules" (not present)
(1 + 0.5) / (2 + 2 * 0.5), # "hello" (present)
]
p_if_spam = math.exp(sum(math.log(p) for p in probs_if_spam))
p_if_ham = math.exp(sum(math.log(p) for p in probs_if_ham))
# Should be about 0.83
assert model.predict(text) == p_if_spam / (p_if_spam + p_if_ham)
此测试通过，因此似乎我们的模型正在做我们认为的事情。如果您查看实际的概率，两个主要驱动因素是我们的消息包含spam（我们唯一的训练垃圾邮件消息包含）和它不包含ham（我们的两个训练非垃圾邮件消息均不包含）。
现在让我们尝试在一些真实数据上运行它。
使用我们的模型
一个流行的（尽管有些陈旧的）数据集是SpamAssassin 公共语料库。我们将查看以20021010为前缀的文件。
这里是一个脚本，将下载并解压它们到您选择的目录（或者您可以手动执行）：
from io import BytesIO # So we can treat bytes as a file.
import requests # To download the files, which
import tarfile # are in .tar.bz format.
BASE_URL = "https://spamassassin.apache.org/old/publiccorpus"
FILES = ["20021010_easy_ham.tar.bz2",
"20021010_hard_ham.tar.bz2",
"20021010_spam.tar.bz2"]
# This is where the data will end up,
# in /spam, /easy_ham, and /hard_ham subdirectories.
# Change this to where you want the data.
OUTPUT_DIR = 'spam_data'
for filename in FILES:
# Use requests to get the file contents at each URL.
content = requests.get(f"{BASE_URL}/{filename}").content
# Wrap the in-memory bytes so we can use them as a "file."
fin = BytesIO(content)
# And extract all the files to the specified output dir.
with tarfile.open(fileobj=fin, mode='r:bz2') as tf:
tf.extractall(OUTPUT_DIR)
文件的位置可能会更改（这在本书的第一版和第二版之间发生过），如果是这样，请相应调整脚本。
下载数据后，您应该有三个文件夹：spam，easy_ham和hard_ham。每个文件夹包含许多电子邮件，每封邮件都包含在单个文件中。为了保持非常简单，我们将仅查看每封电子邮件的主题行。
如何识别主题行？当我们浏览文件时，它们似乎都以“主题：”开头。因此，我们将寻找这个：
import glob, re
# modify the path to wherever you've put the files
path = 'spam_data/*/*'
data: List[Message] = []
# glob.glob returns every filename that matches the wildcarded path
for filename in glob.glob(path):
is_spam = "ham" not in filename
# There are some garbage characters in the emails; the errors='ignore'
# skips them instead of raising an exception.
with open(filename, errors='ignore') as email_file:
for line in email_file:
if line.startswith("Subject:"):
subject = line.lstrip("Subject: ")
data.append(Message(subject, is_spam))
break # done with this file
现在我们可以将数据分为训练数据和测试数据，然后我们就可以构建分类器了：
import random
from scratch.machine_learning import split_data
random.seed(0) # just so you get the same answers as me
train_messages, test_messages = split_data(data, 0.75)
model = NaiveBayesClassifier()
model.train(train_messages)
让我们生成一些预测并检查我们的模型的表现：
from collections import Counter
predictions = [(message, model.predict(message.text))
for message in test_messages]
# Assume that spam_probability > 0.5 corresponds to spam prediction
# and count the combinations of (actual is_spam, predicted is_spam)
confusion_matrix = Counter((message.is_spam, spam_probability > 0.5)
for message, spam_probability in predictions)
print(confusion_matrix)
这给出了 84 个真正的阳性（被分类为“垃圾邮件”的垃圾邮件），25 个假阳性（被分类为“垃圾邮件”的非垃圾邮件），703 个真负（被分类为“非垃圾邮件”的非垃圾邮件）和 44 个假阴性（被分类为“非垃圾邮件”的垃圾邮件）。这意味着我们的精度是 84 /（84 + 25）= 77％，我们的召回率是 84 /（84 + 44）= 65％，对于如此简单的模型来说，这些数字并不差。（假设如果我们查看的不仅仅是主题行，我们可能会做得更好。）
我们还可以检查模型的内部，看看哪些单词最少和最具有指示性的垃圾邮件。
def p_spam_given_token(token: str, model: NaiveBayesClassifier) -> float:
# We probably shouldn't call private methods, but it's for a good cause.
prob_if_spam, prob_if_ham = model._probabilities(token)
return prob_if_spam / (prob_if_spam + prob_if_ham)
words = sorted(model.tokens, key=lambda t: p_spam_given_token(t, model))
print("spammiest_words", words[-10:])
print("hammiest_words", words[:10])
最垃圾的词包括像sale，mortgage，money和rates这样的词，而最不垃圾的词包括像spambayes，users，apt和perl这样的词。因此，这也让我们对我们的模型基本做出正确的判断有了一些直观的信心。
我们如何才能获得更好的性能？一个明显的方法是获得更多的训练数据。还有很多改进模型的方法。以下是您可以尝试的一些可能性：
-
查看消息内容，而不仅仅是主题行。您将需要小心处理消息头。
-
我们的分类器考虑了训练集中出现的每个单词，甚至是仅出现一次的单词。修改分类器以接受可选的
min_count阈值，并忽略不至少出现这么多次的标记。 -
分词器没有类似词汇的概念（例如 cheap 和 cheapest）。修改分类器以接受可选的
stemmer函数，将单词转换为 等价类 的单词。例如，一个非常简单的词干提取函数可以是：def drop_final_s(word): return re.sub("s$", "", word)创建一个良好的词干提取函数很难。人们经常使用 Porter stemmer。
-
尽管我们的特征都是形如“消息包含词 w i”，但这并非必须如此。在我们的实现中，我们可以通过创建虚假标记，如 contains:number，并在适当时修改
tokenizer以发出它们来添加额外的特征，比如“消息包含数字”。
进一步探索
-
保罗·格雷厄姆的文章 “一个垃圾邮件过滤计划” 和 “更好的贝叶斯过滤” 非常有趣，并深入探讨了构建垃圾邮件过滤器背后的思想。
-
scikit-learn 包含一个
BernoulliNB模型，实现了我们这里实现的相同朴素贝叶斯算法，以及模型的其他变体。
第十四章：简单线性回归
艺术，如道德，就在于在某处划出界限。
G. K. Chesterton
在第五章中，我们使用correlation函数来衡量两个变量之间线性关系的强度。对于大多数应用程序，知道存在这样一个线性关系是不够的。我们需要理解关系的本质。这就是我们将使用简单线性回归的地方。
模型
回想一下，我们正在研究 DataSciencester 用户的朋友数量和每天在网站上花费的时间之间的关系。假设您已经确信，拥有更多朋友导致人们在网站上花费更多时间，而不是我们讨论过的其他解释之一。
参与用户参与部长要求您建立描述这种关系的模型。由于您找到了一个相当强的线性关系，线性模型是一个自然的起点。
特别是，您假设存在常数α（alpha）和β（beta），使得：
y i = β x i + α + ε i
其中y i是用户i每天在网站上花费的分钟数，x i是用户i的朋友数，ε是一个（希望很小的）误差项，表示这个简单模型未考虑到的其他因素。
假设我们已经确定了这样的alpha和beta，那么我们可以简单地进行预测：
def predict(alpha: float, beta: float, x_i: float) -> float:
return beta * x_i + alpha
如何选择alpha和beta？嗯，任何alpha和beta的选择都会给我们每个输入x_i预测输出。由于我们知道实际输出y_i，我们可以计算每对的误差：
def error(alpha: float, beta: float, x_i: float, y_i: float) -> float:
"""
The error from predicting beta * x_i + alpha
when the actual value is y_i
"""
return predict(alpha, beta, x_i) - y_i
我们真正想知道的是整个数据集的总误差。但我们不只是想把误差加起来——如果x_1的预测值过高，而x_2的预测值过低，误差可能会抵消掉。
因此，我们将平方误差加起来：
from scratch.linear_algebra import Vector
def sum_of_sqerrors(alpha: float, beta: float, x: Vector, y: Vector) -> float:
return sum(error(alpha, beta, x_i, y_i) ** 2
for x_i, y_i in zip(x, y))
最小二乘解是选择使sum_of_sqerrors尽可能小的alpha和beta。
使用微积分（或繁琐的代数），最小化误差的alpha和beta由以下公式给出：
from typing import Tuple
from scratch.linear_algebra import Vector
from scratch.statistics import correlation, standard_deviation, mean
def least_squares_fit(x: Vector, y: Vector) -> Tuple[float, float]:
"""
Given two vectors x and y,
find the least-squares values of alpha and beta
"""
beta = correlation(x, y) * standard_deviation(y) / standard_deviation(x)
alpha = mean(y) - beta * mean(x)
return alpha, beta
不需要详细进行数学推导，让我们思考为什么这可能是一个合理的解决方案。选择alpha简单地表示，当我们看到自变量x的平均值时，我们预测因变量y的平均值。
选择beta的意义在于，当输入值增加了standard_deviation(x)时，预测值就会增加correlation(x, y) * standard_deviation(y)。如果x和y完全正相关，x增加一个标准差会导致预测值增加一个y的标准差。当它们完全负相关时，x的增加会导致预测值的减少。当相关性为 0 时，beta为 0，这意味着x的变化对预测没有任何影响。
通常情况下，我们来快速测试一下：
x = [i for i in range(-100, 110, 10)]
y = [3 * i - 5 for i in x]
# Should find that y = 3x - 5
assert least_squares_fit(x, y) == (-5, 3)
现在很容易将其应用于第五章中去除异常值的数据：
from scratch.statistics import num_friends_good, daily_minutes_good
alpha, beta = least_squares_fit(num_friends_good, daily_minutes_good)
assert 22.9 < alpha < 23.0
assert 0.9 < beta < 0.905
这给出了alpha = 22.95 和beta = 0.903 的值。因此，我们的模型表明，我们预计一个没有朋友的用户每天在 DataSciencester 上花费大约 23 分钟。对于每个额外的朋友，我们预计用户每天在网站上多花大约一分钟。
在图 14-1 中，我们绘制预测线，以了解模型拟合观察数据的程度。
图 14-1. 我们的简单线性模型
当然，我们需要一种比盯着图表更好的方法来确定我们对数据的拟合程度。一个常见的度量是决定系数（或R 平方），它衡量因变量的总变异中模型所捕获的比例：
from scratch.statistics import de_mean
def total_sum_of_squares(y: Vector) -> float:
"""the total squared variation of y_i's from their mean"""
return sum(v ** 2 for v in de_mean(y))
def r_squared(alpha: float, beta: float, x: Vector, y: Vector) -> float:
"""
the fraction of variation in y captured by the model, which equals
1 - the fraction of variation in y not captured by the model
"""
return 1.0 - (sum_of_sqerrors(alpha, beta, x, y) /
total_sum_of_squares(y))
rsq = r_squared(alpha, beta, num_friends_good, daily_minutes_good)
assert 0.328 < rsq < 0.330
记住，我们选择了使预测误差平方和最小化的alpha和beta。我们可以选择一个线性模型“始终预测mean(y)”（对应于alpha = mean(y)和beta = 0），其预测误差平方和恰好等于总平方和。这意味着 R 平方为 0，表明该模型（在这种情况下显然）的表现不比简单预测均值好。
显然，最小二乘模型至少要与此一样好，这意味着预测误差平方和最多等于总平方和，这意味着 R 平方至少为 0。而预测误差平方和至少为 0，这意味着 R 平方最多为 1。
数字越高，我们的模型拟合数据越好。这里我们计算了一个 R 平方为 0.329，表明我们的模型在拟合数据方面只算是可以接受，显然还有其他因素在起作用。
使用梯度下降法
如果我们写theta = [alpha, beta]，我们也可以用梯度下降法解决这个问题：
import random
import tqdm
from scratch.gradient_descent import gradient_step
num_epochs = 10000
random.seed(0)
guess = [random.random(), random.random()] # choose random value to start
learning_rate = 0.00001
with tqdm.trange(num_epochs) as t:
for _ in t:
alpha, beta = guess
# Partial derivative of loss with respect to alpha
grad_a = sum(2 * error(alpha, beta, x_i, y_i)
for x_i, y_i in zip(num_friends_good,
daily_minutes_good))
# Partial derivative of loss with respect to beta
grad_b = sum(2 * error(alpha, beta, x_i, y_i) * x_i
for x_i, y_i in zip(num_friends_good,
daily_minutes_good))
# Compute loss to stick in the tqdm description
loss = sum_of_sqerrors(alpha, beta,
num_friends_good, daily_minutes_good)
t.set_description(f"loss: {loss:.3f}")
# Finally, update the guess
guess = gradient_step(guess, [grad_a, grad_b], -learning_rate)
# We should get pretty much the same results:
alpha, beta = guess
assert 22.9 < alpha < 23.0
assert 0.9 < beta < 0.905
如果您运行此操作，您将得到与我们使用精确公式相同的alpha和beta值。
最大似然估计
为什么选择最小二乘法？其中一个理由涉及最大似然估计。想象一下，我们有来自依赖于某些未知参数θ（theta）的分布的数据样本v 1 , ... , v n：
p ( v 1 , ... , v n | θ )
如果我们不知道θ，我们可以反过来将这个数量视为给定样本的θ的似然：
L ( θ | v 1 , ... , v n )
在这种方法下，最可能的θ是能够最大化这个似然函数的值——即使得观察数据最有可能出现的值。对于连续分布的情况，我们有一个概率分布函数而不是概率质量函数，我们也可以做同样的事情。
回到回归。关于简单回归模型经常做的一个假设是回归误差服从均值为 0、某个（已知）标准差σ的正态分布。如果是这种情况，那么基于观察到一对(x_i, y_i)的似然是：
L ( α , β | x i , y i , σ ) = 1 2πσ exp ( - (y i -α-βx i ) 2 / 2 σ 2 )
基于整个数据集的似然是每个个体似然的乘积，在alpha和beta被选择以最小化平方误差时最大。也就是说，在这种情况下（在这些假设下），最小化平方误差等价于最大化观察数据的似然。
进一步探索
继续阅读关于多元回归的内容在第十五章！
第十五章：多元回归
我不会看着问题并在里面加入不影响它的变量。
比尔·帕塞尔
虽然副总统对你的预测模型印象深刻，但她认为你可以做得更好。因此，你收集了额外的数据：你知道每个用户每天工作的小时数，以及他们是否拥有博士学位。你希望利用这些额外数据来改进你的模型。
因此，你假设一个包含更多独立变量的线性模型：
minutes = α + β 1 friends + β 2 work hours + β 3 phd + ε
显然，用户是否拥有博士学位不是一个数字——但是，正如我们在第十一章中提到的，我们可以引入一个虚拟变量，对于拥有博士学位的用户设为 1，没有的设为 0，之后它与其他变量一样是数值化的。
模型
回想一下，在第十四章中，我们拟合了一个形式为：
y i = α + β x i + ε i
现在想象每个输入x i不是单个数字，而是一个包含k个数字的向量，x i1 , ... , x ik 。多元回归模型假设：
y i = α + β 1 x i1 + . . . + β k x ik + ε i
在多元回归中，参数向量通常称为β。我们希望这个向量包括常数项，可以通过在数据中添加一列 1 来实现：
beta = [alpha, beta_1, ..., beta_k]
以及：
x_i = [1, x_i1, ..., x_ik]
那么我们的模型就是：
from scratch.linear_algebra import dot, Vector
def predict(x: Vector, beta: Vector) -> float:
"""assumes that the first element of x is 1"""
return dot(x, beta)
在这种特殊情况下，我们的自变量x将是一个向量列表，每个向量如下所示：
[1, # constant term
49, # number of friends
4, # work hours per day
0] # doesn't have PhD
最小二乘模型的进一步假设
为了使这个模型（以及我们的解决方案）有意义，还需要一些进一步的假设。
第一个假设是x的列是线性独立的——没有办法将任何一个写成其他一些的加权和。如果这个假设失败，估计beta是不可能的。在一个极端情况下，想象我们在数据中有一个额外的字段num_acquaintances，对于每个用户都恰好等于num_friends。
然后，从任意beta开始，如果我们将num_friends系数增加任意量，并将相同量从num_acquaintances系数减去，模型的预测将保持不变。这意味着没有办法找到num_friends的系数。（通常这种假设的违反不会那么明显。）
第二个重要假设是x的列与误差ε不相关。如果这一点不成立，我们对beta的估计将会系统错误。
例如，在第十四章中，我们建立了一个模型，预测每增加一个朋友与额外 0.90 分钟的网站使用时间相关。
想象也是这种情况：
-
工作时间更长的人在网站上花费的时间较少。
-
拥有更多朋友的人 tend to work more hours.
换句话说，假设“实际”模型如下：
minutes = α + β 1 friends + β 2 work hours + ε
其中β 2是负数，而且工作时间和朋友数量是正相关的。在这种情况下，当我们最小化单变量模型的误差时：
minutes = α + β 1 friends + ε
我们会低估β 1 。
想象一下，如果我们使用单变量模型并使用“实际”值β 1 进行预测会发生什么。（也就是说，这个值是通过最小化我们称之为“实际”模型的误差得到的。）预测值会倾向于对工作时间较长的用户过大，并且对工作时间较少的用户也稍微偏大，因为β 2 < 0 而我们“忘记”将其包含在内。由于工作时间与朋友数量呈正相关，这意味着对于朋友较多的用户，预测值往往过大，而对于朋友较少的用户，则稍微过大。
这样做的结果是，我们可以通过减少对β 1的估计来减少（单变量模型中的）误差，这意味着误差最小化的β 1小于“实际”值。也就是说，在这种情况下，单变量最小二乘解法会倾向于低估β 1。而且，通常情况下，每当自变量与这些误差相关联时，我们的最小二乘解法都会给我们一个偏倚的β 1估计。
拟合模型
就像我们在简单线性模型中所做的那样，我们会选择beta来最小化平方误差的和。手动找到一个确切的解决方案并不容易，这意味着我们需要使用梯度下降法。同样，我们希望最小化平方误差的和。误差函数与我们在第十四章中使用的几乎完全相同，只是不再期望参数[alpha, beta]，而是会接受任意长度的向量：
from typing import List
def error(x: Vector, y: float, beta: Vector) -> float:
return predict(x, beta) - y
def squared_error(x: Vector, y: float, beta: Vector) -> float:
return error(x, y, beta) ** 2
x = [1, 2, 3]
y = 30
beta = [4, 4, 4] # so prediction = 4 + 8 + 12 = 24
assert error(x, y, beta) == -6
assert squared_error(x, y, beta) == 36
如果你懂得微积分，计算梯度就很容易：
def sqerror_gradient(x: Vector, y: float, beta: Vector) -> Vector:
err = error(x, y, beta)
return [2 * err * x_i for x_i in x]
assert sqerror_gradient(x, y, beta) == [-12, -24, -36]
否则，你需要相信我的话。
此时，我们准备使用梯度下降法找到最优的beta。让我们首先编写一个least_squares_fit函数，可以处理任何数据集：
import random
import tqdm
from scratch.linear_algebra import vector_mean
from scratch.gradient_descent import gradient_step
def least_squares_fit(xs: List[Vector],
ys: List[float],
learning_rate: float = 0.001,
num_steps: int = 1000,
batch_size: int = 1) -> Vector:
"""
Find the beta that minimizes the sum of squared errors
assuming the model y = dot(x, beta).
"""
# Start with a random guess
guess = [random.random() for _ in xs[0]]
for _ in tqdm.trange(num_steps, desc="least squares fit"):
for start in range(0, len(xs), batch_size):
batch_xs = xs[start:start+batch_size]
batch_ys = ys[start:start+batch_size]
gradient = vector_mean([sqerror_gradient(x, y, guess)
for x, y in zip(batch_xs, batch_ys)])
guess = gradient_step(guess, gradient, -learning_rate)
return guess
然后我们可以将其应用到我们的数据中：
from scratch.statistics import daily_minutes_good
from scratch.gradient_descent import gradient_step
random.seed(0)
# I used trial and error to choose num_iters and step_size.
# This will run for a while.
learning_rate = 0.001
beta = least_squares_fit(inputs, daily_minutes_good, learning_rate, 5000, 25)
assert 30.50 < beta[0] < 30.70 # constant
assert 0.96 < beta[1] < 1.00 # num friends
assert -1.89 < beta[2] < -1.85 # work hours per day
assert 0.91 < beta[3] < 0.94 # has PhD
在实践中，你不会使用梯度下降法来估计线性回归；你会使用超出本书范围的线性代数技术来得到精确的系数。如果你这样做，你会得到如下方程：
minutes = 30 . 58 + 0 . 972 friends - 1 . 87 work hours + 0 . 923 phd
这与我们找到的结果非常接近。
解释模型
模型中的系数代表每个因素的其他条件相等估计影响的总和。其他条件相等时，每增加一个朋友，每天会多花一分钟在网站上。其他条件相等时，用户工作日每增加一个小时，每天会少花约两分钟在网站上。其他条件相等时，拥有博士学位与每天在网站上多花一分钟相关联。
这并没有（直接）告诉我们关于变量之间互动的任何信息。有可能工作时间对于朋友多的人和朋友少的人有不同的影响。这个模型没有捕捉到这一点。处理这种情况的一种方法是引入一个新变量，即“朋友”和“工作时间”的乘积。这实际上允许“工作时间”系数随着朋友数量的增加而增加（或减少）。
或者可能是，你有更多的朋友，你在网站上花的时间就越多直到一个点，之后进一步的朋友导致你在网站上花费的时间减少。（也许有太多的朋友经验就太压倒性了？）我们可以尝试通过添加另一个变量，即朋友数量的平方，来捕捉这一点在我们的模型中。
一旦我们开始添加变量，就需要担心它们的系数是否“重要”。我们可以无限制地添加乘积、对数、平方和更高次方。
拟合优度
再次我们可以看一下 R 平方：
from scratch.simple_linear_regression import total_sum_of_squares
def multiple_r_squared(xs: List[Vector], ys: Vector, beta: Vector) -> float:
sum_of_squared_errors = sum(error(x, y, beta) ** 2
for x, y in zip(xs, ys))
return 1.0 - sum_of_squared_errors / total_sum_of_squares(ys)
现在已经增加到 0.68：
assert 0.67 < multiple_r_squared(inputs, daily_minutes_good, beta) < 0.68
请记住，然而，向回归中添加新变量必然会增加 R 平方。毕竟，简单回归模型只是多重回归模型的特殊情况，其中“工作时间”和“博士学位”的系数都等于 0。最佳的多重回归模型将至少有一个与该模型一样小的误差。
因此，在多重回归中，我们还需要看看系数的标准误差，这些标准误差衡量我们对每个β i的估计有多确定。整体回归可能非常适合我们的数据，但如果一些自变量相关（或不相关），它们的系数可能意义不大。
测量这些误差的典型方法始于另一个假设——误差ε i是独立的正态随机变量，均值为 0，有一些共享的（未知）标准偏差σ。在这种情况下，我们（或者更可能是我们的统计软件）可以使用一些线性代数来找出每个系数的标准误差。它越大，我们对该系数的模型越不确定。不幸的是，我们没有设置好可以从头开始执行这种线性代数的工具。
偏离：自助法
想象我们有一个由某个（对我们来说未知的）分布生成的包含n个数据点的样本：
data = get_sample(num_points=n)
在第五章中，我们编写了一个可以计算样本中位数的函数，我们可以将其用作对分布本身中位数的估计。
但我们对我们的估计有多自信呢？如果样本中所有数据点都非常接近 100，那么实际中位数似乎也接近 100。如果样本中大约一半的数据点接近 0，而另一半接近 200，则我们对中位数的估计就不太确定。
如果我们能够重复获得新样本，我们可以计算许多样本的中位数，并查看这些中位数的分布。通常我们做不到这一点。在这种情况下，我们可以通过从我们的数据中有放回地选择n个数据点来bootstrap新数据集。然后我们可以计算这些合成数据集的中位数：
from typing import TypeVar, Callable
X = TypeVar('X') # Generic type for data
Stat = TypeVar('Stat') # Generic type for "statistic"
def bootstrap_sample(data: List[X]) -> List[X]:
"""randomly samples len(data) elements with replacement"""
return [random.choice(data) for _ in data]
def bootstrap_statistic(data: List[X],
stats_fn: Callable[[List[X]], Stat],
num_samples: int) -> List[Stat]:
"""evaluates stats_fn on num_samples bootstrap samples from data"""
return [stats_fn(bootstrap_sample(data)) for _ in range(num_samples)]
例如，考虑以下两个数据集：
# 101 points all very close to 100
close_to_100 = [99.5 + random.random() for _ in range(101)]
# 101 points, 50 of them near 0, 50 of them near 200
far_from_100 = ([99.5 + random.random()] +
[random.random() for _ in range(50)] +
[200 + random.random() for _ in range(50)])
如果计算这两个数据集的median，两者都将非常接近 100。然而，如果你看一下：
from scratch.statistics import median, standard_deviation
medians_close = bootstrap_statistic(close_to_100, median, 100)
你将主要看到数字非常接近 100。但如果你看一下：
medians_far = bootstrap_statistic(far_from_100, median, 100)
你会看到很多接近 0 和很多接近 200 的数字。
第一组中位数的标准偏差接近 0，而第二组中位数的则接近 100：
assert standard_deviation(medians_close) < 1
assert standard_deviation(medians_far) > 90
（这种极端情况下，通过手动检查数据很容易找到答案，但通常情况下这是不成立的。）
回归系数的标准误差
我们可以采取同样的方法来估计回归系数的标准误差。我们重复从我们的数据中取出一个bootstrap_sample，并基于该样本估计beta。如果与一个独立变量（比如num_friends）对应的系数在样本中变化不大，那么我们可以相信我们的估计相对较为精确。如果系数在样本中变化很大，那么我们就不能对我们的估计感到有信心。
唯一的微妙之处在于，在抽样之前，我们需要zip我们的x数据和y数据，以确保独立变量和因变量的相应值一起被抽样。这意味着bootstrap_sample将返回一个成对的列表(x_i, y_i)，我们需要重新组装成一个x_sample和一个y_sample：
from typing import Tuple
import datetime
def estimate_sample_beta(pairs: List[Tuple[Vector, float]]):
x_sample = [x for x, _ in pairs]
y_sample = [y for _, y in pairs]
beta = least_squares_fit(x_sample, y_sample, learning_rate, 5000, 25)
print("bootstrap sample", beta)
return beta
random.seed(0) # so that you get the same results as me
# This will take a couple of minutes!
bootstrap_betas = bootstrap_statistic(list(zip(inputs, daily_minutes_good)),
estimate_sample_beta,
100)
在此之后，我们可以估计每个系数的标准偏差。
bootstrap_standard_errors = [
standard_deviation([beta[i] for beta in bootstrap_betas])
for i in range(4)]
print(bootstrap_standard_errors)
# [1.272, # constant term, actual error = 1.19
# 0.103, # num_friends, actual error = 0.080
# 0.155, # work_hours, actual error = 0.127
# 1.249] # phd, actual error = 0.998
（如果我们收集了超过 100 个样本并使用了超过 5,000 次迭代来估计每个beta，我们可能会得到更好的估计，但我们没有那么多时间。）
我们可以使用这些来测试假设，比如“β i是否等于 0？”在零假设β i = 0（以及我们对ε i分布的其他假设）下，统计量：
t j = β j ^ / σ j ^
其中，我们对β j的估计除以其标准误差的估计值，遵循“ n - k自由度”的学生 t 分布。
如果我们有一个students_t_cdf函数，我们可以为每个最小二乘系数计算p-值，以指示如果实际系数为 0，则观察到这样的值的可能性有多大。不幸的是，我们没有这样的函数。（尽管如果我们不是从头开始工作，我们会有这样的函数。）
然而，随着自由度的增大，t-分布越来越接近于标准正态分布。在像这样的情况下，其中n远大于k，我们可以使用normal_cdf而仍然感觉良好：
from scratch.probability import normal_cdf
def p_value(beta_hat_j: float, sigma_hat_j: float) -> float:
if beta_hat_j > 0:
# if the coefficient is positive, we need to compute twice the
# probability of seeing an even *larger* value
return 2 * (1 - normal_cdf(beta_hat_j / sigma_hat_j))
else:
# otherwise twice the probability of seeing a *smaller* value
return 2 * normal_cdf(beta_hat_j / sigma_hat_j)
assert p_value(30.58, 1.27) < 0.001 # constant term
assert p_value(0.972, 0.103) < 0.001 # num_friends
assert p_value(-1.865, 0.155) < 0.001 # work_hours
assert p_value(0.923, 1.249) > 0.4 # phd
（在不像这样的情况下，我们可能会使用知道如何计算t-分布以及如何计算确切标准误差的统计软件。）
虽然大多数系数的p-值非常小（表明它们确实是非零的），但“PhD”的系数与 0 的差异不“显著”，这使得“PhD”的系数很可能是随机的，而不是有意义的。
在更复杂的回归场景中，有时您可能希望对数据进行更复杂的假设检验，例如“至少一个β j非零”或“β 1等于β 2 和 β 3等于β 4。” 您可以使用F-检验来执行此操作，但遗憾的是，这超出了本书的范围。
正则化
在实践中，您经常希望将线性回归应用于具有大量变量的数据集。这会产生一些额外的复杂性。首先，您使用的变量越多，就越有可能将模型过度拟合到训练集。其次，非零系数越多，就越难以理解它们。如果目标是解释某种现象，那么具有三个因素的稀疏模型可能比稍好的具有数百个因素的模型更有用。
正则化是一种方法，其中我们将惩罚项添加到误差项中，随着beta的增大而增加。然后，我们最小化组合误差和惩罚。我们越重视惩罚项，就越能够阻止大的系数。
例如，在岭回归中，我们添加的惩罚与beta_i的平方和成比例（通常不惩罚beta_0，即常数项）：
# alpha is a *hyperparameter* controlling how harsh the penalty is.
# Sometimes it's called "lambda" but that already means something in Python.
def ridge_penalty(beta: Vector, alpha: float) -> float:
return alpha * dot(beta[1:], beta[1:])
def squared_error_ridge(x: Vector,
y: float,
beta: Vector,
alpha: float) -> float:
"""estimate error plus ridge penalty on beta"""
return error(x, y, beta) ** 2 + ridge_penalty(beta, alpha)
然后我们可以按照通常的方式将其插入梯度下降：
from scratch.linear_algebra import add
def ridge_penalty_gradient(beta: Vector, alpha: float) -> Vector:
"""gradient of just the ridge penalty"""
return [0.] + [2 * alpha * beta_j for beta_j in beta[1:]]
def sqerror_ridge_gradient(x: Vector,
y: float,
beta: Vector,
alpha: float) -> Vector:
"""
the gradient corresponding to the ith squared error term
including the ridge penalty
"""
return add(sqerror_gradient(x, y, beta),
ridge_penalty_gradient(beta, alpha))
然后我们只需修改least_squares_fit函数，以使用sqerror_ridge_gradient而不是sqerror_gradient。（我不会在这里重复代码。）
将alpha设置为 0 后，就没有惩罚了，我们获得了与以前相同的结果：
random.seed(0)
beta_0 = least_squares_fit_ridge(inputs, daily_minutes_good, 0.0, # alpha
learning_rate, 5000, 25)
# [30.51, 0.97, -1.85, 0.91]
assert 5 < dot(beta_0[1:], beta_0[1:]) < 6
assert 0.67 < multiple_r_squared(inputs, daily_minutes_good, beta_0) < 0.69
随着alpha的增加，拟合的好坏变得更差，但beta的大小变小：
beta_0_1 = least_squares_fit_ridge(inputs, daily_minutes_good, 0.1, # alpha
learning_rate, 5000, 25)
# [30.8, 0.95, -1.83, 0.54]
assert 4 < dot(beta_0_1[1:], beta_0_1[1:]) < 5
assert 0.67 < multiple_r_squared(inputs, daily_minutes_good, beta_0_1) < 0.69
beta_1 = least_squares_fit_ridge(inputs, daily_minutes_good, 1, # alpha
learning_rate, 5000, 25)
# [30.6, 0.90, -1.68, 0.10]
assert 3 < dot(beta_1[1:], beta_1[1:]) < 4
assert 0.67 < multiple_r_squared(inputs, daily_minutes_good, beta_1) < 0.69
beta_10 = least_squares_fit_ridge(inputs, daily_minutes_good,10, # alpha
learning_rate, 5000, 25)
# [28.3, 0.67, -0.90, -0.01]
assert 1 < dot(beta_10[1:], beta_10[1:]) < 2
assert 0.5 < multiple_r_squared(inputs, daily_minutes_good, beta_10) < 0.6
特别是，“PhD”的系数在增加惩罚时消失，这与我们先前的结果相符，即其与 0 没有显著不同。
注意
通常在使用这种方法之前，你应该重新缩放你的数据。毕竟，如果你将工作经验从年转换为世纪，其最小二乘系数将增加 100 倍，并且突然受到更严重的惩罚，尽管模型是相同的。
另一种方法是套索回归，它使用惩罚项：
def lasso_penalty(beta, alpha):
return alpha * sum(abs(beta_i) for beta_i in beta[1:])
虽然岭回归的惩罚在整体上缩小了系数，但套索惩罚倾向于强制系数为 0，这使其非常适合学习稀疏模型。不幸的是，它不适合梯度下降，这意味着我们无法从头开始解决它。
进一步探索
-
回归背后有丰富而广泛的理论支持。这是另一个你应该考虑阅读教科书或至少大量维基百科文章的地方。
-
scikit-learn 有一个
linear_model模块，提供类似于我们的LinearRegression模型，以及岭回归、套索回归和其他类型的正则化。 -
Statsmodels是另一个 Python 模块，其中包含（除其他内容外）线性回归模型。
第十六章：logistic 回归
很多人说天才和疯狂之间只有一条细微的界限。我认为这不是细微的界限，实际上是一个巨大的鸿沟。
Bill Bailey
在第一章中，我们简要讨论了试图预测哪些 DataSciencester 用户购买高级帐户的问题。在这里，我们将重新讨论这个问题。
问题
我们有一个约 200 个用户的匿名数据集，包含每个用户的工资、作为数据科学家的经验年数，以及她是否为高级帐户支付费用（图 16-1）。像典型的分类变量一样，我们将依赖变量表示为 0（没有高级帐户）或 1（高级帐户）。
像往常一样，我们的数据是一列行 [experience, salary, paid_account]。让我们把它转换成我们需要的格式：
xs = [[1.0] + row[:2] for row in data] # [1, experience, salary]
ys = [row[2] for row in data] # paid_account
显而易见的第一次尝试是使用线性回归找到最佳模型：
paid account = β 0 + β 1 experience + β 2 salary + ε
图 16-1 付费和非付费用户
当然，我们完全可以用这种方式对问题建模。结果显示在图 16-2 中：
from matplotlib import pyplot as plt
from scratch.working_with_data import rescale
from scratch.multiple_regression import least_squares_fit, predict
from scratch.gradient_descent import gradient_step
learning_rate = 0.001
rescaled_xs = rescale(xs)
beta = least_squares_fit(rescaled_xs, ys, learning_rate, 1000, 1)
# [0.26, 0.43, -0.43]
predictions = [predict(x_i, beta) for x_i in rescaled_xs]
plt.scatter(predictions, ys)
plt.xlabel("predicted")
plt.ylabel("actual")
plt.show()
图 16-2 使用线性回归预测高级账户
但是这种方法会导致一些即时问题：
-
我们希望我们预测的输出是 0 或 1，以表示类别成员资格。如果它们在 0 和 1 之间，我们可以将其解释为概率，例如 0.25 的输出可能表示成为付费会员的 25%的机会。但是线性模型的输出可以是非常大的正数甚至是负数，这样不清楚如何解释。实际上，这里很多我们的预测是负数。
-
线性回归模型假设误差与x的列无关。但是在这里，
experience的回归系数为 0.43，表明经验越多，付费会员的可能性越大。这意味着我们的模型对于经验丰富的人输出非常大的值。但是我们知道实际值最多只能是 1，这意味着非常大的输出（因此非常大的experience值）对应于误差项非常大的负值。由于这种情况，我们对beta的估计是有偏的。
我们希望dot(x_i, beta)的大正值对应接近 1 的概率，大负值对应接近 0 的概率。我们可以通过对结果应用另一个函数来实现这一点。
逻辑函数
在逻辑回归中，我们使用逻辑函数，如图 16-3 所示：
def logistic(x: float) -> float:
return 1.0 / (1 + math.exp(-x))
图 16-3. 逻辑函数
当其输入变得很大且为正时，它逐渐接近 1。当其输入变得很大且为负时，它逐渐接近 0。此外，它具有一个方便的性质，即其导数为：
def logistic_prime(x: float) -> float:
y = logistic(x)
return y * (1 - y)
稍后我们将利用这一点。我们将使用这个来拟合一个模型：
y i = f ( x i β ) + ε i
其中f是logistic函数。
回想一下，对于线性回归，我们通过最小化平方误差来拟合模型，这最终选择了最大化数据的似然的β。
在这里两者并不等价，因此我们将使用梯度下降直接最大化似然。这意味着我们需要计算似然函数及其梯度。
给定一些β，我们的模型表明每个y i应该等于 1 的概率为f ( x i β )，等于 0 的概率为1 - f ( x i β )。
特别地，y i的概率密度函数可以写成：
p ( y i | x i , β ) = f (x i β) y i (1-f(x i β)) 1-y i
因为如果y i为 0，则此等于：
1 - f ( x i β )
如果y i为 1，则它等于：
f ( x i β )
结果表明，最大化对数似然实际上更简单：
log L ( β | x i , y i ) = y i log f ( x i β ) + ( 1 - y i ) log ( 1 - f ( x i β ) )
因为对数是一个严格递增的函数，任何使对数似然最大化的beta也将最大化似然，反之亦然。因为梯度下降是最小化的，所以我们实际上会处理负对数似然，因为最大化似然等同于最小化其负值：
import math
from scratch.linear_algebra import Vector, dot
def _negative_log_likelihood(x: Vector, y: float, beta: Vector) -> float:
"""The negative log likelihood for one data point"""
if y == 1:
return -math.log(logistic(dot(x, beta)))
else:
return -math.log(1 - logistic(dot(x, beta)))
如果我们假设不同数据点彼此独立，则总体似然仅为各个似然的乘积。这意味着总体对数似然是各个对数似然的和：
from typing import List
def negative_log_likelihood(xs: List[Vector],
ys: List[float],
beta: Vector) -> float:
return sum(_negative_log_likelihood(x, y, beta)
for x, y in zip(xs, ys))
一点微积分给了我们梯度：
from scratch.linear_algebra import vector_sum
def _negative_log_partial_j(x: Vector, y: float, beta: Vector, j: int) -> float:
"""
The jth partial derivative for one data point.
Here i is the index of the data point.
"""
return -(y - logistic(dot(x, beta))) * x[j]
def _negative_log_gradient(x: Vector, y: float, beta: Vector) -> Vector:
"""
The gradient for one data point.
"""
return [_negative_log_partial_j(x, y, beta, j)
for j in range(len(beta))]
def negative_log_gradient(xs: List[Vector],
ys: List[float],
beta: Vector) -> Vector:
return vector_sum([_negative_log_gradient(x, y, beta)
for x, y in zip(xs, ys)])
到这一点我们已经拥有所有需要的部分。
应用模型
我们将要将数据分成训练集和测试集：
from scratch.machine_learning import train_test_split
import random
import tqdm
random.seed(0)
x_train, x_test, y_train, y_test = train_test_split(rescaled_xs, ys, 0.33)
learning_rate = 0.01
# pick a random starting point
beta = [random.random() for _ in range(3)]
with tqdm.trange(5000) as t:
for epoch in t:
gradient = negative_log_gradient(x_train, y_train, beta)
beta = gradient_step(beta, gradient, -learning_rate)
loss = negative_log_likelihood(x_train, y_train, beta)
t.set_description(f"loss: {loss:.3f} beta: {beta}")
之后我们发现beta大约是：
[-2.0, 4.7, -4.5]
这些是rescaled 数据的系数，但我们也可以将它们转换回原始数据：
from scratch.working_with_data import scale
means, stdevs = scale(xs)
beta_unscaled = [(beta[0]
- beta[1] * means[1] / stdevs[1]
- beta[2] * means[2] / stdevs[2]),
beta[1] / stdevs[1],
beta[2] / stdevs[2]]
# [8.9, 1.6, -0.000288]
不幸的是，这些并不像线性回归系数那样容易解释。其他条件相同，一年额外的经验将在logistic的输入上增加 1.6。其他条件相同，额外的 10,000 美元薪水将在logistic的输入上减少 2.88。
然而，输出的影响也取决于其他输入。如果dot(beta, x_i)已经很大（对应概率接近 1），即使大幅增加它也不能太大影响概率。如果接近 0，稍微增加可能会大幅增加概率。
我们可以说的是，其他条件相同的情况下，经验丰富的人更有可能支付账户。而其他条件相同的情况下，收入较高的人支付账户的可能性较低。（这在我们绘制数据时也有些明显。）
拟合度
我们还没有使用我们留出的测试数据。让我们看看如果我们在概率超过 0.5 时预测付费账户会发生什么：
true_positives = false_positives = true_negatives = false_negatives = 0
for x_i, y_i in zip(x_test, y_test):
prediction = logistic(dot(beta, x_i))
if y_i == 1 and prediction >= 0.5: # TP: paid and we predict paid
true_positives += 1
elif y_i == 1: # FN: paid and we predict unpaid
false_negatives += 1
elif prediction >= 0.5: # FP: unpaid and we predict paid
false_positives += 1
else: # TN: unpaid and we predict unpaid
true_negatives += 1
precision = true_positives / (true_positives + false_positives)
recall = true_positives / (true_positives + false_negatives)
这给出了 75%的精度（“当我们预测付费账户时，我们有 75%的准确率”）和 80%的召回率（“当用户有付费账户时，我们 80%的时间预测为付费账户”），考虑到我们拥有的数据很少，这并不算糟糕。
我们还可以绘制预测与实际情况对比图（见图 16-4，#logistic_prediction_vs_actual），这也显示出模型表现良好：
predictions = [logistic(dot(beta, x_i)) for x_i in x_test]
plt.scatter(predictions, y_test, marker='+')
plt.xlabel("predicted probability")
plt.ylabel("actual outcome")
plt.title("Logistic Regression Predicted vs. Actual")
plt.show()
图 16-4。逻辑回归预测与实际
支持向量机
dot(beta, x_i)等于 0 的点集是我们类之间的边界。我们可以绘制这个图来精确了解我们的模型在做什么（见图 16-5，#logit_image_part_two）。
图 16-5。付费和未付费用户与决策边界
这个边界是一个超平面，将参数空间分成两个半空间，对应预测付费和预测未付费。我们发现它是在找到最可能的逻辑模型的副作用中发现的。
分类的另一种方法是只需寻找在训练数据中“最佳”分离类别的超平面。这是支持向量机背后的思想，它找到最大化每个类别最近点到超平面距离的超平面（见图 16-6，#separating_hyperplane）。
图 16-6。一个分离超平面
寻找这样的超平面是一个涉及我们过于高级的技术的优化问题。另一个问题是，一个分离超平面可能根本不存在。在我们的“谁付费？”数据集中，简单地没有一条线完全分离付费用户和未付费用户。
我们有时可以通过将数据转换为更高维空间来绕过这个问题。例如，考虑到简单的一维数据集，如图 16-7 所示。
图 16-7。一个不可分离的一维数据集
从明显的来看，没有一个超平面可以将正例和负例分开。然而，当我们通过将点x映射到(x, x**2)的方式将这个数据集映射到二维空间时，看看会发生什么。突然间，可以找到一个可以分割数据的超平面（见图 16-8）。
图 16-8. 数据集在更高维度中变得可分离
这通常被称为核技巧，因为我们不是实际将点映射到更高维度的空间（如果有很多点并且映射很复杂的话，这可能会很昂贵），而是可以使用一个“核”函数来计算在更高维度空间中的点积，并使用这些来找到超平面。
使用支持向量机而不依赖于由具有适当专业知识的人编写的专业优化软件是困难的（也可能不是一个好主意），因此我们将不得不在这里结束我们的讨论。
进一步调查
第十七章：决策树
树是一个难以理解的神秘。
Jim Woodring
DataSciencester 的人才副总裁已经面试了一些来自该网站的求职者，结果各不相同。他收集了一个数据集，其中包含每个候选人的几个（定性）属性，以及该候选人是否面试表现良好或不佳。他问道，你能否利用这些数据建立一个模型，识别哪些候选人会面试表现良好，这样他就不必浪费时间进行面试了？
这似乎非常适合决策树，这是数据科学家工具包中的另一种预测建模工具。
什么是决策树？
决策树使用树结构来表示多条可能的决策路径和每条路径的结果。
如果你玩过Twenty Questions游戏，那么你就熟悉决策树了。例如：
-
“我在想一个动物。”
-
“它有五条腿以上吗？”
-
“不。”
-
“它好吃吗？”
-
“不。”
-
“它出现在澳大利亚五分硬币的背面吗？”
-
“是的。”
-
“它是针鼹吗？”
-
“是的，就是它！”
这对应于路径：
“不超过 5 条腿” → “不好吃” → “在 5 分硬币上” → “针鼹！”
在一个古怪（并不是很全面的）“猜动物”决策树中（Figure 17-1）。
图 17-1. “猜动物”决策树
决策树有很多优点。它们非常容易理解和解释，它们达到预测的过程完全透明。与我们迄今所看到的其他模型不同，决策树可以轻松处理数值型（例如，腿的数量）和分类型（例如，好吃/不好吃）属性的混合数据，甚至可以对缺少属性的数据进行分类。
与此同时，为一组训练数据找到一个“最优”决策树在计算上是一个非常困难的问题。（我们将通过尝试构建一个足够好的树来避开这个问题，尽管对于大型数据集来说，这仍然可能是一项艰巨的工作。）更重要的是，很容易（也很糟糕）构建过度拟合训练数据的决策树，这些树在未见数据上的泛化能力很差。我们将探讨解决这个问题的方法。
大多数人将决策树分为分类树（生成分类输出）和回归树（生成数值输出）。在本章中，我们将专注于分类树，并通过 ID3 算法从一组标记数据中学习决策树，这将帮助我们理解决策树的实际工作原理。为了简化问题，我们将局限于具有二元输出的问题，例如“我应该雇佣这位候选人吗？”或“我应该向这位网站访客展示广告 A 还是广告 B？”或“我在办公室冰箱里找到的这种食物会让我生病吗？”
熵
要构建决策树，我们需要决定提出什么问题以及顺序。在树的每个阶段，我们消除了一些可能性，还有一些没有。学到动物的腿不超过五条后，我们排除了它是蚱蜢的可能性。我们还没排除它是一只鸭子的可能性。每个可能的问题根据其答案将剩余的可能性进行分区。
理想情况下，我们希望选择的问题答案能够提供关于我们的树应该预测什么的大量信息。如果有一个单一的是/否问题，“是”答案总是对应于 True 输出，而“否”答案对应于 False 输出（反之亦然），那这将是一个很棒的问题选择。相反，如果一个是/否问题的任何答案都不能给出关于预测应该是什么的新信息，那可能不是一个好选择。
我们用熵来捕捉这种“信息量”概念。你可能听说过这个术语用来表示无序。我们用它来表示与数据相关的不确定性。
想象我们有一个数据集 S，其中每个成员都被标记为属于有限数量的类别 C 1 , ... , C n 中的一种。如果所有数据点属于同一类，则没有真正的不确定性，这意味着我们希望熵很低。如果数据点均匀分布在各个类别中，就会有很多不确定性，我们希望熵很高。
从数学角度来看，如果 p i 是标记为类别 c i 的数据的比例，我们定义熵如下：
H ( S ) = - p 1 log 2 p 1 - ... - p n log 2 p n
按照（标准）惯例，0 log 0 = 0 。
不必过分担心可怕的细节，每个术语 - p i log 2 p i 都是非负的，当 p i 接近于 0 或接近于 1 时，它接近于 0（图 17-2）。
图 17-2. -p log p 的图示
这意味着当每个 p i 接近于 0 或 1 时（即大多数数据属于单一类别时），熵将很小，当许多 p i 不接近于 0 时（即数据分布在多个类别中时），熵将较大。这正是我们期望的行为。
将所有这些内容整合到一个函数中是相当简单的：
from typing import List
import math
def entropy(class_probabilities: List[float]) -> float:
"""Given a list of class probabilities, compute the entropy"""
return sum(-p * math.log(p, 2)
for p in class_probabilities
if p > 0) # ignore zero probabilities
assert entropy([1.0]) == 0
assert entropy([0.5, 0.5]) == 1
assert 0.81 < entropy([0.25, 0.75]) < 0.82
我们的数据将由成对的(输入，标签)组成，这意味着我们需要自己计算类别概率。注意，我们实际上并不关心每个概率与哪个标签相关联，只关心这些概率是多少：
from typing import Any
from collections import Counter
def class_probabilities(labels: List[Any]) -> List[float]:
total_count = len(labels)
return [count / total_count
for count in Counter(labels).values()]
def data_entropy(labels: List[Any]) -> float:
return entropy(class_probabilities(labels))
assert data_entropy(['a']) == 0
assert data_entropy([True, False]) == 1
assert data_entropy([3, 4, 4, 4]) == entropy([0.25, 0.75])
分区的熵
到目前为止，我们所做的是计算单一标记数据集的熵（想想“不确定性”）。现在，决策树的每个阶段都涉及提出一个问题，其答案将数据分成一个或多个子集。例如，我们的“它有超过五条腿吗？”问题将动物分成有超过五条腿的动物（例如，蜘蛛）和没有超过五条腿的动物（例如，针鼹）。
相应地，我们希望从某种程度上了解通过某种方式对数据集进行分区所产生的熵。如果分区将数据分成具有低熵（即高度确定）的子集，我们希望分区具有低熵；如果包含具有高熵（即高度不确定）的子集（即大且）分区具有高熵。
例如，我的“澳大利亚五分硬币”问题非常愚蠢（尽管非常幸运！），因为它将那时剩下的动物分成S 1 = {针鼹}和S 2 = {其他所有动物}，其中S 2既大又高熵。 (S 1没有熵，但它表示剩余“类别”的小部分。)
从数学上讲，如果我们将我们的数据S分成包含数据比例的子集S 1 , ... , S m，那么我们将分区的熵计算为加权和：
H = q 1 H ( S 1 ) + ... + q m H ( S m )
我们可以实现为：
def partition_entropy(subsets: List[List[Any]]) -> float:
"""Returns the entropy from this partition of data into subsets"""
total_count = sum(len(subset) for subset in subsets)
return sum(data_entropy(subset) * len(subset) / total_count
for subset in subsets)
注意
这种方法的一个问题是，使用具有许多不同值的属性进行分区会由于过拟合而导致熵非常低。例如，假设你在银行工作，试图建立一个决策树来预测哪些客户可能会违约他们的抵押贷款，使用一些历史数据作为你的训练集。进一步假设数据集包含每个客户的社会安全号码。在社会安全号码上进行分区将产生单个人的子集，每个子集的熵必然为零。但是依赖社会安全号码的模型肯定无法超出训练集的范围。因此，在创建决策树时，你应该尽量避免（或适当地分桶）具有大量可能值的属性。
创建决策树
副总裁为您提供了面试者数据，根据您的规定，每位候选人的相关属性是一个NamedTuple——她的级别、她偏爱的语言、她是否在 Twitter 上活跃、她是否有博士学位以及她是否面试表现良好：
from typing import NamedTuple, Optional
class Candidate(NamedTuple):
level: str
lang: str
tweets: bool
phd: bool
did_well: Optional[bool] = None # allow unlabeled data
# level lang tweets phd did_well
inputs = [Candidate('Senior', 'Java', False, False, False),
Candidate('Senior', 'Java', False, True, False),
Candidate('Mid', 'Python', False, False, True),
Candidate('Junior', 'Python', False, False, True),
Candidate('Junior', 'R', True, False, True),
Candidate('Junior', 'R', True, True, False),
Candidate('Mid', 'R', True, True, True),
Candidate('Senior', 'Python', False, False, False),
Candidate('Senior', 'R', True, False, True),
Candidate('Junior', 'Python', True, False, True),
Candidate('Senior', 'Python', True, True, True),
Candidate('Mid', 'Python', False, True, True),
Candidate('Mid', 'Java', True, False, True),
Candidate('Junior', 'Python', False, True, False)
]
我们的树将包含决策节点（提出问题并根据答案引导我们不同路径）和叶节点（给出预测）。我们将使用相对简单的ID3算法构建它，该算法操作如下。假设我们有一些带标签的数据，并且有一个要考虑分支的属性列表：
-
如果所有数据都具有相同的标签，请创建一个预测该标签的叶节点，然后停止。
-
如果属性列表为空（即没有更多可能的问题可问），创建一个预测最常见标签的叶节点，然后停止。
-
否则，尝试按照每个属性对数据进行分割。
-
选择具有最低分区熵的分区。
-
根据选择的属性添加一个决策节点。
-
对剩余属性使用递归对每个分区的子集进行分割。
这被称为“贪婪”算法，因为在每一步中，它选择最即时的最佳选项。给定一个数据集，可能有一个看起来更差的第一步，但却会有一个更好的树。如果确实存在这样的情况，此算法将无法找到它。尽管如此，它相对容易理解和实现，这使得它成为探索决策树的一个很好的起点。
让我们手动按照面试者数据集中的这些步骤进行。数据集具有True和False标签，并且我们有四个可以分割的属性。因此，我们的第一步将是找到熵最小的分区。我们将首先编写一个执行分割的函数：
from typing import Dict, TypeVar
from collections import defaultdict
T = TypeVar('T') # generic type for inputs
def partition_by(inputs: List[T], attribute: str) -> Dict[Any, List[T]]:
"""Partition the inputs into lists based on the specified attribute."""
partitions: Dict[Any, List[T]] = defaultdict(list)
for input in inputs:
key = getattr(input, attribute) # value of the specified attribute
partitions[key].append(input) # add input to the correct partition
return partitions
还有一个使用它计算熵的函数：
def partition_entropy_by(inputs: List[Any],
attribute: str,
label_attribute: str) -> float:
"""Compute the entropy corresponding to the given partition"""
# partitions consist of our inputs
partitions = partition_by(inputs, attribute)
# but partition_entropy needs just the class labels
labels = [[getattr(input, label_attribute) for input in partition]
for partition in partitions.values()]
return partition_entropy(labels)
然后我们只需为整个数据集找到最小熵分区：
for key in ['level','lang','tweets','phd']:
print(key, partition_entropy_by(inputs, key, 'did_well'))
assert 0.69 < partition_entropy_by(inputs, 'level', 'did_well') < 0.70
assert 0.86 < partition_entropy_by(inputs, 'lang', 'did_well') < 0.87
assert 0.78 < partition_entropy_by(inputs, 'tweets', 'did_well') < 0.79
assert 0.89 < partition_entropy_by(inputs, 'phd', 'did_well') < 0.90
最低熵来自于按level分割，因此我们需要为每个可能的level值创建一个子树。每个Mid候选人都标记为True，这意味着Mid子树只是一个预测True的叶节点。对于Senior候选人，我们有True和False的混合，因此我们需要再次分割：
senior_inputs = [input for input in inputs if input.level == 'Senior']
assert 0.4 == partition_entropy_by(senior_inputs, 'lang', 'did_well')
assert 0.0 == partition_entropy_by(senior_inputs, 'tweets', 'did_well')
assert 0.95 < partition_entropy_by(senior_inputs, 'phd', 'did_well') < 0.96
这告诉我们我们接下来应该在tweets上进行分割，这会导致零熵分区。对于这些Senior级别的候选人，“是”推文总是导致True，而“否”推文总是导致False。
最后，如果我们对Junior候选人执行相同的操作，我们最终会在phd上进行分割，之后发现没有博士学位总是导致True，而有博士学位总是导致False。
图 17-3 显示了完整的决策树。
图 17-3. 招聘的决策树
把所有这些整合起来
现在我们已经了解了算法的工作原理，我们希望更普遍地实现它。这意味着我们需要决定如何表示树。我们将使用可能最轻量级的表示。我们将一个树定义为以下内容之一：
-
一个
Leaf（预测单个值），或者 -
一个
Split（包含要拆分的属性、特定属性值的子树，以及在遇到未知值时可能使用的默认值）。from typing import NamedTuple, Union, Any class Leaf(NamedTuple): value: Any class Split(NamedTuple): attribute: str subtrees: dict default_value: Any = None DecisionTree = Union[Leaf, Split]
有了这个表示，我们的招聘树将如下所示：
hiring_tree = Split('level', { # first, consider "level"
'Junior': Split('phd', { # if level is "Junior", next look at "phd"
False: Leaf(True), # if "phd" is False, predict True
True: Leaf(False) # if "phd" is True, predict False
}),
'Mid': Leaf(True), # if level is "Mid", just predict True
'Senior': Split('tweets', { # if level is "Senior", look at "tweets"
False: Leaf(False), # if "tweets" is False, predict False
True: Leaf(True) # if "tweets" is True, predict True
})
})
还有一个问题，即如果我们遇到意外的（或缺失的）属性值该怎么办。如果我们的招聘树遇到level是Intern的候选人会怎么样？我们将通过用最常见的标签填充default_value属性来处理这种情况。
给定这样的表示，我们可以对输入进行分类：
def classify(tree: DecisionTree, input: Any) -> Any:
"""classify the input using the given decision tree"""
# If this is a leaf node, return its value
if isinstance(tree, Leaf):
return tree.value
# Otherwise this tree consists of an attribute to split on
# and a dictionary whose keys are values of that attribute
# and whose values are subtrees to consider next
subtree_key = getattr(input, tree.attribute)
if subtree_key not in tree.subtrees: # If no subtree for key,
return tree.default_value # return the default value.
subtree = tree.subtrees[subtree_key] # Choose the appropriate subtree
return classify(subtree, input) # and use it to classify the input.
剩下的就是从我们的训练数据中构建树表示：
def build_tree_id3(inputs: List[Any],
split_attributes: List[str],
target_attribute: str) -> DecisionTree:
# Count target labels
label_counts = Counter(getattr(input, target_attribute)
for input in inputs)
most_common_label = label_counts.most_common(1)[0][0]
# If there's a unique label, predict it
if len(label_counts) == 1:
return Leaf(most_common_label)
# If no split attributes left, return the majority label
if not split_attributes:
return Leaf(most_common_label)
# Otherwise split by the best attribute
def split_entropy(attribute: str) -> float:
"""Helper function for finding the best attribute"""
return partition_entropy_by(inputs, attribute, target_attribute)
best_attribute = min(split_attributes, key=split_entropy)
partitions = partition_by(inputs, best_attribute)
new_attributes = [a for a in split_attributes if a != best_attribute]
# Recursively build the subtrees
subtrees = {attribute_value : build_tree_id3(subset,
new_attributes,
target_attribute)
for attribute_value, subset in partitions.items()}
return Split(best_attribute, subtrees, default_value=most_common_label)
在我们构建的树中，每个叶子节点完全由True输入或完全由False输入组成。这意味着该树在训练数据集上的预测完全正确。但我们也可以将其应用于训练集中不存在的新数据：
tree = build_tree_id3(inputs,
['level', 'lang', 'tweets', 'phd'],
'did_well')
# Should predict True
assert classify(tree, Candidate("Junior", "Java", True, False))
# Should predict False
assert not classify(tree, Candidate("Junior", "Java", True, True))
也适用于具有意外值的数据：
# Should predict True
assert classify(tree, Candidate("Intern", "Java", True, True))
注意
由于我们的目标主要是演示如何构建树，所以我们使用整个数据集构建了树。一如既往，如果我们真的试图为某事创建一个良好的模型，我们会收集更多数据并将其分割成训练/验证/测试子集。
随机森林
鉴于决策树可以如此紧密地适应其训练数据，他们很容易出现过拟合的倾向并不奇怪。避免这种情况的一种方法是一种称为随机森林的技术，其中我们构建多个决策树并组合它们的输出。如果它们是分类树，我们可以让它们投票；如果它们是回归树，我们可以平均它们的预测。
我们的树构建过程是确定性的，那么我们如何获得随机树呢？
其中一部分涉及引导数据（参见“插曲：自助法”）。我们不是在整个训练集上训练每棵树，而是在bootstrap_sample(inputs)的结果上训练每棵树。由于每棵树都是使用不同的数据构建的，因此每棵树与其他每棵树都不同。（一个副作用是，使用未抽样数据来测试每棵树是完全公平的，这意味着如果在衡量性能时聪明地使用所有数据作为训练集，你就可以获得成功。）这种技术被称为自助聚合或装袋。
第二个随机性源涉及更改选择要拆分的best_attribute的方法。我们不是查看所有剩余属性，而是首先选择它们的随机子集，然后在其中最好的属性上进行拆分：
# if there are already few enough split candidates, look at all of them
if len(split_candidates) <= self.num_split_candidates:
sampled_split_candidates = split_candidates
# otherwise pick a random sample
else:
sampled_split_candidates = random.sample(split_candidates,
self.num_split_candidates)
# now choose the best attribute only from those candidates
best_attribute = min(sampled_split_candidates, key=split_entropy)
partitions = partition_by(inputs, best_attribute)
这是集成学习的一个例子，其中我们结合了多个弱学习器（通常是高偏差、低方差模型），以生成一个总体上强大的模型。
进一步探索
-
scikit-learn 包含许多决策树模型。它还有一个
ensemble模块，其中包括RandomForestClassifier以及其他集成方法。 -
XGBoost 是一个用于训练梯度提升决策树的库，经常在许多 Kaggle 风格的机器学习竞赛中获胜。
-
我们只是浅尝了解决策树及其算法的表面。维基百科是一个更广泛探索的良好起点。
第十八章：神经网络
我喜欢胡说八道；它唤醒了大脑细胞。
博士苏斯
人工神经网络（或简称神经网络）是一种受大脑运作方式启发的预测模型。将大脑视为一组互相连接的神经元。每个神经元查看输入到它的其他神经元的输出，进行计算，然后根据计算结果是否超过某个阈值来“激活”或不激活。
因此，人工神经网络由人工神经元组成，它们对其输入执行类似的计算。神经网络可以解决各种问题，如手写识别和面部检测，在深度学习中广泛使用，这是数据科学中最流行的子领域之一。然而，大多数神经网络是“黑箱”—检查它们的细节并不能让你理解它们如何解决问题。而且，大型神经网络可能难以训练。对于大多数初涉数据科学的问题，它们可能不是正确的选择。当你试图构建一个引发“奇点”的人工智能时，它们可能非常合适。
感知器
最简单的神经网络就是感知器，它模拟了一个具有n个二进制输入的单个神经元。它计算其输入的加权和，如果该加权和大于等于 0，则“激活”：
from scratch.linear_algebra import Vector, dot
def step_function(x: float) -> float:
return 1.0 if x >= 0 else 0.0
def perceptron_output(weights: Vector, bias: float, x: Vector) -> float:
"""Returns 1 if the perceptron 'fires', 0 if not"""
calculation = dot(weights, x) + bias
return step_function(calculation)
感知器只是区分了由点x组成的超平面分隔的半空间：
dot(weights, x) + bias == 0
通过适当选择的权重，感知器可以解决许多简单的问题（图 18-1）。例如，我们可以创建一个AND 门（如果其输入都是 1 则返回 1，但如果其中一个输入是 0 则返回 0）：
and_weights = [2., 2]
and_bias = -3.
assert perceptron_output(and_weights, and_bias, [1, 1]) == 1
assert perceptron_output(and_weights, and_bias, [0, 1]) == 0
assert perceptron_output(and_weights, and_bias, [1, 0]) == 0
assert perceptron_output(and_weights, and_bias, [0, 0]) == 0
如果两个输入都是 1，计算等于 2 + 2 – 3 = 1，并且输出为 1。如果只有一个输入是 1，计算等于 2 + 0 – 3 = –1，并且输出为 0。如果两个输入都是 0，计算等于–3，并且输出为 0。
使用类似的推理，我们可以用以下方式构建一个OR 门：
or_weights = [2., 2]
or_bias = -1.
assert perceptron_output(or_weights, or_bias, [1, 1]) == 1
assert perceptron_output(or_weights, or_bias, [0, 1]) == 1
assert perceptron_output(or_weights, or_bias, [1, 0]) == 1
assert perceptron_output(or_weights, or_bias, [0, 0]) == 0
图 18-1. 二输入感知器的决策空间
我们也可以用以下方式构建一个NOT 门（其只有一个输入，并将 1 转换为 0，将 0 转换为 1）：
not_weights = [-2.]
not_bias = 1.
assert perceptron_output(not_weights, not_bias, [0]) == 1
assert perceptron_output(not_weights, not_bias, [1]) == 0
然而，有一些问题简单地无法通过单个感知器解决。例如，无论你多么努力，你都不能使用一个感知器来构建一个异或门，即当其输入中只有一个是 1 时输出为 1，否则为 0。这时我们需要更复杂的神经网络。
当然，你不需要模拟神经元来构建逻辑门：
and_gate = min
or_gate = max
xor_gate = lambda x, y: 0 if x == y else 1
像真实的神经元一样，人工神经元开始在连接起来时变得更加有趣。
前馈神经网络
大脑的拓扑结构极其复杂，因此常常用由离散的层组成的理想化前馈神经网络来近似它，每个层都与下一层相连接。通常包括一个输入层（接收并不改变输入并向前传递），一个或多个“隐藏层”（每个都由神经元组成，接收前一层的输出，执行某些计算，并将结果传递到下一层），和一个输出层（生成最终的输出）。
就像在感知机中一样，每个（非输入）神经元都有与其每个输入对应的权重和偏置。为了简化我们的表示，我们将偏置添加到权重向量的末尾，并给每个神经元一个偏置输入，其值始终为 1。
与感知机类似，对于每个神经元，我们将计算其输入和权重的乘积之和。但在这里，我们不会输出应用于该乘积的step_function，而是输出其平滑的近似。这里我们将使用sigmoid函数（图 18-2）：
import math
def sigmoid(t: float) -> float:
return 1 / (1 + math.exp(-t))
图 18-2. sigmoid 函数
为什么使用sigmoid而不是更简单的step_function？为了训练神经网络，我们需要使用微积分，而要使用微积分，我们需要平滑的函数。step_function甚至不连续，而sigmoid是它的一个很好的平滑近似。
注意
你可能还记得sigmoid函数，它在第十六章中称为logistic。技术上，sigmoid指的是函数的形状，logistic指的是特定的函数，尽管人们经常将这些术语互换使用。
然后我们计算输出为：
def neuron_output(weights: Vector, inputs: Vector) -> float:
# weights includes the bias term, inputs includes a 1
return sigmoid(dot(weights, inputs))
有了这个函数，我们可以简单地将神经元表示为一个权重向量，其长度比该神经元的输入数量多一个（因为有偏置权重）。然后我们可以将神经网络表示为（非输入）层的列表，其中每一层只是该层中神经元的列表。
那就是，我们将神经网络表示为向量的列表（层），其中每个向量（神经元）是向量（权重）的列表。
给定这样的表示法，使用神经网络非常简单：
from typing import List
def feed_forward(neural_network: List[List[Vector]],
input_vector: Vector) -> List[Vector]:
"""
Feeds the input vector through the neural network.
Returns the outputs of all layers (not just the last one).
"""
outputs: List[Vector] = []
for layer in neural_network:
input_with_bias = input_vector + [1] # Add a constant.
output = [neuron_output(neuron, input_with_bias) # Compute the output
for neuron in layer] # for each neuron.
outputs.append(output) # Add to results.
# Then the input to the next layer is the output of this one
input_vector = output
return outputs
现在很容易构建我们无法用单个感知机构建的 XOR 门。我们只需扩展权重，使得neuron_output要么非常接近 0，要么非常接近 1：
xor_network = [# hidden layer
[[20., 20, -30], # 'and' neuron
[20., 20, -10]], # 'or' neuron
# output layer
[[-60., 60, -30]]] # '2nd input but not 1st input' neuron
# feed_forward returns the outputs of all layers, so the [-1] gets the
# final output, and the [0] gets the value out of the resulting vector
assert 0.000 < feed_forward(xor_network, [0, 0])[-1][0] < 0.001
assert 0.999 < feed_forward(xor_network, [1, 0])[-1][0] < 1.000
assert 0.999 < feed_forward(xor_network, [0, 1])[-1][0] < 1.000
assert 0.000 < feed_forward(xor_network, [1, 1])[-1][0] < 0.001
对于给定的输入（一个二维向量），隐藏层生成一个二维向量，包含两个输入值的“与”和两个输入值的“或”。
输出层接收一个二维向量，并计算“第二个元素但不是第一个元素”。结果是一个执行“或但不是与”的网络，这正是 XOR（图 18-3）。
图 18-3. 一个用于 XOR 的神经网络
一个有启发性的思考方式是，隐藏层正在计算输入数据的特征（在本例中是“and”和“or”），输出层将这些特征组合起来以生成所需的输出。
反向传播
通常我们不会手工构建神经网络。部分原因是因为我们使用它们来解决更大的问题——例如图像识别问题可能涉及数百或数千个神经元。另一部分原因是通常我们无法“推理出”神经元应该是什么。
相反（通常情况下），我们使用数据来训练神经网络。典型的方法是一种称为反向传播的算法，它使用梯度下降或其变种之一。
假设我们有一个训练集，由输入向量和相应的目标输出向量组成。例如，在我们之前的xor_network示例中，输入向量[1, 0]对应于目标输出[1]。假设我们的网络有一些权重。然后，我们使用以下算法调整权重：
-
对输入向量运行
feed_forward以产生网络中所有神经元的输出。 -
我们知道目标输出，所以我们可以计算一个损失，即平方误差的总和。
-
计算这种损失作为输出神经元权重的函数的梯度。
-
“传播”梯度和误差向后计算与隐藏神经元权重相关的梯度。
-
进行梯度下降步骤。
通常我们会对整个训练集运行这个算法多次，直到网络收敛。
首先，让我们编写计算梯度的函数：
def sqerror_gradients(network: List[List[Vector]],
input_vector: Vector,
target_vector: Vector) -> List[List[Vector]]:
"""
Given a neural network, an input vector, and a target vector,
make a prediction and compute the gradient of the squared error
loss with respect to the neuron weights.
"""
# forward pass
hidden_outputs, outputs = feed_forward(network, input_vector)
# gradients with respect to output neuron pre-activation outputs
output_deltas = [output * (1 - output) * (output - target)
for output, target in zip(outputs, target_vector)]
# gradients with respect to output neuron weights
output_grads = [[output_deltas[i] * hidden_output
for hidden_output in hidden_outputs + [1]]
for i, output_neuron in enumerate(network[-1])]
# gradients with respect to hidden neuron pre-activation outputs
hidden_deltas = [hidden_output * (1 - hidden_output) *
dot(output_deltas, [n[i] for n in network[-1]])
for i, hidden_output in enumerate(hidden_outputs)]
# gradients with respect to hidden neuron weights
hidden_grads = [[hidden_deltas[i] * input for input in input_vector + [1]]
for i, hidden_neuron in enumerate(network[0])]
return [hidden_grads, output_grads]
前述计算背后的数学并不是非常困难，但涉及一些繁琐的微积分和仔细的注意细节，所以我会留给你作为练习。
有了计算梯度的能力，我们现在可以训练神经网络。让我们试着通过手工设计的 XOR 网络来学习它。
我们将从生成训练数据开始，并使用随机权重初始化我们的神经网络：
import random
random.seed(0)
# training data
xs = [[0., 0], [0., 1], [1., 0], [1., 1]]
ys = [[0.], [1.], [1.], [0.]]
# start with random weights
network = [ # hidden layer: 2 inputs -> 2 outputs
[[random.random() for _ in range(2 + 1)], # 1st hidden neuron
[random.random() for _ in range(2 + 1)]], # 2nd hidden neuron
# output layer: 2 inputs -> 1 output
[[random.random() for _ in range(2 + 1)]] # 1st output neuron
]
通常情况下，我们可以使用梯度下降来训练它。与我们之前的例子不同之处在于，这里我们有几个参数向量，每个向量都有自己的梯度，这意味着我们需要为每个向量调用gradient_step。
from scratch.gradient_descent import gradient_step
import tqdm
learning_rate = 1.0
for epoch in tqdm.trange(20000, desc="neural net for xor"):
for x, y in zip(xs, ys):
gradients = sqerror_gradients(network, x, y)
# Take a gradient step for each neuron in each layer
network = [[gradient_step(neuron, grad, -learning_rate)
for neuron, grad in zip(layer, layer_grad)]
for layer, layer_grad in zip(network, gradients)]
# check that it learned XOR
assert feed_forward(network, [0, 0])[-1][0] < 0.01
assert feed_forward(network, [0, 1])[-1][0] > 0.99
assert feed_forward(network, [1, 0])[-1][0] > 0.99
assert feed_forward(network, [1, 1])[-1][0] < 0.01
对我来说，得到的网络具有如下权重：
[ # hidden layer
[[7, 7, -3], # computes OR
[5, 5, -8]], # computes AND
# output layer
[[11, -12, -5]] # computes "first but not second"
]
从概念上讲，这与我们之前的定制网络非常相似。
示例：Fizz Buzz
工程副总裁希望通过让技术候选人解决“Fizz Buzz”，来面试技术候选人，这是一个广为人知的编程挑战：
Print the numbers 1 to 100, except that if the number is divisible
by 3, print "fizz"; if the number is divisible by 5, print "buzz";
and if the number is divisible by 15, print "fizzbuzz".
他认为解决这个问题表明了极限编程技能。你认为这个问题如此简单，以至于一个神经网络可以解决它。
神经网络将向量作为输入，并产生向量作为输出。正如所述，编程问题是将整数转换为字符串。因此，第一个挑战是想出一种将其重新定义为向量问题的方法。
对于输出来说并不困难：基本上有四类输出，所以我们可以将输出编码为一个包含四个 0 和 1 的向量：
def fizz_buzz_encode(x: int) -> Vector:
if x % 15 == 0:
return [0, 0, 0, 1]
elif x % 5 == 0:
return [0, 0, 1, 0]
elif x % 3 == 0:
return [0, 1, 0, 0]
else:
return [1, 0, 0, 0]
assert fizz_buzz_encode(2) == [1, 0, 0, 0]
assert fizz_buzz_encode(6) == [0, 1, 0, 0]
assert fizz_buzz_encode(10) == [0, 0, 1, 0]
assert fizz_buzz_encode(30) == [0, 0, 0, 1]
我们将用这个来生成我们的目标向量。输入向量则不那么明显。你不应该只使用一个包含输入数字的一维向量，因为有几个原因。一个单一的输入捕捉到了“强度”，但是 2 是 1 的两倍，4 又是两倍，对这个问题并不感兴趣。此外，只有一个输入，隐藏层将无法计算非常有趣的特征，这意味着它可能无法解决问题。
原来，一个相对有效的方法是将每个数字转换为其二进制表示，由 1 和 0 组成。（别担心，这不明显——至少对我来说不是。）
def binary_encode(x: int) -> Vector:
binary: List[float] = []
for i in range(10):
binary.append(x % 2)
x = x // 2
return binary
# 1 2 4 8 16 32 64 128 256 512
assert binary_encode(0) == [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
assert binary_encode(1) == [1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
assert binary_encode(10) == [0, 1, 0, 1, 0, 0, 0, 0, 0, 0]
assert binary_encode(101) == [1, 0, 1, 0, 0, 1, 1, 0, 0, 0]
assert binary_encode(999) == [1, 1, 1, 0, 0, 1, 1, 1, 1, 1]
因为目标是构建 1 到 100 的输出，所以在这些数字上进行训练是作弊的。因此，我们将在数字 101 到 1,023 上进行训练（这是我们可以用 10 位二进制表示的最大数字）：
xs = [binary_encode(n) for n in range(101, 1024)]
ys = [fizz_buzz_encode(n) for n in range(101, 1024)]
接下来，让我们创建一个具有随机初始权重的神经网络。它将有 10 个输入神经元（因为我们将我们的输入表示为 10 维向量）和 4 个输出神经元（因为我们将我们的目标表示为 4 维向量）。我们将给它 25 个隐藏单元，但我们将使用一个变量，所以这样很容易改变：
NUM_HIDDEN = 25
network = [
# hidden layer: 10 inputs -> NUM_HIDDEN outputs
[[random.random() for _ in range(10 + 1)] for _ in range(NUM_HIDDEN)],
# output_layer: NUM_HIDDEN inputs -> 4 outputs
[[random.random() for _ in range(NUM_HIDDEN + 1)] for _ in range(4)]
]
那就是了。现在我们准备好训练了。因为这是一个更复杂的问题（而且有很多事情可能会出错），我们希望密切监控训练过程。特别是，对于每个时期，我们将跟踪平方误差的总和并将其打印出来。我们希望确保它们会减少：
from scratch.linear_algebra import squared_distance
learning_rate = 1.0
with tqdm.trange(500) as t:
for epoch in t:
epoch_loss = 0.0
for x, y in zip(xs, ys):
predicted = feed_forward(network, x)[-1]
epoch_loss += squared_distance(predicted, y)
gradients = sqerror_gradients(network, x, y)
# Take a gradient step for each neuron in each layer
network = [[gradient_step(neuron, grad, -learning_rate)
for neuron, grad in zip(layer, layer_grad)]
for layer, layer_grad in zip(network, gradients)]
t.set_description(f"fizz buzz (loss: {epoch_loss:.2f})")
这将需要一段时间来训练，但最终损失应该开始稳定下来。
最后，我们准备解决我们最初的问题。我们还有一个问题。我们的网络将生成一个四维向量的数字，但我们想要一个单一的预测。我们将通过取argmax来做到这一点，这是最大值的索引：
def argmax(xs: list) -> int:
"""Returns the index of the largest value"""
return max(range(len(xs)), key=lambda i: xs[i])
assert argmax([0, -1]) == 0 # items[0] is largest
assert argmax([-1, 0]) == 1 # items[1] is largest
assert argmax([-1, 10, 5, 20, -3]) == 3 # items[3] is largest
现在我们终于可以解决“FizzBuzz”了：
num_correct = 0
for n in range(1, 101):
x = binary_encode(n)
predicted = argmax(feed_forward(network, x)[-1])
actual = argmax(fizz_buzz_encode(n))
labels = [str(n), "fizz", "buzz", "fizzbuzz"]
print(n, labels[predicted], labels[actual])
if predicted == actual:
num_correct += 1
print(num_correct, "/", 100)
对我来说，训练后的网络获得了 96/100 的正确率，远高于工程副总裁的招聘门槛。面对这些证据，他屈服了并将面试挑战改为“反转二叉树”。
进一步探索
-
继续阅读：第十九章将更详细地探讨这些主题。
-
我关于“Tensorflow 中的 Fizz Buzz”的博客帖子相当不错。
