二元逻辑回归实现鸢尾花数据分类（python）
2019-03-09 本文已影响0人
马尔代夫Maldives
说明：
- 本文利用python实现二元逻辑回归，没有加正则项。挑选iris数据前100个样本作为训练集，他们分属于两个类别，样本特征选择第1列（花萼长度x1）和第2列（花萼宽度x2）。
- 程序以函数形式实现，附录中给出封装式程序实现。代码均以矩阵计算的思路编写。
- 四部分程序是连续的，拷贝在一起就能完整运行。
- 没考虑过拟合、欠拟合等问题。
1.查看鸢尾花数据集
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt
import numpy as np
#%matplotlib
iris = load_iris() #导入鸢尾花数据
data = iris.data #numpy.ndarray
target = iris.target #numpy.ndarray
#print(data.shape) #查看维度（150,4），即150个样本，每个样本4个属性
#print(target.shape) #（150,）每个样本的类别0,1,2值，三个类型，前100个样本属于0和1类，后50个样本属于2类
#print(data) #查看完整数据
#print(target) #查看类型数据
#根据不同类型0,1,2和第0列（x1），第1列（x2）绘图。
plt.figure(1)
index_0 = np.where(target==0) #找出0所在位置
plt.scatter(data[index_0,0],data[index_0,1],marker='x',color = 'b',label = '0',s = 15)
index_1 =np.where(target==1) #找出1所在位置
plt.scatter(data[index_1,0],data[index_1,1],marker='o',color = 'r',label = '1',s = 15)
index_2 =np.where(target==2) #找出2所在位置
plt.scatter(data[index_2,0],data[index_2,1],marker='s',color = 'k',label = '2',s = 15)
plt.xlabel('sepal length')
plt.ylabel('sepal width')
plt.legend(loc = 'upper left')
上述程序的结果如下：
鸢尾花.png
2.二元逻辑回归程序实现
#定义二元逻辑回归需要的各个函数
def h_sigmod(W,X): #输出1×m
g = np.dot(W,X.T) #1×m 计算h(X)
return sigmod(g)
def sigmod(X):
return 1/(1+np.exp(-X))
def train(X,y,learn_rate=0.001):
#利用已知样本训练模型，得到系数
W = None #存放系数
Loss = [] #存放每次迭代的损失函数值
num_train,num_feature = X.shape
W = 0.001*np.random.randn(num_feature,1).reshape(1,-1) #初始化系数,并转化为1个行向量
#print('W={}'.format(W))
for i in range(5000):
h = h_sigmod(W,X_train) #1×m
#注意loss是一个标量，但是一个ndarray型数据，而非纯粹的int
loss = -(np.log(h).dot(y) + np.log(1-h).dot(1-y))
# print(loss.shape)
loss = loss/num_train
dW = (y.T-h).dot(X_train) #计算(y-h(X))*X，实际上就是J(W)的一阶导数
# print('dW={}'.format(dW))
W += learn_rate*dW #更新系数W=W+a(y-h(X))*X
Loss.append(loss[0])
# print ('i={},Loss={}'.format(i,Loss[i]))
return W,Loss
#利用训练得到的W系数对结果进行预测
def predict(W,X_test):
h = h_sigmod(W,X_test) #利用训练得到的系数W计算h(X)结果就是预测值
y_pred = np.where(h>=0.5,1,0)
return y_pred
3.利用鸢尾花数据对模型进行训练得到系数
##利用鸢尾花部分数据训练模型并查看结果##############################################################
#从原始数据中抽取部分数据作为训练数据
# 前100个样本分别属于0和1类型，本例中用逻辑回归做二分类，不考虑花型数据为2类的数据。
X = data[0:100,[0,1]] #抽出第1和2个特征（x1和x2），分别是花萼的长和宽
y = target[0:100].reshape(-1,1) #每个样本对应的类型，0或1。转化成一个列向量。
# print(X[:5])
# print(y[-5:])
##增加一列“1”，表示x0列，这列只是为了与系数w的第一列w0相匹配（该列不需要与实际特征x1或x2相乘）
#该步很重要
one = np.ones((X.shape[0],1))
X_train = np.hstack((one,X))
W,Loss = train(X_train,y) #利用已知样本训练模型，得到系数和损失函数结果，（）
print('预测的系数：[w0,w1,w2]={}'.format(W)) #输出训练得到的系数
plt.figure(2)
#print('W={}'.format(W))
plt.plot(Loss)
plt.xlabel('number of iteration')
plt.ylabel('Loss')
plt.show()
plt.figure(3)
label = np.array(y)
index_0 = np.where(label==0)
plt.scatter(X[index_0,0],X[index_0,1],marker='x',color = 'b',label = '0',s = 15)
index_1 =np.where(label==1)
plt.scatter(X[index_1,0],X[index_1,1],marker='o',color = 'r',label = '1',s = 15)
#show the decision boundary
x1 = np.arange(4,7.5,0.5)
x2 = (- W[0][0] - W[0][1]*x1) / W[0][2]
plt.plot(x1,x2,color = 'black')
plt.xlabel('X1')
plt.ylabel('X2')
plt.legend(loc = 'upper left')
##利用鸢尾花部分数据训练模型并查看结果
输出结果：
预测的系数：[w0,w1,w2]=[[-2.01401491 5.05730268 -8.19245827]]
损失函数.png
分类边界.png
4.利用前面得到的模型数据预测对其他数据进行分类
#预测测试一
#取X_train的前20个样本，让分类器对其进行预测，如果结果都是0，那么说明分类很准确。
X_test= X_train[:20,[0,1,2]]
result = predict(W,X_test) #根据预测得到的系数W对新的样本进行分类
print("测试一结果:{}".format(result)) #看输出是否为0
#预测测试二
#暂时抽取第2类数据的x1和x2特征来作为测试数据。
#从第一幅途中可以看出，仅从x1和x2两个特征看，第2类鸢尾花是属于第一1类的，即输出应该都是1
Test_data = data[100:,[0,1]] #暂时抽取第2类数据的x1和x2特征来作为测试数据
one = np.ones((Test_data.shape[0],1))
X_test = np.hstack((one,Test_data))
result = predict(W,X_test) #根据预测得到的系数W对新的样本进行分类
print("测试二结果:{}".format(result)) #看输出是否为1
输出分类结果：
测试一结果:[[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]] #全正确
测试二结果:[[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1]] #全正确
附录：上述过程的封装式实现如下
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 7 19:50:46 2019
@author: sean
"""
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt
import numpy as np
# 多cell输出
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
#%matplotlib
iris = load_iris()
data = iris.data #numpy.ndarray
target = iris.target #numpy.ndarray
#print(data.shape) #查看维度（150,4），即150个样本，每个样本4个属性
#print(target.shape) #（150,）每个样本的类别0,1,2值，三个类型，前100个样本属于0和1类，后50个样本属于2类
#print(data) #查看完整数据
#print(target) #查看类型数据
#根据不同类型0,1,2和第0列（x1），第1列（x2）绘图。
plt.figure(1)
index_0 = np.where(target==0) #找出0所在位置
plt.scatter(data[index_0,0],data[index_0,1],marker='x',color = 'b',label = '0',s = 15)
index_1 =np.where(target==1) #找出1所在位置
plt.scatter(data[index_1,0],data[index_1,1],marker='o',color = 'r',label = '1',s = 15)
index_2 =np.where(target==2) #找出2所在位置
plt.scatter(data[index_2,0],data[index_2,1],marker='s',color = 'k',label = '2',s = 15)
plt.xlabel('sepal length (x1)')
plt.ylabel('sepal width (x2)')
plt.legend(loc = 'upper left')
#准备数据############################################################################################
# 前100个样本分别属于0和1类型，本例中用逻辑回归做二分类，不考虑花型数据为2类的数据。
X = data[0:100,[0,1]] #抽出第1和2个特征（x1和x2），分别是花萼的长和宽
y = target[0:100].reshape(-1,1) #每个样本对应的类型，0或1。转化成一个列向量。
# print(X[:5])
# print(y[-5:])
##增加一列“1”，表示x0列，这列只是为了与系数w的第一列w0相匹配（该列不需要与实际特征x1或x2相乘）
one = np.ones((X.shape[0],1))
X_train = np.hstack((one,X))
#构造一个二元逻辑回归分类器########################################################################
class logistic(object):
def __init__(self):
self.W = None
self.Loss = []
def train(self,X,y,learn_rate = 0.001,num_iters = 5000):
num_train,num_feature = X.shape
#init the weight
self.W = 0.001*np.random.randn(num_feature,1).reshape((1,-1))#初始化系数,并转化为1个行向量
for i in range(num_iters):
loss,dW = self.compute_loss(X,y)
# print('i={},loss={}'.format(i,loss))
self.W += learn_rate*dW #更新系数W=W+a(y-h(X))*X
#print(self.W.shape)
self.Loss.append(loss[0][0])
def compute_loss(self,X,y):
num_train = X.shape[0]
h = self.h_sigmod(X) #1×m
#注意loss是一个标量，但是一个ndarray型数据，而非纯粹的int
loss = -(np.log(h).dot(y) + np.log(1-h).dot(1-y))
#print(loss.shape)
loss = loss/num_train
dW = (y.T-h).dot(X_train) #计算(y-h(X))*X，实际上就是J(W)的一阶导数
# print('dW={}'.format(dW))
return loss,dW
def h_sigmod(self,X): #输出1×m
g = np.dot(self.W,X.T) #1×m，计算h(X)
return self.sigmod(g)
def sigmod(self,X):
return 1/(1+np.exp(-X))
def predict(self,X): #利用训练得到的系数W计算h(X)结果就是预测值
h = self.h_sigmod(X)
y_pred = np.where(h>=0.5,1,0)
return y_pred
#产生一个二元逻辑回归分类器对象######################################################################
my_logistic = logistic() #构造一个二元逻辑回归分类器
my_logistic.train(X_train,y) #利用数据X_train 和y对模型进行训练，确定W和Loss值
Loss = my_logistic.Loss
W = my_logistic.W
#print( 'W={}'.format(my_logistic.W) )
#print( 'Loss={}'.format(my_logistic.Loss) )
#查看损失函数
plt.figure(2)
plt.plot(Loss)
plt.xlabel('Iteration number')
plt.ylabel('Loss value')
plt.show()
##绘制根据二元逻辑回归得到的分界线
#将鸢尾花的前100个样本，根据x1和x2特征绘图
plt.figure(3)
label = np.array(y)
index_0 = np.where(label==0)
plt.scatter(X[index_0,0],X[index_0,1],marker='x',color = 'b',label = '0',s = 15)
index_1 =np.where(label==1)
plt.scatter(X[index_1,0],X[index_1,1],marker='o',color = 'r',label = '1',s = 15)
#根据训练得到的系数画分界线
x1 = np.arange(4,7.5,0.5)
x2 = (- W[0][0] - W[0][1]*x1) / W[0][2]
plt.plot(x1,x2,color = 'black')
plt.xlabel('X1')
plt.ylabel('X2')
plt.legend(loc = 'upper left')
#利用训练得到的参数对新的样本进行分类
#预测测试一
#取X_train的前20个样本，让分类器对其进行预测，如果结果都是0，那么说明分类很准确。
X_test= X_train[:20,[0,1,2]]
result = my_logistic.predict(X_test) #对测试数据进行预测
print("测试一结果:{}".format(result)) #看结果有多少是0
#预测测试二
#暂时抽取第2类数据的x1和x2特征来作为测试数据。
#从第一幅途中可以看出，仅从x1和x2两个特征看，第2类鸢尾花是属于第一1类的，即输出应该都是1
Test_data = data[100:,[0,1]]
one = np.ones((Test_data.shape[0],1))
X_test = np.hstack((one,Test_data))
result = my_logistic.predict(X_test)
print("测试二结果:{}".format(result)) #看输出是否为1
