机器学习Python实现之线性模型

本文将详细解释线性分类的几个常用模型：线性回归、对数回归、对数几率回归，并简要介绍其优化方法。文末附有Python代码实现。如果问题，欢迎留言交流～

线性回归（linear regression）

模型

设样本表示为d维列向量x，其标记为y，记x^=(x1)，线性回归的基本模型：y=f(x)=wTx+b=w^Tx^。

令X^=⎛⎝⎜⎜⎜⎜xT0xT1...xTn−1111⎞⎠⎟⎟⎟⎟，即每行是一个x^T。则给定w^，训练集上的损失函数为：∥X^w^−Y∥2，优化目标argminw^∥X^w^−Y∥2。

闭式解方法：w^=(X^TX^)−1X^TY

梯度下降方法：w^t+1=w^t−γt√X^T(X^w^t−Y)

代码参数

在文末的代码中，regression.py文件中的

linear_regression(X,Y,method)

函数用于线性回归。该函数求解了argminw^∥X^w^−Y∥2问题。

参数：

X

：训练样本的矩阵表示末尾添加全1的一列，即上文中的X^；

Y

：标记（label）的矩阵表示（n×1维，n为样本个数）；

method

：方法参数，字符串，若为

cloesd_form

，则用闭式解方法求解。 若为

gd

，则用梯度下降方法求解，如果不想使用梯度下降的默认值，还可以设置

gamma

eps

max_iter

等参数，分别表示梯度下降等步长系数，终止条件，最大迭代次数。

返回：

回归系数w^，为列向量。

对数回归（log regression）

模型

注意！对数回归不是对数几率回归！

对数回归的模型很简单，lny=w^Tx^，求解时，只需将真实标记取对数，即可按照线性回归模型来做。注意，此时要求标记y>0。

代码

在文末的代码中，regression.py文件中的

log_regression(X,Y,method)

函数用于线性回归。该函数求解了argminw^∥X^w^−lnY∥2问题。

参数：

X

：训练样本的矩阵表示末尾添加全1的一列，即上文中的X^；

Y

：标记（label）的矩阵表示（n×1维，n为样本个数）；

method

：方法参数，字符串，若为

cloesd_form

，则用闭式解方法求解。 若为

gd

，则用梯度下降方法求解，如果不想使用梯度下降的默认值，还可以设置

gamma

eps

max_iter

等参数，分别表示梯度下降等步长系数，终止条件，最大迭代次数。

返回：

回归系数w^，为列向量。

对数几率回归（logistic regression）

模型

对数几率回归模型相对复杂，针对二分类问题，标记为1或者0：p(y=1 |x)=11+exp(−w^Tx^)。可以变化为lnp(y=1)1−p(y=1)=w^Tx^。若将y看作x为正例的可能性，1−y看作是负例的可能性，则y1−y是两者的相对可能性，称之为“几率”，lny1−y称之为“对数几率”。

那么，如何优化呢？我们用极大似然法来解出w^。

显然，p(y=1 |x^,w^)=exp(w^Tx^)1+exp(w^Tx^)，p(y=0 |x^,w^)=11+exp(w^Tx^)。给定w^，则模型的似然函数为：

l(w^)=∑i=1nlnp(yi |xi^,w^)=∑i=1nln(yip(1 |xi^,w^)+(1−yi)p(0 |xi^,w^))=∑i=1nln(yiexp(w^Txi^)1+exp(w^Txi^)+(1−yi)11+exp(w^Txi^))=∑i=1nyiw^Txi^−ln(1+exp(w^Txi^)) (考虑到yi取值只能为0或1，故该式成立)

最大化似然函数l(w^)等价于最小化 −l(w^)，幸运的是，−l(w^)是连续可导的凸函数，故有多种方法求解。

代码

在文末的代码中，regression.py文件中的

logistic_regression(X,Y,method)

函数用于线性回归。该函数求解了argmaxw^l(w^)问题。

参数：

X

：训练样本的矩阵表示末尾添加全1的一列，即上文中的X^；

Y

：标记（label）的矩阵表示（n×1维，n为样本个数）；

method

：方法参数，字符串。若为

gd

，则用梯度下降方法求解，如果不想使用梯度下降的默认值，还可以设置

gamma

eps

max_iter

等参数，分别表示梯度下降等步长系数，终止条件，最大迭代次数。

返回：

回归系数w^，为列向量。

Python代码

文件名：regression.py

#!/usr/bin/python
# -*- coding: UTF-8 -*-
import os
import numpy as np

def linear_regression(X, Y, method = "gd", gamma=0.001, eps=0.0001, max_iter=10000):
#求使得 ||Xw-Y||^2最小的w
#X, Y是样本和label,是np.array类型
#如果method=="closed_form"，则用闭式解方法求出w，不需要后面三个参数
#如果method=="gd"，则用梯度下降方法求出w
if method == "closed_form":
return linear_regression_by_closed_form(X, Y)
if method == "gd":
return linear_regression_by_gd(X, Y, gamma, eps, max_iter)
print ("args error")

def linear_regression_by_closed_form(X, Y): #闭式解求线性回归
w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y);
return w

def linear_regression_by_gd(X, Y, gamma=0.001, eps=0.0001, max_iter=10000): #梯度下降求线性回归
pre_w = np.array(np.ones((X.shape[1], 1)))
cur_w = np.array(np.zeros((X.shape[1], 1)))
count = 1
while (cur_w - pre_w).T.dot(cur_w - pre_w) > eps and count < max_iter:
pre_w = cur_w
cur_w = cur_w - gamma / np.sqrt(count) * X.T.dot( X.dot(cur_w) - Y)
count += 1
return cur_w

def log_regression(X, Y, method="gd", gamma=0.001, eps=0.0001, max_iter=10000):
ln_Y = np.log(Y)
return linear_regression(X, ln_Y, method, gamma, eps, max_iter)

def logistic_regression(X, Y, method = "gd", gamma=0.001, eps=0.0001, max_iter=10000):
#X, Y是样本和label,是np.array类型
#要求，label的取值为0或1
#如果method=="gd"，则用梯度下降方法求出w
if method == "gd":
return logistic_regression_by_gd(X, Y, gamma, eps, max_iter)
print ("args error")

def logistic_regression_by_gd(X, Y, gamma=0.001, eps=0.0001, max_iter=10000) :
pre_w = np.array(np.ones((X.shape[1], 1)))
cur_w = np.array(np.zeros((X.shape[1], 1)))
count = 1
while (cur_w - pre_w).T.dot(cur_w - pre_w) > eps and count < max_iter:
pre_w = cur_w
gradient = X.T.dot(sigmoid(X.dot(cur_w)) - Y)
cur_w = cur_w - gamma / np.sqrt(count) * gradient
count += 1
return cur_w

def sigmoid(x):
return 1.0 / (1 + np.exp( - x))

文件名：test_regression.py

其中，

d = np.loadtxt('./data/Iris.data', delimiter=',')

中Iris.data文件来自http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data，取前100行，并且将字符串表示的label改为0，1。

#!/usr/bin/python
# -*- coding: UTF-8 -*-
import os
import numpy as np
from regression import *

if __name__ == "__main__":
##分类数据集Iris.txt来自http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data，取前100行，并且将字符串表示的label改为0，1
##读取矩阵文件 ./data/Iris.data，存放一个矩阵，前m-1列是X，最后一列是Y
d = np.loadtxt('./data/Iris.data', delimiter=',')
np.random.shuffle(d)
n, m = d.shape
train_num = 2 * n / 3
test_num = n - train_num
train_data = d[0:train_num,0: (m-1)]
train_data = np.c_[train_data, np.ones((train_num,1))] #回归的时候会有常数项，故此处加了一列
train_label = d[0:train_num,m-1].reshape(train_num,1) #python中一维数组默认是行向量，需要reshape函数转换
test_data = d[train_num:n,0: (m-1)]
test_data = np.c_[test_data, np.ones((test_num, 1))]
test_label = d[train_num:n,m-1].reshape(test_num,1)

print ("linear regression")
print ("\t training start ...")
threshold = (max(train_label) + min(train_label)) / 2
gamma, eps, max_iter = 0.001, 0.00001, 10000
w = linear_regression(train_data, train_label, 'gd', gamma, eps, max_iter)
print ("\t training done !")
train_y_predict = train_data.dot(w)
test_y_predict = test_data.dot(w)
# print ("\t train mean squre error\t: %f"%((train_y_predict-train_label).T.dot(train_y_predict-train_label)/train_num)[0,0])
print ("\t train predict error\t: %f"%(sum( abs( ((train_y_predict > threshold) + 0) - ((train_label > threshold) + 0) ))[0] / (train_num + 0.0)))
# print ("\t test  mean squre error\t: %f"%((test_data.dot(w)-test_label).T.dot(test_data.dot(w)-test_label)/test_num)[0,0])
print ("\t test predict error \t: %f"%(sum( abs( ((test_y_predict > threshold) + 0) - ((test_label > threshold) + 0) ))[0] / (test_num + 0.0)))

print ("\nlog regression")
print ("\t training start ...")
min_label, max_label = min(train_label), max(train_label)
train_label = train_label - min_label + 1 #保证label>0，才可以取对数
test_label = test_label - min_label + 1 #保证label>0，才可以取对数
threshold = (np.log(max(train_label)) + np.log(min(train_label))) / 2
gamma, eps, max_iter = 0.001, 0.00001, 10000
w = log_regression(train_data, train_label, 'gd', gamma, eps, max_iter)
train_y_predict = train_data.dot(w)
test_y_predict = test_data.dot(w)
print ("\t training done")
# print ("\t train mean squre error\t: %f"%((np.exp(train_y_predict)-train_label).T.dot(np.exp(train_y_predict)-train_label)/train_num)[0,0])
print ("\t train predict error\t: %f"%(sum( abs( ((train_y_predict > threshold) + 0) - ((train_label > threshold) + 0) ))[0] / (train_num + 0.0)))
# print ("\t test  mean squre error\t: %f"%((np.exp(test_data.dot(w))-test_label).T.dot(np.exp(test_data.dot(w))-test_label)/test_num)[0,0])
print ("\t test predict error \t: %f"%(sum( abs( ((test_y_predict > threshold) + 0) - ((test_label > threshold) + 0) ))[0] / (test_num + 0.0)))

print ("\nlogistic regression")
print ("\t training start ...")
min_label, max_label = min(train_label), max(train_label)
train_label = (train_label - min_label) / (max_label - min_label) #将label变为0，1
test_label = (test_label - min_label) / (max_label - min_label) #将label变为0，1
threshold = 0.5
gamma, eps, max_iter = 0.001, 0.00001, 10000
w = logistic_regression(train_data, train_label, 'gd', gamma, eps, max_iter)
train_y_predict = sigmoid(train_data.dot(w))
print ("\t training done")
# print ("\t train mean squre error\t: %f"%((train_y_predict-train_label).T.dot(train_y_predict-train_label)/train_num)[0,0])
print ("\t train predict error \t: %f"%(sum( abs( ((train_y_predict > threshold) + 0) - ((train_label > threshold) + 0) ))[0] / (train_num + 0.0)))

test_y_predict = sigmoid(test_data.dot(w))
# print ("\t test  mean squre error\t: %f"%((test_y_predict-test_label).T.dot(test_y_predict-test_label)/test_num)[0,0])
print ("\t test predict error \t: %f"%(sum( abs( ((test_y_predict > threshold) + 0) - ((test_label > threshold) + 0) ))[0] / (test_num + 0.0)))

运行test_regression.py，结果示例：

linear regression
training start ...
training done !
train predict error    : 0.000000
test predict error     : 0.000000

log regression
training start ...
training done
train predict error    : 0.515152
test predict error     : 0.470588

logistic regression
training start ...
training done
train predict error    : 0.000000
test predict error     : 0.000000

可以看出，在

Iris

数据集上，线性回归和对数几率回归效果比较好，对数回归效果比较差。