20210516 线性回归模型案例

# 模型需要假设，最简单的假设就是线性假设
# 每一个因子，每一个特征与预测变量的关联都是线性的
# 如何求解模型？数据量很大时，唯一比较现实的是不假设能够求得准确的解析解
# 而是用一个代价函数来衡量不准确的解与真实情况的差距

# 用暴力搜索的办法可以找到参数，更好的办法叫做梯度优化
# 不要小看线性假设
# 在什么都不知道的情况下，可以先用线性回归进行尝试，判断是否存在能够使用的因子；
# 什么叫回归？从一种特征预测另一种特征的方法；
# 回归特指预测的特征是一组连续的实数，比如价格
# 是从不正确到正确的过程
# 洞见：从模型中提取和预测相关的内容

# 模型详解
# 引入房价数据集
1-1
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import Image
import warnings
warnings.filterwarnings('ignore')

data = pd.read_csv('train.csv')
print(data.shape,data.shape)
# -->  (1460,81) (1460,81)
print(data.head())
# -->
# 返回一个 5 rows x 81 columns 的数据列表
# 特征有 81 个特征，5个 数据点

1-2
# 数据预处理
# 数太大，归一化，让数据的分布处于同一区间，
# 我们选择一种最简单的数据调整方法，每个数除以其最大值。
# 选择数据集合中的几个重要的特征
data_select = data[['BedroomAbvGr','LotArea','Neighborhood','SalePrice']]
# rename 是 pandas 重命名
data_select = data_select.rename(columns = {'BedroomAbvGr':'room','LotArea':'area'})
# dropna 去掉空值，因为空值会报错
data_select = data_select.dropna(axis = 0)
# col in np.take 是对现有特征进行处理；
# 除以最大值，三个用数字表达的特征，目的是让值在 0-1之间
# 处理起来会比较整齐
# 跳过了 neighborhood 是因为它是文字表达的；
for col in np.take(data_select.columns,[0,1,-1]):
    data_select[col] /= data_select[col].max()
data_select.head()
# SalePrice 是要预测的目标，叫 target
# 而 'BedroomAbvGr','LotArea','Neighborhood' 作为 3 个 用来预测的特征

1-3
# 预测和训练集的分割
# 如何评价算法的参数是好的，通常并不是错误量最小的数据
# 而是将现有数据分成两部分，一部分叫训练数据，一部分叫测试数据
# 训练数据用来得到参数，测试数据用来判断得到的参数好不好
# 为什么这样做？因为不能认为模型在某种数据上好，在其他数据上也好
# 所以用于训练的数据通常与用于测试的数据要求是不同的
# 这就好比考试的练习题和考试真题是不会一样的
from sklearn.cross_validation import train_test_split
train, test = train_test_split(data_select.copy(), test_size = 0.9)
# train.describe 可以看到特征的大概统计特性
train.describe()

2-1-1
# 构建模型
# 线性回归模型
# 从最简单的情况出发，构建一台预测价格的机器，
# 假设函数 h(x) = wx + b 是线性的。
# def linear 构建线性模型；这是一个 numpy 的基础运算
def linear(features, pars):
    price = np.sum(features*pars[:-1], axis = 1) + pars[-1]
    return price

2-2-2
### 模型测试
# 试一试你的模型预测价格， 和真实价格做比较
# par1 = 0.1, par2 = 0.1
# 通过函数 linear 就可以得到预测的数据；用 linear 这样的小小的函数，就可以实现线性回归模型的模样
train['predict'] = linear(train[['room', 'area']].values, np.array([0.1, 0.1, 0.0]))
train.head()
# 可以看到，在该参数下，模型的预测价格和真实价格有较大的差距。那么寻找合适的参数值是我们需要做的事情。

3-1
### 代价函数
# 预测函数为 h(x) = wx + b 度量预测错误，代价函数是误差平方和（RSS）
def mean_squared_error(y_pred, y):
    '''
    这个函数返回 平方和误差； 也就是 y值 和 y_pred预测值的差的平方的求和
    '''
    return sum(np.array(y_pred - y) ** 2)

def Cost(df, features, pars):
    '''
    用前面定义的 linear 函数返回预测值，预测值和真实值做平方误差；得到返回值 cost
    '''
    df['predict'] = linear(df[features].values, pars)
    cost = mean_squared_error(df.predict, df.SalePrice) / len(df)
    return cost

cost = Cost(train, ['room', 'area'], np.array([0.1, 0.1, 0.0]))
print(cost)
# --> 0.0521960417447
# 得到的 cost，是随便用 [0.1, 0.1, 0.0] 算出的，一定是很差的

3-2
### 暴力搜索
# 用暴力搜索求解参数
# 我们需要找到一组参数，使线性回归模型的代价最少！
# 最简单的办法是用格点搜索法。为了减少演算时间，只取两个特征 w1 w2，假设模型的截距项b = 0
Xs = np.linspace(0, 1, 100)
Ys = np.linspace(0, 1, 100)
Zs = np.zeros([100,100])

Xs,Ys = np.meshgrid(Xs,Ys)
print(Xs.shape, Ys.shape)
# --> ((100,100),(100,100))

3-2-1
# 进行一个 double loop 循环，遍历所有 w1 和 w2 的值；用Costs列表记录数值
W1=[]
W2=[]
Costs =[]
for i in range (100):
    for j in range(100):
            W1.append(0.01*i)
            W2.append(0.01*j)
            Costs.append(Cost(train,['room', 'area'], np.array([0.01*i, 0.01*j, 0.])))
index = np.array(Costs).argmin()
print (W1[index], W2[index],Costs[index])
# --> 0.52 0.99 0.00995275656167
# 返回列表中的最小值，和对应的 index索引下的参数值；
# 这样就求得了使cost最小的参数，也就是 0.52 和 0.99

3-2-2 
# 为了验证结果是正确的，可以进行画图
# 把所有代价的值通过格点搜索的方法表达在三维图像上
from mpl_toolkits.mplot3d import Axes3D
fig=plt.figure()
ax=fig.add_subplot(111,projection='3d')
ax.view_init(5,-15)
ax.scatter(W1,W2,Costs,s=10)
ax.scatter(0.58, 0.28, zs = Cost(train,['room', 'area'], np.array([0.58, 0.28, 0.0]) ), s=100,
           color='red')
plt.xlabel('rooms') # 房屋数量
plt.ylabel('lotArea')   # 面积
# 房屋数量和面积组成平面，对于任何上面的点，都有一个 z 值与之对应
# z 值就是代价，组成了一个类似抛物面的曲面，cost是曲面上的最低点
# 从中可以看出我们格点搜索出来的点，确实具备有更小的 Cost.
# (z值就是代价；需要找到z值最小的点)

3-3   ### 梯度下降法
# 格点搜索法计算量通常很大，更多时候，尤其是在参数很多的时候，一般用梯度下降法。
# 代价函数就是一个山谷，梯度就是山谷上升最快的方向，它的反面就是下降最快的方向
# 梯度下降就是逆着梯度，向下走，走到最低点；这个最低点的参数就是 我们想要的参数
# pars 是 w和b的总和

3-3-1 梯度
def gradient(train,features, pars):
    '''
    此函数返回 一个梯度向量，参数的每个分量发生微小改变，代价函数会发生轻微变化
    求解代价函数在每个分量上的变化率，就会组成一个向量，向量的维度等于参数的数量
    向量在每个分量上的值，就是每个方向上参数的变化率
    '''
    Gradient = np.zeros(len(pars))
    for i in range(len(pars)):
        pars_new = pars.copy()
        pars_new[i] += 0.01
        Gradient[i] = (Cost(train, features, pars_new) - Cost(train,features, pars))/0.01
    return Gradient
print(gradient(train, ['room', 'area'],[0.2, 0.1, 0]))
# --> array([-0.12062854,-0.01706893,-0.32515586])
# 这里有 三个参数 w1 w1 和 b，求得任意一个点，
# 比如 [0.2, 0.1, 0] 这个点的梯度的值 就是  array([-0.12062854,-0.01706893,-0.32515586])
# 返回的向量，直接指向山谷的最低点

3-3-2 梯度下降
def GradientDescent(data, epochs, lr, features, pars):
    '''
    梯度下降，这个函数也是一个循环，需要做的就是从原有的点改变一点点
    :param data:  训练的数据
    :param epochs:迭代的次数
    :param lr: 学习速率（每次在梯度方向上下降的数量值）；因为下降的太快或者太慢都不好
    '''
    Costs = []
    for i in range (epochs):
        grad = gradient(data, features, pars)
        if i%50 == 0:   # 迭代 50次，记录代价
            Costs.append(Cost(data, features, pars))
        pars -= grad*lr
    print ('w = ', pars)
    return pars, Costs
# 结果将返回最终的优化参数，w1 w2 和 b
pars,Costs = GradientDescent(train, 500, 0.002,['room', 'area'],[0.1, 0.1, 0] )
# --> w = [0.16047375  0.10910733  0.15654558]
# 返回的是一组最终被优化了的参数 w1 w2 和 b

# 查看梯度下降的过程
plt.plot(Costs)
plt.xlabel('time')
plt.ylabel('cost')

4-1 ## 预测和模型评估
# 用得到的参数预测和评估模型
# 用学习出来的模型做预测
cost = Cost(train,['room', 'area'], pars)
print(cost)
# --> 0.010393847591213893

4-1-1 # 查看预测的均方误差
# 从 sklearn中引入 mean_squared_error 测量误差的工具包
from sklearn.metrics import mean_squared_error
train['predict'] = linear(train[['room', 'area']].values, pars )
print('MSE train: %.3f' % (mean_squared_error(train['SalePrice'], train['predict'])))
# --> MSE train: 0.010
# 经过处理，可以看到MSE的误差，和代价函数是一致的，只是做了平均
train.head()

4-1-2 前面的是训练集，现在看测试集的误差
# 这是测试集的误差，更换数据集检查准确率，留了90%的数据作为测试集
test['predict'] = linear(test[['room', 'area']].values, pars)
print('MSE test: %.3f' % (mean_squared_error(test['SalePrice'], test['predict'])))
# --> MSE test: 0.011
# 之比之前的 差一点，说明参数 取得还是可以的
test.head()

cost = Cost(test,['room', 'area'], pars)
print(cost)
# --> 8.5112068780221547e-06

5-1
## 如何做的更好
# 首先将方法变成一个类，好处是方便扩展，把所有和算法有关的量都整合一起
# 结果是否理想， 能够如何优化，是否存在一个模型，使得误差为0？ 增加特征试一试！
# 我们把之前的模型封装起来 ， 并且使用公式运算梯度， 比较一下和之前的速度差别！
class LinearRegressionGD(object):
    def __init__(self, eta=0.001, n_iter=20):
        self.eta = eta              # eta 是学习速率
        self.n_iter = n_iter     # n_iter 是训练次数

    def fit(self, X, y):
        '''
        梯度下降，使用了 np 简化前面的过程
        '''
        self.w_ = np.zeros(1 + X.shape[1])
        self.cost_ = []

        for i in range(self.n_iter):
            output = self.net_input(X)
            errors = (y - output)
            self.w_[:-1] += self.eta * X.T.dot(errors)
            self.w_[-1] += self.eta * errors.sum()
            cost = (errors**2).sum() / 2.0
            self.cost_.append(cost)
        return self

    # net_input 和 predict 都是在fit后，用于做预测
    def net_input(self, X):
        return np.dot(X, self.w_[:-1]) + self.w_[-1]

    def predict(self, X):
        return self.net_input(X)
# 通常的一个机器学习算法的类，就是包含这些部分
# 第一 初始化，初始化中包含的属性一定有学习速率和循环次数，这是算法训练的基本参数
# 训练过程里，需要使用梯度下降求点参数
# 得到参数后，测试，输入特征，得到结果，将输入变成输出

5-2
# 为了提高准确率，加入地区的特征，这是字符串特征；
# 用 get_dummies 函数，把这些不同类别的地区变成数字，这是哑变量的方法，
data_select = pd.get_dummies(data_select)
train, test = train_test_split(data_select.copy(), test_size = 0.9)

# 分割训练集和测试集，在训练集上进行参数的训练
train_ = train.copy()
train_y = train_.pop('SalePrice')
train_x = train_
test_ = test.copy()
test_y = test_.pop('SalePrice')
test_x = test_

train_x.head()

lr = LinearRegressionGD(n_iter = 1000)
lr.fit(train_x.values, train_y.values)

print(lr.w_)
# 新的特征都有一个w参数与之对应，w相当于这个地区的基准房价
# lr.w_ 这组数据表达了相关性和重要性，数值越大，说明相关性越高；正号是正相关，负号是反相关
# 从正负号可以读取关键信息；得到不同参数的趋势变化后，就可以用这组参数看准确率的变化

5-3 利用参数看准确率的变化
from sklearn.metrics import mean_squared_error
train['predict'] = lr.predict(train_x.values)

print('MSE train: %.3f' % (mean_squared_error(train_y, train['predict'])))
print(train.head())
# --> MSE train: 0.003
# train.head()得到的结果是在训练集上的准确率；前面的结果是 0.01
# 加入地点区域信息后，准确率得到提高

5-4
# 然后在测试集上查看
test['predict'] = lr.predict(test_x.values)
print('MSE test: %.3f' % (mean_squared_error(test_y, test['predict'])))
# --> MSE test: 0.006
test.head()

标签：案例,test,self,20210516,train,pars,线性,np,data
来源： https://blog.51cto.com/u_15149862/2778659