一元线性回归

​由于笔者的数学不太好，而且网上关于一元线性回归的文章有很多，所以相关内容大家可以查找一下，这里我就简单的搬运一下简单概念。

一元线性回归的方程：
$$h(x)={\beta }_0+{\beta }_{1}x$$
其中第一个参数${\beta }_0$为截距，第二个参数${\beta }_1$为斜率。

代价函数

​ 回归分析的主要目的是通过已有的信息，去推测出未知的信息。通过一个例子大家可能会更深刻的理解回归分析的目的。

上图为广告费与销售额的关系图，虚线为我们的一元线性回归方程，通过回归分析，我们可以预测当广告费为14万元时，我们的销售额可能是30万元。

​回归分析属于统计学问题，这就说明在给定自变量x时，我们是无法准确的得出因变量y的。但是我们可以通过一些方法去**“减少误差”**，使得预测的结果尽量的接近真实值。所以我们引入了代价函数，其使用最小二乘法的理论去减少这种误差。代价函数数学公式如下：
$$J({\beta }_0,{\beta }_1)=\frac{1}{2m}\sum _{i=1}^m[y_i-h(x_i){]}^2$$
​ 为了更好的理解代价函数，这里使第一个参数为零，然后观察一下图片（由于是从视频中截取的图片，所以这里的参数又β变成了θ，是我太懒了😂）

当θ=1时，此时回归方程贯穿每一个点，所以误差为零，J(θ)值如右图所示。

当θ=0.5时，我们一顿计算可得，J(θ)≈0.58

当θ=0时，我们同样可以计算出J(θ)的值。

如果我们大量的将θ带入，我们将得到如下的代价函数图

由此我们可以清晰地看到，J(θ)在某一点处是可以取到最小值的，这就是我们的引入代价函数的目的：通过调整参数来减少误差，使得预测的结果尽量的接近真实值。

注：最小二乘法的公式如下：
$$\sum [y_i-({\beta }_0+{\beta }_1{x}_i){]}^2$$

梯度下降法

​梯度下降法就是一个很好的调整参数的方法，它可以通过下面公式不断调整${\beta }_0$和${\beta }_1$的值，从而得到一个全局最小值，或者是一个局部最小值。
r e p e a t     u n t i l     c o n v e r g e n c e { β j = β j − α ∂ ∂ β j J ( β 0 , β 1 ) f o r   j = 1   a n d   j = 0 } repeat\ \ \ until\ \ \ convergence\{ \\ β_j = β_j - α\frac{∂}{∂β_j}J(β_0,β_1) \\ for\ j = 1\ and\ j=0 \\ \} repeat   until   convergence{βj​=βj​−α∂βj​∂​J(β0​,β1​)for j=1 and j=0}
其中α称为学习率(learning rate)，而学习率就是步长，学习率大，一次跨越的距离就远，这样可能会错过全局（局部）最小值点；学习率小，一次跨越的距离就短，这样会花费比较多的时间来处理数据。

求偏导后：
$$\begin{gathered} j=0:{\beta }_0={\beta }_0-\alpha (\frac{1}{m}\sum _{i=1}^m(y_i-{\beta }_0-{\beta }_1{x}_i)\ast -1) \\ ={\beta }_0-\alpha (-\frac{1}{m}\sum _{i=1}^m[y_i-({\beta }_1{x}_i+{\beta }_0)]) \\ ={\beta }_0-\alpha (-\frac{1}{m}\sum _{i=1}^m[y_i-(h(x_i)]) \\ j=1:{\beta }_1={\beta }_1-\alpha (\frac{1}{m}\sum _{i=1}^m(y_i-{\beta }_0-{\beta }_1{x}_i)\ast -x_i) \\ ={\beta }_1-\alpha (-\frac{1}{m}\sum _{i=1}^m{x}_i[y_i-({\beta }_1{x}_i+{\beta }_0)]) \\ ={\beta }_1-\alpha (-\frac{1}{m}\sum _{i=1}^m{x}_i[y_i-(h(x_i)]) \end{gathered}$$

代码实现：

import numpy as np
import matplotlib.pyplot as plt
#载入数据
data = np.genfromtxt("data.csv",delimiter=",")#delimiter分隔符，此分隔符为，
x_data = data[:,0]
y_data = data[:,1]
plt.scatter(x_data,y_data)
plt.show()

#学习率 learning rate
lr = 0.0001
#截距
b = 0
#斜率
k = 0
#最大迭代次数
epochs = 50

#最小二乘法
def compute_error(b, k, x_data, y_data):
    totalError = 0
    for i in range(0, len(x_data)):
        totalError += (y_data[i] - (k * x_data[i] + b)) ** 2
    return totalError / float(len(x_data)) /2.0 #最小二乘法法公式

def gradient_descent_runner(x_data, y_data, b, k, lr, epochs):
    #计算总数据量
    m = float(len(x_data))
    #循环epochs次
    for i in range(epochs):
        b_grad = 0
        k_grad = 0
        #计算梯度的总和再求平均
        for j in range(0, len(x_data)):
            b_grad += -(1/m)*(y_data[j]-(k * x_data[j] + b))
            k_grad += -(1/m)*x_data[j] * (y_data[j] - (k * x_data[j] + b))
        #更新b和k
        b = b - (lr * b_grad)
        k = k - (lr * k_grad)
         # 每迭代5次，输出一次图像
''' 方面查看变化趋势
			if i % 5==0:
            print("epochs:",i)
            plt.plot(x_data, y_data, 'b.')
            plt.plot(x_data, k*x_data + b, 'r')
            plt.show()
            '''
    return b, k

b, k = gradient_descent_runner(x_data, y_data, b, k, lr, epochs)
#画图
plt.scatter(x_data, y_data, c='b')
plt.plot(x_data, k*x_data + b, 'r')
plt.show()

一元线性回归分析

一元回归分析的基本概念
回归模型的建立一般包括：
（1）通过某事物现，转化为具体问题；
（2）确定指标变量，收集整理数据，并构建模型进行参数估计；
（3）模型的检验，当模型检验不通过时，需要重新修改模型；
（4）模型的应用，得出结论，运行给出决策等。

基本概念
通常我们要先收集与研究相关的数据的一组或者多组样本，为直观观察数据的分布规律，我们可以将收集到每组数据绘制二维数据散点图。

一元回归分析的参数估计
一元回归模型的参数估计一般采用极大似然法与最小二乘法，其中最常用的是最小二乘法估计。

相关系数的检验

实际例子：

#给出因变量自变量的值
x <- c(3.4,1.8,4.6,2.3,3.1,5.5,0.7,
       3.0,2.6,4.3,2.1,1.1,6.1,4.8,3.8)
y <- c(26.2,17.8,31.3,23.1,27.5,36.0,14.1,
       22.3,19.6,31.3,24.0,17.3,43.2,36.4,26.1)
plot(x,y)
abline(lm)
#建立回归方程
lm <-lm(y~x)
#输出回归分析的结果
summary(lm)

summary()函数用于结果的展示，其中包括残差的最大值最小值、中位数、回归系数、以及回归系数估计值、显著性检验的值，它包含的结果非常多。其中intercept为截距，也就是回归常数项的估计值，Estimate是回归系数的估计值，std是标准差，后面的为T检验的值与检验P值与显著性结果。判断显著性的方法除了P值外，还可以使用方差分析来判断。
拟合模型

方差分析

#方差分析
anova(lm)

从方差分析的结果来看，自由度为13，sum-sq为回归平方和为SSR=841.77，与残差平方和SSE=69.75，从P值得结果可知，回归方程是显著的。

计算相关系数与显著性检验
回归系数的检验方式有很多，比如kendall检验，spearman检验等，在这里采用常用的pearson相关系数检验。

#计算相关系数并进行显著性检验、皮尔选相关系数检验
cor(x,y,method = "pearson")

#相关系数显著性检验
cor.test(x,y,alternative = "two.sided",
         method = "pearson",conf.level = 0.95)

最后得到相关系数为：0.9609777，95%置信区间为[0.8837722,0.9872459]，其中P值小于0.05，则可以说明两者具有高度显著的线性相关。

残差分析

#残差分析、保留3位小数
e <- resid(lm)
round(e,3)
#残差图的绘制
plot(x,e)

回归系数的置信区间

#计算回归系数的置信区间
confint(lm)

可见常数项的置信区间为[7.209605,13.346252],回归系数的置信区间为[4.070851 5.767811]。

预测值

#模型预测
new_data <- data.frame(x = 3.3)
#计算预测值，给定置信区间0.95
y.pred <- predict(lm,new_data,interval = "prediction",level = 0.95)
#计算因变量平均值的预测区间的置信区间
y.conf <- predict(lm,new_data,interval = "confidence",level = 0.95)

模型诊断

#回归诊断
library("lindia")
gg_diagnose(lm)

其中第一为残差直方图，可以用来判断残差的正态性是否存在极端值；第二张自变量与残差的关系图；第三张为预测值与残差的关系图，可以用来判断方差齐性的评估，每个点随机的分布在参考线的两侧，没有规律；第四张为残差的Q-Q图，用于残差的正态性判断；第五张当参考线水平时证明方差齐性，第六张为杠杆和残差的关系图，第七张为Cook’s Distance图，用于识别是否存在强影响点。

完整代码

#给出因变量自变量的值
x <- c(3.4,1.8,4.6,2.3,3.1,5.5,0.7,
       3.0,2.6,4.3,2.1,1.1,6.1,4.8,3.8)
y <- c(26.2,17.8,31.3,23.1,27.5,36.0,14.1,
       22.3,19.6,31.3,24.0,17.3,43.2,36.4,26.1)
#绘制散点图
plot(x,y)
#添加趋势线
abline(lm)
#建立回归方程
lm <-lm(y~x)
#输出回归分析的结果 
summary(lm)

#方差分析
anova(lm)

#计算相关系数并进行显著性检验、皮尔选相关系数检验
cor(x,y,method = "pearson")

#相关系数显著性检验
cor.test(x,y,alternative = "two.sided",
         method = "pearson",conf.level = 0.95)
#残差分析、保留3位小数
e <- resid(lm)
round(e,3)

#绘制残差图
plot(x,e)

#标准化残差
zre <- e/2.316
#计算学生化残差
sre <- rstandard(lm)

#计算回归系数的置信区间
confint(lm)

#模型预测
new_data <- data.frame(x = 3.3)
#计算预测值，给定置信区间0.95
y.pred <- predict(lm,new_data,interval = "prediction",level = 0.95)
#计算因变量平均值的预测区间的置信区间
y.conf <- predict(lm,new_data,interval = "confidence",level = 0.95)

#回归诊断
library("lindia")
gg_diagnose(lm)

代码如下：

import numpy as np
import matplotlib.pyplot as plt
from numpy import array

""" 使用最小二乘法，拟合出一元线性回归模型：z = wx + b。
一元的意思是样本x通过一个属性描述，原本可能是矢量x_i = (x_i1, x_i2...,x_id)被例如颜色，大小...
属性描述，现在只有一个x_i1描述，则直接把矢量x_i看成标量，w也是标量

计算出使得损失最小的w和b，
画出拟合直线和原始的散点图
点距离拟合直线越远，代表误差越大
"""


# 画出样例的真实分布，输入样本x和真实标记y
def plot_origin(points: array) -> None:
    """
    :param points: array类型的二维数组
    :return:
    """
    arr_x = points[:, 0]  # return list: 所有元素（子数组）中的第一个元素：x
    arr_y = points[:, 1]  # return list: 所有元素（子数组）中的第二个元素：y

    # 画出散点图：照理说学的时间越长，考试分数越高
    plt.scatter(arr_x, arr_y)
    plt.show()


# 2. 策略：求均方误差（损失函数）
def compute_cost(w: float, b: float, points: array) -> float:
    """ 计算均方误差（损失函数）：预测输出和真实标记之间的差距: y - (wx + b)
                z = wx + b 为我们要拟合的线性模型
    :param w: 线性模型参数
    :param b: 线性模型参数
    :param points: array类型的二维数组：所有样例
    :return:  输出E(w，b) 均方误差值
    """
    total_cost = 0
    m = len(points)  # 样本个数m

    # 计算均方误差
    for i in range(m):
        x_i = points[i, 0]  # 第i个样例的第一个元素：x
        y_i = points[i, 1]  # 第i个样例的第二个元素：y
        total_cost += (y_i - w * x_i - b) ** 2
    return total_cost / m  # 均方误差


"""3. 算法：拟合：学得z = wx+b 近似于真实标记y
使用基于均方误差最小化的 最小二乘参数估计
求能使得均方误差最小的w和b：损失函数分别对w，b求偏导=0
以下代码都基于公式推导出来的w，b的表示方法
"""


# 求列表内元素平均值
def avg(lst):
    l = len(lst)
    return sum(lst[i] for i in range(l)) / l


def fit(points: array) -> tuple:
    x_avg = avg(points[:, 0])  # 样本均值
    m = len(points)
    # 求w
    numerator, denominator = 0, -m * x_avg ** 2  # w公式的分子，分母
    for i in range(m):
        x_i, y_i = points[i, 0], points[i, 1]
        numerator += y_i * (x_i - x_avg)
        denominator += x_i ** 2
    w = numerator / denominator

    # 求b
    sum_y = 0
    for i in range(m):
        x_i, y_i = points[i, 0], points[i, 1]
        sum_y += y_i
    b = (sum_y - w * x_avg * m) / m
    return w, b


# 画出拟合函数：一元线性回归模型
def plot_fit(arr_x: array, arr_y: array) -> None:
    plt.scatter(arr_x, arr_y)  # 画散 点 图
    # array类型可以直接对每个元素乘上一个常数，不用for循环慢慢一个个乘
    predict_y = w * arr_x + b  # 拟合的线性模型：预测标记y
    plt.plot(x, predict_y, c='r')  # 画 经过x，y的曲线/直线
    plt.show()


if __name__ == '__main__':
    points = np.genfromtxt('D:\Pycharm\code\MY_machine_learning\mydata\data.csv', delimiter=',')  # array
    x = points[:, 0]  # return array: 所有元素（子数组）中的第一个元素：x
    y = points[:, 1]  # return array: 所有元素（子数组）中的第二个元素：y

    w, b = fit(points)
    print('w, b分别为', w, b)

    print('损失为：', compute_cost(w, b, points))

    plot_fit(x, y)

运行结果如下：

 附data.csv网盘提取

：百度网盘 请输入提取码

提取码：rstx