1. 线性模型

发表信息: by Creative Commons Licence

目录:

线性模型一般用于回归分类

一、线性回归

如图所示,回归主要是完成对数据的拟合。

(一)回归模型

对数据进行线性拟合

这里设定样本参数维度为\(d\),样本数量维度是\(n\)。因而\(X\)的维度是\(n \times d\),\(w\)的维度是\(d \times 1\),\(Y\)的维度是\(n \times 1\)。

(1)针对单个样本

\[y_i = a(W^Tx_i) = w^T x_i = w_0 + w_1 x_{i,1} + \cdots + w_d x_{i,d} = \begin{pmatrix}w_0 & w_1 & \cdots & w_d \end{pmatrix} \begin{pmatrix} 1 \\ x_{i,1} \\ \vdots \\ x_{i,d} \end{pmatrix}\]

(2)针对\(n\)个样本

\[Y = a(Xw) = Xw = \begin{pmatrix} w_0 + w_1 x_{1,1} + \cdots + w_d x_{1,d} \\w_0 + w_1 x_{2,1} + \cdots + w_d x_{2,d} \\ \vdots \\ w_0 + w_1 x_{n,1} + \cdots + w_d x_{n,d}\end{pmatrix} = \begin{pmatrix}1 & x_{11} & \cdots & x_{1d} \\ 1 & x_{21} & \cdots & x_{2d} \\ \vdots \\ 1 & x_{n1} & \cdots & x_{nd} \\ \end{pmatrix} \begin{pmatrix}w_0 \\ w_1 \\ \vdots \\ w_d \\ \end{pmatrix}\]

(二)损失函数

对线性拟合效果进行评估

损失函数(loss function)是用来估量模型的预测值和真实值的不一致程度,不一致程度越小,模型拟合效果越好。

因此目标一般都是最小化损失函数。

换句话说,如何选择损失函数就是如何去选择最优的拟合直线。

线性回归中最常使用的是最小二乘法(OLS),也称为平方损失函数。

(1)最小二乘法

\[\begin{eqnarray} L(w) & = & \frac{1}{2n} \sum_{i=1}^n[a(x_i)-y_i]^2 \\ & = & \frac{1}{2n} \sum_{i=1}^n[w_0 + w_1 x_{i,1} + \cdots + w_d x_{i,d}-y_i]^2\\ & = & \frac{1}{2n} \sum_{i=1}^n[w^T x_i-y_i]^2\\ & = & \frac{1}{2n} \begin{pmatrix}w^T x_1-y_1 & \cdots & w^T x_n-y_n \end{pmatrix} \begin{pmatrix}w^T x_1-y_1 \\ \cdots \\ w^T x_n-y_n \end{pmatrix} \\ & = & \frac{1}{2n} \begin{bmatrix}Xw-Y \end{bmatrix}^T \begin{bmatrix} Xw-Y \end{bmatrix}\\ & = & \frac{1}{2n} \begin{Vmatrix}Xw-Y \end{Vmatrix}^2\\ \end{eqnarray}\]

之所以写成\(\frac{1}{2n}\),而不是\(\frac{1}{n}\),是为了方便求导时,可以与2相消。

(2)代码

# 这里假设X,Y,w都是np.mat格式
def loss(X, Y, w):
	n = Y.shape[0]
	return 1.0/(2*n) * (X*w - Y).T * (X*w - Y)

(三)梯度下降

求解参数w需要最小化损失函数

针对线性回归,最简单粗暴求解方法是,直接求导,令导数等于0。

\[\frac{\partial L(w)}{\partial w} = \frac{1}{n} X^T (Xw-Y) = 0\] \[\begin{eqnarray} X^T (Xw-Y) & = & 0 \\ X^T X w &= & X^T Y\\ w &= & (X^T X)^{-1} X^T Y \end{eqnarray}\]

直接针对矩阵进行运算,与梯度下降法相比,无需多次迭代,但是当数据量很大,即\(X\)过大时,矩阵运算需要很大的计算复杂度。

如果矩阵X是奇异矩阵,求逆则会非常不稳定。因此人们常用梯度下降法作为优化方法,通过不断迭代来得到最小化损失函数的最优解。

(1)梯度下降

第一步:初始化 \(w^0\)

第二步:计算梯度 \(\nabla L(w) = \begin{pmatrix} \frac{\partial L(w)}{\partial w_0} & \cdots & \frac{\partial L(w)}{\partial w_d} \end{pmatrix} ^T\)

第三步:更新参数 \(w^{t} = w^{t-1} - \nabla L(w^{t-1}) \cdot \eta\)

到达终值条件以后,更新停止,最终的\(w\)为线性回归模型的最终参数。

\(\eta\)是学习率。

以本节为例,得到 :

\[\nabla L(w) = \begin{pmatrix} \frac{\partial L(w)}{\partial w_0} & \cdots & \frac{\partial L(w)}{\partial w_d} \end{pmatrix} ^T = \frac{1}{n} X^T (Xw-Y)\]

非矩阵的写法是:

\[\begin{eqnarray} \frac{\partial L(w)}{\partial w_j} & = & \frac{1}{2n} \frac{\partial \sum_{i=1}^n[w_0 + w_1 x_{i,1} + \cdots + w_d x_{i,d}-y_i]^2}{\partial w_j}\\ & = & \frac{1}{n} \sum_{i=1}^{n}(w^T x_i-y_i)x_{i,j} \end{eqnarray}\]

(2)代码

python 中矩阵运行要快于循环运算,因此代码采用上述的矩阵结果。

# 这里假设X,Y,w都是np.mat格式

# 计算梯度值
def OLS_dev(X, Y, w):
	n = Y.shape[0]
	return 1.0/(n) * X.T * (X*w - Y)
	
# 进行梯度下降:
def gradientDescent(X, Y, w, eta, num_iters):
	n = Y.shape[0]  # 样本量
	d = w.shape[0]  # 参数量
	
	L = np.zeros((num_iters, 1))  # 保存每次迭代的损失值
	
	for i in range(num_iters):
		w -= eta * OLS_dev(X, Y, w)   # 前面定义的梯度值
		L[i] = loss(X, Y, w)          # 前面定义的loss值
	
	GD = {'w':w,
		'L':L}
	
	return GD

(四)整合

上面分别实现了损失值、梯度值和梯度下降的实现,这里整合这三部分,以类(class)的形式,实现这一模型:

import numpy as np 

class linearRegression():
	def __init__(self, X, Y, w, n, d):
		self.X = X 
		self.Y = Y
		self.w = w 
		self.n = n
		self.d = d
		
	def loss(self):
		return 1.0/(2*self.n) * (self.X*self.w - self.Y).T * (self.X*self.w - self.Y)
		
	def OLS_dev(self):
		return 1.0/(self.n) * self.X.T * (self.X*self.w - self.Y)
		
	def gradientDescent(self, eta, num_iters):
		L = np.zeros((num_iters, 1))  # 保存每次迭代的损失值
		
		for i in range(num_iters):
			self.w -= eta * OLS_dev()   # 前面定义的梯度值
			L[i] = loss()          # 前面定义的loss值
		
		GD = {'w':self.w,
			'L':L}
		
		return GD
	
	def predict(self):
		y_predict = self.X * GD['w']
		return y_predict

我们获取一些数据来实现这一模型:

from sklearn.datasets import make_regression
import numpy as np
import matplotlib.pyplot as plt

# 从sklearn库生成需要的回归模型
X, Y = make_regression(n_samples=1000, n_features=1, noise=0.5)  # 这里参数即n=1000, d=1,为了便于后续画图展示效果,这里d只取1

# 处理数据,将X,Y转换成矩阵格式
# 先查看X,Y的维度
print("X的维度是:",X.shape, "\nY的维度是:", Y.shape)

X = np.insert(X, 0, values=1, axis=1) # 为了添加截距项
X = np.mat(X)
Y = np.mat(Y).T  #为了矩阵维度与之前设定一致

# 初始化参数w
n = X.shape[0]
d = X.shape[1]
w = np.random.random((d+1, 1))

# 定义线性回归模型
class linearRegression():
	def __init__(self, X, Y, w, n, d):
		self.X = X 
		self.Y = Y
		self.w = w 
		self.n = n
		self.d = d
		
	def loss(self):
		return 1.0/(2*self.n) * (self.X*self.w - self.Y).T * (self.X*self.w - self.Y)
		
	def OLS_dev(self):
		return 1.0/(self.n) * self.X.T * (self.X*self.w - self.Y)
		
	def gradientDescent(self, eta, num_iters):
		L = np.zeros((num_iters, 1))  # 计算每次迭代的损失值
		
		for i in range(num_iters):
			self.w -= eta * self.OLS_dev()   # 前面定义的梯度值
			L[i] = self.loss()          # 前面定义的loss值
		
		GD = {'w':self.w,
			'L':L}
		return GD
		
	def predict(self):
		y_predict = self.X * self.w
		print("train accuracy:{}%".format(100 - 1./n*np.sqrt(np.sum(np.square(y_predict-Y)))*100))
		return y_predict
		
lr = linearRegression(X, Y, n, d)
GD = lr.gradientDescent(eta=0.1, num_iters=1000)
y_predict = lr.predict()


# 绘图,要将矩阵转换成list形式才可以画图
plt.scatter(X[:,1].tolist(), Y.tolist())
plt.plot(X.tolist(), y_predict.tolist(), c='r')
plt.show()

得到结果:

train accuracy:83.580213%

对于训练过程损失值的变化:

L = GD["L"]
plt.plot(L)
plt.show()

(五)scikit-learn库实现

Python提供了很方便的实现方式:

from sklearn.linear_model import LinearRegression
from sklearn.datasets import make_regression
import matplotlib.pyplot as plt
import numpy as np

# 生成数据
X, Y = make_regression(n_samples=1000, n_features=1, noise=5)  
X = np.insert(X, 0, values=1, axis=1) 
X = np.mat(X)
Y = np.mat(Y).T  

# 回归
lr = LinearRegression()
lr.fit(X, Y)   				
Y_predict = lr.predict(X)

# 结果展示
print("train accuracy:{}%".format(100 - 1./n*np.sqrt(np.sum(np.square(Y_predict-Y)))*100))
plt.scatter(X[:,1].tolist(), Y.tolist())
plt.plot(X.tolist(), Y_predict.tolist(), c='r')
plt.show()

train accuracy:84.394867%

接下来介绍线性分类

二、线性分类

分类比较直观,如下图所示:

(一)分类模型

(1)二分类模型

针对单个样本

\[a(x_i) = sign(\sigma(w^T x_i))\]

例如二项逻辑回归

针对单个样本

\[\hat{y_i} = P(y_i=1|x_i,w) = \frac{1}{1+e^{-w^T x_i}}\]

针对多个样本

\[\hat{Y} = P(Y=1|X,w) = \frac{1}{1+e^{-Xw}}\]

X,w,Y维度设定和线性回归部分一致

(2)多分类模型

针对单个样本

\[a(x_i) = arg max(\sigma(w^T_k x_i))\]

例如softmax

\[z_1 = w^T_1 x_i, z_2 = w^T_2 x_i, \cdots, z_k = w^T_k x_i\] \[\sigma(z) = \begin{pmatrix} \frac{e^{z_1}}{\sum_{k=1}^K} & \cdots & \frac{e^{z_K}}{\sum_{k=1}^K} \end{pmatrix}\]

(3)代码

这部分实现二项逻辑回归

def sigmoid(X, w):
	return 1./(1.+np.exp(-X*w))

(二)损失函数

常见损失函数有几种:

(1)squared loss

\[(w^T x_i - 1)^2\]

即:\(\frac{1}{2n}\sum_{i=1}^n(w^T x_i - 1)^2\)

这是针对二分类问题

(2)分类准确度

\[[a(x_i)=y_i]\]

即:\(\frac{1}{n}\sum_{i=1}^n [a(x_i)=y_i]\)

(3)交叉熵(log损失函数)

\[-\sum_{k=1}^K[y=k]log\frac{e^{z_k}}{\sum_{j=1}^k e^{z_j}} = -log(\frac{e^{z_y}}{\sum_{j=1}^k e^{z_j}})\]

即:\(\frac{1}{n}\sum_{i=1}^n log\frac{e^{z_{y_i}}}{\sum_{j=1}^k e^{z_j}} = \frac{1}{n}\sum_{i=1}^n log\frac{e^{w_{y_i}^T x_i}}{\sum_{j=1}^k e^{w_j^T x_i}}\)

针对二项逻辑回归,损失函数是:

令\(\hat{y_i} = \pi(x_i) = P(y_i=1 \mid x_i,w)\),则 \(1 - \pi(x_i) = 1-P(y_i=1 \mid x_i,w) = P(y_i=0 \mid x_i,w)\)

则,

\[\begin{eqnarray} L(w) & = & -\frac{1}{n}\sum_{i=1}^{N}[y_i log \pi(x_i) + (1 - y_i)log(1-\pi(x_i))]\\ & = & -\frac{1}{n}\sum_{i=1}^{N}[y_i \cdot w^T x_i - log(1+ e^{w^T x_i}))] \\ & = & -\frac{1}{n}\sum_{i=1}^{N}[y_i log \hat{y_i} + (1 - y_i) log(1-\hat{y_i})]\\ & = & -\frac{1}{n} [Y^T log(\hat{Y}) + (1 - Y)^T log(1-\hat{Y})] \end{eqnarray}\]

(4)代码

这里实现二项逻辑回归损失函数

def loss(Y, X, w):
	n = Y.shape[0]
	Y_predict = sigmoid(X, w)
	return -(1./n)*(Y.T * np.log(Y_predict) + (1 - Y).T * np.log(1 - Y_predict))

(三)梯度下降

与前文一脉相承,这里计算二项逻辑回归的梯度

(1)梯度

\[\nabla L(w) = \begin{pmatrix} \frac{\partial L(w)}{\partial w_0} & \cdots & \frac{\partial L(w)}{\partial w_d} \end{pmatrix} ^T = -\frac{1}{n} X^T(Y - \hat{Y})\]

(2)代码


def LR_dev(Y, X, w):
	n = Y.shape[0]
	Y_predict = sigmoid(X, w)
	return -(1./n)* X.T *(Y - Y_predict)

def gradientDescent(X, Y, w, eta, num_iters):
	n = Y.shape[0]  # 样本量
	d = w.shape[0]  # 参数量
	
	L = np.zeros((num_iters, 1))  # 计算每次迭代的损失值
	
	Y_predict = sigmoid(X,w)
	for i in range(num_iters):
		w -= eta * LR_dev(Y, Y_predict, X)   # 前面定义的梯度值
		Y_predict = sigmoid(X,w)
		L[i] = loss(Y, Y_predict)          # 前面定义的loss值
	
	GD = {'w':w,
		'L':L}
	return GD

(四)整合

我们获取一些数据来实现这一模型:

from sklearn.datasets import make_classification
import numpy as np
import matplotlib.pyplot as plt

# 从sklearn库生成需要的回归模型
X, Y = make_classification(n_samples=1000, n_features=20,  n_classes=2) # 意思是这里d=20,n=1000,类别是2  

# 处理数据,将X,Y转换成矩阵格式
# 先查看X,Y的维度
print("X的维度是:",X.shape, "\nY的维度是:", Y.shape)

X = np.insert(X, 0 , values=1, axis=1)  # 添加截距项
X = np.mat(X)
Y = np.mat(Y).T  #为了矩阵维度与之前设定一致

# 初始化参数w
n = X.shape[0]
d = X.shape[1]
w = np.random.random((d, 1))

# 定义线性回归模型
class logisticRegression():
	def __init__(self, X, Y, w, n, d):
		self.X = X 
		self.Y = Y
		self.w = w 
		self.n = n
		self.d = d
	
	def sigmoid(self):
		return 1./(1.+np.exp(-self.X*self.w))
		
	def loss(self):
		Y_predict = self.sigmoid()
		return -(1./n)*(self.Y.T * np.log(Y_predict) + (1 - self.Y).T * np.log(1 - Y_predict))
	
	def LR_dev(self):
		Y_predict = self.sigmoid()
		return -(1./n)* self.X.T * (self.Y - Y_predict)

	def gradientDescent(slef, eta, num_iters):
		L = np.zeros((num_iters, 1))  # 计算每次迭代的损失值

		for i in range(num_iters):
			self.w -= eta * self.LR_dev()   # 前面定义的梯度值
			L[i] = self.loss()          # 前面定义的loss值
			
		GD = {'w':self.w,
			'L':L}
		return GD
		
	def predict(self):
		y_predict = self.sigmoid()
		y_predict[y_predict>0.5]=1
		y_predict[y_predict<=0.5]=0
		print("train accuracy:{}%".format(100 - 1./n*np.sqrt(np.sum(np.square(y_predict-Y)))*100))
		return y_predict
		
lr = logisticRegression(X, Y, n, d)
GD = lr.gradientDescent(eta=0.1, num_iters=1000)
y_predict = lr.predict()


# 绘图,要将矩阵转换成list形式才可以画图
plt.scatter(X.tolist(), Y.tolist())
plt.scatter(X.tolist(), y_predict.tolist())
plt.show()

    pos = np.where(y==1)    #找到y==1的坐标位置
    neg = np.where(y==0)    #找到y==0的坐标位置
    #作图
    plt.figure(figsize=(15,12))
    plt.plot(X[pos,0],X[pos,1],'ro')        # red o
    plt.plot(X[neg,0],X[neg,1],'bo')        # blue o
    plt.title(u"两个类别散点图",fontproperties=font)
    plt.show()

结果如下图:

train accuracy:83.580213%

对于训练过程损失值的变化:

L = GD["L"]
plt.plot(L)
plt.show()

(五)scikit-learn库实现

Python提供了很方便的实现方式:

from sklearn.linear_model import LogisticRegression

lr = LogisticRegression()
lr.fit(X_train, Y)

Y_predict = model.predict(X_test)

right = sum(predict==y_test)

predict = np.hstack((predict.reshape(-1,1),y_test.reshape(-1,1)))   # 将预测值和真实值放在一块,好观察       
print ('测试集准确率:%f%%'%(right*100.0/predict.shape[0]))       #计算在测试集上的准确度 

(六)scipy优化

这里使用Python中scipy库对梯度部分进行优化计算。

使用优化算法:拟牛顿法(Broyden-Fletcher-Goldfarb-Shanno)。对应库中的fmin_bfgs。

result = optimize.fmin_bfgs(costFunction, initial_theta, fprime=gradient, args=(X,y,initial_lambda))    

costFunction是自己实现的一个求代价的函数,initial_theta表示初始化的值,fprime指定costFunction的梯度,args是其余测参数,以元组的形式传入,最后会将最小化costFunction的theta返回

三、拓展

(一)二项逻辑回归函数以及交叉熵的由来

这个由来与似然函数有着很大的关系,见这篇(附上链接

在逻辑回归的推导中,它假设样本服从伯努利分布,然后求得满足该分布的似然函数,接着对对数求极值等。

最大化似然函数,即最小化负的似然函数,那么负的似然函数就是逻辑回归的损失函数。

取对数是为了方便计算极大似然估计,因为在MLE中,直接求导比较困难,所以通常都是先取对数再求导找极值点。

损失函数\(L(Y, P(Y \mid X))\)表达的是样本X在分类Y的情况下,使概率\(P(Y \mid X)\)达到最大值。

换言之,就是利用已知的样本分布,找到最有可能,即最大概率导致这种分布的参数值。或者说什么样的参数才能使我们观测到目前这组数据的概率最大。

因为log函数是单调递增的,所以\(logP(Y \mid X)\)也会达到最大值,因此在前面加上负号之后,最大化\(P(Y \mid X)\)就等价于最小化L了

我们已知的是某个事件是否发生(或者关于这个事物的分类),因而可以计算其发生的几率(常见对事几率),根据这个已知信息,去获取我们希望得到的参数w.

写成公式即:

\[log \frac{P(Y=1|x)}{1 - P(Y=1|x)} = wx\]

进行推导可得:

\[\frac{P(Y=1|x)}{1 - P(Y=1|x)} = e^{wx}\] \[\frac{1}{P(Y=1|x)}-1 = e^{-wx}\] \[P(Y=1|x)-1 = \frac{1}{1+e^{-wx}}\]

这是线性函数,\(wx\)的值越接近 \(+\infty\),概率值越接近1;\(wx\)的值越接近 \(-\infty\),概率值越接近0。

与此对应的似然函数可以写成(服从二项分布(伯努利分布(0-1分布))):

\[\prod_{i=1}^N [\pi(x_i)^{y_i} [1 - \pi(x_i)]^{1-{y_i}]\]

其中,\(\pi(x) = P(Y=1 \mid x)\),\(1 - \pi(x) = 1 - P(Y=1 \mid x) = P(Y=0 \mid x)\)

取对数后,得到对数似然函数:

\[L(w) = \sum_{i=1}^{N}[y_i log \pi(x_i) + (1 - y_i)log(1-\pi(x_i))]\]

除了上面似然的解释,使用交叉熵,而不使用线性回归部分的损失函数,还有一个原因是:

因为线性回归的代价函数可能是非凸的,对于分类问题,使用梯度下降很难得到最小值,而交叉熵是凸函数。

对于凸与非凸问题,我们绘制函数图形就可以很直观的看出来:

因为\(y=0或1\),所以交叉熵函数主要是由\(log\pi(x)\)和\(log(1-\pi(x))\)组成:

\(log\pi(x)\)的图像是(即y=1时):

当\(-\pi(x)\)趋于1,loss趋于0,意味着预测值\(\hat{y}\)趋于1,与实际结果一致。

\(-log(1-\pi(x))\)的图像是(即y=0时):

同理,当\(-\pi(x)\)趋于0,loss趋于0,意味着预测值\(\hat{y}\)趋于0,与实际结果一致。

人们也常这样写交叉熵:

\[\begin{cases} J(\theta) = \frac{1}{m}\sum_{i=1}^m cost(h_{\theta}(x^{(i)}), y^{(i)})\\ cost(h_{\theta}(x), y) = \begin{cases} -log(h_{\theta}(x)) & y=1 \\ -log(1 - h_{\theta}(x)) & y=0 \\ \end{cases} \\ n/2, & \text{if $n$ is even} \\ 3n+1, & \text{if $n$ is odd} \\ \end{cases}\]

综合起来为:

\[J(\theta) = - \frac{1}{m}\sum_{i=1}^m [y^{(i)}log(h_{\theta}(x^{(i)})) + (1-y^{(i)})log(1-h_{\theta}(x^{(i)}))]\]

这里的\(h_{\theta}(x^{(i)})\)就是前文中的\(\pi(x_i)\),\(\theta\)即\(w\)即

\[h_{\theta}(x^{(i)}) = \pi(x_i)= \frac{1}{1+e^{-\theta x_i}} = \frac{1}{1+e^{-w x_i}}\]

(二)平方损失函数

平方损失函数可以通过线性回归在假设样本是高斯分布的条件下推导得到。为什么假设成高斯分布呐,原因来自中心极限定理)

通过极大似然估计(MLE)可以推导出最小二乘式子。

最小二乘的基本原则是:最优拟合直线应该是使各点到回归直线的距离和最小的直线,即平方和最小。

换言之,OLS是基于距离的,而这个距离使用最多的是欧几里得距离。

为什么选择欧式距离作为误差度量呢?

简单,计算方便;欧氏距离是一种很好的相似性度量标准;在不同的表示域变换后特征性质不变。

通常说线性有两种情况:

一种是因变量

()标准化

以下是针对此篇的简述,详细点击这里(后面附链接)

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
scaler.fit(X)
X_scaler = scaler.fit_transform(X)

()正则化

以下是针对此篇的简述,详细点击这里(后面附链接)

为了防止过拟合,一般会在损失函数后面加上正则项(或称惩罚项),以交叉熵为例得:

\[L(w) = \sum_{i=1}^{N}[y_i log \pi(x_i) + (1 - y_i)log(1-\pi(x_i))] + \frac{\lambda}{2m} \sum_{j=1}^n w^2_j\]

需要注意的是,这里\(j\)是从1开始的,因为\(w_0\)是截距项,也被视为常数项,一般不做正则化出来。

代码部分将新的损失函数重新求导,替换原先的导数就可以了。

()训练集&测试集划分

from sklearn.cross_validation import train_test_split

X_train , X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

()项目练习

(1)手写数字识别

https://github.com/lawlite19/MachineLearning_Python/tree/master/LogisticRegression

~~~

参考文献:

  1. https://github.com/lawlite19/

  2. https://mp.weixin.qq.com/s/u07yzOfJUVaLbqVmt7IOnQ

  3. https://blog.csdn.net/shenxiaoming77/article/details/51614601