Python 及经典统计模型回顾#
\(\hspace{1.5em}\) 在本章中，我们将简要回顾 Python 的基本命令以及基本的统计模型，包括线性回归模型以及逻辑回归模型。线性回归模型和逻辑回归模型在深度学习中扮演着至关重要的角色，它们不仅是机器学习领域的基础，也是构建和理解深度学习模型的重要基石。在深度学习中，线性回归可以被视为神经网络中的线性层，通过线性变换对输入数据进行初步处理。这种线性变换虽然是简单的，但它是构建复杂神经网络不可或缺的一部分。逻辑回归模型，则是一种专门用于处理二分类问题的模型。它通过线性组合和sigmoid激活函数，将输出结果压缩到0和1之间，表示样本属于正类的概率。通过堆叠多个线性层和非线性激活函数，深度学习模型能够捕捉到数据中的复杂模式和特征，从而实现高精度的预测和分类。逻辑回归模型可被认为是一种最简单的神经网络模型；在某种程度上讲，深度学习模型可被认为是对逻辑回归模型的推广。
Python回顾#
\(\hspace{1.5em}\)Python在人工智能领域发挥着至关重要的作用，其简洁的语法、丰富的库支持以及强大的社区生态为开发者提供了极大的便利。
无论是经典的统计学习、机器学习抑或是以深度学习为代表的人工智能，Python都以其出色的性能和广泛的应用场景成为实现这些功能的首选语言。
此外，Python还具有良好的跨平台性、动态类型与快速迭代的优势，以及面向对象与函数式编程的灵活性，进一步推动了其在人工智能领域的深入应用和发展。
Python还广泛应用于数据分析、Web开发、科学计算等多个领域。掌握Python可以让开发者更加高效地解决问题，提高开发效率，同时也为未来的职业发展和技术提升打下坚实的基础。关于Python基础命令的详细介绍，可参见Python Tutorial等。
在本节中，我们将对本课程中常用的一些Python数据结构及命令进行简单的回顾。
数据类型#
\(\hspace{1.5em}\) 数据在编程中无处不在，无论是进行简单的数学运算，还是处理复杂的数据分析，数据类型都是必须要掌握的基本内容。Python作为一门功能强大的编程语言，提供了多种不同的数据类型，主要包括整数（integer）、浮点数（floating point）、字符串（string）和布尔值（boolean）。
整型变量（integer）
\(\hspace{1.5em}\) 在Python中，整型变量是一种基本且常用的数据类型，专门用于表示整数。整型变量可以存储正整数、负整数以及零，是编程中处理数字数据的基础。Python的整型数据类型具有一个显著的特点，即其精度是任意的。这意味着，与一些具有固定位数限制（如32位或64位）的编程语言不同，Python的整型变量可以表示非常大或非常小的整数，而无需担心溢出或精度损失。Python会根据整型变量的大小自动分配足够的内存空间，以确保数据的准确性和完整性。因此，在使用Python进行编程时，我们可以放心地使用整型变量来处理各种整数运算，无需担心其范围限制。
1 age = 25 #定义整数变量
2 height_in_cm = 175 #定义整数变量
3
4 sum = age + height_in_cm #整数运算
5 print("Sum:", sum) #输出
Sum: 200
浮点型变量（float）
\(\hspace{1.5em}\) 在Python中，浮点型变量是一种用于表示带有小数的数字的数据类型。浮点型变量是Python内置的基本数值类型之一，它能够存储实数，包括正数、负数和零，且这些数可以包含小数点。浮点型变量在科学计算、数据分析以及任何需要精确小数处理的场景中都非常有用。Python程序能够处理各种实数运算，满足开发者在数值计算方面的需求。
1 price = 19.99 #定义浮点数变量
2 weight = 72.5 #定义浮点数变量
3
4 total_cost = price * weight #浮点数运算
5 print("Total Cost:", total_cost) # 输出（输出结果可能因浮点数精度问题略有不同）
Total Cost: 1449.2749999999999
\(\hspace{1.5em}\) 关于整型和浮点型数据的详细介绍，请参见Numeric Types。
字符串（string）
\(\hspace{1.5em}\) 在Python中，字符串是一种用于表示文本数据的数据类型。字符串由一系列字符组成，可以包含字母、数字、符号以及空格等，用于存储和操作文本信息。Python中的字符串是不可变的，这意味着一旦字符串被创建，它的内容就不能被改变。字符串在Python编程中非常常用，无论是处理用户输入、文件操作还是网络数据传输，都离不开字符串的处理。Python提供了丰富的字符串操作方法和函数，使得字符串的处理变得非常方便和高效。
1a = "Hello, " #定义字符串变量
2b = "World!" #定义字符串变量
3
4c = a+b #字符串拼接
5print(c)
Hello, World!
\(\hspace{1.5em}\) 关于字符串的详细介绍，请参见Text Sequence Type。
布尔值（boolean）
\(\hspace{1.5em}\) 在Python中，布尔值是一种用于表示真（True）或假（False）两种状态的数据类型。布尔值是逻辑运算的基础，广泛用于条件判断、循环控制以及逻辑表达式中。在Python里，布尔值不仅可以直接使用True和False来表示，还可以通过逻辑运算符（如and、or、not）对布尔值或其他可转化为布尔值的表达式进行组合和运算。布尔值的简洁明了和逻辑运算的强大功能，使得Python程序在控制流程和决策判断方面更加灵活和高效。
1is_student = True #定义布尔值变量
2has_ticket = False #定义布尔值变量
3
4
5can_enter = is_student or has_ticket #布尔值或运算
6print("Can Enter:", can_enter) # 输出运算结果
7
8if is_student: #在 if 语句中使用布尔值
9 print("The student can enter.")
10else:
11 print("The person needs a ticket to enter.")
Can Enter: True
The student can enter.
\(\hspace{1.5em}\) 关于布尔值的详细介绍，请参见Boolean Type。
基本数据结构#
\(\hspace{1.5em}\)Python提供了多种基本数据结构，以满足不同编程需求。列表（list）是一种有序、可变的数据结构，它可以存储任意类型的对象，并且支持索引、切片等多种操作。列表非常灵活，可以动态地调整大小，允许包含重复的元素，是处理一组有序数据时的首选。在本课程中，我们将基于列表结构存储并更新深度学习模型的参数。元组（tuple）与列表相似，但它是一种不可变的数据结构，一旦创建，其元素就不能被修改。元组常用于存储多个相关的值，由于其不可变性，它为数据提供了一种稳定的表示方式。我们在本课程中，主要利用元组作为函数的输出。集合（set）是另一种重要的数据结构，它无序且不重复，用于保存不重复的元素。集合中的元素必须是可哈希的，这使得集合能够进行高效的成员检测和集合运算，如交集、并集、差集等。集合在处理需要去重或进行集合运算的数据时非常有用。最后，字典（dictionary）是一种无序的由键：值对组成的集合，它用于存储和查找数据。字典中的键是唯一的，而值可以是任意类型。通过键，我们可以快速地查找对应的值，这使得字典在处理需要通过键快速访问数据的场景时非常高效。在本课程中，我们将深度学习模型的参数及其对应的值用字典表示并索引。
列表（list）
\(\hspace{1.5em}\)Python中的列表是一种非常基础且强大的数据结构，它是一种有序、可变的数据集合。列表可以存储任意类型的对象，包括数字、字符串、甚至是其他列表，这使得它非常灵活。列表中的元素按照插入的顺序进行排列，并且每个元素都有一个对应的索引，索引从0开始。这使得我们可以通过索引来快速访问、修改或删除列表中的元素。此外，列表还支持切片操作，我们可以通过指定起始和结束索引来获取列表的一个子序列。需要指出的是，在深度学习中，我们通常利用numpy或者torch产生列表；具体讨论请参见numpy库以及torch.Tensor的介绍。
1numbers = [1, 2, 3, 4, 5] #创建一个包含整数的列表
2print(numbers[2]) #列表索引
3
4fruits = ['apple', 'banana', 'cherry'] #创建一个包含字符串的列表
5print(fruits[:2]) #列表索引
6
7mixed_list = [1, 'hello', 3.14, True] #创建一个包含混合类型的列表
8print(mixed_list[1:3]) #列表索引
3
['apple', 'banana']
['hello', 3.14]
元组（tuple）
\(\hspace{1.5em}\) 元组是Python中的基本数据结构，其是由小括号括起来的若干元素组成。Python规定，元组的元素不能更改。在这里，我们不对元组进行详细的解释，只介绍其打包和拆包的功能。
1a = 1,2,3 #打包
2print(a)
3x1, x2, x3 = a #拆包
4print(x1)
5print(x2)
6print(x3)
(1, 2, 3)
1
2
3
\(\hspace{1.5em}\) 在上面的例子中，我们用命令a = 1,2,3表示一个打包过程，即我们将等式右端用逗号隔开的三个数字打包成一个元组(1,2,3)，再将此元组赋值给变量a。我们用命令x1, x2, x3 = a表示一个拆包过程，即我们将元组a拆成三个元素，分别赋值给变量x1, x2, x3。需要指出的是，有时候我们并不需要对拆包之后的元素都赋值给相应变量。例如，我们只关心元组a中的第三个值，而对前两个数值并不关心。此时，我们可以用符号_代替前两个变量；请参见下例。
1 a = 1,2,3 #打包
2 print(a)
3 _, _, x3 = a #拆包
4 print(x3)
(1, 2, 3)
3
\(\hspace{1.5em}\) 在本课程中，我们将利用元组的打包作为函数的输出，而利用拆包过程对函数的结果进行赋值。在下例中，我们编写一个函数用来计算一个数组x的样本均值和方差。这里，我们用到了Python中一个常见的库numpy，我们将在后面内容对其进行详细介绍。
1 import numpy as np
2 def statistics(x):
3 # x: 一个包含样本数据的数组
4 mean_x = np.mean(x) #均值
5 var_x = np.var(x) #方差
6 return mean_x, var_x #将均值和方差结果打包输出
7
8 x = np.arange(5)
9 print(x)
10 mean_x, var_x = statistics(x) #将函数结果拆包，并进行相应的赋值
11 print(mean_x)
12 print(var_x)
[0 1 2 3 4]
2.0
2.0
\(\hspace{1.5em}\) 关于列表和元组的详细介绍，请参见Sequence Types。
集合（set）
\(\hspace{1.5em}\)Python中的集合是一种无序、不重复的数据结构，它主要用于快速成员关系测试和消除重复元素。集合中的元素没有固定顺序，且每个元素都是唯一的，这保证了数据的互异性。集合支持多种操作，如添加元素、删除元素以及进行集合运算（如交集、并集、差集等）。由于集合内部采用哈希表实现，因此成员关系测试的操作非常高效，时间复杂度通常为\(O(1)\)。此外，集合还可以用于快速去除列表或其他可迭代对象中的重复元素，使得数据处理更加简洁高效。
1 my_set = {1, 2, 3, 2, 3, 4} #重复的元素会被自动去除
2 print(my_set)
3
4 my_set.add(5) #添加元素到集合中
5 print("添加元素", my_set) # 输出可能是: {1, 2, 3, 4, 5}（顺序可能不同）
6
7 my_set.remove(3) #删除元素
8 print("删除元素", my_set) # 输出可能是: {1, 2, 4, 5}
9
10 set1 = {1, 2, 3}
11 set2 = {3, 4, 5}
12 print("并集",set1.union(set2)) # 并集: {1, 2, 3, 4, 5}
13 print("交集",set1.intersection(set2)) # 交集: {3}
14 print("差集",set1.difference(set2)) # 差集: {1, 2}
{1, 2, 3, 4}
添加元素 {1, 2, 3, 4, 5}
删除元素 {1, 2, 4, 5}
并集 {1, 2, 3, 4, 5}
交集 {3}
差集 {1, 2}
\(\hspace{1.5em}\) 关于集合的详细介绍，请参见Set Types。
字典（dictionary）
\(\hspace{1.5em}\) 字典是Python中另一个重要的数据结构。一个字典由键：值对组成。在本节课程中，我们所用到的操作如下例所示。
1 dic = {"b1": 1, "W1": 2, "b2": 3, "W2": 4} #定义字典
2 print(dic)
3 print(dic["b1"]) #根据键"b1"，调取其对应的值
4 dic["b3"] = 5 #对字典新增加一个键：值对
5 print(dic)
6
7 i = 1
8 print(dic["W"+str(i)]) #通过字符串处理生成键值，并调取相应的值。
{'b1': 1, 'W1': 2, 'b2': 3, 'W2': 4}
1
{'b1': 1, 'W1': 2, 'b2': 3, 'W2': 4, 'b3': 5}
2
\(\hspace{1.5em}\) 在上例中，我们首先生成了一个包含四个键：值对的字典，并利用相应的键值调取了键值为b1对应的具体值。我们还通过命令dic["b3"] = 5增加了一个键：值对到字典中。需要指出的是，我们基于字符串处理的基本命令，利用dic["W"+str(i)]调取了键值为W1的具体数值。如果要遍历字典中的所有键值对，我们可以使用如下命令。
1 # 遍历键
2 for key in dic.keys():
3 print(key)
4
5 # 遍历值
6 for value in dic.values():
7 print(value)
8
9 # 遍历键值对
10 for key, value in dic.items():
11 print(key, value)
b1
W1
b2
W2
b3
1
2
3
4
5
b1 1
W1 2
b2 3
W2 4
b3 5
\(\hspace{1.5em}\) 在本课程中，我们将着重利用字典记录模型参数值，并对其进行迭代更新。关于字典的详细介绍，请参见Mapping Types。
习题#
编写函数
initial(n)，该程序的输入为一个正整数\(n\)，用于生成并返回一个字典。在该字典中，包含键b1,...,bn以及键W1,...,Wn；对于\(i=1,\ldots,n\)，键bi对应的值为0，而键Wi对应的值为一个\(2\times 2\)的矩阵，每个元素从标准正态分布随机生成。
1import numpy as np
2
3# 下面的程序对完成本题有帮助
4dict = {}
5dict['b'+str(1)] = 0
6dict['W'+str(1)] = np.random.normal(size = (2,2))
7
8def initial(n):
9 # 完成该函数
假设我们有一个包含多个字典的列表
student_info，每个字典代表一个学生的信息，包含以下键：'name'（名字）、'age'（年龄）和'scores'（一个包含各科成绩的字典）。请编写一个Python程序ave_score(student_info, name)，该程序的输入为字典列表student_info以及一个学生的名字name。若该字典列表中包含这个学生的成绩，则输出该学生成绩的平均分；否则，输出字符串"查无此人"。
1student_info = [
2 {"name": "Alice", "age": 20, "scores": {"Math": 90, "English": 85, "Science": 82}},
3 {"name": "Bob ", "age": 22, "scores": {"Math": 88,"Physics": 80,"Chemistry": 84}},
4 {"name": "Charlie ", "age": 19, "scores": {"Math": 90,"Biology": 96,"History": 80}},
5]
6
7def ave_score(student_info, name):
8 # 完成该函数
numpy 简介#
\(\hspace{1.5em}\)numpy是Python中用于科学计算的最为核心的基础库之一，
其重要性体现在提供了高性能的多维数组对象及对这些数组的操作。它简化了数值计算的复杂度，是数据分析、机器学习、科学计算等领域不可或缺的工具，极大地提高了开发效率和计算速度。在本节中，我们将简单讨论其广播机制（broadcasting）以及在本课程会使用的一些基本命令。
\(\hspace{1.5em}\) 在Python中，加载库（也称为模块）是一个简单而常见的操作，它允许你访问库提供的函数、类和变量等。要加载一个库，我们需要使用import语句。
通常，我们将用np简化numpy。
1 import numpy as np #加载numpy库
numpy列表#
\(\hspace{1.5em}\) 列表是Python中一个重要的数据结构，其可以用来保存一维数组、二维矩阵，还可以用来保存三维图片像素，甚至是多维度张量。numpy为列表操作提供了极大的便利，进而为Python插上了科学计算的翅膀。列表的生成多种多样，我们通过实例介绍本课程中会用到的方式。
1 x1 = np.arange(4) #生成列表 [0，1，2，3]
2 print(x1)
3 print("*"*6)
4 x2 = np.zeros((4,2)) #生成一个4x2的零矩阵
5 print(x2)
6 print("*"*6)
7 x3 = np.ones((4,2)) #生成一个4x2的矩阵，其中矩阵元素都是1
8 print(x3)
9 print("*"*6)
10 x4 = np.ones_like(x2) #生成一个同x2具有相同维度的矩阵，其中矩阵元素都是1
11 print(x4)
12 print("*"*6)
13 x7 = np.random.normal(size = (4,2)) #生成一个4x2的随机矩阵，其中矩阵元素服从标准正态分布
14 print(x7)
15 print("*"*6)
[0 1 2 3]
******
[[0. 0.]
[0. 0.]
[0. 0.]
[0. 0.]]
******
[[1. 1.]
[1. 1.]
[1. 1.]
[1. 1.]]
******
[[1. 1.]
[1. 1.]
[1. 1.]
[1. 1.]]
******
[[-0.19768302 1.6689484 ]
[ 0.6617977 0.78745647]
[-0.0136145 0.66776744]
[-0.63459065 0.50157825]]
******
\(\hspace{1.5em}\) 当给定一个数组时，我们通常需要对其特定位置的元素进行操作。下面的例子中，我们给出了几个常用的命令。
1 x1 = np.random.normal(size = (4,2)) #生成一个4x2的随机矩阵，其中矩阵元素服从标准正态分布
2 print(x1)
3 print("*"*6)
4 print(x1[0,:]) #取出x1的第一行
5 print("*"*6)
6 print(x1[0]) #同样取出x1的第一行
7 print("*"*6)
8 print(x1[:,0]) #取出x1的第一列
9 print("*"*6)
10 print(x1[0,1]) #取出x1中第一行第二列元素
[[ 0.08505252 -1.19209771]
[-0.55964263 -0.90427558]
[-1.18035351 -1.27827524]
[-0.70215563 0.50339898]]
******
[ 0.08505252 -1.19209771]
******
[ 0.08505252 -1.19209771]
******
[ 0.08505252 -0.55964263 -1.18035351 -0.70215563]
******
-1.192097711896629
广播机制#
\(\hspace{1.5em}\) 在Python中，广播机制（Broadcasting）是一种强大的功能，特别是在处理numpy数组时，它能够极大地简化代码并提高计算效率。广播机制允许numpy在进行数组运算时，对形状不同的数组进行操作。通过广播机制，numpy会自动将较小的数组“扩展”到较大数组的形状上，以便进行元素间的运算。这种机制不需要显式地复制数据，而是通过调整数组的形状来实现。本节内容主要参考 Broadcasting 。
\(\hspace{1.5em}\) 对于两个数组a以及 b而言，numpy的代数运算基本是基于这两个数组具有相同形状，即a.shape=b.shape。此时，numpy是基于元素之间的运算进行的，而这种运算正是我们往往在编程中需要的。请参见如下例子。
1a = np.array([1.0, 2.0, 3.0])
2b = np.array([2.0, 2.0, 2.0])
3print(a.shape)
4print(b.shape)
5print(a+b)
(3,)
(3,)
[3. 4. 5.]
\(\hspace{1.5em}\) 然而，要求两个数组具有相同维度是具有明显缺陷的。例如，有时候我们只希望对以上数组a中的每个元素加上相同的因子c。在上例中，我们将数组a的每一个元素加上了2，但利用一个长度为3的数组b实现这个目的，显然是消耗内存的。为了达到这个目的，我们希望利用简化的命令a+c，其中c为一个标量，而不再是与数组a具有相同形状的数组b。numpy中的广播机制可以认为就是对具有相同形状这一假设的放松。例如，我们可以利用numpy轻松实现和上例相同的效果。
1 a = np.array([1.0, 2.0, 3.0])
2 c = 2.0
3 print(a+c)
[3. 4. 5.]
\(\hspace{1.5em}\) 结合以上两个例子，我们可以直观地认为在第二个例子中，标量c的值被拉伸成了一个与数组a具有相同维度的数组；换句话说，我们通过复制c的值若干次，形成了一个与数组a相同维度的新数组，然后进行代数运算。当然，这种基于拉伸或者复制的想法只是理解广播机制的直观，numpy利用了更加高效的方法避免了由于复制同一数值而带来的内存和计算效率的损失。尽管我们用了两种方法实现将数组a的每个元素加上相同的元素，但第二种方法的计算效率更高、对内存要求更低。下图展示了第二种计算的直观理解。
\(\hspace{1.5em}\) 下面，我们将更一般地讨论基于两个数组a以及b间代数运算的广播机制。numpy首先从右向左检验两个数组的维度，即a.shape以及b.shape。当如下两个条件之一满足时，我们说数组 a以及b在某个维度上是兼容的：
两个维数相等，
其中一个数组的维数等于1。
\(\hspace{1.5em}\) 例如，若a.shape = (3,1)以及b.shape = (1,3)，则这两个数组在两个维度上均兼容。具体得讲，在最右边的维度上，数组a的维数为1，而数组b的维数为3；因此，在最右边的维度上，满足上面的第二条。我们可以简单看出，在最左边的维度上，两个数组也是兼容的。再例如，若a.shape = (3,4)但b.shape = (4,3)，则这两个数组在两个维度上均不兼容。如果两个数组存在一个不兼容的维度，则进行相关的二元运算时，Python将会报错：ValueError: operands could not be broadcast together。
\(\hspace{1.5em}\) 通过数组与标量间的二元运算的例子我们可知，numpy的广播机制不一定要求两个数组具有相同的维度。当两个数组a以及b具有不同的维度时，numpy首先将较小维度的数组从左端用1补充维度，直至与较长数组的维度相同。如果补全后，两个数组维度兼容，我们将在每个维度上，比较两个数组的维数，取较大者作为最终运算结果的维度。例如，假设数组a的维度为\((256,256,3)\)数组b的维度为\((3,)\)。此时，数组a的维度为3，而数组b的维度为1。直观上理解，numpy首先将通过左侧补1的方法将数组b的维度变为\((1,1,3)\)。然后，通过对比每个维度上两个数组维数的较大值，运算结果的维度为\((256,256,3)\)。再例如，数组a的维度为\((8，1，6，1)\)数组b的维度为\((7,1,5)\)。当两个数组进行二元运算时，我们首先将较小维度的数组左侧补1，即将数组b的维度变为\((1,7,1,5)\)，通过比较我们可知两个数组维度是兼容的。所以，计算结果的维度为(8,7,6,5)。
\(\hspace{1.5em}\) 当我们已经知道结果的维度时，我们还需要知道numpy是如何通过广播机制得到最终结果的。这个过程与数组与标量间的二元运算例子中展示的过程基本一致。直观上讲，numpy在维度为1的方向上复制相应的结果，当两个数组具有相同形状时，我们即可通过元素间的二元运算得到最终计算结果。我们通过两个例子展示这个过程。
1 a = np.array([[2.0, 3.0, 4.0],[3.0, 4.0, 5.0], [4.0, 5.0, 6.0]])
2 b = np.array([3.0, 4.0, 5.0])
3 print(a.shape)
4 print(b.shape)
5 print(a+b)
(3, 3)
(3,)
[[ 5. 7. 9.]
[ 6. 8. 10.]
[ 7. 9. 11.]]
\(\hspace{1.5em}\) 在以上这个例子中，数组a的维度为\((3,3)\)但数组b的维度为\((3,)\)。我们首先将数组b的维度从左段补1扩充为\((1,3)\)。通过比较扩充后的两个数组可知，他们的维度在各个位置上均是兼容的。通过观察可知，在左边第一个维度上，数组a的维数为3，而数组b的维数为1。因此，直观上理解，numpy将数组b沿着第一个维度，将该数组复制两次，得到一个与数组a相同的数组。此时，我们便可在每个位置上对两个元素进行相加，得到最终结果了。需要指出的是，numpy在运算时，不会复制数组以节省内容和提高计算效率。下图展示了本例子计算过程的直观理解。
1 a = np.array([[1.0],[2.0], [3.0]])
2 b = np.array([3.0, 4.0, 5.0])
3 print(a.shape)
4 print(b.shape)
5 print(a+b)
(3, 1)
(3,)
[[4. 5. 6.]
[5. 6. 7.]
[6. 7. 8.]]
\(\hspace{1.5em}\) 在上例中，数组a的维度为\((3,1)\)但数组b的维度为\((3,)\)。首先，我们将数组b的维度从左段补1扩充为\((1,3)\)。通过比较，我们可知两个数组在各个位置上的维数是兼容的。但通过比较我们可知，我们需要同时复制两个数组，只不过数组a需要沿着第二个维度复制，而数组b需要沿着第一个维度复制。最终结果的维度为\((3,3)\)。下图展示了本例子计算过程的直观理解。
\(\hspace{1.5em}\) 为了避免潜在的错误，当我们利用广播机制对两个数组的二元运算时，我们需要保证这两个二元数组的维度相同。即我们要努力避免对维度不同的两个数组之间进行二元运算。比较下面一段代码中数组b和数组c维度上的差异。在具体的编程中，尽管计算结果相同，但我们应当选择a+c而不是a+b。
1 a = np.array([[1.0],[2.0], [3.0]])
2 b = np.array([3.0, 4.0, 5.0])
3 c = np.array([3.0, 4.0, 5.0]).reshape((1,3))
4 print("数组a的维度为："+str(a.shape))
5 print("数组b的维度为："+str(b.shape))
6 print("数组c的维度为："+str(c.shape))
7 print("a+b的结果为：")
8 print(a+b)
9 print("a+c的结果为：")
10 print(a+c)
数组a的维度为：(3, 1)
数组b的维度为：(3,)
数组c的维度为：(1, 3)
a+b的结果为：
[[4. 5. 6.]
[5. 6. 7.]
[6. 7. 8.]]
a+c的结果为：
[[4. 5. 6.]
[5. 6. 7.]
[6. 7. 8.]]
常用命令#
\(\hspace{1.5em}\)numpy除了以上介绍的基于二元运算的广播机制外，还有很多高效的基于元素的运算。例如，我们如果希望对一个给定数组a的每一个元素求指数，我们可简单用np.exp()命令。下面的例子中，展现了几个常见的numpy命令。
1 a = np.array([1.0, 2.0, 3.0])
2 print(np.exp(a)) #求指数
3 print(a**2) #求平方
4 print(np.sqrt(a)) #开根号
5 print(np.sin(a)) #求正弦
[ 2.71828183 7.3890561 20.08553692]
[1. 4. 9.]
[1. 1.41421356 1.73205081]
[0.84147098 0.90929743 0.14112001]
\(\hspace{1.5em}\) 我们介绍两个命令，包括np.mean()以及np.var()。这两个命令的分别用来求解一个数组a的均值和方差。下例中，我们从标准正态分布中产生了10个随机数，并计算其对应的样本均值和方差。
1 np.random.seed(1)
2 a = np.random.normal(size=10)
3 print("数组a为：")
4 print(a)
5 print("数组a的均值为：")
6 print(np.mean(a)) #求均值
7 print("数组a的方差为：")
8 print(np.var(a)) #求均值
数组a为：
[ 1.62434536 -0.61175641 -0.52817175 -1.07296862 0.86540763 -2.3015387
1.74481176 -0.7612069 0.3190391 -0.24937038]
数组a的均值为：
-0.09714089080609986
数组a的方差为：
1.4182393613078983
\(\hspace{1.5em}\) 在实际应用中，我们通常将训练集表示成矩阵的形式。例如，对于训练集\(\{\bx_i:i=1,\ldots,n\}\)而言，我们通常记\(\bX=(\bx_1,\ldots,\bx_n)\trans\in\mathbb{R}^{n\times d}\)，其中\(n\)表示训练集的规模，\(d\)表示特征的维数。在实际操作中，我们往往需要对训练集中每个特征求均值或者方差。也就是说，我们需要对\(\bX\)中的每一列分别求均值和方差。为了达到这个目的，我们可以使用np.mean()或者 np.var()中的参数axis=0来实现。在下面的例子中，我们生成一个规模为\((1000,6)\)的二维数组X，并介绍如何用np.mean()或者 np.var()求解每一列的均值和方差。
1 np.random.seed(1)
2 X = np.random.normal(size=(1000, 6))
3 print("数组X的列均值为：")
4 print(np.mean(X, axis=0)) #求均值
5 print("数组X的列均值为：")
6 print(np.var(X, axis=0)) #求均值
数组X的列均值为：
[ 0.01679263 -0.01111949 0.01457063 0.02767384 0.03295546 -0.00981674]
数组X的列均值为：
[1.05076426 1.07445554 0.96374936 0.94109343 1.02405612 0.95722715]
\(\hspace{1.5em}\) 接下来，我们对上例中的结果做两方面的具体解释。首先，我们利用了参数axis=0来实现了列均值以及列方差的求解。我们也可以利用axis=1来实现对行均值和方差的求解。axis的数值设定提示了numpy对第几个维度进行计算。例如，当我们使用np.mean(X, axis=0)时，numpy将沿着数组X的第一个维度（即矩阵的行）求均值。类似地，当我们使用np.mean(X, axis=1)时，numpy将沿着数组X的第二个维度（即矩阵的列）求均值。下面例子中，X的元素较为简单，请大家比较两个均值求解的差异。
1 X = np.arange(4).reshape(2,2)
2 print(X)
3 print(np.mean(X, axis=0))
4 print(np.mean(X, axis=1))
[[0 1]
[2 3]]
[1. 2.]
[0.5 2.5]
\(\hspace{1.5em}\) 将以上的讨论进行拓展，当我们有一个3维数组X时，我们可以分别使用np.mean(X, axis=0)、np.mean(X, axis=1)和np.mean(X, axis=2)分别沿着数组X的第一个维度、第二个维度和第三个维度求解均值。
1 X = np.arange(8).reshape(2,2,2)
2 print(X)
3 print("下面是沿着三个维度的均值求解结果：")
4 print(np.mean(X, axis=0))
5 print(np.mean(X, axis=1))
6 print(np.mean(X, axis=2))
[[[0 1]
[2 3]]
[[4 5]
[6 7]]]
下面是沿着三个维度的均值求解结果：
[[2. 3.]
[4. 5.]]
[[1. 2.]
[5. 6.]]
[[0.5 2.5]
[4.5 6.5]]
\(\hspace{1.5em}\) 在上例中，np.mean(X,axis =0)结果中，对应与第一行第一列的元素为X[0,0,0]和
X[1,0,0]的均值。而np.mean(X,axis =1)结果中，对应与第一行第一列的元素为X[0,0,0]和
X[0,1,0]的均值。
\(\hspace{1.5em}\) 通过数据生成过程我们可知数组\(X\)为二维，但计算得到的列均值和列方差均为一维数组；即X.shape = (1000,6)，但求出的列均值和列方差的规模为(6,)。在进行建模时，我们有时候需要对一个二维数组的列进行归一化，这就要求我们将数组X的每一行减掉列均值，再除列标准差（方差的开根号结果）。根据上面的讨论，我们需要将列均值和列方差的规模变为\((1,6)\)。即求解列均值或者列方差后的结果规模也应该为一个二元数组，实现这一点对复杂编程具有非常重要的意义。另一个理解这个重要意义的角度为，在上例中，np.mean(X,axis=0)、np.mean(X,axis=1)以及np.mean(X,axis=2)对应的结果均为一个规模为\((2,2)\)的二元数组。这就为后续造成了困扰，即我们可能不知道该用哪个结果。因此，保证结果也为三维数组对后续编程也很重要。为了达到这个目的，我们可以简单使用参数keepdims=True来实现。下面我们用求解列均值来比较两个不同的结算结果，请注意比较两个结果的维度差异。
1 np.random.seed(1)
2 X = np.random.normal(size=(1000, 6))
3 mean1 = np.mean(X,axis=0)
4 mean2 = np.mean(X,axis=0,keepdims =True)
5 print(mean1)
6 print(mean2)
7 print(mean1.shape)
8 print(mean2.shape)
[ 0.01679263 -0.01111949 0.01457063 0.02767384 0.03295546 -0.00981674]
[[ 0.01679263 -0.01111949 0.01457063 0.02767384 0.03295546 -0.00981674]]
(6,)
(1, 6)
\(\hspace{1.5em}\) 下面，请看三维数组的情况。
1 X = np.arange(8).reshape(2,2,2)
2 print(X)
3 print("下面是沿着三个维度的均值求解结果：")
4 mean0 = np.mean(X,axis=0, keepdims=True)
5 mean1 = np.mean(X,axis=1, keepdims=True)
6 mean2 = np.std(X,axis=2, keepdims=True)
7 print(mean0)
8 print(mean1)
9 print(mean2)
10 print("下面是沿着三个维度的均值求解结果的维度：")
11 print(mean0.shape)
12 print(mean1.shape)
13 print(mean2.shape)
[[[0 1]
[2 3]]
[[4 5]
[6 7]]]
下面是沿着三个维度的均值求解结果：
[[[2. 3.]
[4. 5.]]]
[[[1. 2.]]
[[5. 6.]]]
[[[0.5]
[0.5]]
[[0.5]
[0.5]]]
下面是沿着三个维度的均值求解结果的维度：
(1, 2, 2)
(2, 1, 2)
(2, 2, 1)
\(\hspace{1.5em}\)Python中的矩阵往往是一个二维numpy数组，其运算是通过numpy库来实现的。numpy提供了丰富的矩阵操作函数和方法，使得矩阵运算变得简单而高效。矩阵的基本运算包括加法、减法、乘法、转置、求逆、行列式计算等。例如，使用+和-运算符可以进行矩阵的加法和减法运算；使用numpy.dot()函数可以进行矩阵乘法运算；使用.T属性可以获取矩阵的转置；使用np.linalg.inv()函数可以计算矩阵的逆（如果矩阵可逆）；使用np.linalg.det()函数可以计算矩阵的行列式。我们简单介绍本课程中常用的三种基本的矩阵运算，包括 Hadamard 乘法、标准矩阵乘法以及矩阵转置。
\(\hspace{1.5em}\)numpy中的Hadamard乘法，也称为逐元素乘法或Hadamard乘积，是一种矩阵运算方式。它不同于传统的矩阵乘法，而是对两个矩阵对应位置上的元素进行相乘，得到一个新的矩阵。具体来说，如果矩阵\(\bA=(a_{ij})\)和矩阵\(\bB=(b_{ij})\)的形状相同，即它们都是\(m×n\)的矩阵，那么\(\bA\)和\(\bB\)的Hadamard乘积\(\bC= \bA\circ\bB\)也是一个\(m×n\)的矩阵，其中\(\bC\)的每个元素等于\(\bA\)和\(\bB\)相应位置元素的乘积。具体地，对于\(i=1,\ldots,m\)以及\(j=1,\ldots,n\)，\(c_{ij}=a_{ij}\times b_{ij}\)。在numpy中，实现Hadamard乘法与广播机制相同，我们只需要在保证两个矩阵规模相同时，用A*B即可。
\(\hspace{1.5em}\) 标准矩阵乘法是指按照矩阵乘法的规则，第一个矩阵的行与第二个矩阵的列对应元素相乘后求和，得到结果矩阵的对应元素。在numpy中，这可以通过@操作符或np.matmul()函数来实现。
1 A = np.array([[1, 2], [3, 4]])
2 B = np.array([[5, 6], [7, 8]])
3
4 C = A @ B #使用 @ 操作符进行标准矩阵乘法
5 D = np.matmul(A, B) #或者使用 np.matmul() 函数
6
7 print(C)
8 print(D)
[[19 22]
[43 50]]
[[19 22]
[43 50]]
\(\hspace{1.5em}\) 矩阵转置是指将矩阵的行和列互换，即原矩阵的第\(i\)行变成转置矩阵的第\(i\)列，原矩阵的第\(j\)列变成转置矩阵的第\(j\)行。在numpy中，矩阵转置可以通过数组的.T属性或np.transpose()函数来实现。
1 A = np.array([[1, 2, 3], [4, 5, 6]])
2
3 B = A.T
4 C = np.transpose(A)
5
6 print(B)
7 print(C)
[[1 4]
[2 5]
[3 6]]
[[1 4]
[2 5]
[3 6]]
习题#
利用命令
np.arange以及reshape，写出生成如下二维数组X的命令。
X = np.array([[1,2,3],[4,5,6],[7,8,9]])
print(X)
#在下面用 np.arange 以及 .reshape() 函数写出生成以上数组的命令
X =
编写函数
studentize(X)，该程序的输入为一个二维numpy数组X，用于对输入数组的每一列进行归一化处理，即输出结果的每一列均值为0，方差为1。具体地，我们需将X中每个元素减到对应位置的列均值，再除掉对应位置的列标准差 （可用np.std()实现标准差求解）。该函数输出为一个与X规模相同的归一化的二维数组。
np.random.seed(1)
X = np.random.normal(size=(1000, 6))
def studentize(X):
# 完成该函数
编写函数
BatchNormalization(X, beta, gamma)，该程序的输入为一个二维numpy数组X，beta.shape = (1,d)，gamma.shape = (1,d)，其中d=X.shape[1]。该函数需要实现的目的为：将二维数组X的列均值变为beta，同时将其列标准差变为gamma。该函数输出为一个与X规模相同的归一化的二维数组。
np.random.seed(1)
X = np.random.normal(size=(1000, 6))
def BatchNormalization(X):
# 完成该函数
简单线性回归#
\(\hspace{1.5em}\) 在本节中，我们将依托简单线性回归模型，介绍如何对一个统计模型进行蒙特卡洛模拟，并对模型参数的估计量进行统计推断。在这个过程中，我们将分析步骤拆分成若干独立的函数，并利用numpy的基本命令实现编程；关于本节中使用的部分命令，可参见上一节中的相关介绍。
\(\hspace{1.5em}\) 在正式讨论简单线性回归模型之前，我们需要简单讨论 随机种子 这一基本概念。模拟研究与真实数据分析的主要区别在于，模拟研究的训练集是基于预先设定好的模型产生的。也就是说，我们从该预先设定的模型中随机生成训练集。为了保证模拟结果的可重复性，我们需要需要保证不同研究人员所产生的训练集是一致的，而这一点可通过控制随机种子来实现。换言之，在控制随机种子后，其他人可以生成与我们相同的训练集，进而得到和我们一致的分析结果。因此，通过设定随机种子，我们的分析结果就可以被复现了。当不指定随机种子时，不同研究人员所产生的训练集往往不同，分析结果也有所差异。
\(\hspace{1.5em}\) 下面这段代码中，我们利用np.random.normal()产生两个随机数。由于我们并没有预先设定随机种子，所产生的两个随机数并不相同。
import numpy as np #加载库
a = np.random.normal(1)
b = np.random.normal(1)
print(a,b)
-1.8510695426019397 1.8231234224707915
\(\hspace{1.5em}\) 接下来，我们利用np.random.seed()设定同一个随机种子再用相同的命令产生随机数。对比以上两组结果可知，当我们设定了同一随机种子后，用相同命令所产生的随机数是相同的。设定随机种子对模拟数据分析结果的复现性非常重要。
np.random.seed(1234)
a = np.random.normal(1)
np.random.seed(1234)
b = np.random.normal(1)
print(a,b)
1.471435163732493 1.471435163732493
备注
\(\hspace{1.5em}\) 在使用 Python 进行模拟实验时，我们可能会用到其他的库产生随机数，例如 random 库。对于 random 库而言，我们则需要使用 random.seed() 设定随机种子。
生成训练集#
真实模型假设#
\(\hspace{1.5em}\) 在本节中，我们考虑以下线性回归模型：
其中，\(b_0 = w_0 = 1\) 为模型真实参数，\(x_i \sim N(2, 2^2)\) 为特征（feature），\(\epsilon_i \sim N(0, 1)\) 为白噪声以及 \(n\) 为样本量。在该模拟实验中，我们要根据模型(1)生成训练集。我们称 \(y_i\) 为第 \(i\) 个实例 (example) 的标签（label）。需要指出的是，\(b_0\) 和 \(w_0\) 称为产生训练集的真实模型 (1) 的参数真值，这两个数值只用于产生模拟所需要的训练集。基于所产生的训练集，我们提出带有参数的统计模型，并对所提出模型的参数进行估计。在传统统计分析中，我们通常利用下角标`0`代表产生训练集的模型的真实参数值，而用没有下角标的字母表示基于样本所提出的模型的相关参数。需要指出的是，我们根据训练集所提出的统计模型往往和产生数据的模型不同，但我们在本节不讨论这个问题及其后果。
训练集的生成#
\(\hspace{1.5em}\) 接下来，我们将编写用于产生训练集的Python函数，这里将用到np.random.normal()函数。
def train_data_generation(n, rn):
# n: 样本量
# rn: 随机种子
np.random.seed(rn) #设置随机种子为rn
x = np.random.normal(2, 2, (n,1)) #从均值为2、标准差为2的正态分布中生成大小为nx1的特征向量
epsilon = np.random.normal(0, 1, (n,1)) #从均值为0、标准差为1的正态分布中生成大小为nx1的随机误差
y = 1 + x + epsilon #通过(1)式生成y
return x, y #打包输出相关结果
备注
\(\hspace{1.5em}\) 上面一段代码定义了一个名为 train_data_generation 的函数。该函数的输入由样本量 n 以及随机种子 rn 两部分组成，其输出为一个由 x 和 y 组成的元组。在函数内部，我们利用注释解释了函数的输入以及相应的代码。我们希望各位同学在将来自己编写函数时，也能养成用简单的注释解释代码的习惯。此外，在函数编写时，我们用的是有意义的字符串表示相应的结果。例如，我们用epsilon表示误差项，用x以及y表示特征和标签。我们用到的符号与 真实模型 中的符号相同。我们也希望同学们将来在编写代码时，也能够用具有意义的字符串，以便让代码具有较强的可读性。
\(\hspace{1.5em}\) 下面，我们用 matplotlib.pyplot 库中的命令对训练集数据进行可视化。
import matplotlib.pyplot as plt #我们通常将matplotlib.pyplot简单记为plt
x, y = train_data_generation(1000, 100)
plt.scatter(x,y)
plt.xlabel("Feature")
plt.ylabel("Label")
plt.show()
\(\hspace{1.5em}\) 通过观察，我们发现标签 \(y\) 和特征 \(x\) 之间呈现线性关系。因此，我们选择如下统计模型对训练集进行拟合：
其中，\(b\) 和 \(w\) 为模型参数。
参数估计#
\(\hspace{1.5em}\) 记 \(\tilde{\bx}_i = (1,x_i)\trans\), \(\bX=(\tilde{\bx}_i,\ldots,\tilde{\bx}_n)\trans\in\mathbb{R}^{n\times2}\) 以及 \(\bY=(y_1,\ldots,y_n)\trans\)。利用正则方程可知，模型(2) 中参数的估计为
我们可证明，在真实模型 (1) 的假设下，(3) 中的估计量 \((\hat{b}_0,\hat{w}_0)\) 是无偏的。
\(\hspace{1.5em}\) 第一个等式对应于基于循环运算的“求和法”，而第二个等式对应于向量化。作为本次的编程训练，我们将实现两种方法，并比较两种方法的计算效率。
备注
\(\hspace{1.5em}\) 向量化在数据处理和算法实现中具有重要意义。通过将数据表示为向量或矩阵形式，向量化允许使用线性代数中的矩阵运算，这些运算通常比传统的循环运算更快，特别是在处理大型数据集时。这种表示方式不仅提高了计算效率，还简化了代码实现，因为可以使用统一的数学运算来处理数据，减少了需要编写的代码量。此外，向量化还便于利用现代计算机体系结构的并行计算能力，如GPU，从而进一步加速数据处理和算法执行。因此，向量化是数据科学和机器学习领域中一种基础且关键的数据处理手段。我们在接下来的课程中，将详细介绍该方法在不同深度学习框架下的应用。
基于for循环的求和法#
\(\hspace{1.5em}\) 我们将利用(3)中第一部分编写参数估计函数，并基于此复习np.concatenate()，np.linalg.inv()以及np.flatten()三个函数的用法。此外，我们还将复习矩阵乘法和转置等相关运算。
def estimation_summation(x,y):
# x: 长度为n的特征向量
# y: 长度为n的标签向量
n = len(x) #获取样本量n的大小
aug_x = np.concatenate((np.ones_like(x),x), axis = 1) #将截距项1加入到特征中，得到增广特征矩阵，记为aug_x
xx = np.zeros((2,2)) #用2X2的零矩阵初始化'xx'
xy = np.zeros((2,1)) #用2X1的零矩阵初始化'xy'
for i in range(n): #根据（3）的第一个方程计算'xx'和'xy'的求和
example_i = aug_x[i,:].reshape((2,1))
xx += example_i @ example_i.transpose()
xy += example_i * y[i]
xx_inv = np.linalg.inv(xx) #计算'xx'的逆矩阵
par_est = xx_inv @ xy #根据（3）的第一个方程得到估计量
return par_est.flatten()
向量化#
\(\hspace{1.5em}\) 接下来，我们将基于向量化，即(3)中第二部分，求解参数的估计量。在本部分，我们也将用到基于for循环的求和法所提到的基本运算。
def estimation_vectorization(x,y):
# x: 长度为n的特征向量
# y: 长度为n的标签向量
aug_x = np.concatenate((np.ones_like(x),x), axis = 1) #将截距项1加入到特征中，得到增广特征矩阵，记为aug_x
xx = aug_x.transpose() @ aug_x #通过X^TX获得'xx'，其中X是步骤2得到的增广矩阵
xy = aug_x.transpose() @ y #通过X^TY获得'xy'，其中Y的形状为nX1
xx_inv = np.linalg.inv(xx) #计算'xx'的逆矩阵
par_est = xx_inv @ xy #根据(3)的第二个方程获得估计器
return par_est.flatten()
\(\hspace{1.5em}\) 下面，我们将调用以上两个函数，验证两种方法的估计值是否一致。
print(estimation_summation(x,y))
print(estimation_vectorization(x,y))
print('以上两个结果应该一致')
print('两个参数的真实值为1和1。')
[1.01322026 0.99776934]
[1.01322026 0.99776934]
以上两个结果应该一致
两个参数的真实值为1和1。
备注
\(\hspace{1.5em}\) 尽管两个方法的估计值相同，但其计算效率差异较大。我们将在计算效率一节中对两种方法进行比较。
\(\hspace{1.5em}\) 在传统的统计模型分析中，统计推断是一个重要的内容。但该部分内容涉及到较多的数理统计相关内容，基础较为薄弱的同学可以跳过下面的内容。
统计推断
\(\hspace{1.5em}\) 通常，我们会进行1,000次蒙特卡洛模拟。也就是说，我们根据真实模型 (1) 独立地生成大小为 \(n\) 的训练集1,000次。然后，对于每个生成的训练集，我们获得模型参数估计值 \((\hat{b}, \hat{w})\)。接着，我们将它们汇总在一起，用以计算其估计量的偏差、方差和均方误差（MSE）。如果我们有多种估计方法，我们可以比较它们在偏差、方差和均方误差方面的表现，用以衡量不同模型的优劣。尽管我们有两种不同的方法来估计模型参数，但它们的结果应该是相同的。因此，我们只考虑向量化方法。
\(\hspace{1.5em}\) 我们首先产生1,000个训练集，并将对应的参数估计值保存在一个 \(1000\times2\) 的矩阵中，其中每一行代表基于对应训练集得到的参数估计值。
n = 1000 #样本量
par_result = np.zeros((1000,2)) #生成1000x2的零矩阵，用以保存参数估计值
for i in range(1000):
x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i
par_result[i,:] = estimation_vectorization(x,y)
\(\hspace{1.5em}\) 接下来我们计算1,000次蒙特卡洛模拟估计值的偏差、方差以及均方误差。此处，我们将用到上节所复习的np.mean()和np.var()两个函数。
bias_est = np.mean(par_result,axis = 0)-1 #由于参数真值均为1，故我们减掉1，以得到估计量的蒙特卡洛偏差
var_est = np.var(par_result,axis = 0)
mse_est = bias_est**2 + var_est
print('参数b_0和w_0的估计偏差为：')
print(bias_est)
print("*"*6)
print('参数b_0和w_0的估计方差为：')
print(var_est)
print("*"*6)
print('参数b_0和w_0的MSE为：')
print(mse_est)
参数b_0和w_0的估计偏差为：
[5.59456735e-05 2.83427856e-04]
******
参数b_0和w_0的估计方差为：
[0.00200172 0.0002554 ]
******
参数b_0和w_0的MSE为：
[0.00200173 0.00025548]
\(\hspace{1.5em}\) 请注意，在模型(1)的设定下，给定生成的特征，估计量 \((\hat{b},\hat{w})\)的理论方差应为
其中\(\sigma^2=1\)是模型(1)中误差项 \(\epsilon\) 的方差。根据大数定律，对于每次蒙特卡洛模拟，\(n^{-1}(X\trans X)\) 的变化“不会太大”，尤其是当 \(n\) 足够大时。因此，估计量的方差应该和 \(n^{-1}\) 同阶。为了对模型参数真值 \((b_0,w_0)\) 做统计推断，我们需要编写函数用来得到估计量的方差估计。我们通常用残差项的样本方差估计 (4) 中的参数 \(\sigma^2\)，即利用
的样本方差 \(\hat{\sigma}^2\) 估计 \(\sigma^2\)。
def se_est(x,y):
# x: 长度为n的特征向量
# y: 长度为n的标签向量
est_par = estimation_vectorization(x,y) #获取 (b_0, w_0) 的估计量
error_est = y - est_par[0] - est_par[1] * x #获取误差项
sigma2_est = np.var(error_est) #获取 \hat{\sigma}^2 的值
aug_x = np.concatenate((np.ones_like(x),x), axis = 1) #用截距项增广特征向量
xx = aug_x.transpose() @ aug_x #通过 X^TX 获取 'xx'，其中 X 是大小为 nX2 的增广矩阵
sd_par = np.sqrt(sigma2_est * np.diag(np.linalg.inv(xx))) #根据上述结果获取方差估计量
return sd_par
\(\hspace{1.5em}\) 接下来，我们通过比较计算的方差和1,000次蒙特卡洛模拟所得到估计量的方法，来检验我们的函数是否编写正确。
n = 1000
sd_result = np.zeros((1000,2))
for i in range(1000):
x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i
sd_result[i,:] = se_est(x,y)
print('函数计算结果为:')
print(np.mean(sd_result,axis = 0))
print('蒙特卡洛模拟结果为：')
print(np.array([0.04361731, 0.01530787]))
函数计算结果为:
[0.04472754 0.01581741]
蒙特卡洛模拟结果为：
[0.04361731 0.01530787]
\(\hspace{1.5em}\) 如果我们对方差的估计结果非常接近于蒙特卡洛模拟结果，这就在某种程度上说明了我们所编写的方差估计函数没有问题。基于此，我们可以进一步构造双边置信区间，并检验所置信区间的覆盖率。这一步分析对于传统统计模型来说较为重要。我们可以证明，在模型(1)设定下，当样本量 \(n\) 趋于无穷大时，以下渐近结果成立：
其中 \(\Gamma_n\Gamma_n\trans = (X\trans X)(\hat{\sigma}^2)^{-1}\), \(I_2\) 是一个 \(2\times 2\) 的单位矩阵。利用这个结果，我们还可以检验两侧95%置信区间的覆盖率，即：
其中 \(q_\alpha\) 是标准正态分布的 \(\alpha\) 下分位数，而 \(\hat{\sigma}_b\) 和 \(\hat{\sigma}_w\) 分别是这两个估计量的估计标准误。对于每个模拟的训练集，我们需要检查所构建的置信区间是否包含了参数的真实值，并以此来检验置信区间的覆盖率。在如下的函数中，我们将用到逻辑运算and 。
def cr_indicator(x,y,alpha):
# x: 长度为n的特征向量
# y: 长度为n的标签向量
# alpha: 显著性水平。例如，alpha=0.05对应95%的双侧置信区间。
quan_normal = scipy.stats.norm.ppf(1-alpha/2) #获取标准正态分布的第 (1-alpha/2) 分位数。
est_par = estimation_vectorization(x,y) #获取参数估计值。
est_se = se_est(x,y) #获取标准误估计值。
ind_b = est_par[0] - quan_normal * est_se[0] < 1 and est_par[0] + quan_normal * est_se[0] > 1 #生成一个指标，表示b_0的真实值是否位于置信区间内。
ind_w = est_par[1] - quan_normal * est_se[1] < 1 and est_par[1] + quan_normal * est_se[1] > 1 #生成一个指标，表示w_0的真实值是否位于其置信区间内。
return np.array([ind_b, ind_w])
\(\hspace{1.5em}\) 下面，我们计算对应置信水平下的覆盖率。
import scipy.stats # 为了得到标准正态分布的分位数
cr_est = np.ones_like(par_result)-10 # 通常来讲，我们需要空矩阵记录每一次蒙特卡洛模拟结果。
for i in range(1000):
x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i
cr_est[i,:] = cr_indicator(x, y, alpha = 0.05)
print('在置信水平为0.95时，双侧置信区间的覆盖率为:')
print(np.mean(cr_est, axis = 0))
在置信水平为0.95时，双侧置信区间的覆盖率为:
[0.955 0.952]
\(\hspace{1.5em}\) 以上结果非常接近于0.95，则说明我们的理论和程序是正确的。
计算效率#
\(\hspace{1.5em}\) 尽管基于for循环的求和法与向量化方法得出的估计结果相同，但求和法的计算效率较低，特别是在样本量较大时。在本节中，我们将对这两种方法的计算效率进行比较。我们依然以样本量 \(n=1\,000\) 的情况为例进行说明。
1import time # 用来记录算法的计算时间
2T1 = time.time()
3for i in range(1000):
4 x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i
5 par_est = estimation_summation(x,y)
6T2 = time.time()
7print('基于for循环的求和法耗费了大约%4.4s秒' % ((T2 - T1)))
8
9T1 = time.time()
10for i in range(1000):
11 x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i
12 par_est = estimation_vectorization(x,y)
13T2 = time.time()
14print('向量化方法耗费了大约%4.4s秒' % (T2 - T1))
基于for循环的求和法耗费了大约1.32秒
向量化方法耗费了大约0.02秒
相比之下，我们可知向量化方法的计算速度远快于求和法。
习题#
在对一个统计模型进行蒙特卡洛模拟时，我们通常会考虑多种不同的设置。为了简化讨论，我们仅考虑样本量变化的情况。设样本量 \(n\) 分别取 \(\{100, 300, 500, 1\,000\}\)，并使用图表来汇报以下结果。请注意，在所有图表中，样本量应作为 \(x\) 轴。
1.1. 偏差盒装图：以 \(y\) 轴表示偏差，分别绘制并讨论两个参数估计量在1000次蒙特卡洛模拟中的偏差。
1.2. 计算效率图表：展示两种估计方法的计算效率，不同方法应使用不同的颜色。绘制图表并对结果进行评论。
广义线性回归#
\(\hspace{1.5em}\) 在本节中，我们将简单回顾我们已经学过的统计模型，包括基本的线性回归以及逻辑逻辑回归。需要指出的是，全连接神经网络与逻辑回归模型具有较强的相似性。
线性回归#
\(\hspace{1.5em}\) 我们已经在简单线性回归部分中，结合Python基本命令简要介绍了如何分析一个特殊的、具有一元特征的线性回归模型。在本节中，我们首先将在更一般的情况下，介绍线性模型的基本统计结果。我们考虑以下线性回归模型：
其中，\(y_i\in\mathbb{R}\) 为第 \(i\) 个实例的标签，\(b_0\) 为模型偏置项(bias)2在不产生混淆的情况下，我们用下角标``0``代表模型真实值，用没有下角标的对应参数作为损失函数的自变量。，\(\bx_i\in\mathbb{R}^d\) 为长度为 \(d\) 的特征向量，\(\bw_0\in\mathbb{R}^d\) 为权重(weight)向量，\(\epsilon_i\) 为满足 \(E(\epsilon_i\mid\bx_i)=0\) 的白噪声，\(n\) 为样本量。通常，我们还假设噪声的同质性(homogeneity)，即 \(\var(\epsilon_i\mid\bx_i)=\sigma^2\,(i=1,\ldots,n)\)，其中 \(\sigma^2>0\) 为固定常数。
备注
\(\hspace{1.5em}\) 本部分所讨论的模型(5)是对之前线性模型的一般化推广。需要指出的是，“偏置项”和“权重项”两个概念也将在神经网络中有所涉及，这两个概念在传统的统计学模型中被称为“截距项”和“斜率项”。请不要将这里的“偏置项(bias)”与统计学中的基本概念“偏差(bias)”相混淆。
\(\hspace{1.5em}\) 记 \(\tilde{\bx}_i=(1,\bx_i\trans)\trans\)，对应地记模型参数为 \({\btheta}_0=(b_0,\bw_0\trans)\trans\)。在新的符号下，模型(5)等价于：
\(\hspace{1.5em}\) 在本节中，我们所提出的模型形式和(6)相同，只不过我们用 \(\btheta\) 作为所提出模型的模型参数。为了估计所提出模型的参数，我们考虑如下 \(l_2\) 损失函数：
其中，\(\bY = (y_1,\ldots,y_n)\trans\)，\(\tilde{\bX}=(\tilde{\bx}_1,\ldots,\tilde{\bx}_n)\trans\)。
备注
\(\hspace{1.5em}\) 在进一步假设噪声项 \(\epsilon_i\) 服从均值为0，方差为 \(\sigma^2\) 的正态分布时，最小化损失函数 \(\mathcal{J}(\tilde{\bw})\) 等价于最大化关于 \(\btheta\) 的（对数）似然函数。
\(\hspace{1.5em}\) 在一般条件下，我们通过计算 \(\partial^2\mathcal{J}(\btheta)/(\partial\btheta\partial\btheta\trans)\) 可知，\(\mathcal{J}(\btheta)\) 是关于参数 \(\btheta\) 的凸函数，故模型参数估计值 \(\hat{\btheta}\) 满足
\(\hspace{1.5em}\) 通过计算(8)，我们可知参数估计量 \(\hat{\btheta}\) 满足如下正则方程
备注
\(\hspace{1.5em}\) 通过正则方程(9)，我们可知对应于估计量 \(\hat{\btheta}\) 对应的残差 \(\bY-\tilde{\bX}\hat{\btheta}\) 与设计矩阵 \(\tilde{\bX}\) 的列正交。特别地，我们有 \(\sum_{i=1}^n(y_i - \tilde{\bx}_i\trans\hat{\btheta})\tilde{\bx}_i=0\)。
\(\hspace{1.5em}\) 通过(9)可知，线性模型的参数估计量为
逻辑回归#
\(\hspace{1.5em}\) 逻辑回归的模型为
其中 \(\mbox{Bernoulli}(p)\) 为成功概率为 \(p\) 的伯努利二项分布，\(\tilde{\bx}_i = (1,\bx_i\trans)\trans\)，\(\btheta_0=(b_0,\bw_0\trans)\trans\) 为模型参数真值以及 \(\sigma(z) = \{1+\exp(-z)\}^{-1}\)。在（早期）深度学习中，\(\sigma(z)\) 通常被用作非线性变换所需要的函数（激活函数）。由于特征向量 \(\bx_i\) 以线性形式存在于成功概率 \(\sigma(\tilde{\bx}_i\trans\btheta_0)\) 中，逻辑回归模型也是一种常见的广义线性模型。显然，随着 \(\tilde{\bx}_i\trans\btheta_0\) 的增加，\(Y_i=1\) 的概率也在增加3我们通常用小写字母表示观测，而用相应的大写字母表示随机变量。例如，\(y_i\) 表示第 \(i\) 个样本（观测到）的标签，而 \(Y_i\) 表示其对应的随机变量。。
备注
\(\hspace{1.5em}\) 关于（广义）线性模型更加深入的探讨，请参见Hastie等（2009）。
\(\hspace{1.5em}\) 在给定特征向量 \(\bx_i\) 时，由于我们在逻辑回归模型(11)中假设了关于 \(y_i\) 的条件分布，故我们可以通过最大化关于 \(\btheta\) 的对数似然函数对该参数进行估计。因此，我们将其负对数似然函数作为估计参数 \(\btheta\) 的损失函数，即
其中，\(a_i(\btheta) = \sigma\{𝑧_i(\btheta)\}\)，\(z_i(\btheta) = \tilde{\bx}_i\trans\btheta\)。注意，此处我们用了符号a和z。在随后的深度学习部分内容中，a代表了激活后的结果，而z将代表线性运算后的结果；其具体定义我们将在后面的章节中进行详细讨论。
\(\hspace{1.5em}\) 与线性回归结果类似，在一般条件下，我们可证明 \(\mathcal{J}(\btheta)\) 是关于参数 \(\btheta\) 的凸函数，故损失函数 \(\mathcal{J}(\btheta)\) 具有唯一最小值，其可通过求解
得到。然而，对于逻辑回归的参数估计，我们没有对应的正则方程，即我们无法得到参数估计的显示表达式。在传统的统计学习中，我们通常利用Newton-Raphson算法4参见Hastie等（2009）第4.4.1节。对方程(13)进行求解。
Newton-Raphson算法#
\(\hspace{1.5em}\) 该算法需要计算损失函数的二阶混合偏导，即黑森矩阵(Hessian matrix)，
Newton-Raphon算法
初始化权重向量 \(\btheta^{(0)}\)。
给定 \(\btheta^{(t)}\)，计算
重复第二步直至收敛。
\(\hspace{1.5em}\) 得益于黑森矩阵，Newton-Raphson算法的收敛速度非常快。然而，当权重向量的维度很大时，计算黑森矩阵及其逆矩阵的计算复杂度较大。在深度神经网络模型中，其参数规模往往较大。在参数求解时，我们往往不用Newton-Raphson算法，取而代之的为梯度下降算法。
（批量）梯度下降算法#
\(\hspace{1.5em}\) 对于具有较高维度的模型，例如神经网络模型等，计算黑森矩阵的复杂度非常高。因此，Newton-Raphson算法对这种模型的训练而言计算效率比较低下。在训练神经网络等复杂模型时，我们通常使用梯度下降算法。
（批量）梯度下降法
初始化权重向量 \(\btheta^{(0)}\)。
给定 \(\btheta^{(t)}\)，计算
重复第二步直至收敛。
\(\hspace{1.5em}\) 与Newton-Raphson算法不同，（批量）梯度下降法5批量(Batch)是指我们在每步更新时，使用所有的训练样本。如果每次更新只使用部分样本，对应的算法叫做（小批量）梯度下降法或者随机梯度下降法。利用学习率 \(\alpha\) 替代了黑森矩阵，这就大大提高了该算法的计算效率。我们将在接下来的章节中详细讲解该算法及其变种。
交叉熵#
\(\hspace{1.5em}\) 逻辑回归模型通常被用于分析二分类问题，而交叉熵是该类问题中十分重要的概念，其与损失函数(12)的形式非常类似。在本节中，我们简单讨论交叉熵的基本概念。需要指出的是，本节所用到的数学符号与其他章节略有不同。
\(\hspace{1.5em}\) 首先，基于Jessen不等式6可参见概率论基本教材。，我们有如下结果。
极大似然
\(\hspace{1.5em}\) 在特定条件下，我们有
其中，\(X \sim P\) 表示随机变量 \(X\) 的分布为 \(P\)，该分布的密度函数为 \(p(x)\) ，\(q(x)\) 表示任意分布。
\(\hspace{1.5em}\) 根据以上定理，我们可知当随机变量 \(X\) 的真实分布为 \(P\)，且该分布具有密度函数 \(p(x)\) 时，该真实密度函数 \(p(x)\) 最小化如下关于密度函数 \(q(x)\) 的目标函数：
我们将(14)定义为“分布” \(p\) 和 \(q\) 的交叉熵，记为
其中，\(p\) 为随机变量 X 的真实密度函数，\(q\) 为任意密度函数。交叉熵在衡量分布 \(q\) 对随机变量 \(X\) 分布近期的效果。
\(\hspace{1.5em}\) 在实际问题中，我们往往无法观测到随机变量 \(X\) 的总体分布，而往往只能观测到一个样本量为 \(n\) 的随机样本 \(S=\{x_1,\ldots,x_n\}\)。此时，我们可以利用大数定理，得到对(15)的近似，即
我们希望找到一个密度函数 \(q\) 使得 \(\hat{H}(p,q)\) 达到最小。
\(\hspace{1.5em}\) 例如，如果我们想要找一个正态密度函数 \(q(x;\mu,\sigma^2) = (2\pi\sigma^2)^{-1/2}\exp\{-(x-\mu)^2/(2\sigma^2)\}\) 对真实密度函数 \(p(x)\) 进行近似，那么我们需要寻找参数 \(\mu\) 和 \(\sigma^2\) 最小化
我们可以观察到，(17)正是关于参数 \(\mu\) 和 \(\sigma^2\) 的负对数似然函数。因此，最小化 (17) 等价于最大化对应的对数似然函数。从本例中可以看出，最小化交叉熵与最大化对数似然在一定的条件下是等价的。
备注
\(\hspace{1.5em}\) 在以上分析中，我们是通过大数定理揭示(15)和(16)的关系的。理解(16)的另一个角度为 \(\hat{H}(p,q) =H(p_n,q)\)， 其中 \(p_n(x)=n^{-1}\delta_S(x)\) 为样本 \(S\) 的经验密度函数，\(\delta_S(x)\) 为示性函数，若 \(x\in S\)，则 \(\delta_S(x)=1\)，否则为0。即 \(\hat{H}(p,q)\) 所衡量的是密度函数 \(q(x)\) 对经验密度函数 \(p_n(x)\) 的近似效果。
\(\hspace{1.5em}\) 关于逻辑回归模型的编程，请参见 与逻辑回归模型的比较。
逻辑回归与交叉熵#
\(\hspace{1.5em}\) 在本节的最后，我们简单讨论下逻辑回归模型的参数估计过程与交叉熵的关系。不失一般性，我们假设 \(S=\{(\bx_i,y_i):i=1,\ldots,n\}\) 为一个独立同分布的随机样本。假设 \(p(\bx)\) 与 \(p(y\mid \bx)\) 分别为特征向量 \(\bx\) 的真实边际密度函数和随机变量 \(y\) 在给定特征向量 \(\bx\) 下的真实条件密度函数7在讨论密度函数时，我们应用了Radon-Nikodym定理。。根据逻辑回归的模型假设(11)，\(p(y\mid \bx)=\sigma(\tilde{\bx}\trans\btheta_0)^y\{1-\sigma(\tilde{\bx}\trans\btheta_0)\}^{1-y}\)，而我们不对边际密度函数 \(p(\bx)\) 做任何参数化假设。
\(\hspace{1.5em}\) 根据(16)，权重向量 \(\btheta_0\) 的估计量 \(\hat{\btheta}\) 最小化
其中，\(q(y\mid \bx;\btheta) = \sigma(\tilde{\bx}\trans\btheta)^y\{1-\sigma(\tilde{\bx}\trans\btheta)\}^{1-y}\)。从(18)可知，其第一项与参数无关，第二项即为逻辑回归模型中的损失函数(12)。因此，最小化(12)对应的目标函数 \(\mathcal{J}(\btheta)\) 等价于最小化交叉熵(18)。
习题#
记 \(\bx=(x_1,\ldots,x_d)\trans\in\mathbb{R}^d\)，\(\ba=(a_1,\ldots,a_d)\trans\in\mathbb{R}^d\)，\(\bA=(a_{ij})\in\mathbb{R}^{d\times d}\)，回答如下问题：
- 1.1 假设向量 \(\ba\) 为定值，考虑关于 \(\bx\) 的函数 \(f(\bx)=\bx\trans\ba=\sum_{i=1}^dx_ia_i\)。对于 \(i=1,\ldots,d\)，计算 \(\partial f(\bx)/\partial x_i\)，并证明如下结论：
- 1.2 假设矩阵 \(\bA\) 为定值，考虑关于 \(\bx\) 的函数 \(g(\bx)=\bx\trans \bA\bx=\sum_{i=1}^d\sum_{j=1}^da_{ij}x_ix_j\)。对于 \(i=1,\ldots,d\) 以及 \(j=1,\ldots,d\) 计算 \(\partial g(\bx)/\partial x_i$以及$\partial^2g(\bx)/(\partial x_i\partial x_j)\)，并证明如下两个结论：