手写一个简单的线性回归、岭回归

自定义类，实现线性回归和岭回归，完成主体部分，一个字--清爽简洁。

内容：

1.简单复制粘贴下数学原理

2.代码及测试

3.思考

一、数学原理

由于数学公式太多，直接粘贴了；

图片复制自：Python手写实现线性回归（从原理到代码） - 知乎 (zhihu.com)

python机器学习手写算法系列——线性回归_python手写代码实现一元线性回归设计思想-CSDN博客

 

 解法1.最小二乘法：

在数据量大、样本稍微有点难度下，很难算出来

解法2：梯度下降

万物皆可迭代，梯度下降是个不错的选择

二、数据及代码

2.1 数据

import numpy as np
import pandas as pd 
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
np.set_printoptions(suppress=True)
# 导入数据
diabetes_x,diabetes_y = datasets.load_diabetes(return_X_y=True,as_frame=True)
diabetes_x

自带的糖尿病数据集，其实该数据集用单纯的线性回归效果很一般。

 首先我们瞄一眼sklearn中自带的线性回归的参数，对比分析下要实现什么功能：

1是否要带截距，肯定要啊！2是否复制数据 ，3是用几个核心，4是强制让所有系数都为正，看起来不是很聪明。

2.2 代码

比较简单，如下，但是悄悄埋下一个坑，后面分解：

# 避免重名
class MyLinearRegression(object):
    def __init__(self, alpha=0.01, n_iters=1000,norm=False):
        self.alpha = alpha # 学习率
        self.n_iters = n_iters # 迭代次数
        self.norm = norm # 是否标准化
        self.losslist=[] # 记录损失函数变化好画图

    def fit(self, x, y):
        num_samples, num_features = x.shape
        if self.norm:
            # 标准化
            x = (x-x.mean(axis=0))/x.std(axis=0)

        # 初始化模型参数
        self.theta = np.zeros(num_features) # 系数
        self.b = 0 # 截距
        for i in range(self.n_iters):
            # 矩阵预算 (m,n)dot(n,1)
            y_pred = np.dot(x, self.theta) + self.b
            # 计算损失 (均方误差, Mean Squared Error)
            loss = 1/2 * np.mean((y - y_pred) ** 2) # 乘个1/2
            self.losslist.append(loss)
            # 计算梯度
            dw = -(1 / num_samples) * np.dot(x.T,(y - y_pred))
            db = -(1 / num_samples) * np.sum(y - y_pred)
            # 更新参数
            self.theta -= self.alpha * dw
            self.b -= self.alpha * db
        
    def predict(self, x):
        if self.norm:
            x = (x-x.mean(axis=0))/x.std(axis=0)
        return np.dot(x, self.theta) + self.b
    
    def plotloss(self,title='Loss_descend'):
        plt.plot(range(self.n_iters),self.losslist,color='lightblue',alpha=0.8)
        plt.title(label=title)
        plt.show()

2.3对比分析

那我们来看看跟sklearn自带的有什么区别：

model_test = MyLinearRegression(n_iters=5000,norm=False,alpha=0.05)
model_test.fit(x_train,y_train)
y_ = model_test.predict(x_test)
from sklearn.metrics import mean_squared_error
error1 = mean_squared_error(y_test,y_)
error1

-> 3312.5459764918764

用sklearn的看看：

sklearn看起来要强一点，但强得不多，应该是学习率和迭代次数的区别，这个数据集用单纯的线性回归还是不行，不过自带好找的数据集并不是那么多。

画图看看：

 画了又没有完全画，把索引调一下：

# 将y_test从小到大排列
# 要先重塑索引
y_test_reindex= y_test.reset_index(drop=True)
y_test_sort_index = np.argsort(y_test_reindex)
# 画图
plt.plot(range(len(y_test)),y_[y_test_sort_index],color='green',label='pred',alpha=.5)
plt.plot(range(len(y_test)),y_test_reindex[y_test_sort_index],color='blue',label='true',alpha=.5)
plt.plot(range(len(y_test)),y_2[y_test_sort_index],color='orange',label='sklearn',alpha=.5)
plt.legend()

看起来好一点。

注意看代码的话，我们悄悄加入了是否归一化的选项，来看看归一化是否好点。

这个数字大到离谱，但自带的数据集，用于回归的好像就只有这一个，波士顿房价因为零元购而被正确掉了。

三、思考

要不要归一化呢？

如果要归一化，不要在算法里面搞，在数据处理阶段就搞！！！

Z-score标准化，z = x_i - mu / sigma

只有保证训练、测试是先标准化，再切分的，才有可比性。

假如在算法里面标准化，几十几百也会被弄到0附近，几千几万也会被弄到0附近，比如你一家店日活几百跟几万，差距肯定很大啊！

3.2 岭回归

顺手把岭回归也搞一下：

删掉偷偷藏下的雷--norm

也是梯度下降法，主要在损失函数和梯度计算上有区别：

其实还缺点东西，比如tolerance,不过这样讲还有solver呢，会越写越多，就到这吧

class MyRidge(object):
    def __init__(self, learning_rate=0.01, n_iters=1000,lambda_=1):
        self.learning_rate = learning_rate # 学习率改下避免和L1重复
        self.n_iters = n_iters # 迭代次数
        self.losslist=[] # 记录损失函数变化好画图
        self.lambda_ = lambda_ # L2正则，加_避免重名

    def fit(self, x, y):
        num_samples, num_features = x.shape
        # 初始化模型参数
        self.theta = np.zeros(num_features) # 系数
        self.b = 0 # 截距
        for i in range(self.n_iters):
            # 矩阵预算 (m,n)dot(n,1)
            y_pred = np.dot(x, self.theta) + self.b
            # 都乘1/2比较方便
            loss = 1/2 * np.mean((y - y_pred) ** 2) +1/2*self.lambda_*np.sum(self.theta**2)
            self.losslist.append(loss)
            # 系数的梯度&更新
            gradient = np.dot(x.T,(y_pred-y))+self.lambda_*self.theta
            gradient = gradient/num_samples
            self.theta -= self.learning_rate*gradient
            # 截距偏导跟正则项没关系
            db = -(1 / num_samples) * np.sum(y - y_pred)
            self.b -= self.learning_rate * db
        
    def predict(self, x):
        return np.dot(x, self.theta) + self.b
    
    def plotloss(self,title='Loss_descend'):
        plt.plot(range(self.n_iters),self.losslist,color='lightblue',alpha=0.8)
        plt.title(label=title)
        plt.show()

 因为代码比较简单，所以直接全部重写了，实际上聪明一点应该继承的写法：

# 使用继承的方法
class MyRidge(MyLinearRegression): 
    def __init__(self,alpha=0.05,n_iters=100,lambda_=1): # 必须附上默认值
        super().__init__(alpha, n_iters) # 调用父类的构造方法
        self.lambda_ = lambda_ # 
    def fit(self, x, y):
        num_samples, num_features = x.shape
        self.theta = np.zeros(num_features) #
        self.b = 0 
        for i in range(self.n_iters):
            y_pred = np.dot(x, self.theta) + self.b
            loss = 1/2 * np.mean((y - y_pred) ** 2) +1/2*self.lambda_*np.sum(self.theta**2)
            self.losslist.append(loss)
            gradient = np.dot(x.T,(y_pred-y))+self.lambda_*self.theta
            gradient = gradient/num_samples
            self.theta -= self.alpha*gradient
            db = -(1 / num_samples) * np.sum(y - y_pred)
            self.b -= self.alpha * db
    # 这样少写下面两个方法

 结果差不多