ML-线性回归

线性回归

对于每一个样本数据(x=(x_{1},x_{2},...,x_{n}))，希望拟合出一个模型 f(x)。当有新的数据输入时，可以给出误差最小的估计值。假设函数如下：

[y = f(x)= heta_{0}+ heta_{1}x_{1}+ heta_{2}x_{2}+...+ heta_{n}x_{n} ]

可以写成向量的形式，其中( heta ={ heta_{0}, heta_{1}, heta_{2}, ..., heta_{n}})，(x^{(i)} ={x_{0},x_{1},x_{2},...,x_{n}}), ( heta x^{(i)})表示向量的内积,即 (y^{(i)}) ，下面的式子表示的是向量的运算

[ heta =( heta_{0}, heta_{1}, heta_{2}, ..., heta_{n})^T ]

[x^{(i)}=(x_{0},x_{1},x_{2},...,x_{n}) ]

[y^{(i)} = x^{(i)} heta ]

最后，对于所有的训练数据，将训练数据表示成矩阵形式 (X_b), ( heta) 是权重矩阵

[X_b = egin{bmatrix} 1 & X_{1}^{(1)} & X_{2}^{(1)} & ... & X_{n}^{(1)} \\ 1 & X_{1}^{(2)} & X_{2}^{(2)} & ... & X_{n}^{(2)} \\ ... & & & ...& \\ 1 & X_{1}^{(m)} & X_{2}^{(m)} & ... & X_{n}^{(m)} end{bmatrix} ]

[ heta =( heta_{0}, heta_{1}, heta_{2}, ..., heta_{n})^T ]

所有的预测结果可以表示为列矩阵(hat{y}),注意此时的(hat{y})是预测的结果，计算误差时还需要与给定的训练集中的结果y比较:

[hat{y} = X_b heta ]

对于数据集 (D={(x_{1}, y_{1}), (x_{2}, y_{2}), ...,(x_{m}, y_{m})}) ，算法希望得到一个模型，是的预测的值(f(x))与(y_{i})的误差最小。预测的误差可以用代价函数来表示,以下式子也称为均方误差：

[J( heta ) = frac{1}{m}sum_{i=1}^{m}(f(x_{i})-y_{i})^{2} ]

[J( heta ) = frac{1}{m}sum_{i=1}^{m}(x^{(i)} heta-y_{i})^{2} ]

注意：这里的(J( heta)) 中只有 ( heta ^ {(i)}) 是变量，变量一共有m个，x都是已知的，对每一个变量求偏导数，更新每一个 ( heta^{(i)})的值。一直到出现了某组 ( heta) 可以使(J( heta))取最小值（梯度下降法）

那么接下来要做的，就是求出这个代价函数的最小值。当代价函数取最小值时 $ heta $的值就是拟合的结果。

交叉验证

将数据随机划分为训练集和测试集。感觉这里写的还不是很好，没有考虑到正反情况的均匀分布。先copy来参考一下：

def train_test_split(X, y, test_ratio=0.2, seed=None):
    """将数据 X 和 y 按照test_ratio分割成X_train, X_test, y_train, y_test"""
    assert X.shape[0] == y.shape[0], 
        "the size of X must be equal to the size of y"
    assert 0.0 <= test_ratio <= 1.0, 
        "test_ration must be valid"

    if seed:
        np.random.seed(seed)

    shuffled_indexes = np.random.permutation(len(X))

    test_size = int(len(X) * test_ratio)
    test_indexes = shuffled_indexes[:test_size]
    train_indexes = shuffled_indexes[test_size:]

    X_train = X[train_indexes]
    y_train = y[train_indexes]

    X_test = X[test_indexes]
    y_test = y[test_indexes]

    return X_train, X_test, y_train, y_test

正规方程法

• 不需要迭代，直接算出 ( heta)
• 当规模n很大时，会影响效率
• 数据不需要作归一化处理

为了求 (J( heta ) = sum_{i=1}^{m}(f(x_{i})-y_{i})^{2})的最小值，对每个 ( heta)分量求导数，并令分量上的导数为0。最后会得到如下结果：

[ heta =(X^{T}X)^{-1}X^{T}y ]

但是，不是所有X都能算出 ((X^{T}X)^{-1}),如果该矩阵不可逆，就会算出多个结果。但是，一般情况下，都是可逆的

使用python实现很方便：

class LinearRegression:
    """具有通用性的多元线性回归的实现"""

    def __init__(self):
        self.coef_ = None       # 表示参数，theta_[1:]
        self.intercept_ = None   # 表示截距 ==>theta[0]
        self._thera = None       # 表示完整的theta==> theta[:]

    def fit_normal(self, x_train, y_train):
        """拟合"""
        # 构造数据集矩阵
        X_b = np.hstack([np.ones((len(x_train), 1)), x_train])
        # 计算系数
        self._thera = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)
        self.coef_ = self._thera[1:]
        self.intercept_ = self._thera[0]
        return self

    def predict(self, x_test):
        X_b = np.hstack(np.ones(len(x_test), 1), x_test)
        return X_b.dot(self._thera)

使用 scikit-learn 解决回归问题

使用sklearn的数据集

    boston = datasets.load_boston()
    X = boston.data     # 特征数据集
    y = boston.target   # 目标变量

去除一些噪音，这种语法只能用在ndarray中

    plt.scatter(np.arange(len(y)), y, color="red")  # 有一些y顶在50的上限上
    X = X[y < 50]
    y = y[y < 50]

使用sklearn的LinearRegression计算参数

    from LinerRegression.model_selection import train_test_split
    from sklearn.linear_model import LinearRegression
    from LinerRegression.metrics import r2_score

    X_train, x_test, y_train, y_test = train_test_split(X, y, seed=666)

    liner = LinearRegression()
    liner.fit(X_train, y_train)
    print("参数：", liner.coef_)
    print("偏置：", liner.intercept_)

    score = r2_score(y_test, liner.predict(x_test))
    print(score)            # R方： 0.8129794056212811

使用前面编写的LinearRegression可以得到一样的答案

    reg = LinearRegression()
    reg.fit(X_train, y_train)
        
    score = r2_score(y_test, liner.predict(x_test)) # 0.8129794056212811
    score2 = r2_score(y_test, reg.predict(x_test))  # 0.8129794056212811

评估

有几个公式可以评价算法的效果

均方误差 MSE(Mean Squared Error)

[frac{1}{m}sum_{i=1}^{m}(hat{y}-y_{i})^{2} ]

# x_test:测试数据的特征，y_test:测试数据的目标变量，y_predict:根据x_test预测出的y
# 就是计算(x_test, y_test) 与 (x_test, y_predict)之间的误差
y_predict = reg.predict(x_test)
np.sum((y_predict - y_test)**2) / len(y_test)

均方根误差 RMSE(Root Mean Squared Error)

[sqrt{frac{1}{m}sum_{i=1}^{m}(hat{y}-y_{i})^{2}} = sqrt {MSE} ]

from math import sqrt
rmse = sqrt(mse_test)

平均绝对误差 MAE(Mean Absolute Error)

MAE不可导，在拟合时不方便求极值,所以没法拿来作代价函数。但是在此处却可以用来衡量算法的效果。因为此处不需要求导，只是做一个衡量的标准

[frac{1}{m}sum_{i=1}^{m}left | hat{y}-y_{i} ight | ]

mae = np.sum(np.absolute(y_predict - y_test))/len(y_test)

R Square(R方)

[R^2 = 1 - frac{SS_{residual}}{SS_{total}} = 1 - frac{sum (hat{y}^{(i)} - y^{(i)})^2}{sum (overline{y} - y^{(i)})^2}= 1 - frac{sum (hat{y}^{(i)} - y^{(i)})^2/m}{sum (overline{y} - y^{(i)})^2/m} = 1-frac{MSE(hat{y}, y)}{Var(y)} ]

• (R^2<=1)
• (R^2)越大越好。当预测模型没有犯错时，得到最大值1
• 当预测模型的效果等于基准模型时，为0
• 如果(R^2<0)，说明学习到的预测模型还不如基准模型。此时很有可能我们的数据不存在线性关系

def r2_score(y_true, y_predict):
    """R方"""
    return 1 - mean_squared_error(y_true, y_predict) / np.var(y_true)

scikit-learn 中的评估函数

• 没有 RMSE，可以用MSE开根号
• score(X, Y) 传递的参数是X和Y，并不是(hat{y}) 和 (y)
• MSE 和 MAE 放在 sklearn.metrics模块下

    from sklearn.metrics import mean_squared_error
    from sklearn.metrics import mean_absolute_error

    # MSE
    mean_squared_error(y_test, liner.predict(x_test)
    # MAE
    mean_absolute_error(y_test, liner.predict(x_test))
    # R方
    liner.score(x_test, y_test)