机器学习Chapter-1（线性模型）

最新推荐文章于 2024-07-22 16:16:12 发布

DFann

最新推荐文章于 2024-07-22 16:16:12 发布

阅读量5.8k

点赞数 13

分类专栏：机器学习文章标签：机器学习 python 线性模型逻辑回归线性判别分析

本文链接：https://blog.csdn.net/u011974639/article/details/77102663

版权

机器学习专栏收录该内容

5 篇文章 2 订阅

订阅专栏

引言
- 先验条件
- 好了下面开始正文本节的主角是线性模型
线性模型理论
Python实战线性模型
参考资料

hello world!

引言

前段时间一直在看深度学习和TensorFlow的教程，在学习的过程中，深感自己的数学基础较差，尤其是在对数据做预处理等操作时，看着书上那一堆公式，整个人都是懵逼的。看着桌上三本机器学习相关的书：西瓜书《机器学习》、《Python大战机器学习》、《机器学习实战》，思索再三，我决定接下来三个月时间努力啃完它们，顺带写blog，算是一边记笔记一边督促自己坚持下去。

先验条件

在啃这三本书之前，我已经看了一段时间的Andrew Ng在Coursera上的机器学习视频，也看了一点台大的《机器学习基石》.也建议大家先看这些视频，了解一下机器学习中的基本概念.

好了，下面开始正文，本节的主角是线性模型.

线性模型理论

概述

相信看过Andrew Ng在Coursera上关于机器学习的教程的同学都知道（没看过也没关系）,该教程讲的第一个模型就是线性模型(预测房价)。这里我们也以房价预测做例子。打个比方：我们现在要预测某地方的房价price，我们已有的信息是住房的面积size和离市区中心的距离s，我们可以列出一个定性的关系式(关系式根据是生活经验,当然这只是一个假设式):

预 测 房 价 : p r i c e = α * s i z e + β * s . (其 中 α 和 β 为 预 测 模 型 的 参 数)

$预测房价: price = \alpha*size + \beta*s .(其中\alpha和\beta为预测模型的参数)$

上式是基于size和s对于房价的预测，在实际问题中，我们可能多种数据，例如本地房产均价、周围教育设施、体育馆、住房楼层、住房受光面积等，这些因素都会影响房价。

现在，我们拓展一下，如果我们给定了d个属性(size和s就是两种属性)描述的样本如下:

x = (x 1; x 2; . . .; x d); 其 中 x i 是 x 在 第 i 个 属 性 上 的 取 值

$\boldsymbol{x} = (x_1;x_2;...;x_d) ;其中x_i是\boldsymbol{x} 在第i个属性上的取值$

线性模型试图学习到一个通过属性的线性组合来进行预测的函数，即

f (x) = ω 1 x 1 + ω 2 x 2 + . . . + ω d x d + b; (b 为 偏 置 单 元)

$f(\boldsymbol{x}) = \omega_1x_1+ \omega_2x_2+...+\omega_dx_d+b; (b为偏置单元)$

我们可以用向量形式表示

f (x) = ω T x + b;

$f(\boldsymbol{x}) = \boldsymbol{\omega^Tx} +b;$

其 中 ω = (ω 1; ω 2; . . .; ω d), ω 和 b 是 学 习 得 到 的 ， 确 定 后 模 型 就 确 定 了

$其中\boldsymbol{\omega}=(\omega_1;\omega_2;...;\omega_d) ,\boldsymbol{\omega}和b是学习得到的，确定后模型就确定了$

上面的式子就是线性模型的一般表达式了，可以看到线性模型的形式很简洁，而且易于建模。线性模型蕴含着机器学习中一些重要的基本思想.许多功能更为强大的非线性模型可在线性模型的基础上通过引入层次结构或高维映射而得。线性模型还有多种推广形式,常见的有广义线性模型：逻辑回归、岭回归等。接下来我们将介绍多种线性回归的基本思想、优缺点以及如何使用Python实现。

线性回归

1. 什么是回归？

回归分析本质上是一个函数估计问题，通俗的来说，就是找出因变量和自变量之间的因果关系。回归分析的因变量应该是连续变量(房价是连续的，或者预测某些病症发病的概率值);若因变量为离散变量(预测图片分类等),则问题就是分类问题，回归分析是一个有监督学习问题。

2.什么是线性回归?

线性回归(linear regression)是一种回归分析技术，例如在给定数据集 $D$ 。线性回归试图学习到一个线性模型以尽可能准确地预测实值输出标记。通过在数据集上建立线性模型，建立代价函数(loss function),最终以优化代价函数为目标确定模型参数 $w和b$ ,从而得到模型用以后续的预测。整个线性回归算法流程如下:

这里写图片描述

3.如何确定参数w和b?

根据已知的数据集D

D = ((x 1, y 1); (x 2, y 2); . . .; (x m, y m)) 其 中 x i = (x i 1; x i 2; . . .; x i d), y i \in R

$\boldsymbol{D} = ((x_1,y_1);(x_2,y_2);...;(x_m,y_m)) 其中x_i=(x_{i1};x_{i2};...;x_{id}) ,y_i∈R$

(x i 为 对 于 的 属 性 数 据 ， y i 为 对 于 的 标 签 值)

$(x_i为对于的属性数据，y_i为对于的标签值)$
来计算参数w和b，对于给定的样本数据

xi $x_i$ 其预测值为

yi~=wxi+b $\tilde{y_i}=\boldsymbol{w}x_i+b$ ，这时候我们需要一个可以衡量预测值与数据集真实值之间差距的函数，我们称之为代价函数，可想而知如果代价函数在整个数据集上的值越小，就代表的预测值与真实值越接近，而这就是我们想要的，所以我们的优化目标就是让代价函数的值最小。

这里代价函数我们选用均方误差(square loss)，即在样本 $x_i$ 上预测值为 $\tilde{y_i}$ ，故均方误差为 $(\tilde{y_i}-y_i)^2$ ，则在整个数据集D上，模型的代价函数为

L (f) = \sum i = 1 D (y i ~ - y i) 2 = \sum i = 1 D (w x i + b - y i) 2

$L(f) = \sum_{i=1}^D(\tilde{y_i}-y_i)^2=\sum_{i=1}^D(\boldsymbol{w}x_i+b-y_i)^2$
我们的要确定的参数为:

(w *, b *) = a r g min w, b (\sum i = 1 D (w x i + b - y i) 2)

$(\boldsymbol{w^*},b^*)=arg\min_{\boldsymbol{w},b}(\sum_{i=1}^D(\boldsymbol{w}x_i+b-y_i)^2)$

这样的基于均方误差最小化来进行模型求解的方法称为“最小二乘法”(least square method,LSM)。有趣的是当年数学王子高斯在推测行星运动时用的就是这个最小二乘法。

到此，我们的优化目标已经确定了，接下来就是求解问题了。

4.如何求解w和b?

求解 $w$ 和 $b$ 使 $L（f）$ 最小化的的过程，称为线性回归模型的最小二乘“参数估计”(parameter estimation).由凸优化知识可知， $L_{(w,b)}$ 是关于 $w$ 和 $b$ 的凸函数(证明略过)，在数学上我们可知，在凸函数上当函数在导数为0处取极值，同时极值就是函数的最值。故当 $L_{(w,b)}$ 关于 $w$ 和 $b$ 的导数均为零时，得到 $w$ 和 $b$ 的最优解。

我们考虑更为一般的情况，在数据集D上，样本由d个属性描述，此时线性回归又称为“多元线性回归”。为了便于讨论，我们把 $w 和 b统一记为\hat{\boldsymbol{w}}=(\boldsymbol{w};b)$ 与之对应的，数据集D表示为一个 $m\cdot (d+1)$ 大小的矩阵X,表示为

X = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ x T 1 x T 2 ⋮ x T m 11 ⋮ 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟, y = (y 1, y 2, . . ., y m)

$X=\begin{pmatrix} x_1^T &1 \\ x_2^T & 1\\ \vdots & \vdots\\ x_m^T & 1 \end{pmatrix},\boldsymbol{y}=(y_1,y_2,...,y_m)$
可以得到新的优化目标表达式为:

w^* = a r g min w (y - X w^) T (y - X w^)

$\hat{\boldsymbol{w}}^*=arg\min_{\boldsymbol{w}}(\boldsymbol{y-X\hat{w}})^T(\boldsymbol{y-X\hat{w}})$
新的代价函数为

L w^= (y - X w^) T (y - X w^)

$L_{\hat{w}}=(\boldsymbol{y-X\hat{w}})^T(\boldsymbol{y-X\hat{w}})$

我们可以对 $L_{\hat{w}} 关于\hat{w}$ 求导，得到

\partial L w ^ \partial w ^= 2 X T (X w^- y)

$\frac{\partial L_{\hat{w}}}{\partial \hat{w}} = 2\boldsymbol{X^T}(\boldsymbol{X\hat{w}-y})$

令上式为零，可以得到 $\hat{w}$ 的最优闭式解，此时需要分类讨论：

如果 $\boldsymbol{X^TX}$ 为满秩矩阵或正定矩阵时，可得：

$w^* = (X T X) - 1 X T y$ $\hat{\boldsymbol{w}}^*=(\boldsymbol{X^TX})^{-1}\boldsymbol{X^Ty}$
线性回归模型为:
$f (x i^) = x i^T w^* = x i^T (X T X) - 1 X T y$ $f(\hat{x_i})=\hat{x_i}^T\hat{\boldsymbol{w}}^*=\hat{x_i}^T(\boldsymbol{X^TX})^{-1}\boldsymbol{X^Ty}$
如果 $\boldsymbol{X^TX}$ 不是满秩矩阵时(样本数量小于属性数量)， $\boldsymbol{X}$ 的行数大于列数，必然不满秩，存在多个解析解，它们都能使得损失函数最小化，选择哪一个解析解作为输出，由学习算法的归纳偏好决定，这时候常常引入正则化(regularization)项.常见的正则化项如 $L_1正则化或L_2正则化，我们以L_2正则化为例$ :

$w^* = a r g min w [(y - X w^) T (y - X w^) + λ | | w^22 | |]; (λ > 0 为正则比例)$ $\hat{\boldsymbol{w}}^*=arg\min_{\boldsymbol{w}}[(\boldsymbol{y-X\hat{w}})^T(\boldsymbol{y-X\hat{w}})+\lambda||\boldsymbol{\hat{w}}_2^2||];( \lambda>0为正则比例 )$
则线性回归模型为:
$f (x i^) = x i^T w^* = a r g min w [(y - X w^) T (y - X w^) + λ | | w^22 | |]$ $f(\hat{x_i})=\hat{x_i}^T\hat{\boldsymbol{w}}^*=arg\min_{\boldsymbol{w}}[(\boldsymbol{y-X\hat{w}})^T(\boldsymbol{y-X\hat{w}})+\lambda||\boldsymbol{\hat{w}}_2^2||]$

5.线性回归小结

这里写图片描述

广义线性模型

概述

线性模型虽然简单，却有丰富的变化.当我们希望线性模型的预测值逼近真实标记值 $y$ 时，就得到了线性回归模型.为了便于观察，我们简写线性回归模型为 $y=w^Tx+b$ ，此时我们可以让线性回归模型去逼近 $y$ 的其他目标，例如 $lny$ .
照着这个想法，我们可以更一般地，考虑单调可导函数 $h(·),令h(y)=w^Tx+b$ ，这样 $w^Tx+b$ 可以逼近一个 $h(·)$ ，如果我们选好的 $h(·)$ 就可以使用线性模型去逼近不同的函数了，我们称这样的模型为广义线性模型(generalized linear model). 上述表达式可以记为:

y = h - 1 (w T x + b)

$y=h^{-1}(w^Tx+b)$
有一个典型的例子，当

h(⋅)=ln(⋅) $h(·)=ln(·)$ 时广义线性模型就是对数线性回归，即

y = e (w T x + b)

$y=e^{(w^Tx+b)}$

可以看到，虽然是叫广义线性模型，其实大部分时候变换函数 $h^{-1}$ 都是非线性变换，即实质上已经是非线性的。

逻辑回归

1.怎么用线性模型做分类？

在一开始我们说到，回归分析本质上是一个函数估计问题，就是找出因变量和自变量之间的因果关系。如果因变量是连续值那就叫回归，如果因变量是离散值那就叫分类。前面说的都是使用线性模型进行回归学习，下面我们就说说怎么用线性模型做分类。

其实用线性模型做分类的办法已经在上一节介绍，就是广义线性模型,考虑到在广义线性模型中，如果我们找到一个单调可微函数 $h(·)$ 可以将分类任务和线性回归联系到一起，那就可以用线性模型完成分类任务了。

2.什么样的函数满足分类要求？

这里，我们先简化分类任务，考虑二分类问题，即预测值为 $y=0或y=1$ ,考虑到 $z=w^Tx+b的取值z为(-∞,+∞)$ ,如果我们有这样一个理想函数(信号分析里面的阶跃函数):

这里写图片描述

当 $z<0 取值为0时即取负例，z>0时取值为1时即取正例$ ，可以看到，这样的函数能够很好的解决分类问题。这样的函数能作为广义线性模型里面的 $h(·)$ 吗？答案是不能。这是因为这样的函数有一个问题：函数不连续故不可导。那我们能不能找到一个类似这样的函数，并且连续可导?答案当然是有的:例如对数概率函数(logistic function)，下图为对数概率函数的示意图:

这里写图片描述

由图可以看出对数概率函数将自变量 $z$ 转化为一个接近0或者1的 $y$ 值，且输出值在 $z=0$ 处变化剧烈，重要的是对数概率函数满足连续可微的条件，故我们可将对数概率函数和广义线性模型结合到一起:

y = 1 1 + e - z, z = w T x + b

$y=\frac{1}{1+e^{-z}},z=w^Tx+b$

看到这里，或许你和我有一样的问题，为什么我们非要选用这样的一个函数，难道就是因为它连续可微和阶跃函数很像就足够了吗? 那么下面这些函数也有类似的性质：
这里写图片描述

为啥不选它们，这个其实是考虑到对数概率函数的导函数，我们看一下它的导函数形式:

d y d z = e - z ( 1 + e - z ) 2 = y e - z 1 + e - z = y (1 - y)

$\frac{\mathrm{d} y}{\mathrm{d} z}=\frac{e^{-z}}{(1+e^{-z})^2}=y\frac{e^{-z}}{1+e^{-z}}=y(1-y)$
对数概率函数的导函数可以直接通过原函数很方便的计算出来，我们做参数估计时是需要作求导操作的，函数有这样的特性是非常方便的。也就是它在众多函数中脱颖而出的原因之一

我们将上面的线性模型称为逻辑回归，为啥叫逻辑回归。这是因为有

l n y 1 - y = w T + b

$ln\frac{y}{1-y}=w^T+b$
如果我们将 $y$ 的输出值认为是 $P(y=1|x)$ ,那么 $1-y就是P(y=0|x)$ ，而 $\frac{y}{1-y}=\frac{P(y=1|x)}{P(y=0|x)}$ 称为概率(又称几率,odds)，概率反映了 $x$ 作为正例的相对可能性，对概率取对数则到了“对数概率”，所以这个函数又叫对数概率函数(logistic function),而对应的模型称为“逻辑回归”(logistic regression)。虽然叫做逻辑回归，实际上却是一种分类学习方法，这种方法有许多优点，例如对数概率函数的任意阶可导的凸函数，这有这很多很好的数学性质，有许多数值优化算法都可以直接应用求解。

3.怎么计算逻辑回归模型的参数估计？

确定好逻辑回归的模型后，下面该计算模型中的参数 $w和b$ ,这里我们依旧是把 $w 和 b统一记为\hat{\boldsymbol{w}}=(\boldsymbol{w};b)$ ,此时有

P (y = 1 | x) = e w ^ T x 1 + e w ^ T x

$P(y=1|\boldsymbol{x})=\frac{e^{\hat{\boldsymbol{w}}^Tx}}{1+e^{\hat{\boldsymbol{w}}^Tx}}$

P (y = 0 | x) = 1 1 + e w ^ T x

$P(y=0|\boldsymbol{x})=\frac{1}{1+e^{\hat{\boldsymbol{w}}^Tx}}$

我们可以通过极大似然法来估计 $\hat{w}$ ，学过数理统计的同志们应该对极大似然估计法不陌生，我当年是背着公式勉强没挂科，没怎么理解这到底是干啥的，现在用到了，就好好的回顾一下极大似然法。

极大似然估计

事实上，我们认为概率模型的训练过程就是参数估计过程。对于参数估计，统计学界有两个学派提供不同的解决方法：

频率主义学派(Frequentist):认为参数虽然未知，但却是客观存在的固定值，因此，可以通过优化似然函数等准则来确定参数值
贝叶斯学派(Bayesian)：认为参数是未观察到的随机变量，其本身也可有分布，因此，可假设参数服从一个先验分布，然后基于观测到的数据来计算参数的后验分布.

极大似然估计(Maximum Likelihood Estimation,MLE)是源于频率主义学派的根据数据采样来估计概率分布参数的经典方法。

从感性上认识

我们尽量的感性的去理解极大似然估计，假设我们现在已经拿到了很多样本(即因变量)，极大似然估计就是去找未知的参数估计值，使得样本发生的概率最大(因为我已经有了很多样本，所有我认为能使得样本发生概率最大的参数才是合逻辑的)。这个时候是在求样本所有观测的联合概率最大化，而这样的计算时一个连乘，我们可以通过取对数，变换为连加。然后在通过对未知参数求导，令导数为零，就得到了最大似然估计值。

数学表示

讲完了感性的认识，下面我们看一下怎样用数学表示：
令 $D_c$ 表示训练集D中第c类样本组成的集合，假设这样的样本是独立同分布的，参数 $\theta_c对于数据集D_c的似然是$

P (D c | θ c) = \prod x \in D c P (x | θ c)

$P(D_c|\theta_c) = \prod _{x∈D_c}P(x|\theta_c)$

对 $\theta_c$ 进行极大似然估计，就是去寻找能最大化似然 $P(D_c|\theta_c)的参数值\hat{\theta}_c$ ,直观上看，极大似然估计是试图在 $\theta_c$ 所有可能的取值中，找到一个能使数据出现的“可能性”最大的值。而这个找最大值的我们常用的方法就是求导。
上式的连成不易计算，且容易造成下溢，通常使用对数似然将连乘转换为连加，这不改变数据极值点。

L （ θ c ） = l o g P (D c | θ c) = \prod x \in D c l o g P (x | θ c)

$L（\theta_c）=logP(D_c|\theta_c) = \prod _{x∈D_c}logP(x|\theta_c)$
此时参数

θc的极大似然估计θc^为 $\theta_c的极大似然估计\hat{\theta_c}为$ :

θ c^= a r g θ c m a x L (θ c)

$\hat{\theta_c}=arg_{\theta_c}maxL(\theta_c)$
一般到这步就是对参数求偏导取零求参数值了。

需要注意这种参数化的方法虽能使类条件概率估计变得相对简单，单估计结果的准确性严重依赖于所假设的概率分布形式是否符合潜在的真实数据分布，在应用时想要做出能较好地接近潜在真实分布的假设，需要在一定程度上利用关于应用任务本身的经验知识。

回顾了下极大似然估计，现在我们继续计算 $\hat{w}$ ，上面我们说了可以通过极大似然法来估计 $\hat{w}$ ，在给定数据集 $\{(x_i,y_i)\}_{i=1}^m$ ，逻辑回归模型的对数似然函数为:

L (w^) = l n (\prod i = 1 m [p (y = 1 | x)] y i [p (y = 0 | x)] 1 - y i) = \sum i = 1 m [y i p (y = 1 | x) + (1 - y i) p (y = 0 | x)]

$L(\hat{w}) = ln(\prod _{i=1}^m[p(y=1|x)]^{y_i}[p(y=0|x)]^{1-y_i})=\sum _{i=1}^m[y_ip(y=1|x)+(1-y_i)p(y=0|x)]$
因为有

P (y = 1 | x) = e w ^ T x 1 + e w ^ T x

$P(y=1|\boldsymbol{x})=\frac{e^{\hat{\boldsymbol{w}}^Tx}}{1+e^{\hat{\boldsymbol{w}}^Tx}}$

P (y = 0 | x) = 1 1 + e w ^ T x

$P(y=0|\boldsymbol{x})=\frac{1}{1+e^{\hat{\boldsymbol{w}}^Tx}}$

故带入可得

L (w^) = \sum i = 1 m (- y i w^T x i + l n (1 + e w^T x i))

$L(\hat{w})=\sum _{i=1}^m(-y_i\hat{w}^Tx_i+ln(1+e^{\hat{w}^Tx_i}))$

上式是关于 $\hat{w}$ 的高阶可导连续凸函数，可用梯度下降法或者拟牛顿法求解。

4.多分类回归模型该怎么办？

以上的讨论都是二分类的逻辑回归模型，我们可以顺势推广到多分类逻辑回归模型，这里考虑到一对多等可能，我们可以设因变量Y的离散值取值集合为{1,2,…,K},则多分类逻辑回归模型为:

P (Y = k | x) = e w k ^ T x 1 + \sum K - 1 k = 1 e w k ^ T x; k = 1, 2, . . ., K - 1

$P(Y=k|\boldsymbol{x})=\frac{e^{\hat{\boldsymbol{w_k}}^Tx}}{1+\sum _{k=1}^{K-1}e^{\hat{\boldsymbol{w_k}}^Tx}};k=1,2,...,K-1$

P (Y = K | x) = 1 1 + \sum K - 1 k = 1 e w ^ T x

$P(Y=K|\boldsymbol{x})=\frac{1}{1+\sum _{k=1}^{K-1}e^{\hat{\boldsymbol{w}}^Tx}}$
参数估计方法类似二分类逻辑回归模型的参数估计的方法

5.逻辑回归总结

这里写图片描述

线性判别分析

1.什么线性判别分析?

线性判别分析(Linear Discriminant Analysis,LDA)是一种经典的线性学习方法，思想是:

训练时:设法将训练样本投影到一条直线上，使得同类样本的投影点尽可能地接近、异类的样本投影点尽可能的远
预测时:将待预测样本投影到学到的直线上，根据它的投影点的位置判断类别

两类样本的线性判别分析示意图如下:
(图中直线方程 $y=w^Tx省略了参数b$ ，是因为我们计算的是样本在直线上的投影，而我们总可以平移直线过原点且保持投影不变，故可以省略b。)

这里写图片描述

这里你可能会想，为什么非要投影搞来搞去的这么麻烦，为啥不是在两类样本中找到一条直线分开两类样本就好，对啊，我也是这么想的，这样的方法就是支持向量机(support vector machine,SVM)的思想。这里我们先不介绍SVM，继续探讨LDA.

2.怎样找到线性判别分析的优化目标?

给定数据集 $D=\{(x_i,y_i) \}_{i=1}^m,y_i∈\{0,1\},令X_i、\mu_i、\sum_i$ 分别表示第 $i∈\{0,1\}$ 类样本的集合、均值向量、协方差矩阵.若将数据投影到直线上，则：

两类样本的中心在直线上的投影分别是 $w^T\mu_0和w^T\mu_1$
两类样本的协方差(多变量的方差)分别为 $w^T\sum_0w和w^T\sum_1w$

我的优化目标:

同类样本投影点尽可能的近，故 $w^T\sum_0w+w^T\sum_1w$ 尽可能的小
异类样本投影点尽可能的远，故是 $||w^T\mu_0-w^T\mu_1||_2^2$ 尽可能的大

故联合起来，我们的优化目标为:

J = | | w T μ 0 - w T μ 1 | | 2 2 w T \sum 0 w + w T \sum 1 w = w T ( μ 0 - μ 1 ) ( μ 0 - μ 1 ) T w w T ( \sum 0 + \sum 1 ) w

$J=\frac{{||w^T\mu_0-w^T\mu_1||_2^2}}{w^T\sum_0w+w^T\sum_1w}=\frac{w^T(\mu_0-\mu_1)(\mu_0-\mu_1)^Tw}{w^T(\sum_0+\sum_1)w}$

这里我们定义类内散度矩阵(within-class scatter matrix)：

S w = \sum 0 + \sum 1 = \sum x \in X 0 (x - μ 0) (x - μ 0) T + \sum x \in X 1 (x - μ 1) (x - μ 1) T

$S_w=\sum_0+\sum_1=\sum_{x∈X_0}(x-\mu_0)(x-\mu_0)^T + \sum_{x∈X_1}(x-\mu_1)(x-\mu_1)^T$

定义类间散度矩阵(between-class scatter matrix)：

S b = (μ 0 - μ 1) (μ 0 - μ 1) T

$S_b=(\mu_0-\mu_1)(\mu_0-\mu_1)^T$

带入优化目标，我们可得到LDA最终的优化目标:

J = w T S b w w T S w w

$J = \frac{w^TS_bw} {w^TS_ww}$

3.怎样计算LDA模型中参数?

由上面的优化目标确定 $w$ 的求解问题为:

w^= a r g max w^w T S b w w T S w w

$\hat{w}=arg\max_{\hat{w}} \frac{w^TS_bw} {w^TS_ww}$

注意观察上式，无论是分母还是分子，都是 $w^T\_w$ 的结构，这说明了分子分母都是关于 $w$ 的二次项，故对 $\hat{w}$ 的求解与 $w$ 的长度无关。只与其方向有关。既然与长度无关，我们就简化计算，令 $w^TS_ww=1$ ，则新的优化目标写为:

w^= a r g max w^w T S b w = a r g - min w^w T S b w

$\hat{w}=arg\max_{\hat{w}} {w^TS_bw} =arg-\min_{\hat{w}} {w^TS_bw}$

s . t . w T S w w = 1

$s.t.w^TS_ww=1$

可以看到优化目标与限制条件之间的关系，很自然的想到了拉格朗日乘子法(看过SVM的证明很容易想到，想不到拉倒)，由拉格朗日乘子法，得：

- w T S b w + λ w T S w w = 0

$-w^TS_bw+\lambda w^TS_ww=0$

化 简 为 ： S b w = λ S w w; λ 是 拉 格 朗 日 乘 子

$化简为：S_bw= \lambda S_ww ;\lambda 是拉格朗日乘子$

注意到

∵ 类 间 散 度 矩 阵 定 义 为 : S b = (μ 0 - μ 1) (μ 0 - μ 1) T, 这 是 向 量 (μ 0 - μ 1) 与 自 身 的 外 积

$∵类间散度矩阵定义为:S_b=(\mu_0-\mu_1)(\mu_0-\mu_1)^T,这是向量(\mu_0-\mu_1)与自身的外积$

∵ S b w = (μ 0 - μ 1) (μ 0 - μ 1) T w

$∵S_bw=(\mu_0-\mu_1)(\mu_0-\mu_1)^Tw$

令 λ w = (μ 0 - μ 1) T w

$令\lambda_w=(\mu_0-\mu_1)^Tw$

那 么 S b w = λ w (μ 0 - μ 1) = λ S w w ， 因 为 与 w 的 长 度 无 关 ， 我 们 可 令 λ w = λ ， 则 有

$那么S_bw=\lambda_w(\mu_0-\mu_1)=\lambda S_ww，因为与w的长度无关，我们可令\lambda_w=\lambda，则有$

S w w = (μ 0 - μ 1) \Rightarrow w = S - 1 w (μ 0 - μ 1)

$S_ww=(\mu_0-\mu_1) \Rightarrow w=S_w^{-1}(\mu_0-\mu_1)$

上述讨论的就是 $\hat{w}$ 的求解,考虑到数值解的稳定性，在实践中通常是对 $S_w$ 进行奇异值分解(数据量大时候).

4.怎样在多分类任务中使用LDA？

我们可以将LDA推广到多分类任务中，假定存在M个类，属于第i个类的样本集合为 $T_i$ ，则样本的样例集合为 $i=1,2，...,M$ ; $T_i$ 中的样例个数为 $m_i$ ，则有 $\sum_{i=1}^Mm_i=N，N$ 为样本总数，我们可以得到样例的均值向量为：

μ i = (μ (1) i, μ (2) i, . . ., μ (n) i) T = 1 m i \sum x i \in T i x i = 1 N \sum x = 1 N x i

$\mu_i =(\mu_i^{(1)},\mu_i^{(2)},...,\mu_i^{(n)})^T=\frac{1} {m_i}\sum_{x_i∈T_i}x_i= \frac{1} {N}\sum_{x=1}^Nx_i$

这些样例的特征之间协方差矩阵为 $\sum_i$ .

要使得同类样例的投影点尽可能的近，主要到所有类内散度矩阵 $S_w$ 定义为每个类别的散度矩阵之和,故

$S w = \sum i = 1 M S w i; 其中 S w i = \sum x \in T i (x i - μ i) (x i - μ i) T$ $S_w=\sum_{i=1}^MS_{wi} ;其中S_{wi}=\sum_{x∈T_i}(x_i-\mu_i)(x_i-\mu_i)^T$
要使得异类间的投影距离尽可能的远,则可以使异类样例的中心点投影距离尽可能的远，这里有多个中心点。这里用每一类样例的中心点和总样例中心点的距离作为度量，考虑到每一类样例的样本集大小可能不同，故我们队这个距离作加权处理，因此定义类间散度矩阵为:

$S b = \sum i = 1 M m i (μ i - μ) (μ i - μ) T$ $S_b = \sum_{i=1}^Mm_i(\mu_i-\mu)(\mu_i-\mu)^T$
$(\mu_i-\mu)(\mu_i-\mu)^T$ 也是一个协方差矩阵，刻画的是第i类与总体之间的关系。

我们的优化目标为:

J = max W t r ( W T S b W ) t r ( W T S w W )

$J = \max_W\frac{tr(W^TS_bW)} {tr(W^TS_wW)}$

$W$ 的闭式解是 $S_w^{-1}S_b的d'个最大非零广义特征值所对应的特征向量组成的矩阵，d'≤M-1$

若将 $W$ 视为一个投影矩阵，则多分类LDA将样本投影到一个 $d'$ 维空间，而 $d'常小于数据原有的属性数M$ ,于是，可以通过这个投影来减少样本点的维数，且投影过程中使用了训练样本的类别信息，故这是一个经典的监督降维技术.

线性判别分析小结

这里写图片描述

Python实战线性模型

说完了理论部分，接下来就是Python代码实现了

使用模块与数据集介绍

1.使用的模块

导入线性模型需要用到的模块

    import numpy as np   
    import matplotlib.pyplot as plt
    from sklearn import datasets, linear_model,cross_validation

上述用到的numpy和matplotlib都是老相识了，最近也在系统的学习numpy模块.
sklearn是机器学习的常用的工具包，sklearn包含了大部分常见的机器学习模型，且还有一些自带的数据集。
上述的开发包如果一个一个安装还是比较麻烦的，建议直接安装一个Python的发行版，我安装是Anaconda，配合的是Pycharm IDE，感觉还是很方便的。

下图是我的开发环境:
这里写图片描述

2.数据集

在线性回归问题中，使用的数据集是scikit-learn自带的一个糖尿病人的数据集，该数据集特点如下：

数据集有442个样本
每个样本有10个特征
每个特征都是浮点数，数据都在-0.2~0.2之间
样本的目标在整数25~346之间

给出加载数据集的函数：

def load_data():
    '''
    加载用于回归问题的数据集

    :return: 一个元组，用于回归问题。
    元组元素依次为：训练样本集、测试样本集、训练样本集对应的值、测试样本集对应的值
    '''
    diabetes = datasets.load_diabetes()#使用 scikit-learn 自带的一个糖尿病病人的数据集
    return cross_validation.train_test_split(diabetes.data,diabetes.target,
        test_size=0.25,random_state=0) # 拆分成训练集和测试集，测试集大小为原始数据集大小的 1/4

数据归一化的作用

数据集上有一点值得我们注意:数据都在-0.2~0.2之间，这是原数据集帮我们做了特征归一化工作，在机器学习中，我们常使用梯度下降法求解最优化问题，在使用梯度下降法时，要注意特征归一化(Feature Scaling).特征归一化有两个明显的好处:

1.提升模型的收敛速度，例如模型只有两个特征 $x_1和x_2,x_1$ 的取值为0~2000，而 $x_2$ 的取值为1~5，如果对其优化，会得到一个窄长的椭圆形，导致在梯度下降时，梯度的方向在垂直等高线的方向而走之字形路线，这样迭代速度很降低。如果归一化，会使一个椭圆，梯度的方向为直接指向圆心，迭代速度就会很快。
2.提升模型的精度，这在涉及一些距离计算的算法时效果显著，比如算法要计算欧式距离，上面的 $x_2$ 的取值范围较小，计算时对结果的影响远比 $x_1$ 带来的小，可能会造成精度的损失.

在逻辑回归问题中，使用的数据集是鸢尾花数据集，该数据集特点如下：

数据集有150个
数据分为3类(分别为setora、versicolor、virginica)
每类有50个数据
每个数据包含4个属性:萼片(sepal)长度、萼片宽度、花瓣(petal)长度、花瓣宽度

给出加载数据集的函数：

    def load_data():
        '''
        加载用于分类问题的数据集

        :return: 一个元组，用于分类问题
        元组元素依次为：训练样本集、测试样本集、训练样本集对应的标记、测试样本集对应的标记
        '''
        iris=datasets.load_iris() # 使用 scikit-learn 自带的 iris 数据集
        X_train=iris.data
        y_train=iris.target
        return cross_validation.train_test_split(X_train, y_train,test_size=0.25,
            random_state=0,stratify=y_train)# stratify=y_train指定分层采样，拆分成训练集和测试集，测试集大小为原始数据集大小的 1/4

这里有一点要注意的是，原数据中，前50个数据都是类别0，中间50是类别1，后50数据是类别2。如果我们按顺序切分采样，那么切分的测试数据就不是无偏的了(都是类型2了),故我们指定了stratify=y_train，采取分层采样.

线性回归模型

1.线性回归

linear_model是sklearn包下提供的一个有关于线性模型的工具包。本次使用的是linear_model下提供了LinearRegression线性回归模型

LinearRegression类介绍

LinearRegression类原型为:

 class LinearRegression(LinearModel, RegressorMixin): 
            '''Ordinary least squares Linear Regression.'''

父类

LinearModel:
            '''Linear Models的基类'''
        方法：
            predict(self, X):
                Predict class labels for samples in X

        RegressorMixin：
            '''Mixin class for all regression estimators in scikit-learn'''
        方法：
            score(self, X, y, sample_weight=None)
                '''返回预测性能得分。'''

这里要讲一下score函数的原理，设预测集 $T_{test}$ ，真实值为 $y_i$ ，真实值的均值为 $\bar{y}$ ，预测值为 $\hat{y_i}$ ，则:

s o c r e = 1 - \sum T t e s t ( y i - y i ^ ) 2 ( y i - y ¯ ) 2

$socre = 1-\frac{\sum_{T_{test}}(y_i-\hat{y_i})^2} {(y_i-\bar{y})^2}$

score不超过1，可能为负值(预测效果太差)
score越大，预测性能越好

参数

def __init__(self, fit_intercept=True, normalize=False, copy_X=True,
                 n_jobs=1):

        ----------
        fit_intercept : boolean, optional
            是否计算b值,为false不计算

        normalize : boolean, optional, default False
            为True，训练样本在使用前会被归一化，
            This parameter is ignored when `fit_intercept` is set to False.
            When the regressors are normalized, note that this makes the
            hyperparameters learnt more robust and almost independent of the number
            of samples. 

        copy_X : boolean, optional, default True
            If True, 会复制X

        n_jobs : int, optional, default 1
            使用任务并行时指定的CPU数量
            If -1 all CPUs are used.

属性

        ----------
        coef_ : array, shape (n_features, ) or (n_targets, n_features)
            权重系数

        intercept_ : array
           b值

方法

        fit(self, X, y, sample_weight=None) : 训练模型

        predict(self, X):  ：用模型进行预测 (从父类继承)

        score(self, X, y, sample_weight=None)  :返回预测性能得分。(从父类继承)

使用LinearRegression的示例函数如下

def test_LinearRegression(*data):
        '''
        测试 LinearRegression 的用法

        param:*data: 可变参数。它是一个元组
        要求其元素依次为：训练样本集、测试样本集、训练样本的值、测试样本的值
        :return: None
        '''
        X_train,X_test,y_train,y_test=data
        regr = linear_model.LinearRegression()
        regr.fit(X_train, y_train)
        print('Coefficients:%s, intercept %.2f'%(regr.coef_,regr.intercept_))
        # 求方差：  先求残差平方和  在使用np.mean为求均值 
        print("Residual sum of squares: %.2f"% np.mean((regr.predict(X_test) - y_test) ** 2))
        # 预测分数
        print('Score: %.2f' % regr.score(X_test, y_test))

本节程序小结：

    # -*- coding: utf-8 -*-
    """
        广义线性模型
        ~~~~~~~~~~~~~~~~~~~~~~~~~~
        LinearRegression
    """
    import matplotlib.pyplot as plt
    import numpy as np
    from sklearn import datasets, linear_model,cross_validation

    def load_data():
        '''
        加载用于回归问题的数据集
        '''
        diabetes = datasets.load_diabetes()#使用 scikit-learn 自带的一个糖尿病病人的数据集
        return cross_validation.train_test_split(diabetes.data,diabetes.target,
            test_size=0.25,random_state=0) # 拆分成训练集和测试集，测试集大小为原始数据集大小的 1/4


    def test_LinearRegression(*data):
        '''
        测试 LinearRegression 的用法
        :return: None
        '''
        X_train,X_test,y_train,y_test=data
        regr = linear_model.LinearRegression()
        regr.fit(X_train, y_train)
        print('Coefficients:%s, intercept %.2f'%(regr.coef_,regr.intercept_))
        print("Residual sum of squares: %.2f"% np.mean((regr.predict(X_test) - y_test) ** 2))
        print('Score: %.2f' % regr.score(X_test, y_test))

    if __name__=='__main__':
        X_train,X_test,y_train,y_test=load_data() # 产生用于回归问题的数据集
        test_LinearRegression(X_train,X_test,y_train,y_test) # 调用 test_LinearRegression

输出:

   Coefficients:[ -43.26774487 -208.67053951  593.39797213  302.89814903 -560.27689824
      261.47657106   -8.83343952  135.93715156  703.22658427   28.34844354], intercept 153.07
    Residual sum of squares: 3180.20
    Score: 0.36

均方误差为3180.20,预测性能得分为0.36(效果不咋的)

带正则化的线性回归模型

1.概述

在前面的理论部分，我们讲到了对于线性回归问题，如果 $X^TX$ 不是满秩矩阵时存在多个解析解，都能使得均方误差最小，在求解的过程中，常引入正则项来求解。正则项说白了就是对模型的参数添加一些先验假设，控制模型空间，达到使得模型复杂度较小的目的.根据不同的正则化项，有不同的方法：

岭回归：Ridge Regression. 正则化项为 $\alpha ||w||_2^2,\alpha ≥0$ ，即加入L2范数惩罚项
Lasso回归： Lasso Regression. 正则化项为 $\alpha ||w||_1,\alpha ≥0$ ，即加入L1范数惩罚项
ElasticNet回归：Elastic Net.正则化项为 $\alpha\rho ||w||_1+\frac{\alpha(1-\rho)} {2}||w||_2^2,\alpha ≥0,1≥\rho≥0$ ，混合L1范数和L2范数惩罚项

正则化系数 $\alpha$ 的选择很关键，一般在选择 $\alpha$ 时，会先调节 $\alpha$ 为0，确定好一个learning rate，在调节 $\alpha$ 的值并观察validation accuracy,经过粗调到微调的过程确定好的 $\alpha$ .(下面使用的都是集成好的函数，没有涉及到learning rate的调整)

2.岭回归

Ridge Regression通过在代价函数中加入L2范数惩罚项，从而控制线性模型的复杂程度，使得模型更为robust.inear_model下提供了Ridge类实现岭回归模型。

Ridge类介绍

Ridge类原型为:

 class Ridge(_BaseRidge, RegressorMixin):  
            '''Linear least squares with l2 regularizationn.'''

父类

    _BaseRidge:
            ''' Ridge Models的基类'''            
        方法：
            fit(self, X, y, sample_weight=None):
                训练模型

    RegressorMixin：
            介绍见前面

参数

def __init__(self, alpha=1.0, fit_intercept=True, normalize=False,
                 copy_X=True, max_iter=None, tol=1e-3, solver="auto",
                 random_state=None):
        ----------
        alpha : {float, array-like}, shape (n_targets)
            a值，越大代表正则化占比越大

        copy_X : boolean, optional, default True
            If True, 会复制X

        fit_intercept : boolean
            为False就不计算b
            (e.g. 模型会假设你的数据已经中心化了).

        max_iter : int, optional
            指定最大迭代次数，If None，则为默认值(不同的solver默认值不同)

        normalize : boolean, optional, default False
            If True, 数据在使用前会归一化

        solver : {'auto', 'svd', 'cholesky', 'lsqr', 'sparse_cg', 'sag'}
            选择不同的求解最优化问题的算法.
            - 'auto' 依据数据集自动选择.
            - 'svd' 使用奇异值分解来计算回归系数.比'cholesky'稳定.
            - 'cholesky' 使用scipy.linalg.solve求解
            - 'sparse_cg' 使用scipy.sparse.linalg.cg求解.对于大型数据，比'cholesky'合适
            - 'lsqr' uses 使用scipy.sparse.linalg.lsqr求解，速度快但是可能不支持低版本(0.17以下)
            - 'sag' 使用 Stochastic Average Gradient descent求解。当n_samples and n_feature比较大时候，这个方法速度更快.

        tol : float
            指定判别迭代收敛的阈值.

        random_state : int seed, RandomState实例, or None (default)，在slover=sag时使用
            - 如果为整数，指定了随机数生成器的种子
            - 如果为RandomState实例，指定了随机数生成器
            - 如果为None，使用默认的随机数生成器

属性

        ----------
        coef_ : array, shape (n_features, ) or (n_targets, n_features)
            权重系数

        intercept_ : array
           b值   

        n_iter_ : array or None, shape (n_targets,)
           实际的迭代次数

方法

        fit(self, X, y, sample_weight=None) : 训练模型

        predict(self, X):  ：用模型进行预测 (从父类继承)

        score(self, X, y, sample_weight=None)  :返回预测性能得分。(从父类继承)

使用Ridge的函数如下(观察函数预测性能随着正则系数变化趋势)

def test_Ridge_alpha(*data):
    '''
    测试 Ridge 的预测性能随 alpha 参数的影响
    :param data: 可变参数。它是一个元组
    这里要求其元素依次为：训练样本集、测试样本集、训练样本的值、测试样本的值
    :return: None
    '''
    X_train,X_test,y_train,y_test=data
    alphas=[0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50,100,200,500,1000]  #使用不同的正则系数
    scores=[]
    for i,alpha in enumerate(alphas):
        regr = linear_model.Ridge(alpha=alpha)
        regr.fit(X_train, y_train)
        scores.append(regr.score(X_test, y_test))
    ## 绘图  
    fig=plt.figure()
    ax=fig.add_subplot(1,1,1)
    ax.plot(alphas,scores)
    ax.set_xlabel(r"$\alpha$")
    ax.set_ylabel(r"score")
    ax.set_xscale('log')
    ax.set_title("Ridge")
    plt.show()

本节程序小结：

    # -*- coding: utf-8 -*-
    """
        广义线性模型
        ~~~~~~~~~~~~~~~~~~~~~~~~~~
        Ridge回归
    """
    import matplotlib.pyplot as plt
    import numpy as np
    from sklearn import datasets, linear_model,cross_validation

    def load_data():
        '''
        加载用于回归问题的数据集
        '''
        diabetes = datasets.load_diabetes()
        return cross_validation.train_test_split(diabetes.data,diabetes.target,
            test_size=0.25,random_state=0) 

    def test_Ridge_alpha(*data):
        '''
        测试 Ridge 的预测性能随 alpha 参数的影响

        :param data: 可变参数。它是一个元组
        这里要求其元素依次为：训练样本集、测试样本集、训练样本的值、测试样本的值
        :return: None
        '''
        X_train,X_test,y_train,y_test=data
        alphas=[0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50,100,200,500,1000]
        scores=[]
        for i,alpha in enumerate(alphas):
            regr = linear_model.Ridge(alpha=alpha)
            regr.fit(X_train, y_train)
            scores.append(regr.score(X_test, y_test))
        ## 绘图
        fig=plt.figure()
        ax=fig.add_subplot(1,1,1)
        ax.plot(alphas,scores)
        ax.set_xlabel(r"$\alpha$")
        ax.set_ylabel(r"score")
        ax.set_xscale('log')
        ax.set_title("Ridge")
        plt.show()



    if __name__=='__main__':
        X_train,X_test,y_train,y_test=load_data() # 产生用于回归问题的数据集
        test_Ridge_alpha(X_train,X_test,y_train,y_test) # 调用 test_Ridge_alpha

输出:

这里写图片描述

可以看到，当 $\alpha$ 超过1时，随着 $\alpha$ 增长，模型预测性能急剧下降。这是因为 $\alpha$ 较大，正则化项 $\alpha ||w||_2^2$ 影响变大，模型趋于简单。极端的来讲当 $\alpha \rightarrow ∞$ ，而正则化项 $\alpha ||w||_2^2$ 对损失函数的影响非常大，我们的目标是损失函数尽量小，故此时 $||w||_2^2$ 会尽量的小，也就是 $w$ 会趋于0，这时候模型基本就剩下一个b了(输入啥都预测b)。这样就会导致模型变得简单，从而预测性能大大降低。

3.Lasso回归

Lasso Regression和岭回归的区别在于引入的是入L1范数惩罚项，可以将系数控制收缩到0，从而达到变量选择的效果(参数稀疏的功能)。linear_model下提供了Lasso类实现Lasso回归模型

Lasso类介绍

Lasso类原型为:

 class Lasso(ElasticNet):  
            '''Linear Model trained with L1 prior as regularizer (aka the Lasso).'''

父类

    _ElasticNet:
            ''' 后面会讲解'''

参数

def __init__(self, alpha=1.0, fit_intercept=True, normalize=False,
             precompute=False, copy_X=True, max_iter=1000,
             tol=1e-4, warm_start=False, positive=False,
             random_state=None, selection='cyclic'):

    ----------
    precompute : True | False | array-like, default=False
        是否提前计算Gram矩阵来加速计算.

    warm_start : bool, optional
        When set to True, 使用前一次训练的结果训练

    positive : bool, optional
        When set to ``True``, 强制要求权值向量都是正数

    selection : {'random','cyclic'} default 'cyclic'  每次迭代的时候，选择权重向量的哪个分量来更新
        - 'random':随机选择权重向量的一个分量来更新
        - 'cyclic':从前向后依次选择权重向量的一个分量来更新

    alpha : {float, array-like}, shape (n_targets)
        a值，越大代表正则化占比越大

    copy_X : boolean, optional, default True
        If True, 会复制X

    fit_intercept : boolean
        为False就不计算b
        (e.g. 模型会假设你的数据已经中心化了).

    max_iter : int, optional
        指定最大迭代次数，If None，则为默认值(不同的solver默认值不同)

    normalize : boolean, optional, default False
        If True, 数据在使用前会归一化

    tol : float
        指定判别迭代收敛的阈值.

    random_state : int seed, RandomState实例, or None (default)，在slover=sag时使用
        - 如果为整数，指定了随机数生成器的种子
        - 如果为RandomState实例，指定了随机数生成器
        - 如果为None，使用默认的随机数生成器

属性

        ----------
        coef_ : array, shape (n_features, ) or (n_targets, n_features)
            权重系数

        intercept_ : array
           b值   

        n_iter_ : array or None, shape (n_targets,)
           实际的迭代次数

方法

        fit(self, X, y, sample_weight=None) : 训练模型

        predict(self, X):  ：用模型进行预测 (从父类继承)

        score(self, X, y, sample_weight=None)  :返回预测性能得分。(从父类继承)

使用Lasso的函数如下(观察函数预测性能随着正则系数变化趋势)

def test_Lasso_alpha(*data):
    '''
    测试 Lasso 的预测性能随 alpha 参数的影响

    :param data: 可变参数。它是一个元组
    要求其元素依次为：训练样本集、测试样本集、训练样本的值、测试样本的值
    :return: None
    '''
    X_train,X_test,y_train,y_test=data
    alphas=[0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50,100,200,500,1000]
    scores=[]
    for i,alpha in enumerate(alphas):
        regr = linear_model.Lasso(alpha=alpha)
        regr.fit(X_train, y_train)
        scores.append(regr.score(X_test, y_test))
    ## 绘图
    fig=plt.figure()
    ax=fig.add_subplot(1,1,1)
    ax.plot(alphas,scores)
    ax.set_xlabel(r"$\alpha$")
    ax.set_ylabel(r"score")
    ax.set_xscale('log')
    ax.set_title("Lasso")
    plt.show()

本节程序小结：

    # -*- coding: utf-8 -*-
    """
    广义线性模型
    ~~~~~~~~~~~~~~~~~~~~~~~~~~
    Lasso
    """
    import matplotlib.pyplot as plt
    import numpy as np
    from sklearn import datasets, linear_model,cross_validation
    def load_data():
        '''
        加载用于回归问题的数据集
        '''
        diabetes = datasets.load_diabetes()
        return cross_validation.train_test_split(diabetes.data,diabetes.target,
            test_size=0.25,random_state=0) 

    def test_Lasso_alpha(*data):
        '''
        测试 Lasso 的预测性能随 alpha 参数的影响

        :param data: 可变参数。它是一个元组
        这里要求其元素依次为：训练样本集、测试样本集、训练样本的值、测试样本的值
        :return: None
        '''
        X_train,X_test,y_train,y_test=data
        alphas=[0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50,100,200,500,1000]
        scores=[]
        for i,alpha in enumerate(alphas):
            regr = linear_model.Lasso(alpha=alpha)
            regr.fit(X_train, y_train)
            scores.append(regr.score(X_test, y_test))
        ## 绘图
        fig=plt.figure()
        ax=fig.add_subplot(1,1,1)
        ax.plot(alphas,scores)
        ax.set_xlabel(r"$\alpha$")
        ax.set_ylabel(r"score")
        ax.set_xscale('log')
        ax.set_title("Lasso")
        plt.show()

    if __name__=='__main__':
        X_train,X_test,y_train,y_test=load_data() # 产生用于回归问题的数据集
        test_Lasso_alpha(X_train,X_test,y_train,y_test) # 调用 test_Lasso_alpha

输出:

这里写图片描述

和Ridge差不多。

4.ElasticNet回归

ElasticNet是对岭回归和Lasso回归的综合，惩罚项是对L1和L2的均衡。linear_model下提供了ElasticNet类实现ElasticNet回归模型.

ElasticNet类介绍

ElasticNet类原型为:

 class ElasticNet(LinearModel, RegressorMixin):  
            '''Linear regression with combined L1 and L2 priors as regularizer.'''

父类

    LinearModel:
        """Base class for Linear Models"""

    RegressorMixin:
        """Mixin class for all regression estimators in scikit-learn."""

参数

def __init__(self, alpha=1.0, l1_ratio=0.5, fit_intercept=True,
             normalize=False, precompute=False, max_iter=1000,
             copy_X=True, tol=1e-4, warm_start=False, positive=False,
             random_state=None, selection='cyclic'):

    ----------
    l1_ratio : float
        p值.

    alpha : {float, array-like}, shape (n_targets)
        a值，越大代表正则化占比越大

    copy_X : boolean, optional, default True
        If True, 会复制X

    fit_intercept : boolean
        为False就不计算b
        (e.g. 模型会假设你的数据已经中心化了).

    max_iter : int, optional
        指定最大迭代次数，If None，则为默认值(不同的solver默认值不同)

    normalize : boolean, optional, default False
        If True, 数据在使用前会归一化

    precompute : True | False | array-like, default=False
        是否提前计算Gram矩阵来加速计算.

    warm_start : bool, optional
        When set to True, 使用前一次训练的结果训练

    positive : bool, optional
        When set to ``True``, 强制要求权值向量都是正数

    selection : {'random','cyclic'} default 'cyclic'  每次迭代的时候，选择权重向量的哪个分量来更新
        - 'random':随机选择权重向量的一个分量来更新
        - 'cyclic':从前向后依次选择权重向量的一个分量来更新

    tol : float
        指定判别迭代收敛的阈值.

    random_state : int seed, RandomState实例, or None (default)，在slover=sag时使用
        - 如果为整数，指定了随机数生成器的种子
        - 如果为RandomState实例，指定了随机数生成器
        - 如果为None，使用默认的随机数生成器

属性

        ----------
        coef_ : array, shape (n_features, ) or (n_targets, n_features)
            权重系数

        intercept_ : array
           b值   

        n_iter_ : array or None, shape (n_targets,)
           实际的迭代次数

方法

        fit(self, X, y, sample_weight=None) : 训练模型

        predict(self, X):  ：用模型进行预测 (从父类继承)

        score(self, X, y, sample_weight=None)  :返回预测性能得分。(从父类继承)

使用ElasticNet的函数如下(观察函数预测性能随着正则系数变化趋势)

def test_ElasticNet_alpha_rho(*data):
   '''
    测试 ElasticNet 的预测性能随 alpha 和 l1_ratio 的影响

    :param data: 可变参数。它是一个元组
    要求其元素依次为：训练样本集、测试样本集、训练样本的值、测试样本的值
    :return: None
    '''
    X_train,X_test,y_train,y_test=data
    alphas=np.logspace(-2,2)
    rhos=np.linspace(0.01,1)
    scores=[]
    for alpha in alphas:
            for rho in rhos:
                regr = linear_model.ElasticNet(alpha=alpha,l1_ratio=rho)
                regr.fit(X_train, y_train)
                scores.append(regr.score(X_test, y_test))
    ## 绘图
    alphas, rhos = np.meshgrid(alphas, rhos) #  np.meshgrid接收两个一维数组，并产生两个二维矩阵 alphas(50,50) rhos(50,50)
    scores=np.array(scores).reshape(alphas.shape) #转换为array并转换为二维矩阵 scores(50,50)
    from mpl_toolkits.mplot3d import Axes3D  #两个超参数，为三维图像
    from matplotlib import cm
    fig=plt.figure()
    ax=Axes3D(fig)  
    surf = ax.plot_surface(alphas, rhos, scores, rstride=1, cstride=1, cmap=cm.jet,
        linewidth=0, antialiased=False)
    fig.colorbar(surf, shrink=0.5, aspect=5)
    ax.set_xlabel(r"$\alpha$")
    ax.set_ylabel(r"$\rho$")
    ax.set_zlabel("score")
    ax.set_title("ElasticNet")
    plt.show()

本节程序小结：

     # -*- coding: utf-8 -*-
    """
        广义线性模型
        ~~~~~~~~~~~~~~~~~~~~~~~~~~
        ElasticNet

    """
    import matplotlib.pyplot as plt
    import numpy as np
    from sklearn import datasets, linear_model,cross_validation

    def load_data():
        '''
        加载用于回归问题的数据集
        '''
        diabetes = datasets.load_diabetes()
        return cross_validation.train_test_split(diabetes.data,diabetes.target,
            test_size=0.25,random_state=0) 

    def test_ElasticNet_alpha_rho(*data):
        '''
        测试 ElasticNet 的预测性能随 alpha 和 l1_ratio 的影响

        :param data: 可变参数。它是一个元组
        要求其元素依次为：训练样本集、测试样本集、训练样本的值、测试样本的值
        :return: None
        '''
        X_train,X_test,y_train,y_test=data
        alphas=np.logspace(-2,2) # 默认50个点 从10^-2到10^2取等分向量
        rhos=np.linspace(0.01,1) #默认100个点 从0.01到1取等分向量
        scores=[]
        for alpha in alphas:
                for rho in rhos:
                    regr = linear_model.ElasticNet(alpha=alpha,l1_ratio=rho)
                    regr.fit(X_train, y_train)
                    scores.append(regr.score(X_test, y_test))
        ## 绘图
        alphas, rhos = np.meshgrid(alphas, rhos) #  np.meshgrid接收两个一维数组，并产生两个二维矩阵 alphas(50,50) rhos(50,50)
        scores=np.array(scores).reshape(alphas.shape) #转换为array并转换为二维矩阵 scores(50,50)
        from mpl_toolkits.mplot3d import Axes3D  #两个超参数，为三维图像
        from matplotlib import cm
        fig=plt.figure()
        ax=Axes3D(fig)  
        surf = ax.plot_surface(alphas, rhos, scores, rstride=1, cstride=1, cmap=cm.jet,
            linewidth=0, antialiased=False)
        fig.colorbar(surf, shrink=0.5, aspect=5)
        ax.set_xlabel(r"$\alpha$")
        ax.set_ylabel(r"$\rho$")
        ax.set_zlabel("score")
        ax.set_title("ElasticNet")
        plt.show()

    if __name__=='__main__':
        X_train,X_test,y_train,y_test=load_data() # 产生用于回归问题的数据集
        test_ElasticNet_alpha_rho(X_train,X_test,y_train,y_test) # 调用 test_ElasticNet_alpha_rho

输出:

这里写图片描述

随着 $\alpha$ 的增大，预测性能下降。而 $\rho$ 主要影响性能下降的速度，因为这个参数控制L1和L2的比例.

逻辑回归模型

逻辑回归模型

LogisticRegression回归

linear_model下提供了LogisticRegression类实现逻辑回归模型。

LogisticRegression类介绍

LogisticRegression类原型为:

class LogisticRegression(BaseEstimator, LinearClassifierMixin, _LearntSelectorMixin, SparseCoefMixin): '''Logistic Regression (aka logit, MaxEnt) classifier.'''
父类

BaseEstimator: """Base class for all estimators in scikit-learn""" LinearClassifierMixin """Mixin for linear classifiers. _LearntSelectorMixin '''Transformer mixin selecting features based on importance weights.''' SparseCoefMixin: """Mixin for converting coef_ to and from CSR format."""
参数

def __init__(self, penalty='l2', dual=False, tol=1e-4, C=1.0, fit_intercept=True, intercept_scaling=1, class_weight=None, random_state=None, solver='liblinear', max_iter=100, multi_class='ovr', verbose=0, warm_start=False, n_jobs=1) ---------- penalty : str, 'l1' or 'l2', default: 'l2' 指定了正则化策略. The 'newton-cg','sag' and 'lbfgs' solvers support only l2 penalties.

penalty=’l2’，则优化目标函数为: $0.5||w||_2^2+CL(w),C>0；L(w)$ 为极大似然函数
penalty=’l1’，则优化目标函数为: $||w||_1+CL(w),C>0；L(w)$ 为极大似然函数

dual : bool, default: False 如果为True，则求解对偶形式(penalty='l2'且solver=’liblinear‘有对偶形式) False求解原形式 C : float, default: 1.0 指定正则化惩罚系数的倒数，值越小，正则化项越大 fit_intercept : boolean 为False就不计算b (e.g. 模型会假设你的数据已经中心化了). intercept_scaling : float, default 1. 只有当solver='liblinear'才有意义，采用intercept_scaling时，相当于人造一个特征出来，特征恒为1，权重为b.计算正则化项的时候，这个项也被考虑进去了，为了降低这个人造特征的影响，需要提供intercept_scaling class_weight : dict or 'balanced', default: None - dict: 字典给出了每个分类的权重,如{class_label:weight} - 'balanced'：每个分类的权重与该分类在样本集中出现频率成反比 - 未指定:每个分类的权重都为1 max_iter : int, default: 100 Useful only for the newton-cg, sag and lbfgs solvers. 最大迭代次数 random_state : int seed, RandomState instance, default: None Used only in solvers 'sag' and 'liblinear'. - 如果为整数，指定了随机数生成器的种子 - 如果为RandomState实例，指定了随机数生成器 - 如果为None，使用默认的随机数生成器 solver : {'newton-cg', 'lbfgs', 'liblinear', 'sag'}, default: 'liblinear' 指定求解最优化问题的算法. - 对于小数据集, 'liblinear' is a good choice - 对于大数据集,whereas 'sag' is faster for large ones. - 对于多分类问题, only 'newton-cg', 'sag' and 'lbfgs' handle multinomial loss; 'liblinear' is limited to one-versus-rest schemes. - 'newton-cg', 'lbfgs' and 'sag' only 针对 L2 正则化. tol : float 指定判别迭代收敛的阈值. multi_class : str, {'ovr', 'multinomial'}, default: 'ovr' 指定对于多分类问题的策略 - 'ovr', 采用one-vs-rest策略 - 'multinomial':直接采用多分类了逻辑回归策略 verbose : int, default: 0 用于开启/关闭迭代中间输出日志功能 warm_start : bool, optional When set to True, 使用前一次训练的结果训练 n_jobs : int, default: 1 指定任务并行时的CPU数量如果为-1, all cores are used.
属性

---------- coef_ : array, shape (n_features, ) or (n_targets, n_features) 权重系数 intercept_ : array b值 n_iter_ : array or None, shape (n_targets,) 实际的迭代次数
方法

fit(self, X, y, sample_weight=None) : 训练模型 predict(self, X): ：用模型进行预测 (从父类继承) predict_log_proba(self, X): 返回一个数组，数组的元素依次是X预测为各个类别的概率的对数值 predict_proba(self, X):返回一个数组，数组元素依次是X预测为各个类别的概率值 score(self, X, y, sample_weight=None) :返回预测性能得分。(从父类继承)
使用LogisticRegression的函数如下(观察多分类下影响、观察函数预测性能随着正则系数变化趋势)

def test_LogisticRegression(*data): ''' 测试 LogisticRegression 的用法简单地从训练集中学习，参数都使用默认的 :param data: 可变参数。它是一个元组这里要求其元素依次为：训练样本集、测试样本集、训练样本的标记、测试样本的标记 :return: None ''' X_train,X_test,y_train,y_test=data regr = linear_model.LogisticRegression() regr.fit(X_train, y_train) print('Coefficients:%s, intercept %s'%(regr.coef_,regr.intercept_)) print('Score: %.2f' % regr.score(X_test, y_test)) def test_LogisticRegression_multinomial(*data): ''' 测试 LogisticRegression 的预测性能随 multi_class 参数的影响 :param data: 可变参数。它是一个元组要求其元素依次为：训练样本集、测试样本集、训练样本的标记、测试样本的标记 :return: None ''' X_train,X_test,y_train,y_test=data regr = linear_model.LogisticRegression(multi_class='multinomial',solver='lbfgs') regr.fit(X_train, y_train) print('Coefficients:%s, intercept %s'%(regr.coef_,regr.intercept_)) print('Score: %.2f' % regr.score(X_test, y_test)) def test_LogisticRegression_C(*data): ''' 测试 LogisticRegression 的预测性能随 C 参数的影响 :param data: 可变参数。它是一个元组要求其元素依次为：训练样本集、测试样本集、训练样本的标记、测试样本的标记 :return: None ''' X_train,X_test,y_train,y_test=data Cs=np.logspace(-2,4,num=100) scores=[] for C in Cs: regr = linear_model.LogisticRegression(C=C) regr.fit(X_train, y_train) scores.append(regr.score(X_test, y_test)) ## 绘图 fig=plt.figure() ax=fig.add_subplot(1,1,1) ax.plot(Cs,scores) ax.set_xlabel(r"C") ax.set_ylabel(r"score") ax.set_xscale('log') ax.set_title("LogisticRegression") plt.show()

本节程序小结：

# -*- coding: utf-8 -*- """ 广义线性模型 ~~~~~~~~~~~~~~~~~~~~~~~~~~ Logistic 回归 """ import matplotlib.pyplot as plt import numpy as np from sklearn import datasets, linear_model,cross_validation def load_data(): ''' 加载用于分类问题的数据集 ''' iris=datasets.load_iris() # 使用 scikit-learn 自带的 iris 数据集 X_train=iris.data y_train=iris.target return cross_validation.train_test_split(X_train, y_train,test_size=0.25, random_state=0,stratify=y_train)# 分层采样拆分成训练集和测试集，测试集大小为原始数据集大小的 1/4 def test_LogisticRegression(*data): ''' 测试 LogisticRegression 的用法 ''' X_train,X_test,y_train,y_test=data regr = linear_model.LogisticRegression() regr.fit(X_train, y_train) print('Coefficients:%s, intercept %s'%(regr.coef_,regr.intercept_)) print('Score: %.2f' % regr.score(X_test, y_test)) def test_LogisticRegression_multinomial(*data): ''' 测试 LogisticRegression 的预测性能随 multi_class 参数的影响 ''' X_train,X_test,y_train,y_test=data regr = linear_model.LogisticRegression(multi_class='multinomial',solver='lbfgs') regr.fit(X_train, y_train) print('Coefficients:%s, intercept %s'%(regr.coef_,regr.intercept_)) print('Score: %.2f' % regr.score(X_test, y_test)) def test_LogisticRegression_C(*data): ''' 测试 LogisticRegression 的预测性能随 C 参数的影响 :return: None ''' X_train,X_test,y_train,y_test=data Cs=np.logspace(-2,4,num=100) scores=[] for C in Cs: regr = linear_model.LogisticRegression(C=C) regr.fit(X_train, y_train) scores.append(regr.score(X_test, y_test)) ## 绘图 fig=plt.figure() ax=fig.add_subplot(1,1,1) ax.plot(Cs,scores) ax.set_xlabel(r"C") ax.set_ylabel(r"score") ax.set_xscale('log') ax.set_title("LogisticRegression") plt.show() if __name__=='__main__': X_train,X_test,y_train,y_test=load_data() # 加载用于分类的数据集 print('调用 test_LogisticRegression') print('=========================================================================') test_LogisticRegression(X_train,X_test,y_train,y_test) # 调用 test_LogisticRegression print('调用 test_LogisticRegression_multinomial') print('=========================================================================') test_LogisticRegression_multinomial(X_train,X_test,y_train,y_test) # 调用 test_LogisticRegression_multinomial print('调用 test_LogisticRegression_C') print('=========================================================================') test_LogisticRegression_C(X_train,X_test,y_train,y_test) # 调用 test_LogisticRegression_C

输出:

调用 test_LogisticRegression ========================================================================= Coefficients:[[ 0.39310895 1.35470406 -2.12308303 -0.96477916] [ 0.22462128 -1.34888898 0.60067997 -1.24122398] [-1.50918214 -1.29436177 2.14150484 2.2961458 ]], intercept [ 0.24122458 1.13775782 -1.09418724] Score: 0.97 调用 test_LogisticRegression_multinomial ========================================================================= Coefficients:[[-0.38336889 0.85483782 -2.27270458 -0.98438566] [ 0.34335964 -0.37376524 -0.03024797 -0.86146323] [ 0.04000925 -0.48107258 2.30295255 1.84584889]], intercept [ 8.79968016 2.46967258 -11.26935274] Score: 1.00 调用 test_LogisticRegression_C =========================================================================

注意到C是正则化系数的倒数，越小则正则化的权重越大。可以看到随着C的增大，准确率在上升，增大到一定程度后，准确率保持在一个高水平.

线性判别分析

线性判别分析

LinearDiscriminantAnalysis

discriminant_analysis下提供了LinearDiscriminantAnalysis类实现线性判别分析(这是使用的包不一样了)

LinearDiscriminantAnalysis类介绍

LinearDiscriminantAnalysis类原型为:

class LinearDiscriminantAnalysis(BaseEstimator, LinearClassifierMixin,TransformerMixin): '''Linear Discriminant Analysis. '''
父类

BaseEstimator: """Base class for all estimators in scikit-learn""" LinearClassifierMixin """Mixin for linear classifiers. TransformerMixin '''Mixin class for all transformers in scikit-learn.'''
参数

def __init__(self, solver='svd', shrinkage=None, priors=None, n_components=None, store_covariance=False, tol=1e-4): ---------- solver : string, {'svd','lsqr','eigen'}optional , default:'svd' Solver to use, possible values: - 'svd': 奇异值分解，对于大规模特征数据，推荐用这个 - 'lsqr':最小平方差算法, can be combined with shrinkage. - 'eigen': 特征值分解算法, can be combined with shrinkage. shrinkage : string('auto') or float, optional - None: no shrinkage (default). - 'auto': 根据Ledoit-Wolf引理自动决定shrinkage参数 - float between 0 and 1: 指定shrinkage parameter. Note that shrinkage works only with 'lsqr' and 'eigen' solvers. priors : array, optional, shape (n_classes,) 指定每个类别的先验概率，为None则说明每个类别等可能 n_components : int, optional 指定数据降维后的维度(该值必须小于n_classes-1). store_covariance : bool, optional ,default False If Ture，需要额外计算每个类别的协方差矩阵 . tol : float, optional 指定用 SVD solver判别迭代收敛的阈值
属性

---------- coef_ : array, shape (n_features,) or (n_classes, n_features) 权重系数 intercept_ : array, shape (n_features,) b值 covariance_ : array-like, shape (n_features, n_features) 依次给出每个类别的协方差矩阵. means_ : array-like, shape (n_classes, n_features) 给出每个类别的均值向量. priors_ : array-like, shape (n_classes,) Class priors (sum to 1). scalings_ : array-like, shape (rank, n_classes - 1) Scaling of the features in the space spanned by the class centroids. xbar_ : array-like, shape (n_features,) 整体样本的均值向量. classes_ : array-like, shape (n_classes,) Unique class labels. n_iter_ : array or None, shape (n_targets,) 实际的迭代次数
方法

fit(self, X, y, sample_weight=None) : 训练模型 predict(self, X): ：用模型进行预测 (从父类继承) predict_log_proba(self, X): 返回一个数组，数组的元素依次是X预测为各个类别的概率的对数值 predict_proba(self, X):返回一个数组，数组元素依次是X预测为各个类别的概率值 score(self, X, y, sample_weight=None) :返回预测性能得分。(从父类继承)
使用LinearDiscriminantAnalysis的函数如下

def test_LinearDiscriminantAnalysis(*data): ''' 测试 LinearDiscriminantAnalysis 的用法 :param data: 可变参数。它是一个元组要求其元素依次为：训练样本集、测试样本集、训练样本的标记、测试样本的标记 :return: None ''' X_train,X_test,y_train,y_test=data lda = discriminant_analysis.LinearDiscriminantAnalysis() lda.fit(X_train, y_train) print('Coefficients:%s, intercept %s'%(lda.coef_,lda.intercept_)) print('Score: %.2f' % lda.score(X_test, y_test))

输出:

Coefficients:[[ 6.575853 9.75807593 -14.34026669 -21.39076537] [ -1.98385061 -3.49791089 4.21495042 2.60304299] [ -4.47116022 -6.09542385 9.85886057 18.29330864]], intercept [-15.33097142 0.46730077 -30.53297367] Score: 1.00

还算可以.
检查LinearDiscriminantAnalysis分析后的降维数据集:

def plot_LDA(converted_X,y): ''' 绘制经过 LDA 转换后的数据 :param converted_X: 经过 LDA转换后的样本集(150,3)降维后的表示方法 :param y: 样本集的标记(150,1) :return: None ''' from mpl_toolkits.mplot3d import Axes3D fig=plt.figure() ax=Axes3D(fig) colors='rgb' markers='o*s' for target,color,marker in zip([0,1,2],colors,markers): #zip()方法用在for循环中，支持并行迭代 pos=(y==target).ravel() #由标签转换为bool矩阵，分成不同标签类型的矩阵(150,) 用于拆分数据 X=converted_X[pos,:] # 取出对应pos中True的元素，目的是将降维后的数据拆分为三类X(50,3) ax.scatter(X[:,0], X[:,1], X[:,2],color=color,marker=marker, # X[:,0]取出第一行的数据 label="Label %d"%target) ax.legend(loc="best") fig.suptitle("Iris After LDA") plt.show() def run_plot_LDA(): ''' 执行 plot_LDA 。其中数据集来自于 load_data() 函数 :return: None ''' # X_train(112,4) X_test(38,4) y_train(112,) y_test(38,) X_train,X_test,y_train,y_test=load_data() # X(150,4) X=np.vstack((X_train,X_test)) # 考虑到lda的fit接受的参数问题,这里需要组合X_Train和X_test # Y(150,1) 原本维度为1,所以要先reshape Y=np.vstack((y_train.reshape(y_train.size,1),y_test.reshape(y_test.size,1))) # 组合y_train和y_test lda = discriminant_analysis.LinearDiscriminantAnalysis() lda.fit(X, Y) # converted_X(150,3) 对数据降维了权值lda.coef_(3, 4) lda.intercept_(3,) converted_X=np.dot(X,np.transpose(lda.coef_))+lda.intercept_ # X*权值+偏置b 就是输出值 plot_LDA(converted_X,Y)

输出:

可以看到不同的鸢尾花之间间隔较远，已经聚类到一起了.
检查LinearDiscriminantAnalysis中不同solver对预测性能的影响:

def test_LinearDiscriminantAnalysis_solver(*data): ''' 测试 LinearDiscriminantAnalysis 的预测性能随 solver 参数的影响 :param data: 可变参数。它是一个元组要求其元素依次为：训练样本集、测试样本集、训练样本的标记、测试样本的标记 :return: None ''' X_train,X_test,y_train,y_test=data solvers=['svd','lsqr','eigen'] for solver in solvers: if(solver=='svd'): lda = discriminant_analysis.LinearDiscriminantAnalysis(solver=solver) else: lda = discriminant_analysis.LinearDiscriminantAnalysis(solver=solver, shrinkage=None) lda.fit(X_train, y_train) print('Score at solver=%s: %.2f' %(solver, lda.score(X_test, y_test)))

输出:

Score at solver=svd: 1.00 Score at solver=lsqr: 1.00 Score at solver=eigen: 1.00

数据集不大，效果都差不多.
检查LinearDiscriminantAnalysis中solver=lsqr中引入抖动(相当于正则化):

def test_LinearDiscriminantAnalysis_shrinkage(*data): ''' 测试 LinearDiscriminantAnalysis 的预测性能随 shrinkage 参数的影响 :param data: 可变参数。它是一个元组，这里要求其元素依次为：训练样本集、测试样本集、训练样本的标记、测试样本的标记 :return: None ''' X_train,X_test,y_train,y_test=data shrinkages=np.linspace(0.0,1.0,num=20) scores=[] for shrinkage in shrinkages: lda = discriminant_analysis.LinearDiscriminantAnalysis(solver='lsqr', shrinkage=shrinkage) lda.fit(X_train, y_train) scores.append(lda.score(X_test, y_test)) ## 绘图 fig=plt.figure() ax=fig.add_subplot(1,1,1) ax.plot(shrinkages,scores) ax.set_xlabel(r"shrinkage") ax.set_ylabel(r"score") ax.set_ylim(0,1.05) ax.set_title("LinearDiscriminantAnalysis") plt.show()

输出:

还算可以.

本节程序小结：

# -*- coding: utf-8 -*- """ 广义线性模型 ~~~~~~~~~~~~~~~~~~~~~~~~~~ 线性判别分析 """ import matplotlib.pyplot as plt import numpy as np from sklearn import datasets, discriminant_analysis,cross_validation #注意这里导入模型的包变了 def load_data(): ''' 加载用于分类问题的数据集 ''' iris=datasets.load_iris() X_train=iris.data y_train=iris.target return cross_validation.train_test_split(X_train, y_train,test_size=0.25, random_state=0,stratify=y_train) def test_LinearDiscriminantAnalysis(*data): ''' 测试 LinearDiscriminantAnalysis 的用法 :return: None ''' X_train,X_test,y_train,y_test=data lda = discriminant_analysis.LinearDiscriminantAnalysis() lda.fit(X_train, y_train) print('Coefficients:%s, intercept %s'%(lda.coef_,lda.intercept_)) print('Score: %.2f' % lda.score(X_test, y_test)) def plot_LDA(converted_X,y): ''' 绘制经过 LDA 转换后的数据 :param converted_X: 经过 LDA转换后的样本集(150,3)降维后的表示方法 :param y: 样本集的标记(150,1) :return: None ''' from mpl_toolkits.mplot3d import Axes3D fig=plt.figure() ax=Axes3D(fig) colors='rgb' markers='o*s' for target,color,marker in zip([0,1,2],colors,markers): #zip()方法用在for循环中，支持并行迭代 pos=(y==target).ravel() #由标签转换为bool矩阵，分成不同标签类型的矩阵(150,) 用于拆分数据 X=converted_X[pos,:] # 取出对应pos中True的元素，目的是将降维后的数据拆分为三类X(50,3) ax.scatter(X[:,0], X[:,1], X[:,2],color=color,marker=marker, # X[:,0]取出第一行的数据 label="Label %d"%target) ax.legend(loc="best") fig.suptitle("Iris After LDA") plt.show() def run_plot_LDA(): ''' 执行 plot_LDA 。其中数据集来自于 load_data() 函数 :return: None ''' # X_train(112,4) X_test(38,4) y_train(112,) y_test(38,) X_train,X_test,y_train,y_test=load_data() # X(150,4) X=np.vstack((X_train,X_test)) # 考虑到lda的fit接受的参数问题,这里需要组合X_Train和X_test # Y(150,1) 原本维度为1,所以要先reshape Y=np.vstack((y_train.reshape(y_train.size,1),y_test.reshape(y_test.size,1))) # 组合y_train和y_test lda = discriminant_analysis.LinearDiscriminantAnalysis() lda.fit(X, Y) # converted_X(150,3) 对数据降维了权值lda.coef_(3, 4) lda.intercept_(3,) converted_X=np.dot(X,np.transpose(lda.coef_))+lda.intercept_ # X*权值+偏置b 就是输出值 plot_LDA(converted_X,Y) def test_LinearDiscriminantAnalysis_solver(*data): ''' 测试 LinearDiscriminantAnalysis 的预测性能随 solver 参数的影响 :return: None ''' X_train,X_test,y_train,y_test=data solvers=['svd','lsqr','eigen'] for solver in solvers: if(solver=='svd'): lda = discriminant_analysis.LinearDiscriminantAnalysis(solver=solver) else: lda = discriminant_analysis.LinearDiscriminantAnalysis(solver=solver, shrinkage=None) lda.fit(X_train, y_train) print('Score at solver=%s: %.2f' %(solver, lda.score(X_test, y_test))) def test_LinearDiscriminantAnalysis_shrinkage(*data): ''' 测试 LinearDiscriminantAnalysis 的预测性能随 shrinkage 参数的影响 :return: None ''' X_train,X_test,y_train,y_test=data shrinkages=np.linspace(0.0,1.0,num=20) scores=[] for shrinkage in shrinkages: lda = discriminant_analysis.LinearDiscriminantAnalysis(solver='lsqr', shrinkage=shrinkage) lda.fit(X_train, y_train) scores.append(lda.score(X_test, y_test)) ## 绘图 fig=plt.figure() ax=fig.add_subplot(1,1,1) ax.plot(shrinkages,scores) ax.set_xlabel(r"shrinkage") ax.set_ylabel(r"score") ax.set_ylim(0,1.05) ax.set_title("LinearDiscriminantAnalysis") plt.show() if __name__=='__main__': X_train,X_test,y_train,y_test=load_data() # 产生用于分类的数据集 #test_LinearDiscriminantAnalysis(X_train,X_test,y_train,y_test) # 调用 test_LinearDiscriminantAnalysis # run_plot_LDA() # 调用 run_plot_LDA # test_LinearDiscriminantAnalysis_solver(X_train,X_test,y_train,y_test) # 调用 test_LinearDiscriminantAnalysis_solver test_LinearDiscriminantAnalysis_shrinkage(X_train,X_test,y_train,y_test) # 调用 test_LinearDiscriminantAnalysis_shrinkage

参考资料

《机器学习》周志华
《Python大战机器学习》华校专
https://www.zhihu.com/question/20447622 关于最大似然估计的通俗理解