5.4多元逻辑回归实战
5.4.1多元逻辑回归实战之预测病马的死亡率
本次实战内容，将使用Logistic回归来预测患疝气病的马的存活问题。原始数据集下载地址：
archive.ics.uci.edu/ml/datasets…
这里的数据包含368个样本和28个特征。我并非育马专家，从一些文献中了解到，疝病是描述马胃肠痛的术语。然而，这种病不一定源自马的胃肠问题，其他问题也可能引发马疝病。该数据集中包含了医院检测马疝病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。
另外需要说明的是，除了部分指标主观和难以测量外，该数据还存在一个问题，数据集中有30%的值是缺失的。下面将首先介绍如何处理数据集中的数据缺失问题，然后再利用Logistic回归和随机梯度上升算法来预测病马的生死。
数据中的缺失值是个非常棘手的问题，有很多文献都致力于解决这个问题。那么，数据缺失究竟带来了什么问题？假设有100个样本和20个特征，这些数据都是机器收集回来的。若机器上；的某个传感器损坏导致一个特征无效时该怎么办？此时是否要扔掉整个数据？这种情况下，另外19个特征怎么办？它们是否还可用？答案是肯定的。因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以必须采用一些方法来解决这个问题。
下面给出了一些可选的做法：
使用可用特征的均值来填补缺失值；
使用特殊值来填补缺失值，如1；
忽略有缺失值的样本；
使用相似样本的均值添补缺失值；
使用另外的机器学习算法预测缺失值。
原始的数据集经过处理，保存为两个文件：horseColicTest.txt和horseColicTraining.txt。参见实例。现在我们有一个“干净”可用的数据集和一个不错的优化算法，下面将把这些部分融合在一起训练出一个分类器，然后利用该分类器来预测病马的生死问题。
有了这些数据集，我们只需要一个Logistic分类器，就可以利用该分类器来预测病马的生死问题了。看代码吧。
# -*- coding:utf-8 -*-
import numpy as np
import time
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['axes.unicode_minus'] = False
"""
函数说明:读取数据
Parameters:
filename - 文件名
Returns:
dataSet_Arr - 数据集
Labels_Arr - 标签
"""
def loadDataSet_old(filename):
fr = open(filename) #打开训练集
dataSet = []
Labels = []
for line in fr.readlines():
currLine = line.strip().split('\t')
lineArr = []
for i in range(len(currLine)-1):
lineArr.append(float(currLine[i]))
dataSet.append(lineArr)
Labels.append(float(currLine[-1]))
# 将X,Y列表转化成矩阵
dataSet_Arr = np.array(dataSet)
Labels_Arr = np.array(Labels)
return dataSet_Arr,Labels_Arr
"""
函数说明:读取数据
Parameters:
filename - 文件名
Returns:
xArr - x数据集
yArr - y数据集
"""
def loadDataSet(filename):
X = []
Y = []
with open(filename, 'rb') as f:
for idx, line in enumerate(f):
line = line.decode('utf-8').strip()
if not line:
continue
eles = line.split()
eles = list(map(float, eles))
if idx == 0:
numFea = len(eles)
#去掉每行的最后一个数据再放入X中
X.append(eles[:-1])
#获取每行的租后一个数据
Y.append([eles[-1]])
# 将X,Y列表转化成矩阵
xArr = np.array(X)
yArr = np.array(Y)
return xArr,yArr
"""
函数说明:定义sigmoid函数
Parameters:
z
Returns:
返回
"""
def sigmoid(z):
return 1.0 / (1.0 + np.exp(-z))
"""
函数说明:定义代价函数
Parameters:
theta, X, Y, theLambda
Returns:
返回
"""
def J(theta, X, Y, theLambda=0):
m, n = X.shape
h = sigmoid(np.dot(X, theta))
J = (-1.0/m)*(np.dot(np.log(h).T,Y)+np.dot(np.log(1-h).T,1-Y)) + (theLambda/(2.0*m))*np.sum(np.square(theta[1:]))
if np.isnan(J[0]):
return np.inf
# 其实J里面只有一个数值，需要取出该数值
return J.flatten()[0]
"""
函数说明:梯度上升求解参数
Parameters:
X, Y, options
Returns:
返回参数
"""
def gradient(X, Y, options):
'''
options.alpha 学习率
options.theLambda 正则参数λ
options.maxLoop 最大迭代轮次
options.epsilon 判断收敛的阈值
options.method
- 'sgd' 随机梯度上升
- 'bgd' 批量梯度上升
'''
m, n = X.shape
# 初始化模型参数，n个特征对应n个参数
theta = np.zeros((n,1))
cost = J(theta, X, Y) # 当前误差
costs = [cost]
thetas = [theta]
# Python 字典dict.get(key, default=None)函数返回指定键的值，如果值不在字典中返回默认值
alpha = options.get('alpha', 0.01)
epsilon = options.get('epsilon', 0.00001)
maxloop = options.get('maxloop', 1000)
theLambda = float(options.get('theLambda', 0)) # 后面有 theLambda/m 的计算，如果这里不转成float，后面这个就全是0
method = options.get('method', 'bgd')
# 定义随机梯度上升
def _sgd(theta):
count = 0
converged = False
while count < maxloop:
if converged :
break
# 随机梯度上升，每一个样本都更新
for i in range(m):
h =sigmoid(np.dot(X[i].reshape((1,n)), theta))
theta = theta + alpha*((1.0/m)*X[i].reshape((n,1))*(Y[i]-h) + (theLambda/m)*np.r_[[[0]], theta[1:]])
thetas.append(theta)
cost = J(theta, X, Y, theLambda)
costs.append(cost)
if abs(costs[-1] - costs[-2]) < epsilon:
converged = True
break
count += 1
return thetas, costs, count
# 定义批量梯度上升
def _bgd(theta):
count = 0
converged = False
while count < maxloop:
if converged :
break
h = sigmoid(np.dot(X, theta))
theta = theta + alpha*((1.0/m)*np.dot(X.T, (Y-h)) + (theLambda/m)*np.r_[[[0]],theta[1:]])
thetas.append(theta)
cost = J(theta, X, Y, theLambda)
costs.append(cost)
if abs(costs[-1] - costs[-2]) < epsilon:
converged = True
break
count += 1
return thetas, costs, count
methods = {'sgd': _sgd, 'bgd': _bgd}
return methods[method](theta)
"""
函数说明:绘制代价函数图像
Parameters:
costs
Returns:
返回参数
"""
def plotCost(costs):
# 绘制误差曲线
plt.figure(figsize=(8,4))
plt.plot(range(len(costs)), costs)
plt.xlabel(u'迭代次数')
plt.ylabel(u'代价J')
"""
函数说明:绘制参数曲线
Parameters:
thetass - 参数
Returns:
返回参数
"""
def plotthetas(thetas):
# 绘制参数theta变化
thetasFig, ax = plt.subplots(len(thetas[0]))
thetas = np.asarray(thetas)
for idx, sp in enumerate(ax):
thetaList = thetas[:, idx]
sp.plot(range(len(thetaList)), thetaList)
sp.set_xlabel('Number of iteration')
sp.set_ylabel(r'$\theta_%d$'%idx)
"""
函数说明:分类函数
Parameters:
inX - 特征向量
weights - 回归系数
Returns:
分类结果
"""
def classifyVector(inX, weights):
prob = sigmoid(sum(inX*weights))
if prob > 0.5:
return 1.0
else:
return 0.0
"""
函数说明:测试数据
Parameters:
Weights - 回归系数
dataSet - 数据集
Labels - 标签
Returns:
predict - 预测结果
errorRate - 错误率
"""
def test(Weights,dataSet,Labels):
errorCount = 0
numTestVec = 0.0
predict = []
for i in range(len(test_Y)):
numTestVec += 1.0
result = classifyVector(dataSet[i], Weights)
predict.append(result)
if int(result)!= int(Labels[i]):
errorCount += 1
errorRate = (float(errorCount)/numTestVec) * 100 #错误率计算
return predict , errorRate
"""
函数说明:绘制预测值与真实值
Parameters:
X_test,Y_test,y_predict
Returns:
无
"""
def plotfit(X_test,Y_test,y_predict):
#画图对预测值和实际值进行比较
p_x = range(len(X_test))
plt.figure(figsize=(18,4),facecolor="w")
plt.ylim(0,6)
plt.plot(p_x,Y_test,"ro",markersize=8,zorder=3,label=u"真实值")
plt.plot(p_x,y_predict,"go",markersize=14,zorder=2,label=u"预测值")
plt.legend(loc="upper left")
plt.xlabel(u"编号",fontsize=12)
plt.ylabel(u"类型",fontsize=12)
plt.title(u"Logistic算法对数据进行分类",fontsize=12)
#plt.savefig("Logistic算法对数据进行分类.png")
plt.show()
if __name__ == '__main__':
## Step 1: load data...
print("Step 1: load data...")
origin_train_X, train_Y = loadDataSet('horseColicTraining.txt')
origin_test_X, test_Y = loadDataSet('horseColicTest.txt')
train_m,train_n = origin_train_X.shape
test_m,test_n = origin_test_X.shape
train_X = np.concatenate((np.ones((train_m,1)), origin_train_X), axis=1)
test_X = np.concatenate((np.ones((test_m,1)), origin_test_X), axis=1)
#print(train_X)
#print(train_Y)
## Step 2: training sgd...
print("Step 2: training sgd...")
# sgd随机梯度上升
options = {
'alpha':0.5,
'epsilon':0.00000001,
'maxloop':100000,
'method':'sgd'
}
# 训练模型
start = time.time()
sgd_thetas, sgd_costs, sgd_iterationCount = gradient(train_X, train_Y, options)
print("sgd_thetas:")
print(sgd_thetas[-1])
end = time.time()
sgd_time = end - start
print("sgd_time:"+str(sgd_time))
## Step 3: testing...
print("Step 3: testing...")
predict_Y , errorRate = test(sgd_thetas[-1].all(),test_X,test_Y)
## Step 4: show the result...
print("Step 4: show the result...")
print('预测值:',predict_Y)
print("测试集错误率为: %.2f%%" % errorRate)
## Step 5: show the plotfit...
print("Step 5: show the plotfit...")
plotfit(test_X,test_Y,predict_Y)
## Step 6: show costs sgd...
print("Step 6: show costs sgd...")
plotCost(sgd_costs)
## Step 7: show thetas sgd...
print("Step 7: show thetas sgd...")
plotthetas(sgd_thetas)
结果如下所示。
【完整代码参考3.Multi_Logistic_Regression_Classify\Multi_Logistic_Regression_Classify_v1\ Multi_Logistic_Regression_Classify_v1.0\Multi_Logistic_Regression_Classify_v1.0.py】
5.4.2多元逻辑回归实战之预测病马的死亡率调用sklearn
和前面一样还是要使用sklearn来实现的，直接看代码吧，注意和前文的例子比较。Python代码如下。
# @Language : Python3.6
"""
# -*- coding:UTF-8 -*-
import numpy as np
import random
import matplotlib.pyplot as plt
import matplotlib
import time
from sklearn.linear_model import LogisticRegression
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['axes.unicode_minus'] = False
"""
函数说明:读取数据
Parameters:
filename - 文件名
Returns:
dataSet_Arr - 数据集
Labels_Arr - 标签
"""
def loadDataSet_old(filename):
fr = open(filename) #打开训练集
dataSet = []
Labels = []
for line in fr.readlines():
currLine = line.strip().split('\t')
lineArr = []
for i in range(len(currLine)-1):
lineArr.append(float(currLine[i]))
dataSet.append(lineArr)
Labels.append(float(currLine[-1]))
# 将X,Y列表转化成矩阵
dataSet_Arr = np.array(dataSet)
Labels_Arr = np.array(Labels)
return dataSet_Arr,Labels_Arr
"""
函数说明:读取数据
Parameters:
filename - 文件名
Returns:
xArr - x数据集
yArr - y数据集
"""
def loadDataSet(filename):
X = []
Y = []
with open(filename, 'rb') as f:
for idx, line in enumerate(f):
line = line.decode('utf-8').strip()
if not line:
continue
eles = line.split()
eles = list(map(float, eles))
if idx == 0:
numFea = len(eles)
#去掉每行的最后一个数据再放入X中
X.append(eles[:-1])
#获取每行的租后一个数据
Y.append([eles[-1]])
# 将X,Y列表转化成矩阵
xArr = np.array(X)
yArr = np.array(Y)
return xArr,yArr
"""
函数说明:绘制预测值与真实值
Parameters:
X_test,Y_test,y_predict
Returns:
无
"""
def plotfit(X_test,Y_test,y_predict):
#画图对预测值和实际值进行比较
p_x = range(len(X_test))
plt.figure(figsize=(18,4),facecolor="w")
plt.ylim(0,6)
plt.plot(p_x,Y_test,"ro",markersize=8,zorder=3,label=u"真实值")
plt.plot(p_x,y_predict,"go",markersize=14,zorder=2,label=u"预测值")
plt.legend(loc="upper left")
plt.xlabel(u"编号",fontsize=12)
plt.ylabel(u"类型",fontsize=12)
plt.title(u"Logistic算法对数据进行分类",fontsize=12)
#plt.savefig("Logistic算法对数据进行分类.png")
plt.show()
if __name__ == '__main__':
## Step 1: load data...
print("Step 1: load data...")
origin_train_X, train_Y = loadDataSet('horseColicTraining.txt')
origin_test_X, test_Y = loadDataSet('horseColicTest.txt')
train_m,train_n = origin_train_X.shape
test_m,test_n = origin_test_X.shape
train_X = np.concatenate((np.ones((train_m,1)), origin_train_X), axis=1)
test_X = np.concatenate((np.ones((test_m,1)), origin_test_X), axis=1)
#print(train_X)
#print(train_Y)
## Step 2: init lr...
print("Step 2: init lr...")
mlr = LogisticRegression(solver = 'sag',max_iter = 5000)
## Step 3: training...
print("Step 3: training...")
start = time.time()
#使用改进的随即上升梯度训练
#输出y做出ravel()转换。不然会有警告
mlr.fit(train_X,train_Y.ravel())
end = time.time()
mlr_time = end - start
print('mlr_time:',mlr_time)
## Step 4: show weights...
print("Step 4: show weights...")
#模型效果获取
r = mlr.score(train_X,train_Y)
print("R值(准确率):",r)
#print("参数:",mlr.coef_)
#print("截距:",mlr.intercept_)
w0 = np.array(mlr.intercept_)
w1 = np.array(mlr.coef_.T[1:22])
#合并为一个矩阵
weights = np.vstack((w0,w1))
#print('weights:')
#print(weights)
## Step 5: predict ...
print("Step 5: predict...")
predict_Y = mlr.predict(test_X)
print('预测值:',predict_Y)
## Step 6: show the plotfit...
print("Step 6: show the plotfit...")
plotfit(test_X,test_Y,predict_Y)
结果如下所示。
【完整代码参考3.Multi_Logistic_Regression_Classify\Multi_Logistic_Regression_Classify-sklearn_v2\Multi_Logistic_Regression_Classify-sklearn_v2.0.py】
本章参考代码
点击进入
