图解Python逻辑回归：调包实现+徒手实现

1、sklearn实现

from sklearn.linear_model import LogisticRegression
from sklearn.datasets.samples_generator import make_blobs
import numpy as np, matplotlib.pyplot as mp
from mpl_toolkits import mplot3d

# 创建数据
X, y = make_blobs(centers=[[2, 2, 2], [0, 0, 0]], random_state=3)

# 训练模型
model = LogisticRegression()
model.fit(X, y)

# 系数（coefficient）和截距
k = model.coef_[0]
b = model.intercept_[0]

# 可视化
fig = mp.figure()
ax = mplot3d.Axes3D(fig)  # 获取三维坐标轴
# 散点图
ax.scatter(X[:, 0], X[:, 1], X[:, 2], s=99, c=y)
# 决策边界
x1 = np.arange(X[:, 0].min()-1, X[:, 0].max()+1, 2)
x2 = np.arange(X[:, 1].min()-1, X[:, 1].max()+1, 2)
x1, x2 = np.meshgrid(x1, x2)
x3 = (b + k[0] * x1 + k[1] * x2) / - k[2]  # 边界函数
ax.plot_surface(x1, x2, x3, alpha=0.2, color='g')
mp.show()

2、徒手写一个！

import numpy as np, matplotlib.pyplot as mp
from sklearn.datasets import make_blobs

"""创建随机样本"""
X, Y = make_blobs(centers=2, cluster_std=5)

"""数据处理"""
X = np.insert(X, 0, 1, axis=1)  # 增加一列，用于矩阵相乘
Y = Y.reshape(-1, 1)

"""sigmoid函数"""
sigmoid = lambda x: 1 / (1 + np.exp(-x))

"""梯度上升"""
d = X.shape[1]  # 维数
theta = np.mat([[1]] * d)  # 初始化回归系数
for i in range(5999, 8999):
    alpha = 1 / i  # 步长（先大后小）
    h = sigmoid(X * theta)
    theta = theta + alpha * X.T * (Y - h)  # 最终梯度上升迭代公式

"""数据可视化"""
x1, x2 = X[:, 1], X[:, 2]
mp.axis([x1.min() + 1, x1.max() + 1, x2.min() + 1, x2.max() + 1])
mp.scatter(x1, x2, c=Y.reshape(-1))  # 原始样本点
x = np.array([x1.min(), x1.max()])
y = (-theta[0, 0] - theta[1, 0] * x) / theta[2, 0]  # 决策边界
mp.plot(x, y)
mp.show()

2.1、创建样本

X, y = make_blobs(centers=2, cluster_std=6)
X = np.insert(X, 0, 1, axis=1)  # 增加一列，用于矩阵相乘
y = y.reshape(-1, 1)  # 转换成n*1维的矩阵

$$X={\left(\begin{array}{ccc}1& x_{12}& x_{13}\\ 1& x_{22}& x_{33}\\ ⋮& ⋮& ⋮\\ 1& x_{n2}& x_{n3}\end{array}\right)}_{n行3列}y={\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right]}_{n行1列}$$

2.2、Sigmoid函数

亦称S型生长曲线。在信息科学中，常被用作神经网络的阈值函数，将变量映射到(0,1)

import numpy as np, matplotlib.pyplot as mp
sigmoid = lambda x: 1 / (1 + np.exp(-x))  # sigmoid函数
x = np.linspace(-8, 8, 65)
mp.plot(x, sigmoid(x), label=r'$S(x)=\frac{1}{1+e^{-x}}$')
mp.legend()  # 图例
mp.show()

2.3、梯度上升

d = X.shape[1]  # dimensionality
theta = np.mat([[1]] * d)  # 初始化回归系数
for i in range(9, 9999):
    alpha = 1 / i  # 学习率（前期大步迭代，后期小步迭代）
    h = sigmoid(X * theta)
    theta = theta + alpha * X.transpose() * (y - h)  # 最终梯度上升迭代公式

3、梯度下降原理

以 $f(x)=x^2-2x+1$ 为例，求函数的极值点

1、求导数： $f^{\prime }(x)=2x-2$
2、计算机使用迭代的方法，一步一步逼近极值点
3、迭代公式： $x_{i+1}=x_i-\alpha \frac{\partial f(x_i)}{x_i}=x_i-\alpha (2x_i-2)$

import matplotlib.pyplot as mp, numpy as np

# Δx
delta = 1e-9
# 步长，即学习速率
alpha = 0.1
# 允许误差范围
precision = delta

def derivative(fx, x0):
    """求f'(x0)"""
    return (fx(x0 + delta) - fx(x0)) / delta

# 原函数
f = lambda x: x ** 2 - 2 * x + 1
x = np.linspace(0, 2, 9999)
mp.plot(x, f(x), c='black')  # 绘制函数曲线

# 梯度下降求极值点
extreme_point = 0  # 初始值
error = np.inf  # 初始化误差
while abs(error) >= precision:  # 误差小于精度时退出迭代
    mp.scatter(extreme_point, f(extreme_point), c='b', alpha=.3)  # 绘制散点
    error = alpha * derivative(f, extreme_point)  # 步伐
    extreme_point -= error  # 梯度下降
    mp.pause(.1)
    print(extreme_point, error)

4、附录

En	Cn
logistic regression	逻辑回归
alpha	希腊字母的第1个字母：$\alpha$
theta	希腊字母的第8个字母： $\theta$
presision	精度
gradient ascent	梯度上升
linspace	abbr. 线性等分向量（linear space）
abbreviation	缩写词
shape	n. 形状；身材；v 形成
extreme point	极值点
derivative	导数
partial derivatives	偏导

动态图

import numpy as np, matplotlib.pyplot as mp
from sklearn.datasets import make_blobs

"""创建随机样本"""
X, Y = make_blobs(centers=2, cluster_std=5)

"""数据处理"""
X = np.insert(X, 0, 1, axis=1)  # 增加一列，用于矩阵相乘
Y = Y.reshape(-1, 1)
x1, x2 = X[:, 1], X[:, 2]
fig, ax = mp.subplots()  # 创建绘图对象

"""sigmoid函数"""
sigmoid = lambda x: 1 / (1 + np.exp(-x))

"""梯度上升"""
d = X.shape[1]  # 维数
theta = np.mat([[1]] * d)  # 初始化回归系数
for i in range(2000, 2200):
    alpha = 1 / i  # 步长（先大后小）
    h = sigmoid(X * theta)
    theta = theta + alpha * X.T * (Y - h)  # 最终梯度上升迭代公式
    """数据可视化"""
    ax.cla()  # 清除
    ax.axis([x1.min() + 1, x1.max() + 1, x2.min() + 1, x2.max() + 1])
    ax.scatter(x1, x2, c=Y.reshape(-1))  # 原始样本点
    x = np.array([x1.min(), x1.max()])
    y = (-theta[0, 0] - theta[1, 0] * x) / theta[2, 0]  # 决策边界
    ax.plot(x, y)
    mp.pause(.1)