现在我们每天几乎都会通过各种渠道看到几则有关人工智能(AI)的新闻,介绍AI的各式研发成果,很多人基于好奇也许会想一窥究竟,了解背后运用的技术与原理,就会发现一堆数学符号及统计公式,也许会产生疑问:要从事AI系统开发,非要搞定数学、统计不可吗?答案是肯定的,我们都知道机器学习是从数据中学习到知识(Knowledge Discovery from Data, KDD),而算法就是从数据中萃取出知识的“榨汁机”,它必须以数学及统计为理论基础,才能证明其解法具有公信力与精准度,然而数学/统计理论都有局限,只有在假设成立的情况下,算法才是有效的,因此,如果不了解算法中的各个假设,随意套用公式,就好像无视交通规则,在马路上任意飚车一样危险。
(图片来源: Decision makers need more math [1] )
因此,以深度学习而言,我们至少需要熟悉下列学科:①线性代数(Linear Algebra);②微积分(Calculus);③概率论与数理统计(Probability and Statistics);④线性规划(Linear Programming),如图2.1所示。
图2.1 必备的数学与统计知识
以神经网络优化求解的过程为例,上述四门学科就全部用上了,如图2.2所示。
图2.2 神经网络权重求解过程
(1)正向传导:借由线性代数计算误差及损失函数。
(2)反向传导:透过偏微分计算梯度,同时,利用线性规划优化技巧寻找最佳解。
(3)统计串联整个环节,如数据的探索与分析、损失函数定义、效果衡量指标等,都是基于统计的理论架构而成的。
(4)深度学习的推论以概率为基础,预测目标值。
四项学科相互为用,贯穿整个求解过程,因此,要通晓深度学习的运作原理,并且正确选用各种算法,甚至进而能够修改或创新算法,都必须对其背后的数学和统计有一定基础的认识,以免误用或滥用。
四门学科在大学至少都是两学期的课程,对已经离开大学殿堂很久的工程师而言,在上班之余,还要重修上述课程,相信大部分的人都会萌生退意了!那么我们是否有速成的快捷方式呢?
笔者在这里借用一个概念Statistical Programming,原意是“以程序解决统计问题”,换个角度想,我们是不是也能 以程序设计的方式学统计 ,以缩短学习历程呢?
通常我们按部就班地学习数学及统计,都是从“假设”→“定义/定理”→“证明”→“应用”,一步一步学习。
(1)“假设”是“定义/定理”成立的前提。
(2)“证明”是“定理”的验证。
(3)“应用”是“定义/定理”的实践。
由于“证明”都会有一堆的数学符号及公式推导,经常会让人头晕脑涨,降低学习的效率,因此,笔者大胆建议,工程师将心力着重在假设、定义、定理的理解与应用,并利用程序进行大量个案的验证,虽然忽略了证明的做法,会让学习无法彻底地融会贯通,但是对已进入职场的工程师会是一种较为可行的快捷方式。
接下来我们就以上述做法,对四项学科进行重点介绍,除了说明深度学习需要理解的知识外,更强调如何 以撰写程序实现相关理论及解题方法 。
张量(Tensor)是描述向量空间(Vector Space)中物体的特征,包括零维的纯量、一维的向量(Vector)、二维的矩阵(Matrix)或更多维度的张量,线性代数则是说明张量如何进行各种运算,它被广泛应用于各种数值分析的领域。以下就以实例说明张量的概念与运算。
向量是一维的张量,它与线段的差别是除了长度(Magnitude)以外,还有方向(Direction),其数学表示法为
图形表示如图2.3所示。
图2.3 向量长度与方向
下面使用程序计算向量的长度与方向, 请参阅02_01_线性代数_向量.ipynb 。
1. 长度
计算公式为欧几里得距离(Euclidean Distance),即
程序代码如下:
也可以使用np.linalg.norm()计算向量长度,程序代码如下:
2. 方向
使用tan -1 ()函数计算
移项为
程序代码如下:
3. 向量四则运算规则
(1)加减乘除一个常数:常数直接对每个元素作加减乘除。
(2)加减乘除另一个向量:两个向量相同位置的元素作加减乘除,所以两个向量的元素个数必须相等。
4. 向量加减法
向量加减一个常数,长度、方向均改变。
(1)程序代码如下:
(2)执行结果:如图2.4所示。
图2.4 向量加减一个常数,长度、方向均改变
5. 向量乘除法
向量乘除一个常数,长度改变、方向不改变。
(1)程序代码如下:
(2)执行结果:如图2.5所示。
图2.5 向量乘除一个常数,长度改变、方向不改变
6. 向量加减乘除另一个向量
两个向量的相同位置的元素作加减乘除。
(1)程序代码如下:
(2)执行结果:如图2.6所示。
图2.6 向量加另一个向量
7. “内积”或称“点积”(Dot Product)
公式为
numpy是以@作为内积的运算符号,而非*。程序代码如下:
8. 计算两个向量的夹角
公式为
移项,得
再利用cos -1 ()计算夹角 θ 。程序代码如下:
矩阵是二维的张量,拥有行(Row)与列(Column),可用于表达一个平面 N 个点( N ×2)、或一个3D空间 N 个点( N ×3),例如:
矩阵加法/减法与向量相似,相同位置的元素作运算即可,但乘法运算通常是指内积,使用@。
以下程序请参考02_02_线性代数_矩阵.ipynb。
1. 两个矩阵相加
程序代码如下:
2. 两个矩阵相乘
解题:左边矩阵的第二维须等于右边矩阵的第一维,即(m, k)×(k, n)=(m, n),则有
其中左上角的计算过程为(1,2,3)×(9,7,5)=(1×9)+(2×7)+(3×5)=38,右上角的计算过程为(1,2,3)×(8,6,4)=(1×8)+(2×6)+(3×4)=32,以此类推,如图2.7所示。
图2.7 矩阵相乘
程序代码如下:
3. 矩阵( A 、 B )相乘
A × B 是否等于 B × A ?程序代码如下:
执行结果: A × B 不等于 B × A 。
4. 特殊矩阵
矩阵在运算时,除了一般的加减乘除外,还有一些特殊的矩阵,包括转置矩阵(Transpose)、反矩阵(Inverse)、对角矩阵(Diagonal Matrix)、单位矩阵(Identity Matrix)等。
(1)转置矩阵:列与行互换。例如:
( A T ) T = A :进行两次转置,会恢复成原来的矩阵。
对上述矩阵作转置。程序代码如下:
也可以使用np.transpose(A)。
(2)反矩阵( A -1 ): A 必须为方阵,即列数与行数须相等,且必须是非奇异方阵(Non-singular),即每一列或行之间不可以相异于其他列或行。程序代码如下:
执行结果如下:
(3)单位矩阵:若 A 为非奇异(Non-singular)矩阵,则 A @ A -1 =单位矩阵( I )。所谓的非奇异矩阵是指任一行不能为其他行的倍数或多行的组合,包括各种四则运算,矩阵的列也须符合相同的规则。
试对下列矩阵验算 A @ A -1 是否等于单位矩阵( I )。
程序代码如下:
执行结果如下:
结果为单位矩阵,表示 A 为非奇异矩阵。
试对下列矩阵验算 A @ A -1 是否等于单位矩阵( I )。
程序代码如下:
执行结果如下:
A 为奇异(Singular)矩阵,因为
第二列=第一列+ 1
第三列=第一列+ 2
故 A @ A -1 不等于单位矩阵。
在中学阶段,我们通常会以高斯消去法(Gaussian Elimination)解联立方程式。以下列方程式为例:
x + y =16
10 x +25 y =250
将第一个方程式两边乘以-10,加上第二个方程式,即可消去 x ,变成
−10( x + y )=−10(16)
10 x +25 y =250
简化为
15 y =90
得到 y =6,再代入任一方程式,得到 x =10。
以上过程,如果以线性代数求解就简单多了。
(1)以矩阵表示: A 为方程式中未知数( x , y )的系数, B 为等号右边的常数且 AX = B 。
(2)其中 ,则 X = A -1 B 。证明如下:
①两边各乘 A -1 得
A -1 AX = A -1 B
② A -1 A 等于单位矩阵,且任一矩阵乘以单位矩阵,还是等于原矩阵,故
X = A -1 B
③以上式求得( x , y )。注意:前提是 A 须为非奇异矩阵。
以下程序均收录在02_03_联立方程式求解.ipynb。
(1)以NumPy库求解上述联立方程式。程序代码如下:
①inv(A): A 的反矩阵。
②执行结果: x =10, y =6。
③也可以直接使用np.linalg.solve()函数求解。程序代码如下:
(2)画图,交叉点即联立方程式的解。程序代码如下:
执行结果:如图2.8所示。
图2.8 交叉点即联立方程式的解
(3)以NumPy库求解下列联立方程式。
−1 x +3 y =−72
3 x +4 y −4 z =−4
−20 x −12 y +5 z =−50
程序代码如下:
①执行结果:( x , y , z )=(12, -20, -10)。
②也可以使用SymPy库求解,直接将联立方程式整理在等号左边,使用solve()函数,参数内的多项式均假设等号(=)右边为0。
微积分包括微分(Differentiation)与积分(Integration),微分是描述函数某一点的变化率,借由微分可以得到特定点的斜率(Slope)或梯度(Gradient),即变化的速度,二次微分可以得到加速度。积分则是微分的逆运算,积分可以计算长度、面积、体积,也可以用来计算累积的概率密度函数(Probability Density Function)。
在机器学习中,微分在求解的过程中占有举足轻重的地位,因此TensorFlow框架就提供自动微分(Automatic Differentiation)的功能,用以计算各特征变量的梯度,进而求得模型的权重(Weight)参数。
微分用于描述函数的变化率(Rate of Change),如 y =2 x +5,表示 x 每增加一单位, y 会增加2。因此,变化率就等于2,也称为斜率;5为截距(Intercept),或称偏差(Bias),如图2.9所示。
图2.9 斜率与截距
我们先不管截距,只看斜率,算法为:取非常相近的两个点(距离 h 趋近于0), y 坐标值之差(Δ y )除以 x 坐标值之差(Δ x ),有
这就是微分的定义,但上述极限值(limit)不一定存在,其存在的要素如下:
(1) h 为正值时的极限值等于 h 为负值时的极限值,亦即函数在该点时是连续的。
(2)上述极限值不等于无穷大(∞)或负无穷大(-∞)。
如图2.10所示的函数在 x =5的地方是连续的,由上方(5.25)逼近,或由下方(4.75)逼近是相等的, 相关彩色图形可参考02_04_微分.ipynb 。
图2.10 连续函数
相反地,图2.11所示函数在 x =0时是不连续的,逼近 x =0时有两个解。
图2.11 不连续函数
接着来看看几个应用实例。
以下程序请参考02_04_微分.ipynb。
(1)试绘制一次方函数 f ( x )=2 x +5。程序代码如下:
执行结果:如图2.12所示。
图2.12 一次方函数执行结果
由执行结果可以看出,一次方函数每一点的斜率均相同。
(2)试绘制二次方曲线 f ( x )=−10 x 2 +100 x +5,求最大值。程序代码如下:
执行结果:如图2.13所示。由执行结果可以得到以下结论。
图2.13 二次方曲线执行结果
①一次方函数整条在线的每一个点的斜率都相同,但是二次方曲线上的每一个点的斜率就都不一样了,如图2.13所示, 相关彩色图形可参考02_04_微分.ipynb 。
● 绿线(细拋物线):二次曲线,是一条对称的拋物线。
● 紫线(斜线):拋物线的一阶导数。
● 红线(拋物线的切线):三个点(2, 5, 8)的斜率。
②每一个点的斜率即该点与二次曲线的切线(红线),均不相同,斜率值可通过微分求得一阶导数(图中的斜线),随着 x 变大,斜率越来越小,二次曲线的最大值就发生在斜率等于0的地方,当 x =5时, f ( x )=255。
(3)试绘制二次方曲线 f ( x )= x 2 +2 x +7,求最小值。程序代码如下:
执行结果:如图2.14所示。
图2.14 二次方曲线执行结果
由执行结果可得:斜率值可通过微分求得一阶导数(图中的斜线),随着 x 变大,斜率越来越大,二次曲线的最小值就发生在斜率等于0的地方,当 x =-1时, f ( x )=6。
综合范例(2)(3),可以得知微分两次的 二阶导数( f "( x ))为常数 ,且为正值时,函数有最小值,反之,为负值时,函数有最大值。但若 f ( x )为三次方(以上)的函数 ,一阶导数等于0的点,可能只是区域的最佳解(Local Minimum/Maximum),而不是全局最佳解(Global Minimum/Maximum)。
(4)试绘制三次方曲线 f ( x )= x 3 −2 x +100,求最小值。程序代码如下:
执行结果:如图2.15所示。由执行结果可以得到以下结论。
图2.15 三次方曲线执行结果
①三次方曲线 f ( x )= x 3 −2 x +100在斜率等于0的点只是区域的最佳解。
②三次方曲线一般为凸函数时才有全局最佳解。
上述程序中的fd()函数为一阶导数,它们是如何求得的?只要运用以下微分的定理,就可以轻易解出上述范例的一阶导数,相关定理整理如下。
(1) f ( x )一阶导数的表示法为 f ′( x )或 。
(2) f ( x )为常数( C )→ f ′( x )=0。
(3) f ( x )= Cg ( x )→ f ′( x )= Cg ′( x )。
(4) f ( x )= g ( x )+ h ( x )→ f ′( x )= g ′( x )+ h ′( x )。
(5)次方的规则: f ( x )= x n → f ′( x )= nx n −1 。
(6)乘积的规则: 。
(7)商的规则:若 r ( x )= s ( x )/ t ( x ),则 。
(8)链式法则(Chain Rule):
以上一节范例(2) f ( x )=−10 x 2 +100 x +5为例,针对多项式的每一项个别微分再相加,就得到 f ( x )的一阶导数为
f ′( x )= -20 x + 100
SymPy库直接支持微积分函数的计算,可以验证定理,接下来我们就写一些程序来练习一下。
以下程序请参考02_04_微分.ipynb。
(1) f ( x )为常数 。程序代码如下:
执行结果:0。
(2) 。程序代码如下:
执行结果均为10 x 。
(3)乘积的规则: 。程序代码如下:
执行结果:5 x 4 。
(4)链式法则: 。程序代码如下:
执行结果:6 x 5 。
(5)验证 f ( x )=−10 x 2 +100 x +5。程序代码如下:
执行结果:100 -20 x 。
接着,利用一阶导数等于0,求最大值。程序代码如下:
执行结果: x =5时,最大值为255。
(6)验证 f ( x )= x 2 +2 x +7。程序代码如下:
执行结果:2 x + 2。
接着,利用一阶导数等于0,求最小值。程序代码如下:
执行结果: x =-1时,最小值为6。
在机器学习中,偏微分(Partial Differentiation)在求解的过程中占有举足轻重的地位,常用来计算各特征变量的梯度,进而求得最佳权重。梯度与斜率相似,单一变量的变化率称为斜率,多变量的斜率称梯度。
在上一节中, f ( x )只有单变量,如果 x 是多个变量 x 1 、 x 2 、 x 3 …,要如何找到最小值或最大值呢?这时,我们就可以使用偏微分求取每个变量的梯度,让函数沿着特定方向寻找最佳解,如图2.16所示即沿着等高线(Contour),逐步向圆心逼近,这就是所谓的梯度下降法(Gradient Descent)。
图2.16 梯度下降法图解
以下程序请参考02_05_偏微分.ipynb。
(1)假设 f ( x , y )= x 2 + y 2 ,请分别对 x 、 y 作偏微分。
解题:先针对 x 作偏微分,有
将其他变量视为常数,对每一项个别微分,得
加总,得
再针对 y 作偏微分,同样将其他变量视为常数,有
总结得到
程序代码如下:
其中第8、12行的y1.diff(x)、y1.diff(y)分别表示对 x 、 y 作偏微分。
(2)假设 f ( x )= x 2 ,请使用梯度下降法找最小值。
解题:程序逻辑如下。
①任意设定一起始点(x_start)。
②计算该点的梯度fd( x )。
③沿着梯度更新 x ,逐步逼近最佳解,幅度大小以学习率控制。新的 x =0-学习率(Learning Rate)*梯度。
④重复步骤②和③,判断梯度是否接近于0,若已很逼近于0,即可找到最佳解。
程序代码如下:
执行结果:如图2.17所示。
图2.17 梯度下降法实践
得到 x 每一点的坐标为: 。
我们可以改变第24~26列的参数,观察执行结果。
①改变起始点x_start = -5,依然可以找到最小值,如图2.18所示。
图2.18 改变x_start参数
②设定学习率lr = 0.9:如果函数较复杂,可能会跳过最小值,如图2.19所示。
图2.19 设定lr = 0.9
③设定学习率lr = 0.01:还未逼近最小值,就提早停止了,可以增加执行周期数来解决问题,如图2.20所示。
图2.20 设定lr = 0.01
上述程序是神经网络优化器求解的简化版,在后续的章节会进行详细的探讨,目前只聚焦在说明偏微分在深度学习的应用。
在优化求解中,也经常利用一阶导数等于0的特性,求取最佳解。例如,以普通最小二乘法(Ordinary Least Square, OLS)对简单线性回归 y = wx + b 求解。首先定义目标函数(Object Function)或称损失函数(Loss Function)为均方误差(MSE),即预测值与实际值差距的平方和,MSE当然越小越好,所以它是一个最小化的问题,所以,我们可以利用偏微分推导出公式,过程如下。
(1)
式中: ε 为误差,即实际值( y )与预测值( )之差; n 为为样本个数。
(2)MSE = SSE / n ,我们忽略常数 n ,可以只考虑SSE,即
(3)分别对 w 及 b 作偏微分,并且令一阶导数等于0,可以得到两个联立方程式,进而求得 w 及 b 的解。
(4)对 b 偏微分,又因
两边同除以-2,得
分解得
除以 n , 为 x 、 y 的平均数,有
移项得
(5)对 w 偏微分,得
两边同除以-2,得
分解得
代入(4)的计算结果 ,有
化简得
(6)结论。
以下程序请参考02_06_线性回归.ipynb。
现有一个世界人口统计数据集,以年度(year)为 x ,人口数为 y ,依上述公式计算回归系数 w 、 b 。程序代码如下:
执行结果: w =0.061159358661557375, b =-116.35631056117687。
使用NumPy的现成函数polyfit()验算。程序代码如下:
执行结果:
w =0.061159358661554586, b =-116.35631056117121,答案相差不多。
积分则是微分的相反运算,与微分互为逆运算,微分用于斜率或梯度的计算,积分则可以计算长度、面积、体积,也可以用来计算累积的概率密度函数。
(1)积分一般数学表示式为
(2)若 f ( x )= x ,则积分结果为 。
(3)限定范围的积分:先求积分,再将上、下限代入多项式,相减可得结果。例如
接下来撰写程序,看看积分如何计算。
以下程序请参考02_07_积分.ipynb。
1. 进行下列积分计算并作图
(1)积分运算:SciPy库支持积分运算。程序代码如下:
执行结果:积分值为4.5;误差为4.9960036108132044e-14。
程序说明如下。
①第3行载入SciPy库。
②第9行呼叫integrate.quad(),作限定范围的积分,参数如下。
● 函数 f ( x );
● 范围下限:设为负无穷大(-∞);
● 范围上限:设为正无穷大(∞);
● 输出:含积分结果及误差值。
(2)作图。程序代码如下:
执行结果:如图2.21所示。
图2.21 积分运算执行结果
程序说明:第11行fill_between()会填满整个积分区域。
2. 计算正态分布(Normal Distribution)的概率并作图
(1)正态分布的概率密度函数(Probability Density Function, PDF)即
(2)计算正态分布(-∞, ∞)的概率。程序代码如下:
执行结果:累积概率为1.0;误差为1.0178191437091558e-08。
程序说明如下。
①第7、8行定义平均数(mean)为0、标准差(std)为1,也就是“标准”正态分布,又称Z分布。
②第11列定义正态分布的概率密度函数(PDF)。
③注意:任何概率分布的总和必然为1,即所有事件发生的概率总和必然为100%。
3. 正态分布常见的置信区间
包括±1、±2、±3倍标准差,我们可以计算其概率。有关置信区间的定义会在概率与统计的章节进行说明。
①±1倍标准差的置信区间概率为68.3%。程序代码如下:
● 执行结果:累积概率为0.683;误差为7.579375928402476e-15。
● integrate.quad()参数设为-1及1。
②±2倍标准差的置信区间概率为95.4%。程序代码如下:
● 执行结果:累积概率为0.954;误差为1.8403560456416134e-11。
● integrate.quad()参数设为-2及2。
③±3倍标准差的置信区间概率为99.7%。程序代码如下:
● 执行结果:累积概率为0.997;误差为1.1072256503105314e-14。
● integrate.quad()参数设为-3及3。
④我们常用±1.96倍标准差的置信区间,概率为95%,概率刚好为整数。程序代码如下:
● 执行结果:累积概率为0.95;误差为1.0474096492701325e-11。
● integrate.quad()参数设为-1.96及1.96。
⑤另外,可以使用随机数及直方图,绘制标准正态分布图。程序代码如下:
执行结果:如图2.22所示。
图2.22 正态分布±1、±2、±3倍标准差的置信区间
程序说明如下。
● 第9行sns.distplot(x, hist=False):使用Seaborn画直方图,参数hist=False表示不画阶梯直方图,只画平滑曲线。
● 第17行axvline():画垂直线,标示±1、±2、±3倍标准差。
统计是从数据所推衍出来的信息,包括一些描述统计量、概率分布等,如平均数、标准差等。而我们搜集到的数据统称为数据集(Dataset),它包括一些观察值(Observations)或案例(Cases),即数据集的行,以及一些特征(Features)或称属性(Attributes),即数据集的字段。例如,要预测某市下一任市长,我们做一次问卷调查,相关定义如下。
(1)数据集:全部问卷调查数据。
(2)观察值:每张问卷。
(3)特征或属性:每一个问题,通常以 X 表示。
(4)属性值或特征值:每一个问题的回答。
(5)目标(Target)字段:市长候选人,通常以 y 表示。
下面我们会依序介绍下列内容:
(1)抽样;
(2)描述统计量;
(3)概率;
(4)概率分布函数;
(5)假设检定。具体如图2.23所示。
图2.23 基础统计介绍的范围及其关联
依照特征的数据类型,分为定性(Qualitative)及定量(Quantitative)的字段。定性字段为非数值型的字段,通常包含有限类别,又可以分为名义数据(Nominal Data)及有序数据(Ordinal Data);定量字段为数值型字段,又可以分为离散型数据(Discrete Data)及连续型数据(Continuous Data)。预测的目标字段如为离散型变量,则适用分类(Classification)的算法;反之,目标字段为连续型变量,则适用回归(Regression)的算法。
(1)名义数据:域值并没有顺序或大小的隐含意义,如颜色,红、蓝、绿并没有谁大谁小的隐含意义。转换为代码时,应以One-Hot Encoding处理,将每一类别转为个别的哑变量(Dummy Variable),每个哑变数只有True/False或1/0两种值,如图2.24所示。
图2.24 哑变量
Color特征有三种类别,经过One-Hot Encoding处理,会转换为三个哑变量。
以下程序请参考02_08_特征转换.ipynb。程序代码如下:
数据框(Data Frame)包含要转换的字段、前置符号(prefix)、分隔符(prefix_sep)。
(2)有序数据(Ordinal Data):域值有顺序或大小的隐含意义,如衣服尺寸XL >L > M > S,如图2.25所示。
图2.25 有序数据示例
size字段依尺寸大小,分别编码为XL(4)、L(3)、M(2)、S(1),转换后如图2.26所示。
图2.26 编码转换结果
程序代码如下:
程序说明如下。
①第2行:以字典定义转换规则,key为原值,value为转换后的值。
②第8行map():以字典转换size字段的每一个值。
再以预测某市下一任市长的选举为例,由于人力、时间及经费的限制,不太可能调查所有市民的投票倾向,通常我们只会随机抽样1000份或更多的样本进行调查,这种方式称之为抽样(Sampling),相关名词定义如下。
● 母体(Population):全体有投票权的市民。
● 样本(Sample):被抽中调查的市民。
● 分层抽样(Stratified Sampling):依母体某些属性的比例,进行相同比例的抽样,希望能充分代表母体的特性,如政党、年龄、性别比例等。
照例我们看看程序如何撰写, 请参阅02_09_抽样.ipynb 。
(1)简单抽样,从一个集合中随机抽出 n 个样本。程序代码如下:
执行结果:[1, 6, 2, 8, 5]。
程序说明如下。
①第8行random.sample():自集合中随机抽样。
②每次执行结果均不相同,且每个项目不重复,此种抽样方法称为不放回式抽样(Sampling Without Replacement)。
(2)放回式抽样(Sampling With Replacement)。程序代码如下:
执行结果:[8, 4, 7, 4, 6]。
程序说明如下。
①第8行random.choices():自集合中随机抽样。
②每次执行结果均不相同,但项目会重复,如上,4被重复抽出两次,此种抽样方法称为放回式抽样。
(3)以Pandas库进行抽样。程序代码如下:
执行结果:如图2.27所示。
图2.27 以Pandas库进行抽样执行结果
程序说明如下。
①第12行df.sample(5):自集合中随机抽样5个。
②此抽样为不放回式抽样。
(4)以Pandas库进行分层抽样。程序代码如下:
执行结果如下:
程序说明如下。
①第13行StratifiedShuffleSplit(test_size=6):将集合重新洗牌,并从中随机抽取6个数据。
②第14行stratified.split(df, df['y']):以y字段为key,分层抽样,得到的结果是generator数据类型,需转为list,才能一次取出。
③第一个输出是重新洗牌的全部数据,第二个是抽出的索引值。
④利用df.iloc[x[0][1]]指令,取得数据如图2.28所示。观察y列,每个类别各有两个与母体比例相同。
图2.28 利用df.iloc[x[0][1]]指令取得的数据
母体比例可以df['y'].value_counts()指令取得数据如下,0/1/2类别各50个。
(5)以Pandas库进行不分层抽样。程序代码如下:
执行结果:如图2.29所示。
图2.29 不分层抽样执行结果
程序说明如下。
①第13行train_test_split(test_size=6):将数据集重新洗牌,并从中随机抽样6个为测试数据,其余为训练数据。
②取得数据如图2.29所示。观察y列,类别0有1个,类别1有3个,类别2有2个,与母体比例不相同。
通过各种方式搜集到一组样本后,我们会先对样本进行探索,通常会先衡量其集中趋势(Central Tendency)及数据离散的程度(Measures of Variance),这些指标称之为描述统计量(Descriptive Statistics),较为常见的描述统计量如下。
1. 集中趋势
(1)平均数(Mean):母体以 μ 表示,样本以 x 表示。有
(2)中位数(Median):将所有样本由小排到大,以中间的样本值为准。若样本数为偶数,则取中间样本值的平均数。中位数可以避免离群值(Outlier)的影响。例如,统计平均年收入,如果加入一位超级富豪,则平均数将大幅提高,但中位数不受影响。假设有一组样本:[100, 200, 300, 400, 500],其平均数=中位数=300。现在把500改为50000,则样本为:[100, 200, 300, 400, 50000],平均数为100200,被50000影响,大幅提升,无法反映大多数数据的类型;而中位数为300,仍然不变。
(3)众数(Mode):频率发生最高的数值,以大多数的数据为主。
2. 数据离散的程度
(1)级距(Range):级距=最大值-最小值。
(2)百分位数(Percentiles)与四分位数(Quartiles):例如100为75百分位数表示有75%的样本小于100。
(3)变异数(Variance):母体以 δ 表示,样本以 s 表示。则有
3. 箱形图
可以使用箱形图(Box Plot,或称盒须图)直接显示上述相关的统计量。
下面进行实践, 请参阅02_10_基础统计.ipynb 。
(1)以美国历届总统的身高数据,计算各式描述统计量。程序代码如下:
执行结果如下。
①平均数=179.73809523809524。
②中位数=182.0。
③众数=183。
④级距=30。
⑤182cm的百分位数=47.61904761904761。
⑥变异数=7.02。
(2)以美国历届总统的身高数据,计算四分位数。程序代码如下:
执行结果如下。
①四分位数
②0.25 174.25
③0.50 182.00
④0.75 183.00
⑤1.00 193.00
执行结果如下。
①count 42.000000
②mean 179.738095
③std 7.015869
④min 163.000000
⑤25% 174.250000
⑥50% 182.000000
⑦75% 183.000000
⑧max 193.000000
(3)绘制箱形图,定义如图2.30所示。
图2.30 箱形图定义
由箱形图可观察以下特性。
①中间的线:中位数。
②中间的箱子:Q1~Q3共50%的数据集中在箱子内。
③是否存在离群值。注意:离群值发生的可能原因包括记录或输入错误、传感器失常造成量测错误或是重要的影响样本(Influential Sample),应仔细探究为何发生,不要直接剔除。一般而言,会更进一步收集更多样本,来观察离群值是否再次出现,另一方面,更多的样本也可以稀释离群值的影响力。
程序代码如下:
执行结果:如图2.31所示。
图2.31 执行结果
除了单变量的统计量,我们也会做多变量的分析,了解变量之间的关联度(Correlation),在数据探索与分析(Exploratory Data Analysis, EDA)的阶段,我们还会依据数据字段的属性进行不同面向的观察,如图2.32所示。
图2.32 数据探索与分析不同面向的观察
常用的统计图功能简要说明如下,同时也使用程序实践相关图表。
(1)直方图(Histogram):观察数据集中趋势、离散的程度、概率分布及偏态(Skewness)/峰态(Kurtosis),如图2.33所示。
图2.33 直方图
程序示例如下:
(2)饼图(Pie Chart):显示单一变量的各类别所占比例,如图2.34所示。
图2.34 饼图
程序示例如下:
(3)折线图(Line Chart):观察趋势,尤其是时间的趋势,如股价、气象温度、营收情况、运力分析等,如图2.35所示。
图2.35 折线图
程序示例如下:
(4)散点图(Scatter Chart):显示两个变量的关联,亦用于观察是否有离群值,如图2.36所示。
图2.36 散点图
程序示例如下:
(5)气泡图(Bubble Chart):散点图再加一个变量,以该变量值作为点的大小,并且做三个维度的分析,如图2.37所示。
图2.37 气泡图
程序示例如下:
(6)热图(Heatmap):显示各变量之间的关联度,便于轻易辨识出较高的关联度,如图2.38所示。输入数据须为各变量之间的关联系数(Correlation Coefficient),可使用Pandas的corr()函数计算。公式为
图2.38 热图
程序示例如下:
(7)另外,Seaborn库还提供了许多指令,使用一个指令就可以产生多张的图表,如pairplot(见图2.39所示)、facegrid等,读者可参阅Seaborn官网 [2] 。
图2.39 指令生成图表
以某市市长选举为例,当我们收集到问卷调查后,接着就可以建立模型来预测谁会当选,预测结果通常是以概率表示各候选人当选的可能性。
1. 概率的定义与定理
接下来我们先了解概率的相关术语(Terminology)及定义。
(1)实验(Experiment)或称试验(Trial):针对计划目标进行一连串的行动,如掷硬币、掷骰子、抽扑克牌、买彩票等。
(2)样本空间(Sample Space):一个实验会出现的所有可能结果,如掷硬币一次的样本空间为{正面、反面},掷硬币两次的样本空间为{正正、正反、反正、反反}。
(3)样本点(Sample Point):样本空间内任一可能的结果,如掷硬币两次会有四种样本点,抽一张扑克牌有52个样本点。
(4)事件(Event):某次实验发生的结果,一个事件可能含多个样本点。例如,抽一张扑克牌出现红色的事件,该事件就包含26个样本点。
(5)概率:某次事件发生的可能性,因此
例如,掷硬币两次出现两个正面为{正正} 1个样本点,样本空间所有样本点为{正正、正反、反正、反反}共4个,故概率等于1/4。同样地,掷硬币两次出现一正一反的概率={正反、反正} / {正正、正反、反正、反反}=1/2。
(6)条件概率(Conditional Probability)与相依性(Dependence)。
①事件独立(Independent): A 事件发生,不会影响 B 事件出现的概率,则称 A 、 B 两事件独立。例如,掷硬币两次,第一次出现正面或反面,都不会影响第二次的掷硬币结果。
②事件相依(Dependent):抽一张扑克牌,不放回,再抽第二张,概率会被第一张抽出结果所影响。例如:
● 第一张抽到红色的牌,第二张再抽到红色的概率:第一张抽到红色 红色牌剩25张,黑色牌剩26张,则第二张抽到红色概率=25 /(25+ 26)。
● 第一张抽到红色的牌,第二张再抽到黑色的概率:第二张抽到黑色概率= 26 /(25 + 26)。
上述两个概率是不相等的,故两次抽牌的事件是相依的。
③事件互斥(Mutually Exclusive): A 事件发生,就不会出现 B 事件,两事件不可能同时发生。例如,次日天气出现晴天的概率为1/2,阴天的概率为1/4,雨天的概率为1/4,次日出现晴天,就不可能出现阴天,故次日不是雨天的概率为1/2 + 1/4 = 3/4,即互斥事件的概率等于个别事件的概率直接相加。
依照上述定义,衍生的相关定理如下。
① A 、 B 事件独立,则同时发生 A 及 B 事件的概率 P ( A ∩ B )= P ( A )× P ( B )。
② A 、 B 事件互斥,则发生 A 或 B 事件的概率 P ( A ∪ B )= P ( A )+ P ( B )。
③ A 、 B 相依事件,则发生 A 或 B 事件的概率 P ( A ∪ B )= P ( A )+ P ( B )- P ( A ∩ B )。
④样本空间内所有互斥事件的概率总和等于1。
以上条件概率可以使用集合论(Set Theory)表示,如图2.40所示。
图2.40 交集(Intersection)及并集(Union)
2. 排列与组合
有些事件我们会关心发生的顺序,相反地,有些时候则不关心事件发生的顺序,计算的公式会有所不同,结果也不一样,关心事件发生的顺序称为排列(Permutation),反之,不关心事件发生的顺序称为组合(Combination)。例如,掷硬币两次的排列样本空间为{正正、正反、反正、反反},但组合样本空间为{正正、一正一反、反反},因为“正反”和“反正”都是“一正一反”。
以下列举各种案例说明排列与组合的相关计算公式及程序撰写, 请参阅02_11_概率.ipynb 。
范例1. 有三颗球,标号各为1、2、3,它们排列的事件共有几种?
程序代码如下:
执行结果:所有事件如下。
(1, 2, 3)
(1, 3, 2)
(2, 1, 3)
(2, 3, 1)
(3, 1, 2)
(3, 2, 1)
排列数=6
范例2. 有三颗球,抽出两颗球排列的事件共有几种?
程序代码如下:
(1)执行结果:所有事件如下。
(1, 2)
(1, 3)
(2, 1)
(2, 3)
(3, 1)
(3, 2)
排列数=6
(2)出现(1, 2)的概率=[(1, 2),(2, 1)] / 6 = 2/6 = 1/3。
范例3. 有三颗球,抽出两颗球组合的事件共有几种?
程序代码如下:
(1)执行结果:所有事件如下。
(1, 2)
(1, 3)
(2, 3)
组合数=3
(2)出现(1, 2)的概率=[(1, 2)] / 3 = 1/3。
范例4. 再进一步深化,袋子中有10颗不同颜色的球,抽出3颗球共有几种排列方式?
解题:排列计算公式为
式中: n 为样本点数10;
k 为3。
程序代码如下:
(1)执行结果:排列方式共720种。
(2)公式由来:抽第一次有10种选择,抽第二次有9种选择,抽第三次有8种选择,故
范例5. 掷硬币三次,出现0、1、2、3次正面的组合共有几种?
解题:因为只考虑出现正面的次数,不管出现的顺序,故属于组合的问题,因为出现的顺序不同,仍视为同一事件,故要除以 k !,公式为
式中 n 为实验次数,4次;
k 为出现正面次数,0、1、2、3次。
程序代码如下:
执行结果:
出现0次正面的组合共1种
出现1次正面的组合共3种
出现2次正面的组合共3种
出现3次正面的组合共1种
3. 二项分布
掷硬币、掷骰子、抽扑克牌一次,均假设每个单一样本点发生的概率都均等,如果不是均等,则上述的排列/组合的公式需要再考虑概率 p ,如掷硬币出现正面的概率=0.4, p k 代表出现 k 次正面的概率,(1- p ) ( n - k ) 代表出现 n - k 次反面的概率。
范例6. 掷硬币出现正面的概率为0.4,掷硬币三次,出现0、1、2、3次正面的概率为何?
程序代码如下:
(1)执行结果如下。
出现0次正面的概率=0.2160
出现1次正面的概率=0.4320
出现2次正面的概率=0.2880
出现3次正面的概率=0.0640
(2)使用SciPy comb组合函数,执行结果同上。程序代码如下:
(3)直接使用SciPy的统计模块binom.pmf(k,n,p),执行结果同上。程序代码如下:
范例7. 试算今彩539平均回报率 [3] 。
基本玩法:从39个号码中选择5个号码。各奖项的中奖方式见表2.1。
表2.1 中奖方式
各奖项金额见表2.2。
表2.2 各奖项金额
解题:
(1)假定每个号码出现的概率是相等的(1/39),固定选5个号码,故计算概率时不需考虑 p 。
(2)计算从39个号码选5个号码,总共有几种组合。
程序代码如下:
执行结果如下:
从39个号码选5个号码,总共有575757种
头等奖号码的个数:1
贰等奖号码的个数:170
叁等奖号码的个数:5610
肆等奖号码的个数:59840
(3)计算平均中奖金额。程序代码如下:
执行结果如下:
平均中奖金额:27.92
今彩539平均回报率=-44.16%
注意:排列与组合的理论看似简单,然而实践上的应用千变万化,建议读者多找一些案例实践,才能运用自如。
概率分布(Distribution)结合了描述统计量及概率的概念,对数据做进一步的分析,希望推测母体是呈现何种形状的分布,如正态分布、均匀分布、泊松(Poisson)分布或二项分布等,并且依据样本推估概率分布相关的母数,如平均数、变异数等。有了完整概率分布信息后,就可以进行预测、区间估计、假设检定等。
概率分布的种类非常多,如图2.41所示。这里我们仅介绍几种常见的概率分布。
图2.41 各种概率分布及其关联 [11]
首先介绍几个专有名词。
(1)概率密度函数:发生各种事件的概率。
(2)累积分布函数(Cumulative Distribution Function, CDF):等于或低于某一观察值的概率。
(3)概率质量函数(Probability Mass Function, PMF):如果是离散型的概率分布,则PDF改称为PMF。
接下来看几个常见的概率分布。
1. 正态分布
因为正态分布(Normal Distribution)是由Carl Friedrich Gauss提出的,因此又称为高斯分布(Gauss Distribution),世上大部分的事件都属于正态分布,如考试的成绩,考低分及高分的学生人数会比较少,中等分数的人数会占大部分。再如,全年的温度、业务员的业绩、一堆红豆的重量等。正态分布的概率密度函数定义为
两个重要参数介绍如下。
(1) µ :平均数,是全体样本的平均数。
(2) δ :标准差,描述全体样本的离散的程度。
概率密度函数可简写成 N ( µ , δ )。
我们在“2-3-5积分”节已经介绍过相关内容,并撰写了一些程序,这里直接使用SciPy的统计模块做一些实操。
以下程序请参考02_12_正态分布.ipynb。
范例1. 使用SciPy绘制概率密度函数(PDF)。
程序代码如下:
执行结果:如图2.42所示。
图2.42 概率密度函数
范例2. 绘制正态分布的±1、±2、±3倍标准差区间及其概率密度分布。
程序代码如下:
(1)执行结果:如图2.43所示。
图2.43 正态分布±1倍标准差区间及概率密度分布
(2)将第5行改为(-2, 2)、(-3, 3),结果如图2.44所示。
图2.44 正态分布±2、±3倍标准差区间及概率密度分布
范例3. 使用SciPy绘制累积分布函数(CDF)。
程序代码如下:
第12行从PDF改为CDF。
执行结果:如图2.45所示。
图2.45 累积分布函数
有一个Student's t分布,很像正态分布,常被用来做假设检定(Hypothesis Test),SciPy也有资源。
范例4. 使用SciPy绘制Student ' s t分布的概率密度函数(PDF),并与正态分布比较。
第13行为t为Student's t分布。
执行结果:如图2.46所示。
图2.46 Student's t分布的概率密度函数
图2.46中,较粗的线(红色)t为Student's t分布。
2. 均匀分布
均匀分布也是一种常见的分布,通常是属于离散型的数据,所有样本点发生的概率都相同,如掷硬币、掷骰子,出现每一面的概率都一样。均匀分布的概率密度函数为
以下程序请参考02_13_其他分布.ipynb。
范例5. 绘制掷骰子的概率密度函数。
程序代码如下:
(1)执行结果:如图2.47所示。
图2.47 范例5执行结果
(2)可以使用随机数模拟,程序代码为np.random.randint(1, 6)。
3. 二项分布
顾名思义,二项分布就是只有两种观察值,如成功/失败、正面/反面、有/没有等,类似的分布有以下三种。
(1)伯努利(Bernoulli)分布:做 一次二分类 的实验,概率密度函数为
f ( x )= p x (1− p ) 1− x , x =0或1
式中: p 为成功的概率。
(2)二项(Binomial)分布:做 多次二分类 的实验。概率密度函数为
式中: n 为实验次数; k 为成功次数; p 为成功的概率。
(3)多项(Multinomial)分布:做 多次多分类 的实验。概率密度函数为
范例6. 绘制掷硬币的概率密度函数。
程序代码如下:
执行结果:如图2.48所示。
图2.48 伯努利分布执行结果
范例7. 绘制掷硬币10次的概率密度函数。
程序代码如下:
执行结果:如图2.49所示。
图2.49 二项分布执行结果
范例8. 绘制掷骰子10次的概率。
程序代码如下:
(1)执行结果:0.011520000000000013。
(2)若掷骰子60次,1~6点各出现10次的概率,执行结果为4.006789054168284e-06。概率会小很多,这是由于掷骰子60次出现的样本空间会有更多的组合。程序代码如下:
(3)若掷骰子1000次,给定1~6点各出现的次数。程序代码如下:
执行结果:如图2.50所示。
图2.50 执行结果
4. 泊松分布
泊松分布经常运用在“给定的时间内发生 k 次事件”的概率分布函数,我们得以建立一套顾客服务的模型,用来估计一个服务柜台的等候人数,进而计算出需要安排几个人服务柜台,来达成特定的服务水平(Service Level Agreement, SLA),其应用场景非常广,如售票柜台、便利商店结账柜台、客服中心的电话线数、依个人发生车祸的次数决定车险定价等。
泊松分布的概率密度函数为
式中: λ 为平均发生次数;e为自然对数; k 为发生次数。
范例9. 绘制各种 λ 值的泊松分布概率密度函数。
程序代码如下:
执行结果:如图2.51所示。
图2.51 泊松分布执行结果
根据维基百科的定义 [5] ,假设检定是推论统计中用于检定统计假设的一种方法。而统计假设可以通过观察一组随机变量的模型进行检定。通常判定一种新药是否有效,我们并不希望因为随机抽样的误差,而造成判定错误,所以,治愈人数的比例必须超过某一显著水平,才能够认定为具有医疗效果。在谈到假设检定之前,我们先来了解什么是置信区间(Confidence Interval)。
1. 置信区间
之前我们谈到平均数、中位数都是以单一值表达样本的集中趋势,这称为点估计,不过,这种表达方式并不精确。以正态分布而言,如图2.52所示的概率密度函数,样本刚好等于平均数的概率也只有0.4,因此,我们改以区间估计会是比较稳健的做法。例如,估计95%的样本会落在特定区间内,这个区间我们便称为置信区间(Confidence Interval)。
图2.52 置信区间示例
以正态分布而言,1倍的标准差置信水平约为68.3%,2倍的标准差置信水平约为95.4%,3倍的标准差置信水平约为99.7%,如图2.53所示。
图2.53 正态分布的置信水平
以业界常用的95%置信水平为例,就是“我们确信95%的样本会落在(-1.96 δ , 1.96 δ )区间内”,一般医药有效性的检定也是以此概念表达的,如疫苗等。
以下程序请参考02_14_置信区间.ipynb。
范例1. 以美国历届总统的身高数据,计算各式描述统计量及95%的置信区间。
程序代码如下:
(1)执行结果:平均数=179.74,标准差=7.02,置信区间=(165.71, 193.77)。
(2)如上所述,如果单纯以平均数179.74说明美国总统的身高并不完整,因为其中有多位总统的身高低于170cm,因此若再加上“95%美国总统的身高介于(165.71,193.77)”,则会让信息更加完整。
范例2. 利用随机数生成正态分布的样本,再使用SciPy的统计模块(stats)计算置信区间。
解题:
(1)利用随机数生成10000个样本,样本来自正态分布 N (5, 2),即平均数为5,标准差为2,使用norm.interval()计算置信区间,参数分别为置信水平、平均数、标准差。程序代码如下:
执行结果:置信区间=(1.0487, 8.9326),约为( μ –1.96 δ , μ +1.96 δ )。
(2)绘图。程序代码如下:
执行结果:如图2.54所示。
图2.54 95%置信区间执行结果
范例3. 利用随机数生成二项分布的样本,再使用SciPy的统计模块(stats)计算置信区间。
程序代码如下:
(1)执行结果:平均数=0.5039,标准差=0.4999847897686489,置信区间=(4941.0,5137.0),约为( μ –1.96 δ , μ +1.96 δ )。
(2)二项分布标准差公式= ,可以验算( p *(1- p ))**.5 = 0.5,随机数的样本与理论值相差不远。
2. 中心极限定理
中心极限定理(Central Limit Theorem, CLT)是指当样本数越大时,不管任何概率分布,每批样本的平均数之概率分布会近似于正态分布。因此,我们就可以使用正态分布估计任何样本平均数之概率分布的置信区间或做假设检定。
然而,因为每一个样本是一批的数据,假设个数为 n ,则每批样本的平均数所形成的标准差不是原样本标准差,而是
例如,我们有1000批的样本,每批的样本有10个数据,则每批样本的平均数所形成的标准差为
又如,样本来自二项分布,标准差= ,则每批样本的平均数所形成的标准差为
平均数的标准差称为标准误差(Standard Error),它有别于原样本的标准差。
以下程序请参考02_15_中心极限定理.ipynb。
范例4. 利用随机数生成二项分布的10000批样本,每批含100个数据,请计算平均数、标准误差,并绘图验证是否近似于正态分布。
程序代码如下:
(1)执行结果:平均数=0.5003,标准差=0.04982,置信区间=(4906.0, 5102.0),如图2.55所示。
图2.55 中心极限定理执行结果
(2)平均数与原样本的平均数几乎相同(0.5003≈0.5)。
(3)标准差与理论值也很接近(0.04982≈0.05),理论值=(p*(1-p)/ trials)**.5=0.04982。
(4)直方图近似于正态分布。
3. 假设检定
我们设定特定的显著水平(Significance Level),以 α 表示,例如5%,或者相对的置信水平95%,如果样本的平均数落在置信区间之外,我们就说新药是显著有效,这样的判定只有5%概率是错误的,这种检定方法就称为假设检定(Hypothesis Testing)。当然,我们也可以用较严谨的显著水平(0.3%),即三倍标准差,检定新药是否显著有效。
假设检定依检定范围可分为单边(One-side)、双边(Two-side),依样本的设计可分为单样本(Single-sample)、双样本(Two-sample),依检定的统计量,可分为平均数检定或标准差的检定。
单边检定只关心概率分布的一边,例如,A班的成绩是否优于B班( µ A > µ B ),为右尾检定(Right-tail Test),病毒核酸检测(PCR),ct值30以下即视为确诊( µ <30),为左尾检定(Left-tail Test);双边检定则关心两边,如A班的成绩是否与B班有显著差异( µ A ≠ µ B )。
单样本检定只有一份样本,如抽样一群顾客,调查是否喜欢公司特定的产品;双样本则有两份样本互相作比较,如A班的成绩是否优于B班,或者新药有效性测试,通常会将实验对象分为两组,一组为实验组(Treatment Group),让他们服用新药,另一组为对照组,又称控制组(Control Group),让他们服用安慰剂,这种设计又称为A/B Test。
检定时会有两个假设,原假设(Null Hypothesis)及备择假设(Alternative Hypothesis)。
(1)原假设H 0 : µ =0。
(2)备择假设H 1 : µ ≠0,H 1 也可以写成H A 。
备择假设是原假设反面的条件,通常备择假设是我们希望的结果,因此,检定后,如果原假设为真时,我们会使用 不能拒绝 (fail to reject)原假设,而不会说原假设成立。
以下就以实例说明各式的检定。
以下程序请参考02_16_假设检定.ipynb。
范例5. 特朗普的身高190cm,是否与历届的美国总统平均身高有显著的不同?假设历届的美国总统身高为正态分布,以显著水平5%检定。
解题:显著水平5%时,置信区间为( µ –1.96 δ , µ +1.96 δ )。程序代码如下:
(1)执行结果:平均数=179.74,标准差=7.02,置信区间=(165.99, 193.49),特朗普的身高190cm在置信区间内,表示特朗普身高并不显著比历届的美国总统高。
(2)从图2.56可以看出,特朗普的身高190cm在左右两条虚线(置信区间)之间。
图2.56 特朗普身高的假设检验结果
范例6. 要调查顾客是否喜欢公司新上市的食品,公司进行问卷调查,取得客户的评价,范围介于0~10分,已知母体平均分数为5分,标准差为2分,请使用假设检定确认顾客是否喜欢公司新上市的食品。
解题:
原假设H 0 : µ ≤5
备择假设H 1 : µ > 5
检验平均数,我们会使用 t 检定,当样本大于30时,会近似于正态分布。 t 统计量为
(1)使用随机数模拟问卷调查结果。程序代码如下:
执行结果:最小值:1,最大值:9,平均数:5.46。
(2)使用stats.ttest_1samp()做假设检定,参数alternative='greater'为右尾检定,SciPy需v1.6.0以上才支持alternative参数。程序代码如下:
①执行结果: t 统计量2.0250, p 值0.024。90%的置信区间约在( µ –1. 645 δ , µ +1. 645 δ )之间,因 t 统计量(2.0250)>1.645,故备择假设为真,确认顾客喜欢公司新上市的食品。使用 t 统计量比较,需考虑单/双尾检定,并比较不同的值,有点麻烦,专家会改用 p 值与显著水平比较,若小于显著水平,则备择假设为真,反之,不能拒绝原假设,因为 p 值的计算公式=1-CDF(累积分布函数),如图2.57所示。
图2.57 问卷调查结果的假设检验结果
②图2.57实线为“母体平均数+ t 统计量 x 母体标准差”,旁边虚线为置信区间的右界,前者大于后者,表示效果显著,有明显差异。
③若使用左尾检定,也是呼叫stats.ttest_1samp(),参数alternative='less'。
④若使用双尾检定,也是呼叫stats.ttest_1samp(),删除参数alternative即可。
范例7. 要检定新药有效性,将实验对象分为两组,一组为实验组,让他们服用新药,另一组为对照组,又称控制组,让他们服用安慰剂,检验两组疾病复原状况是否有明显差异。
解题:
原假设H 0 : µ 1 = µ 2
备择假设H 1 : µ 1 ≠ µ 2
(1)使用随机数模拟复原状况。程序代码如下:
执行结果:控制组平均数:66.04,实验组平均数:66.46。
(2)使用stats.ttest_ind()作假设检定,参数放入两组数据及alternative='greater'表示右尾检定。程序代码如下:
①执行结果: t 统计量2.2390, p 值0.0129, p 值若小于显著水平(0.05),则备择假设为真,表示新药有显著疗效,如图2.58所示。
图2.58 新药有效性的假设检验结果
②图2.58实线为“母体平均数+ t 统计量 x 母体标准差”,旁边虚线为置信区间的右界,前者大于后者,表示效果显著,有明显差异。
范例8. 另外有一种配对检定(Paired Tests),如学生同时参加期中及期末考试,我们希望检验期末时学生是否有显著进步,因两次考试都是以同一组学生做实验,故称为配对检定。
解题:
原假设H 0 : µ 1 = µ 2
备择假设H 1 : µ 1 ≠ µ 2
(1)使用随机数模拟学生两次考试成绩。程序代码如下:
执行结果:期中考:60.14,期末考:60.90。
(2)使用stats.ttest_rel()做假设检定,参数放入两组数据及alternative='greater'表右尾检定。程序代码如下:
①执行结果: t 统计量1.0159, p 值0.1561, p 值大于显著水平(0.05),备择假设为假,表示学生成绩没有显著进步,如图2.59所示。
图2.59 学生成绩的配对检验结果
②图2.59实线为“母体平均数+ t 统计量 x 母体标准差”,右侧虚线为置信区间的右界,前者小于后者,表示效果不显著。
以上基础统计的介绍,环环相扣,以某市市长选举为例介绍如下,详见图2-60。
(1)抽样:使用抽样方法,从有投票权的市民中随机抽出一批市民意见作为样本。
(2)计算描述统计量:推估母体的平均数、标准差、变异数等描述统计量,描绘出大部分市民的想法与差异。
(3)概率:推测候选人当选的概率,再进而推估概率分布函数。
(4)假设检定:依据概率分布函数,进行假设检定,推论候选人当选的可信度或确定性,这就是一般古典统计的基础流程,与机器学习以数据为主的预测相辅相成,可彼此验证对方。
图2.60 上述介绍的基础统计范围及其关联
线性规划是运筹学(Operations Research, OR)中非常重要的领域,它是给定一些线性的限制条件,求取目标函数的最大值或最小值。
以下程序请参考02_17_线性规划.ipynb。
范例1. 最大化目标函数 z = 3 x +2 y
限制条件:2 x + y ≤100
x + y ≤80
x ≤40
x ≥0, y ≥0
解题:
(1)先画图,涂色区域(黄色)为可行解(Feasible Solutions),即符合限制条件的区域,如图2.61所示。
图2.61 确定可行解
程序代码如下:
(2)上述线性规划求解可使用单形法(Simplex Method)。上述问题比较简单,凸集合的最佳解发生在可行解的顶点(Vertex),所以依上图,我们只要求每一个顶点对应的目标函数值,比较并找到最大的数值即可。不过,深度学习的变数动辄几百个,且神经层及神经元又很多,无法使用单形法求解,而是采取优化的方式,逐渐逼近找到近似解。因此,我们只运用程序来解题,不介绍单形法的原理。首先安装pulp包。
(3)以pulp包求解,定义目标函数及限制条件。程序代码如下:
执行结果如下,显示定义内容:
(4)调用model.solve()求解。程序代码如下:
执行结果如下,当 x =20、 y =60时,目标函数最大值=180。
范例2. 实例应用,运用线性规划来安排客服中心各时段的人力配置,请参考文献 [6] 。
以上是简单的线性规划,还有其他的模型,如整数规划(Integer Programming)、二次规划(Quadratic Programming)、非线性规划(Non-linear Programming)。
普通最小二乘法(Ordinary Least Squares, OLS)、最大似然估计法(Maximum Likelihood Estimation, MLE)是常见的优化和估算参数值的方法,如回归系数、正态分布的平均数/变异数。笔者戏称两者是优化器求解的倚天剑与屠龙刀,许多算法问题凭借这两把利器便可迎刃而解。
在“2-3-4简单线性回归求解”章节,利用一阶导数等于0的特性,对简单线性回归( y = wx + b )求取最佳解,使用的就是普通最小二乘法,这里再强化一下,使用矩阵运算,让解法可应用到多元回归、深度学习。
y = w 1 x 1 + w 2 x 2 + w 3 x 3 +L+ w n x n + b
首先定义线性回归的目标函数(Objective Function)或称损失函数(Loss Function)为均方误差(MSE),公式为误差平方和(SSE),除以样本数( n ),其中
MSE当然越小越好,所以它是一个最小化的问题。
(1)目标函数 。
①MSE = SSE / n ,我们忽略常数 n ,可以只考虑SSE,有
②以矩阵表示线性回归 y = wx 。其中 b = b * x 0 ,故将 b 也可视为 w 的一个项目。
③以矩阵表示SSE:
( y - wx ) T *( y - wx )
y T y - 2 w T x T y + w T x T wx
(2)对 w 偏微分得
结合矩阵、微分,我们就能够轻松求出线性回归的系数,此原理是神经网络优化求解的基石。
以下程序请参考02_18_普通最小二乘法.ipynb。
范例1. 以普通最小二乘法建立线性回归模型,预测波士顿(Boston)房价。
解题:
(1)依上述公式计算回归系数。程序代码如下:
执行结果:W=[-1.08011358e-01 4.64204584e-02 2.05586264e-02 2.68673382e+00
-1.77666112e+01 3.80986521e+00 6.92224640e-04 -1.47556685e+00
3.06049479e-01 -1.23345939e-02 -9.52747232e-01 9.31168327e-03
-5.24758378e-01 3.64594884e+01]
(2)计算相关效果衡量指标。程序代码如下:
(3)以Sklearn库验证。程序代码如下:
执行结果:与公式计算一致,验证无误。
范例2. 使用SciPy以普通最小二乘法计算函数 x 2 +5的最小值。
(1)先对函数绘图。程序代码如下:
执行结果:如图2.62所示。
图2.62 函数绘图
(2)调用scipy.optimizea模块的leastsq()函数进行优化求解。
①在leastsq()函数中,第一个参数是求解的函数;第二个参数是起始点;leastsq是采逼近法,而非纯数学公式求解,nfev显示它经过22次执行周期,才找到最小值5(fvec),当时x=1.72892379e-05≈0。
②执行结果如下:
③起始点可设为任意值,通常采用随机数或是直接给0。指定值设定不佳的话,仍然可以找到最佳解,不过,需要较多次的执行周期,也就是所谓的较慢收敛(Convergence)。
④当面对较复杂的函数或较多的变量时,我们很难单纯运用数学去求解,因此,逼近法会是一个比较实用的方法,深度学习的梯度下降法就是一个典型的例子,后面章节我们将使用TensorFlow程序来进行说明。
最大似然估计法(MLE)也是估算参数值的方法,目标是找到一个参数值,使出现目前事件的概率最大。如图2.63所示,圆点是样本点,曲线是四个可能的正态概率分布(平均数/变异数不同),我们希望利用最大似然估计法找到最适配(Fit)的一个正态概率分布。
图2.63 最大似然估计法示意图
下面使用最大似然估计法求算正态分布的参数 µ 及 δ 。
(1)正态分布
(2)假设来自正态分布的多个样本互相独立,联合概率就等于个别的概率相乘。即
(3)目标是求取参数值 µ 及 δ ,最大化概率值,即最大化目标函数。即
(4) N 个样本概率全部相乘,不易计算,通常我们会取对数,变成一次方,所有的概率值取对数后,大小顺序并不会改变。即
(5)取对数后,化简为
(6)对 µ 及 δ 分别偏微分,得
(7)一阶导数为0时,有最大值,得到目标函数最大值下的 µ 及 δ 。所以,使用最大似然估计法算出的 µ 及 δ 就如同我们常见的公式。有
以下程序请参考02_19_最大似然估计法.ipynb。
范例1. 如果样本点 x =1,计算来自正态分布 N (0,1)的概率。
(1)使用NumPy计算概率密度函数(PDF)。程序代码如下:
执行结果:0.24。
(2)也可以使用scipy.stats模块计算。程序代码如下:
执行结果:0.24。
(3)绘制概率密度函数。程序代码如下:
执行结果:如图2.64所示。
图2.64 绘制概率密度函数
范例2. 如果有两个样本点 x =1、3,来自正态分布 N (1,1)及 N (2,3)的可能性,哪一个比较大?
假设两个样本是独立的,故联合概率为两个样本概率相乘,使用scipy.stats模块计算概率。程序代码如下:
执行结果如下,表明来自 N (1,1)可能性比较大。
范例3. 如果有一组样本,计算来自哪一个正态分布 N ( μ , δ )的概率最大,请依上面证明计算 μ 、 δ 。
(1)对正态分布的PDF取对数(log)。程序代码如下:
执行结果如下,Jupyter Notebook可显示LaTeX数学式。有
(2)带入样本数据。程序代码如下:
(3)上述函数使用diff()各对平均数、变异数偏微分。程序代码如下:
(4)使用simplify()简化偏导数。程序代码如下:
(5)令一阶导数为0,有最大值,可得到联立方程式。程序代码如下:
执行结果:[(3.62, -1.57),(3.62, 1.57)]
(6)使用NumPy验证。程序代码如下:
执行结果:(3.62, 1.57),与MLE求解的结果相符。
有了以上的基础后,我们将综合微积分、矩阵、概率等数学知识,对神经网络求解,这是进入深度学习领域非常重要的概念。
神经网络是深度学习最重要的算法,它主要是仿真生物神经网络的传导系统,希望通过层层解析,归纳出预测的结果,如图2.65所示。
图2.65 生物神经网络的传导系统
生物神经网络中表层的神经元接收到外界信号,归纳分析后,再通过神经末梢,将分析结果传给下一层的每个神经元,下一层神经元进行相同的动作,再往后传导,最后传至大脑,大脑做出最后的判断与反应,如图2.66所示。
图2.66 单一神经元结构
于是,AI科学家将上述生物神经网络简化成如图2.67所示的网络结构:
图2.67 AI神经网络
AI神经网络最简单的连接方式称为完全连接(Full Connected, FC),即每一神经元均连接至下一层的每个神经元,因此,我们可以把第二层以后的神经元均视为一条回归线的 y ,它的特征变量 x 就是前一层的每一个神经元。如图2.68所示, y 1 、 z 1 两条回归线表示为
y 1 = w 1 x 1 + w 2 x 2 + w 3 x 3 + b
z 1 = w 1 y 1 + w 2 y 2 + w 3 y 3 + b
所以,简单地讲,一个神经网络可视为多条回归线组合而成的模型。
图2.68 一个神经网络模型
以上的回归线是线性的,为了支持更通用性的解决方案(Generic Solution),模型还会乘上一个非线性的函数,称为激励函数(Activation Function),期望也能解决非线性的问题,如图2.69所示。因激励函数并不能明确表达原意,因此下面直接以Activation Function表示。
图2.69 激励函数
如果不考虑激励函数,每一条线性回归线的权重及偏差可以通过普通最小二乘法求解,请参阅2-6-1说明,但如果乘上非线性的Activation Function,就比较难用单纯的数学求解了,因此,学者就利用优化(Optimization)理论,针对权重、偏差各参数分别偏微分,沿着切线(即梯度)逐步逼近,找到最佳解,这种算法就称为梯度下降法(Gradient Descent)。
有人用了一个很好的比喻,“当我们在山顶时,不知道下山的路,于是,就沿路往下走,遇到叉路时,就选择坡度最大的叉路走,直到抵达平地为止”,所以梯度下降法利用偏微分(Partial Differential)求算斜率,依斜率的方向,一步步地往下走,逼近最佳解,直到损失函数没有显著改善为止,这时我们就认为已经找到最佳解了,如图2.70所示。
图2.70 梯度下降法示意图
上图及以下程序请参考02_20_梯度下降法.ipynb。
梯度其实就是斜率,单变量回归线的权重称为斜率,多变数回归线时,需分别作偏微分求取权重值,此时就称为梯度。下面我们先针对单变量求解,示范如何使用梯度下降法(Gradient Descent)求取最小值。
范例1. 假定损失函数 f ( x )= x 2 ,而非MSE,请使用梯度下降法求取最小值。
注意:损失函数又称为目标函数或成本函数,在神经网络相关文献中大多称为损失函数,本书从善如流,以下将统一以损失函数取代目标函数 。
(1)定义函数(func)及其导数(dfunc)。程序代码如下:
(2)定义梯度下降法函数,反复更新 x ,更新的公式如下,后面章节我们会推算公式的由来。
新的 x =目前的 x -学习率(learning_rate)*梯度(gradient)。程序代码如下:
(3)设定起始点、学习率(lr)、执行周期数(epochs)等参数后,呼叫梯度下降法求解。程序代码如下:
①执行结果:如图2.71所示。
图2.71 梯度下降法执行结果
②每一执行周期的损失函数如图2.72所示。随着 x 变化,损失函数逐渐收敛,即前后周期的损失函数差异逐渐缩小,最后当 x =0时,损失函数 f ( x )等于0,为函数的最小值,与普通最小二乘法(OLS)的计算结果相同。
[5. 2, 0.8, 0.32, 0.13, 0.05, 0.02, 0.01, 0, 0, 0, 0, 0, 0, 0, 0]
图2.72 每一执行周期的损失函数
③如果改变起始点(x_start)为其他值,如-5,依然可以找到相同的最小值。
范例2. 假定损失函数 f ( x )= 2 x 4 - 3 x - 20,请使用梯度下降法求取最小值。
(1)定义函数及其微分。程序代码如下:
(2)绘制损失函数。程序代码如下:
执行结果:如图2.73所示。
图2.73 梯度下降法损失函数
(3)梯度下降法函数(GD)不变,执行程序,如果学习率不变(lr = 0.3),会出现错误信息“Result too large”,原因是学习率过大,梯度下降过程错过最小值,往函数左方逼近,造成损失函数值越来越大,最后导致溢位。程序代码如下:
(4)修改学习率(lr = 0.001),同时增加执行周期数(epochs = 15000),避免还未逼近到最小值程序就先提早结束。程序代码如下:
执行结果:当 x =0.51时,函数有最小值,如图2.74所示。
图2.74 修改学习率后执行结果
观察上述范例,不管函数为何,我们以相同的梯度下降法(GD函数)都能够找到函数最小值,最重要的关键是 x 的更新公式
新的 x =目前的 x -学习率(learning_rate)×梯度(gradient)
接着我们会说明此公式的由来,也就是神经网络求解的精华所在。
神经网络求解是一个正向传导与反向传导反复执行的过程,如图2.75所示。
图2.75 神经网络权重求解过程
(1)由于神经网络是多条回归线的组合,建立模型的主要任务就是计算出每条回归线的权重( w )与偏差( b )。
(2)依上述范例的逻辑,一开始我们指定 w 、 b 为任意值,建立回归方程式 y = wx + b ,将特征值( x )带入方程式,可以求得预测值( ),进而计算出损失函数,例如 MSE ,这个过程称为正向传导(Forward Propagation)。
(3)透过最小化 MSE 的目标和偏微分,可以找到更好的 w 、 b ,并依学习率来更新每一层神经网络的 w 、 b ,此过程称之为反向传导(Backpropagation)。这部分可以借由微分的链法则(Chain Rule),依次逆算出每一层神经元对应的 w 、 b ,公式为
W t+1 = W t -学习率×梯度
式中: [后续证明];学习率为优化器事先设定的固定值或动能函数。
(4)重复(2)(3)步骤,一直到损失函数不再有明显改善为止。
梯度公式证明如下。
(1)损失函数 ,因 n 为常数,故仅考虑分子,即 SSE 。
(2)
(3)以矩阵表示,SSE= y 2 −2 ywx + w 2 x 2
(4)
(5)同理,
(6)为了简化公式,常把系数2拿掉。
(7)最后公式为
调整后权重=原权重+(学习率*梯度)
(8)有些文章将梯度负号拿掉,公式就修正为
调整后权重=原权重-(学习率*梯度)
以上是以MSE为损失函数时的梯度计算公式,若是其他损失函数,梯度计算结果也会有所不同,如果再加上Activation Function,梯度公式计算就更加复杂了,好在深度学习框架均提供自动微分、计算梯度的功能,我们就不用烦恼了。后续有些算法会自定义损失函数,会因此产生意想不到的功能,如风格转换(Style Transfer)可以合成两张图像,生成对抗网络(Generative Adversarial Network, GAN),可以产生几可乱真的图像。也因为如此关键,我们才花费了这么多的篇幅介绍梯度下降法。
基础的数学与统计就介绍到此,告一段落,下一章起,我们将开始以TensorFlow实现各种算法,并介绍相关的应用。