❗写在前面 ❗：
 
… … … … … … … … . … … …
1️⃣ 提供两种语言实现：本文提供两种语言实现，分别为Python和C++，分别依赖的核心库为Numpy和Eigen.
 
2️⃣ 提供关键步的公式推导：写出详尽的公式推导太花时间，所以只提供了关键步骤的公式推导。更详尽的推导，可移步其他博文或相关书籍.
 
3️⃣ 重在逻辑实现：文中的代码主要以实现算法的逻辑为主，在细节优化上面，并没有做更多的考虑，如果发现有不完善的地方，欢迎指出.
 
 
🌞 欢迎关注专栏，相应的内容会持续更新，祝大家变得更强！！
… … … … … … … … . … … …

 

1、原理推导

 

对数几率回归通过Sigmoid函数将线性回归的预测值映射到0-1之间，通过在0-1之间的设定划分阈值，从而实现分类。

Sigmoid函数：

    $y=\frac{1}{1+e^{-z}}$         （1）

导数：
    $f^{{}^{\prime }}(z)=f(z)(1-f(z))$   （2）
 
其中线性回归模型$z=w^{T}x+b$，代入（1）式中可得：

    $y=\frac{1}{1+e^{-(w^T+b)}}$        （3）
 
两边取对数：

    $ln\frac{y}{1-y}=w^{T}x+b$      （4）
 
其中$y$看做是正例的概率，$1-y$看做是反例的概率。则$\frac{y}{1-y}$为几率，对几率求对数就是对数几率。

为了求（4）中模型的$w$和$b$，需要知道对数几率回归的损失函数，然后对损失函数最小化，得到$w$和$b$的估计值。这里省略推导过程。

交叉熵损失函数：
 
    $-lnp(y|x)=-\frac{1}{m}\sum _{i=1}^m(yln\hat{y}+(1-y)ln(1-\hat{y}))$  （5）
 
其中（5）中的$\hat{y}=\frac{1}{1+e^{-(w^{T}x+b)}}$。令$L=lnp(y|x)$，基于$L$分别对$w$和$b$求偏导，有：
 
    $\frac{\partial L}{\partial w}=\frac{1}{m}(\hat{y}-y)\mathbf{x}$     （6）

    $\frac{\partial L}{\partial b}=\frac{1}{m}\sum _{i=1}^m(\hat{y}-y)$    (7)

由（6）（7）式可知，基于$w$和$b$的梯度下降对交叉熵损失最小化，相应的参数就是模型的最有参数。
  

2、算法实现

  

🌈Python实现

  

# 基于梯度下降的对数几率回归
class LogisticRegression:

    def __init__(self, threshold=0.5):
        self.W = None
        self.b = None
        self.threshold = threshold  # 划分阈值

    # sigmoid函数
    def sigmoid(self, z):
        
        return 1 / (1 + np.exp(-z))

    # 计算当前的交叉熵损失和关于W和b的梯度
    def calculateLossGrident(self, X, y):
        # 数据样本总量和特征维度
        m, n = X.shape
        # 线性回归预测结果
        z = np.dot(X, self.W) + self.b
        # 通过sigmoid映射为概率
        p = self.sigmoid(z)
        # 计算当前参数下的损失
        cost = - (1 / m) * np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))
        # 计算当前参数下关于W的梯度
        dW = np.dot(X.T, (p - y)) / m
        # 计算当前参数下关于b的梯度
        db = np.sum(p - y) / m

        cost = np.squeeze(cost)

        return p, cost, dW, db

    # 模型训练
    def fit(self, X, y, alpha=0.01, epochs=1000, tol=0.01):
        """
        alpha: 学习率
        epochs: 迭代次数
        tol: 停止迭代的损失最大值
        """
        # 初始化参数
        m, n = X.shape
        self.W = np.zeros((n, 1))
        self.b = 0.0

        # 记录损失的列表
        costList = []

        # 梯度下降迭代求解
        for i in range(epochs):
            
            p, cost, dW, db = self.calculateLossGrident(X, y)

            # 当前损失符合要求时，停止迭代
            if (cost < tol):
                break
            
            # 更新参数
            self.W -= alpha * dW
            self.b -= alpha * db

            costList.append(cost)
            
            if (i % 500 == 0):
                print(f"epochs {i}: cost = {cost}")
    
    # 预测方法
    def predict(self, X):

        yPred = self.sigmoid(np.dot(X, self.W)) + self.b

        return [1 if y_ > self.threshold else 0 for y_ in yPred]

  

🌟算法验证

  
生成测试的数据集：

import matplotlib.pyplot as plt
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

# 生成数据
X, y = make_classification(n_samples=1000, 
						   n_features=2, 
						   n_informative=2, 
						   n_redundant=0,
                           n_clusters_per_class=1, 
                           random_state=43, 
                           class_sep=0.5)
# 将标签调整为 0 和 1
y = y.reshape(-1, 1)

# 数据集划分
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 可视化生成的数据集
plt.figure(figsize=(10, 6))
plt.scatter(X[:, 0], X[:, 1], c=y.flatten(), cmap='coolwarm', marker='o', s=50, edgecolors='k')
plt.title("Generated Challenging Classification Data")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
plt.colorbar()
plt.grid(True)
plt.show()

  
Numpy实现的对数几率回归算法测试：

# 创建并训练模型
model = LogisticRegression()
model.fit(X_train, y_train)

# 模型预测
y_pred = model.predict(X_test)

# 计算模型精度
accuracy = accuracy_score(y_test, y_pred)
print(f"Accuracy: {accuracy:.2f}")

测试的结果:

epochs 0: cost = 0.6931471805599452
epochs 500: cost = 0.37211354564747723
epochs 1000: cost = 0.270250459992998
epochs 1500: cost = 0.22207604120765112
epochs 2000: cost = 0.19404827550549517
epochs 2500: cost = 0.17568359002576364
epochs 3000: cost = 0.16269779425357095
epochs 3500: cost = 0.153018072888011
epochs 4000: cost = 0.14551888410245217
epochs 4500: cost = 0.13953565971259974
epochs 5000: cost = 0.13465026362003016
epochs 5500: cost = 0.1305861620774047
epochs 6000: cost = 0.12715297797229222
epochs 6500: cost = 0.12421524834328823
epochs 7000: cost = 0.12167388441889265
epochs 7500: cost = 0.11945468451353636
epochs 8000: cost = 0.11750095206458656
epochs 8500: cost = 0.11576860086513328
epochs 9000: cost = 0.11422282016419599
epochs 9500: cost = 0.11283574788999319
Accuracy: 0.95

  
从生成的数据图片的分布来看，数据并不是完全线性可分的，所以从结果来看，手写的对数几率回归的效果还是可以的。
  

⚡C++实现

  

#include <iostream>
#include <Eigen/Dense>
#include <algorithm>

using namespace Eigen;
using namespace std;


// 声明参数的结构体
struct Theta
{
    VectorXd W;
    double b;
};

class LogisticRegression
{
private:
    Theta theta;

private:
    // Sigmoid函数
    double sigmoid(double z)
    {
        return 1.0 / (1.0 + exp(-z));
    }

    // 初始化模型参数
    void initParams(int dim)
    {
        theta.W = VectorXd::Zero(dim);
        theta.b = 0.0;
    }

    // 计算损失和梯度
    struct LossGradient
    {
        VectorXd z;
        VectorXd yHat;
        double loss;
        VectorXd dW;
        double db;
    };

    LossGradient logisticLoss(const MatrixXd& X, const VectorXd& y)
    {
        int m = X.rows();  // 样本数
        int n = X.cols();  // 特征数

        VectorXd z = (X * theta.W).array() + theta.b;  // 计算z
        VectorXd yHat(m);  // 模型输出

        for (int i = 0; i < m; ++i)
        {
            yHat(i) = sigmoid(z(i));  // 计算每个样本的输出
        }

        double loss = (-1.0 / m) * (y.array() * yHat.array().log() + (1.0 - y.array()) * ((1.0 - yHat.array()).log())).sum();  // 计算损失

        VectorXd dW = (1.0 / m) * X.transpose() * (yHat - y);  // 计算权重W的梯度
        double db = (yHat - y).mean();  // 计算偏置b的梯度

        return { z, yHat, loss, dW, db };
    };

public:
    LogisticRegression() {}
    // 模型训练
    void fit(const MatrixXd& X, const VectorXd& y, double alpha = 0.01, int epochs = 1000)
    {
        int numFeatures = X.cols();

        // 初始化参数
        initParams(numFeatures);

        // 使用梯度下降求解最佳参数
        for (int i = 0; i < epochs; i++)
        {
            auto result = logisticLoss(X, y);  // 计算损失和梯度
            theta.W -= alpha * result.dW;  // 更新权重W
            theta.b -= alpha * result.db;  // 更新偏置b

            if (i % 100 == 0)
            {
                cout << "Epoch " << i << ": loss = " << result.loss << endl;
            }
        }
    }

    // 获取训练后的参数
    Theta getParams() const
    {
        return theta;
    }

    // 使用训练好的模型进行预测
    VectorXd predict(const MatrixXd& X) 
    {
        VectorXd z = (X * theta.W).array() + theta.b;
        VectorXd yHat(X.rows());

        for (int i = 0; i < X.rows(); ++i)
        {
            yHat(i) = sigmoid(z(i));  // 计算每个样本的输出
        }

        return yHat;
    }
};

  

3、总结

对数几率回归是一个经典的分类算法，简单粗暴，还是感知机模型、神经网络和支持向量机等模型的基础。在数学原理推导上有些部分没有详尽，感兴趣的可以参考其他书籍或者博文。