线性回归与逻辑回归及其实现

线性回归与逻辑回归及其实现

回归与分类

• 预测值定性分析，即离散变量预测时，称之为分类；
• 预测值定量分析，即连续变量预测时，称之为回归。

如预测一张图片是猫还是狗，是分类问题；预测明年的房价，是回归问题。

线性回归

线性回归模型

• 一元线性回归模型
  $$h_{\theta }(x)={\theta }_0{x}_0+{\theta }_1{x}_1,(x_0=1)$$
  即简单的 $y=kx+b$

• 多元线性回归模型
  $$h_{\theta }(x)={\theta }_0{x}_0+{\theta }_1{x}_1+\cdots +{\theta }_n{x}_n=\sum _{i=0}^n{\theta }_i{x}_i={\theta }^{T}X,(x_0=1)$$

这里的 $x_0=1$ 用来表示偏置 $b$。

一元线性回归实例——房价预测

在房价问题中我们以房子面积、离地铁站的距离等属性作为自变量，房价是因变量

• $x_0$：房子面积

• $x_1$：离最近地铁站距离

• $x_2$：绿化规模

线性回归可以理解为将这些特征线性组合起来，即可得到房价。而我们的目标就是求出一组合理的权重参数 $\theta$ ，能够较为准确的预测出真实的房价 $h_{\theta }(x)$。

这里的一元自变量我们使用年份，$y=kx+b$ ，这里 $x$ 是年份，$y$​ 是该年的房价。我们采用一元线性回归模型来拟合上海市房价，通过已知历史数据拟合出k和b的值（模型的参数）然后我们就可以通过 $k$、$b$ 来估算2019年、2020年房价。

定义问题

如果我们有 $n$ 组样本：$(x_1,y_1),(x_2,y_2),\dots ,(x_3,y_3)$ ，我们希望设计一个损失函数 $L(x)$ 来拟合这一组样本，当 $\frac{1}{n}(L(x_1)+L(x_2)+\cdots +L(x_n)$ 最小时，我们得到了想要的模型，这个最小化 $L$ 的思路，也称为经验风险最小化（Empirical Risk Minimization**,** ERM）。

在这里，我们使用 MSE/L2 损失：
$$L=\frac{1}{2}(y-h_{\theta }(x){)}^2$$
这与我们在之前 梯度下降法和牛顿法计算开根号中介绍的梯度下降法的损失函数基本一致。

Why MSE ?

我们为什么要选用 MSE 呢？

实际上我们知道，按照解方程的思想，对于只有两个样本的时候，$y=kx+b$ 有唯一解 $k,b$，但是当有大于两个样本时，$k$ 和 $b$ 只有近似解，也就是说，不可能满足所有样本方程，即不可能通过某一组 $k$ 和 $b$ 完全拟合所有数据，总会带有一些差异，我们要最小化这个差异。

真实值和预测值的关系可以表示为：
$$y=h_{\theta }(x)+ϵ$$
这里的 $y$ 是真实值，$h_{\theta }(x)$ 是我们的预测值，而 $ϵ$ 则是误差值。而我们总希望预测值足够的接近真实值，也意味着希望不可控的误差值足够接近0，因此需要对误差项进行分析。

MSE的推导

假设误差是服从均值为0方差为 ${\sigma }^2$ 的正态分布。

• 为什么这里均值为0：因为均值是预测值

• 为什么是正态分布，因为现实生活中绝大部分事物都很容易的服从正态分布。这是统计出来的规律做的假设，使得对误差的假设能够尽可能满足更多的场景。若数据不满足正态分布，则假设失败，模型容易出现非预期结果。

使用最大似然估计法，最大化误差出现的概率。这里是最大化概率密度函数：
$$argma{x}_{\theta }L(\theta )=\prod _{i=1}^n\frac{1}{\sqrt{2\pi }\sigma }\cdot e^{-\frac{{ϵ}_i^2}{2{\sigma }^2}}$$
其中 $ϵ=y-h_{\theta }(x)$ 。当 $ϵ=0$ ，即误差值为 0 时，整式取到最大值。

给定n个样本，每个样本出现误差的可能性乘在一起，表示了整体出现误差的可能性（联合概率），所以是连乘。当所有误差都接近均值时，其发生的概率最大，最大似然其实就是表示了大家都别出现误差好了，用数学上的概率密度函数表示这件事。

然后就是对数似然，属于是本科阶段概率论与数理统计的基操：
$$\begin{array}{rl}lnL(\theta )& =ln\prod _{i=1}^n\frac{1}{\sqrt{2\pi }\sigma }\cdot e^{-\frac{{ϵ}_i^2}{2{\sigma }^2}}\\ & =\sum _{i=1}^{n}ln\frac{1}{\sqrt{2\pi }\sigma }\cdot e^{-\frac{{ϵ}_i^2}{2{\sigma }^2}}\\ & =\sum _{i=1}^n(ln\frac{1}{\sqrt{2\pi }\sigma }+ln{e}^{-\frac{{ϵ}_i^2}{2{\sigma }^2}})\\ & =nln\frac{1}{\sqrt{2\pi }\sigma }-\frac{1}{2{\sigma }^2}\sum _{i=1}^n{ϵ}_i^2\end{array}$$
前面一项是常数，${\sigma }^2$ 是数据集的方差，也是常数，因此这里的最大似然， 我们只需要最大化下面这个 $J(\theta )$：
$$J(\theta )=-\frac{1}{2}\sum _{i=1}^n{ϵ}_i^2=-\frac{1}{2}(y-h_{\theta }(x){)}^2$$
即最小化 $\frac{1}{2}(y-h_{\theta }(x){)}^2$，就是我们的 MSE 损失函数了。

MSE可以用“对数据误差采用正态分布假设”推导而得来，这里也体现了正太分布在模型预测中的重要性，如果不服从假设，容易出现训练异常、不收敛、效率低等各种问题。

欠拟合与过拟合

现实中，我们的要拟合的数据很可能不止是一个简单的线性分布。当我们试图用简单的线性模型去拟合非线性的数据，无疑误差是会非常大的，这就是欠拟合现象。为了增强模型的表达能力，我们可以引入更多的（非线性）参数。比如下式中，我们共有三个参数 $k_{identity},k_{sin},k_{cos}$ ：
$$y=k_{identity}x+k_{sin}sin(x)+k_{cos}cos(x)+b$$
但是，如果我们加的参数过多，反而会导致模型表达能力过强，将数据的偏差也学习了进来，这就会使得模型在面对新的数据时预测不准，这种情况称为过拟合。对于过拟合，我们有很多种缓解的方式，如数据增强、正则化（权重衰减）、神经网络中的dropout、BN等。这里我们介绍一下正则化的方法。

正则化

正则化的思路是限制模型的复杂度，最好能够得到较为稀疏的模型（多个参数 $k$ 为 0）。具体来说，正则化会在模型原本的损失函数中加入一项正则项，对于模型复杂度过高的情况，加以惩罚。在我们之前的例子中，可以是：
$$L=\frac{1}{2}(k_{identity}x+k_{sin}sin(x)+k_{cos}cos(x)+b{)}^2+\lambda (k_{identity}^2+k_{sin}^2+k_{cos}^2+b^2)$$
这里后面一项 $\lambda (k_{identity}^2+k_{sin}^2+k_{cos}^2+b^2)$ 就是正则项，而 $\lambda$ 称为正则项系数。

我们之前提到正则化相当于对模型的参数做了约束，要求其越简单越好，最小化模型参数，这也称为：结构风险最小化(Structural Risk Minimization，SRM)。

正则与范数

正则化也有许多种不同的形式，主要是区别于正则项所选用的范数，比如上面的正则项中，就是选用的 $L_2$ 范数，故也称为 $L_2$ 正则化，岭回归等。

类似的，还有 $L_1$ 正则项 $\lambda (|k_{identity}|+|k_{sin}|+|k_{cos}|+|b|)$ ，这也称为 lasso 回归。

更一般地，对于 $P$ 型范数 $||\theta |{|}_p=\sqrt[p]{{\sum }_{i=1}^n|{\theta }_i{|}^p}$ ，有带正则化的损失函数：$L=L_{\theta }(x)+\lambda ||\theta |{|}_p^p$ 。

代码

#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <cmath>
#include <tuple>
#include <iomanip>

using namespace std;

struct Item {
    float year;
    float price;
};

vector<Item> load_data(const string& file) {

    vector<Item> output;
    fstream ifile(file, ios::binary | ios::in);

    string line;
    getline(ifile, line);

    while(getline(ifile, line)){
        int p = line.find(",");
        Item item;
        item.year  = atof(line.substr(0, p).c_str());
        item.price = atof(line.substr(p + 1).c_str());

        output.emplace_back(item);
    }
    return output;
}

tuple<Item, Item> compute_mean_std(const vector<Item>& items) {

    Item mean{0, 0}, stdval{0, 0};

    for(auto& item : items){
        mean.year  += item.year;
        mean.price += item.price;
    }
    mean.year  /= items.size();
    mean.price /= items.size();

    for(auto& item : items){
        stdval.year  += std::pow(item.year - mean.year, 2.0f);
        stdval.price += std::pow(item.price - mean.price, 2.0f);;
    }
    stdval.year  = std::sqrt(stdval.year / items.size());
    stdval.price = std::sqrt(stdval.price / items.size());
    return make_tuple(mean, stdval);
}

int main() {

    auto datas = load_data("shanghai.csv");

    Item mean, stdval;
    tie(mean, stdval) = compute_mean_std(datas);

    for(auto& item : datas){
        item.year  = (item.year - mean.year) / stdval.year;
        item.price = (item.price - mean.price) / stdval.price;
    }

    float k_identity = 0.1;
    float k_sin      = 0.1;
    float k_cos      = 0.1;
    float lambda     = 1e-5;
    float b = 0;
    float lr = 0.01;

    for(int iter = 0; iter < 1000; ++iter){

        float loss = 0;
        float delta_k_identity = 0;
        float delta_k_sin      = 0;
        float delta_k_cos      = 0;
        float delta_b          = 0;

        for(auto& item : datas){
            float predict = k_identity * item.year + k_sin * std::sin(item.year) + k_cos * std::cos(item.year) + b;
            float L = 0.5 * std::pow(predict - item.price, 2.0f) + lambda * (k_identity*k_identity + k_sin*k_sin + k_cos*k_cos + b*b);

            delta_k_identity += (predict - item.price) * item.year + k_identity * lambda;
            delta_k_sin      += (predict - item.price) * std::sin(item.year) + k_sin * lambda;
            delta_k_cos      += (predict - item.price) * std::cos(item.year) + k_cos * lambda;
            delta_b          += (predict - item.price) + b * lambda;

            loss += L;
        }

        if(iter % 100 == 0)
            cout << "Iter " << iter << ", Loss: " << setprecision(3) << loss << endl;

        k_identity = k_identity - lr * delta_k_identity;
        k_sin      = k_sin      - lr * delta_k_sin;
        k_cos      = k_cos      - lr * delta_k_cos;
        b          = b          - lr * delta_b;
    }
    printf(
        "模型参数：k_identity = %f, k_sin = %f, k_cos = %f, b = %f\n"\
        "数据集：xm = %f, ym = %f, xs = %f, ys = %f\n",
        k_identity,  k_sin, k_cos, b, mean.year, mean.price, stdval.year, stdval.price
    );

    float year          = 2023;
    float x             = (year - mean.year) / stdval.year;
    float predict       = x * k_identity + std::sin(x) * k_sin + std::cos(x) * k_cos + b;
    float predict_price = predict * stdval.price + mean.price;

    printf("预计%d年，房价将是：%.3f 元\n", (int)year, predict_price);
    return 0;
}

逻辑回归

逻辑回归虽然名字里带着“回归”二字，但实际上是进行（二）分类。

同样是以房价为例，如果我们有两个特征，分别是房屋面积和离地铁站的距离，而我们想要将所有房屋分为两类：让人幸福和让人不幸福。我们该怎么做呢？

基于线性回归的思考

如果继续按照线性回归的思路，我们或许会这样考虑：线性回归能够根据输入，输出一个标量预测值这个值的区间没有限制，从负无穷到正无穷都可以，而我们现在想要做的事情是二分类。那么只要把线性回归的输出值压缩在0到1之间，就可以将它看做是一种概率，再设定一个阈值（比如0.5），概率大于等于 0.5 则分类为幸福，小于 0.5 则分类为不幸福，就可以实现分类了。即我们希望用一个函数 $h_{\theta }(x)$ 将线性回归模型得到的输出值，映射在 0-1 之间，表示为：
$$f(x)=\{\begin{array}{l}1,h_{\theta }(x)\ge 0.5\\ 0,h_{\theta }(x)<0.5\end{array}$$

sigmoid函数

我们选择的是这样一个函数：
$$g(z)=\frac{1}{1+e^{-z}}$$
这就是著名的 Sigmoid 函数，其函数图像如下：

为什么选择 Sigmoid 函数？

• sigmoid函数的定义域是(-INF, +INF),而值域为(0, 1)。
• sigmoid函数是可导的。
• sigmoid函数对于给定的输入变量，会根据选择的参数计算输出变量=1的可能性，也就是说它的输出表示概率，都是0到1之间。

Sigmoid函数的推导

根据之前的内容，我们可以这么理解：逻辑回归=线性回归+Sigmoid函数：

• 线性回归：$z=kx+b$
• Sigmoid 函数：$y=\frac{1}{1+e^{-z}}$

则逻辑回归：$y=\frac{1}{1+e^{-(kx+b)}}$

接下来我们开始进行推导：

由 $h_{\theta }(x)=\frac{1}{1+e^{-(kx+b)}}$ 得（注意我们这里 $kx+b$ 是专门针对一元的，如果是多元的，常用 ${\theta }^{T}x$ 表示）：
$$ln\frac{h_{\theta }(x)}{1-h_{\theta }(x)}=kx+b$$
由上式，我们可以这样说：模型预测的，是正反例的对数概率

如果我们把 $h_{\theta }(x)$ 视为样本 $x$ 为正例的可能性，那么，自然在二分类中，$1-h_{\theta }(x)$ 即为负例可能性，两者的比值的对数，称为对数几率，这也是为什么逻辑回归也称为对数几率回归。同样，反过来，我们通过对样本为正反例可能性的伯努利分布推导（和指数族定义），也能得到上式的结果，因此，选择 Sigmoid 的原因在于他符合描述样本 $x$ 为正反例的可能性，而逻辑回归则是对该可能性建模，求解参数极大似然估计的过程

损失函数与伯努利分布

如果我们这里依然用之前线性回归中使用的 MSE 作为损失函数，则得到：
$$L(\theta )=\frac{1}{2}\sum _{i=1}^m(y^{(i)}-h_{\theta }(x^{(i)}){)}^2$$
其中 $h_{\theta }(x)=\frac{1}{1+e^{-(kx+b)}}$

我们观察这个函数机器导数的曲线发现，函数的导数值在多个区域接近 0，使得在部分位置时难以迭代，不利于求解全局最优解，因此，直接采用线性回归的方法，可能在某些时候具备效果，但肯定不是最合适的方案。

对于二分类问题，$y$ 的取值 0, 1服从伯努利分布，则有：

$y$ 为 1 时：
$$P(y=1|x;\theta )=h_{\theta }(x)$$
$y$ 为 0 时：
$$P(y=0|x;\theta )=1-h_{\theta }(x)$$
写成统一的形式：
$$P(y|x;\theta )=h_{\theta }(x{)}^y(1-h_{\theta }(x){)}^{1-y},y=0,1$$
为了估计参数 $\theta$ ，我们有似然函数：
$$L(\theta )=P(y|x;\theta )=\prod _{i=1}^m(y^{(i)}|x^{(i)};\theta )$$
代入后，得：
$$L(\theta )=\prod _{i=1}^m{h}_{\theta }(x^{(i)}{)}^{y^{(i)}}(1-h_{\theta }(x^{(i)}){)}^{1-y^{(i)}}$$
我们要求 $\theta$ 的极大似然估计，找到参数 $\theta$ 使得似然函数 $L$ 的值最大，即找到一组参数 $\theta$ ，使得：

• 若 $y$ 为正类，模型预测 $y=1$ 的概率最大；
• 若 $y$ 为负类，模型预测 $y=0$ 的概率最大。

同样还是取对数：
$$lnL(\theta )=\sum _{i=1}^m{y}^{(i)}ln{h}_{\theta }(x^{(i)})+\ast (1-y^{(i)})ln(1-h_{\theta }(x^{(i)}))$$
我们希望得到的是最大化值的参数估计，相当于是求 $L$ 函数的最大值。对于求最大值，相比线性回归，我们唯一区别是这里使用了梯度上升法，也可以给 $L$ 加个负号，然后仍然用梯度下降法，可能熟悉的朋友已经发现了，这里的 $-lnL(\theta )$就是我们熟悉的交叉熵损失：
$${\theta }^{+}=\theta -\alpha \frac{-\partial lnL(\theta )}{\partial \theta }$$
求导得：
$$\frac{\partial lnL(\theta )}{\partial {\theta }_j}=(y-h_{\theta }(x))x_j$$
即可得到第 $j$ 个参数的更新公式：
$${\theta }_j^{+}={\theta }_j+\alpha (y^{(i)}-h_{\theta }(x^{(i)}))x_j^{(i)}$$

代码

给出上海房价与幸福关系的预测代码

#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <cmath>
#include <tuple>
#include <iomanip>
#include <stdarg.h>

using namespace std;

struct Item{
	float area;
	float distance;
	float label;
};

#define INFO(...)  ::__printf(__FILE__, __LINE__, __VA_ARGS__)

void __printf(const char* file, int line, const char* fmt, ...){

	va_list vl;
	va_start(vl, fmt);

	// None   = 0,     // 无颜色配置
	// Black  = 30,    // 黑色
	// Red    = 31,    // 红色
	// Green  = 32,    // 绿色
	// Yellow = 33,    // 黄色
	// Blue   = 34,    // 蓝色
	// Rosein = 35,    // 品红
	// Cyan   = 36,    // 青色
	// White  = 37     // 白色
	/* 格式是： \e[颜色号m文字\e[0m   */
	printf("\e[32m[%s:%d]:\e[0m ", file, line);
	vprintf(fmt, vl);
	printf("\n");
}

/* 通过csv文件加载数据 */
vector<Item> load_data(const string& file){

	vector<Item> output;
	fstream ifile(file, ios::binary | ios::in);

	string line;
	getline(ifile, line);

	while(getline(ifile, line)){
		int p0 = line.find(",");
		int p1 = line.find(",", p0 + 1);
		Item item;
		item.area     = atof(line.substr(0, p0).c_str());
		item.distance = atof(line.substr(p0+1, p1).c_str());
		item.label    = atof(line.substr(p1+1).c_str());

		// cout << item.area << ", " << item.distance << ", " << item.label << endl;
		output.emplace_back(item);
	}
	return output;
}

/* 计算数据的均值和标准差 */
tuple<Item, Item> compute_mean_std(const vector<Item>& items){

	Item mean{0, 0}, stdval{0, 0};

	for(auto& item : items){
		mean.area  += item.area;
		mean.distance += item.distance;
	}
	mean.area  /= items.size();
	mean.distance /= items.size();

	for(auto& item : items){
		stdval.area  += std::pow(item.area - mean.area, 2.0f);
		stdval.distance += std::pow(item.distance - mean.distance, 2.0f);;
	}
	stdval.area  = std::sqrt(stdval.area / items.size());
	stdval.distance = std::sqrt(stdval.distance / items.size());
	return make_tuple(mean, stdval);
}

double sigmoid(double x){
	return 1 / (1 + std::exp(-x));
}
    
int run(){

	auto datas = load_data("shanghai.csv");

	Item mean, stdval;
	tie(mean, stdval) = compute_mean_std(datas);

	/* 对数据进行减去均值除以标准差，使得均值为0，标准差为1 */
	for(auto& item : datas){
		item.area     = (item.area - mean.area) / stdval.area;
		item.distance = (item.distance - mean.distance) / stdval.distance;
	}
	
	float k_area = 0.1;
	float k_dist = 0.1;
	float b = 0;
	float lr = 0.1;

	for(int iter = 0; iter < 1000; ++iter){

		float loss = 0;
		float delta_k_area = 0;
		float delta_k_dist = 0;
		float delta_b      = 0;

		for(auto& item : datas){
			float predict = k_area * item.area + k_dist * item.distance + b;
			double logistic = sigmoid(predict);
			float L = -(std::log(logistic) * item.label + std::log(1 - logistic) * (1 - item.label));

			delta_k_area += (logistic - item.label) * item.area;
			delta_k_dist += (logistic - item.label) * item.distance;
			delta_b      += (logistic - item.label);
			loss += L;
		}

		if(iter % 100 == 0)
			cout << "Iter " << iter << ", Loss: " << setprecision(3) << loss << endl;

		k_area = k_area - lr * delta_k_area;
		k_dist = k_dist - lr * delta_k_dist;
		b      = b      - lr * delta_b;
	}

	INFO("模型参数：k_area = %f, k_dist = %f, b = %f", k_area, k_dist, b);
	INFO("数据集：area_mean = %f, dist_mean = %f, area_std = %f, dist_std = %f", 
		mean.area, mean.distance, stdval.area, stdval.distance
	);

	INFO(""
		"k_area = %f\n"
		"k_distance = %f\n"
		"b = %f\n"
		"area_mean = %f\n"
		"area_std = %f\n"
		"distance_mean = %f\n"
		"distance_std = %f"
	, k_area, k_dist, b, mean.area, stdval.area, mean.distance, stdval.distance);

	float area          = 100;
	float distance      = 2000;
	float norm_area     = (area - mean.area) / stdval.area;
	float norm_dist     = (distance - mean.distance) / stdval.distance;
	float predict       = k_area * norm_area + k_dist * norm_dist + b;
	float logistic      = sigmoid(predict);

	INFO("在上海，对于房屋面积 %.0f 平方米，距离地铁 %.0f 米的住户而言，他们觉得生活是 【%s】 的.",
		area, distance, logistic > 0.5 ? "幸福" : "并不幸福"
	);
	return 0;
}
int main(){
    return run();
}

Ref

1. 【机器学习】逻辑回归十分钟学会，通俗易懂（内含spark求解过程）