逻辑回归算法动手实现
在《逻辑回归算法》一问中，我们介绍了逻辑回归算法的原理，这一篇文章，我们将通过逻辑回归算法的原理，动手用python实现一遍，以加深对逻辑回归算法的理解。同时跟scikit-learn的逻辑回归工具进行效果对比。
一、问题描述
这一次，我们还是通过一个实际的问题案例来引入。我们使用iris数据集，这是一个比较官方的数据来源。这个数据集市关于鸢尾花的数据，该数据有5个特征，分别是花瓣(petal)的长和宽、萼片(sepal)的长和宽、花的种类。其中三类鸢尾花分别是：setosa, versicolor, virginica。数据集可以参考如下链接：
http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data
为了简化问题，我们使用50-150之间的数据，也就是只使用 versicolor, virginica这两类的鸢尾花数据，共100个样本。同时选取萼片(sepal)的长和宽作为输入变量（特征）。也就是，我们要解决的是一个而分类的问题。
下面，我们先看下这两个鸢尾花的数据分布：
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
# 加载Iris数据集
df = pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', header=None)
df.columns = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'category']
print df.head()
输出：
# 鸢尾花的分布图
print '----------------鸢尾花的分布图-------------------'
plt.scatter(df.iloc[50:100, 2], df.iloc[50:100, 3],color='red', marker='o', label='versicolor') # 中间50个样本的散点图
plt.scatter(df.iloc[100:150, 2], df.iloc[100:150, 3],color='blue', marker='x', label='virginica') # 后50个样本的散点图
plt.xlabel('petal length')
plt.ylabel('petal width')
plt.legend(loc=2)
plt.show()
二、用python实现逻辑回归算法
1、sigmod函数
\[
g(z) = \frac{1}{1 + e ^ {-z}}
\]
def sigmod(z):
"""
sigmod函数
:param z: mx1矩阵
:return: mx1矩阵
"""
return 1.0 / (1 + np.exp(-z))
2、模型函数
\[
h(X)=g(X \theta)
\]
def model(theta, X):
"""
预测函数
:param theta: nx1矩阵
:param X: mxn矩阵
:return: mx1矩阵
"""
return sigmod(X.dot(theta))
3、代价函数
def cosst_function(theta, X, y):
"""
代价函数
:param theta: nx1矩阵
:param X: mxn矩阵
:param y: mx1矩阵
:return: 实数
"""
m = y.size
h = sigmod(X.dot(theta))
J = -(1.0/m)*((np.log(h)).T.dot(y) + (np.log(1-h)).T.dot(1-y))
return J[0, 0]
4、梯度下降
def derivative(theta, X, y):
"""
对theta求偏导
:param theta: nx1矩阵
:param X: mxn矩阵
:param y: mx1矩阵
:return: nx1矩阵
"""
m = y.size
h = sigmod(X.dot(theta))
grad = (1.0/m) * X.T.dot(h - y)
return grad
def gradient(grade_theta, theta, alpha):
"""
𝜃更新
:param grade_theta:求导后的𝜃矩阵(nx1)
:param theta: 更新后的𝜃矩阵(nx1)
:param alpha: 学习率
:return: 返回更新后的𝜃矩阵(nx1)
"""
return theta - alpha * grade_theta
5、逻辑回归算法
def logistic_regression(X, y):
"""
逻辑回归算法
:param X: mxn矩阵
:param y: mx1矩阵
:return: nx1矩阵
"""
start = time.clock()
m = X.shape[0] # 样本个数
n = X.shape[1] # 特征个数
y = y.reshape(m, 1)
cost_record = [] # 记录代价函数的值
alpha = 0.1 # 学习率
maxiters = 1000 # 最大迭代次数
theta = np.zeros(n).reshape(n, 1) # 设置权重参数的初始值
cost_val = cosst_function(theta, X, y)
cost_record.append(cost_val)
iters = 0
while True:
if iters >= maxiters:
break
grad = derivative(theta, X, y)
# 权重参数更新
new_theta = gradient(grad, theta, alpha)
theta = new_theta
cost_update = cosst_function(new_theta, X, y)
cost_val = cost_update
cost_record.append(cost_val)
iters += 1
end = time.clock()
print("cost time: %f s" % (end - start))
return cost_record, theta, iters
6、预测函数
def predict(theta, X, threshold=0.5):
"""
预测函数
:param theta: 已求解的theta函数
:param X: 输入变量
:param threshold: 阈值
:return: 输出变量（预测）
"""
p = model(theta, X) >= threshold
return p.astype('int')
最后，我们通过给定的数据集，对逻辑回归模型进行训练（80%训练集，20%测试集）。
"""
通过判定花萼长度，花萼宽度，花瓣长度，花瓣宽度的尺寸大小来识别鸢尾花的类别
前四列分别为：花萼长度（sepal），花萼宽度，花瓣长度（petal），花瓣宽度
第5列为鸢尾花的类别（包括Setosa，Versicolour，Virginica三类）
"""
# 加载Iris数据集
df = pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', header=None)
df.columns = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'category']
print df.head()
# 鸢尾花的分布图
print '----------------鸢尾花的分布图-------------------'
plt.scatter(df.iloc[50:100, 2], df.iloc[50:100, 3],color='red', marker='o', label='versicolor') # 中间50个样本的散点图
plt.scatter(df.iloc[100:150, 2], df.iloc[100:150, 3],color='blue', marker='x', label='virginica') # 后50个样本的散点图
plt.xlabel('petal length')
plt.ylabel('petal width')
plt.legend(loc=2)
plt.show()
print '----------------数据预处理-------------------'
df = df.iloc[50:150, :]
df.loc[df['category'] == 'Iris-versicolor', 'category'] = 0
df.loc[df['category'] == 'Iris-virginica', 'category'] = 1
df['category'] = df['category'].astype('int')
print df.head()
print df.shape
print df.dtypes
# 取其中两个特征，即：花瓣长度，花瓣宽度，即第2列和第3列
X = df[['petal_length', 'petal_width']].values
y = df[['category']].values
X_std = np.c_[np.ones((df.shape[0], 1)), min_max(X)]
X_train, X_test, y_train, y_test = train_test_split(X_std, y, test_size=0.2, random_state=0)
print X_train.shape
print X_test.shape
cost_record, theta, iters = logistic_regression(X_train, y_train)
print '----------------训练结果-------------------'
print "代价函数数组：", cost_record
print "theta参数：", theta
print "迭代次数：", iters
print '-----------------测试-----------------'
train_right_counts = 0
for i in range(X_train.shape[0]):
original_val = y_train[i, :]
predict_val = predict(theta, X_train[i, :])
# print "Original:", original_val, " Predict:", predict_val, ("o" if original_val == predict_val else "x")
if original_val == predict_val:
train_right_counts += 1
print "训练集准确率：", ((train_right_counts * 1.0) / X_train.shape[0])
test_right_counts = 0
for i in range(X_test.shape[0]):
original_val = y_test[i, :]
predict_val = predict(theta, X_test[i, :])
# print "Original:", original_val, " Predict:", predict_val, ("o" if original_val == predict_val else "x")
if original_val == predict_val:
test_right_counts += 1
print "测试集准确率：", ((test_right_counts * 1.0) / X_test.shape[0])
# 代价函数值收敛趋势图
plt.plot(np.arange(len(cost_record)-1), cost_record[1:], 'r')
plt.show()
训练的模型，给出的测试效果还是不错的，测试集准确度达到了95%。
代价函数收敛趋势图：
三、scikit-learn的逻辑回归算法工具
# -*- coding: UTF-8 -*-
import pandas as pd
from sklearn.model_selection import train_test_split
"""
通过判定花萼长度，花萼宽度，花瓣长度，花瓣宽度的尺寸大小来识别鸢尾花的类别
前四列分别为：花萼长度（sepal），花萼宽度，花瓣长度（petal），花瓣宽度
第5列为鸢尾花的类别（包括Setosa，Versicolour，Virginica三类）
"""
from sklearn import datasets
iris = datasets.load_iris()
# 加载Iris数据集作为DataFrame对象
df = pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', header=None)
df.columns = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'category']
df = df.iloc[50:150, :]
df.loc[df['category'] == 'Iris-versicolor', 'category'] = 0
df.loc[df['category'] == 'Iris-virginica', 'category'] = 1
df['category'] = df['category'].astype('int')
# 取其中两个特征，即：花瓣长度，花瓣宽度，即第2列和第3列
X = df[['petal_length', 'petal_width']].values
y = df[['category']].values
# plt.scatter(X[:50, 0], X[:50, 1],color='red', marker='o', label='setosa') # 前50个样本的散点图
# plt.scatter(X[50:100, 0], X[50:100, 1],color='blue', marker='x', label='versicolor') # 中间50个样本的散点图
# plt.scatter(X[100:, 0], X[100:, 1],color='green', marker='+', label='Virginica') # 后50个样本的散点图
# plt.xlabel('petal length')
# plt.ylabel('sepal length')
# plt.legend(loc=2) # 把说明放在左上角，具体请参考官方文档
# plt.show()
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
#对特征进行标准化处理（特征缩放）
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
sc.fit(X_train)
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(C=1000.0, random_state=0)
lr.fit(X_train_std, y_train.ravel())
print lr.coef_
print '-----------------测试-----------------'
train_right_counts = 0
for i in range(X_train_std.shape[0]):
original_val = y_train[i, :]
predict_val = lr.predict(X_train_std[i, :].reshape(1, -1))
# print "Original:", original_val, " Predict:", predict_val, ("o" if original_val == predict_val else "x")
if original_val == predict_val:
train_right_counts += 1
print "训练集准确率：", ((train_right_counts * 1.0) / X_train_std.shape[0])
test_right_counts = 0
for i in range(X_test_std.shape[0]):
original_val = y_test[i, :]
predict_val = lr.predict(X_test_std[i, :].reshape(1, -1))
# print "Original:", original_val, " Predict:", predict_val, ("o" if original_val == predict_val else "x")
if original_val == predict_val:
test_right_counts += 1
print "测试集准确率：", ((test_right_counts * 1.0) / X_test_std.shape[0])
# # 化边界函数
# X_combined_std = np.vstack((X_train_std, X_test_std))
# y_combined = np.hstack((y_train, y_test))
#
# plot_decision_regions(X_combined_std, y_combined, classifier=lr, test_idx=range(105,150))
# plt.xlabel('petal length [standardized]')
# plt.ylabel('petal width [standardized]')
# plt.legend(loc='upper left')
# plt.show()
测试结果如下，测试集的准确率跟我们手动实现的逻辑回归算法测试结果是一样的，只是训练集的准确率略高一些。
微信扫一扫，打赏作者吧～相关文章...
