线性回归

线性模型

目标值是特征的线性组合就称为线性回归。

普通的最小二乘法

寻找$w=(w_1,w_2,...,w_p)$使得$mi{n}_w||Xw-y|{|}_2^2$。

from sklearn import linear_model
reg = linear_model.LinearRegression()
reg.fit([[0, 0], [1, 1], [2, 2]], [0, 1, 2])
##fit函数中的参数是X和y，训练得到的系数存储在coef_中
print(reg.coef_)

当矩阵X中的特征相关并且X中的列有线性关系时，得到的系数极大地依赖于观测值得随机误差，可能导致较大的方差。

import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model
from sklearn.metrics import mean_squared_error, r2_score

diabetes_X, diabetes_y = datasets.load_diabetes(return_X_y=True)
print(diabetes_X)
##442*10
##加载糖尿病数据集，并且根据参数return_X_y设置返回X和y的二元组
diabetes_X = diabetes_X[:, np.newaxis, 2]
##仅使用一个特征
##np.newaxis是None的别名
##这里相当于将数据集转换成442*1*10之后再第三个维度的10个系列中选择第三个系列,得到的还是一个二维的数组
##因为fit函数得到的X需要是二维的，所以做相关处理
diabetes_X_train = diabetes_X[:-20]
diabetes_X_test = diabetes_X[-20:]
diabetes_y_train = diabetes_y[:-20]
diabetes_y_test = diabetes_y[-20:]
##将数据分成训练集和验证集
regr = linear_model.LinearRegression()
regr.fit(diabetes_X_train, diabetes_y_train)
diabetes_y_pred = regr.predict(diabetes_X_test)
plt.scatter(diabetes_X_test, diabetes_y_test, color='black')
plt.plot(diabetes_X_test, diabetes_y_pred, color='blue', linewidth=3)

plt.xticks(())##不显示坐标刻度
plt.yticks(())
plt.show()

假设X的大小是$(n_{s}amples,n_{f}eatures)$则最小二乘法的时间复杂度是$O(n_{samples}{n}_{features}^2)$

正则化的代价函数

$$mi{n}_w||Xw-y|{|}_2^2+\alpha ||w|{|}_2^2$$
其中的系数$\alpha$控制着正则向和残差平方项的比重。用法如下：

from sklearn import linear_model

reg = linear_model.Ridge(alpha=.5)
reg.fit([[0, 0], [0, 0], [1, 1]], [0, .1, 1])
print(reg.coef_) ##存储的系数是w1,...,
print(reg.intercept_)##存储w0

Lasso讲解

Lasso是一个估计稀疏系数的线性模型。它在某些情况下是有用的，因为它倾向于选择具有较少非零系数的解，有效地减少了给定解所依赖的特征的数量。因此，Lasso和它的变体是压缩感知领域的基础。
最小化目标函数是：$$mi{n}_w\frac{1}{2n_{samples}}||Xw-y|{|}_2^2+\alpha ||w|{|}_1$$
即Lasso回归是在损失函数中加上L1正则化，从而可以选择系数绝对值比较小的作为解。Lasso的复杂程度由$\alpha$控制，$\alpha$越大，对变量越多的线性模型惩罚力度就越大。
Lasso方法可以达到变量选择的效果，将不显著的变量系数压缩为0。Ridge方法虽然也对系数进行压缩，但是任一系数都不会被压缩为0，最终模型保留了所有的变量。

Lasso参数求解

因为Lasso的代价函数中的正则化项使用的是绝对值的和，所以损失函数在一些点是不可导的。因为梯度下降法等优化算法对Lasso失效。但是可以采用坐标轴下降法和最小角回归法。博客传送门

设置参数$\alpha$

scikit-learn通过交叉验证暴露设置Lasso alpha参数的对象:LassoCV和LassoLarsCV。LassoLarsCV基于下面介绍的最小角度回归算法。
LassoCV:适用于具有共线特征的高维数据。
LassoLarsCV:的优势在于探索更多的$\alpha$参数相关值，当样本数量与特征数量相比非常小的时候，LassoLarsCV往往比LassoCV更快。

弹性网络

优化目标：

LARS Lasso

是分段线性拟合
Examples：

正交匹配追踪OMP

贝叶斯回归

贝叶斯回归技术可以用于在估计过程中包含正则化参数:正则化参数不是硬性设定的，而是根据手头的数据进行调整。
优点：
可以适应现有的数据
可以用于在估计过程中包含正则化参数
缺点:
模型的推断可能会花费大量的时间

贝叶斯脊回归

贝叶斯脊回归对不适性问题具有更强的鲁棒性。
用法：

from sklearn import linear_model

X = [[0., 0.], [1., 1.], [2., 2.], [3., 3.]]
Y = [0., 1., 2., 3.]
reg = linear_model.BayesianRidge()
reg.fit(X, Y)
reg.predict([[1, 0.]])

#估计拉索和弹性网回归模型的一个人工产生的稀疏信号与加性噪声的破坏。估计系数与真实值进行比较。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score

##首先生成一些待处理的数据
np.random.seed(42)
n_samples, n_features = 50, 100
X = np.random.randn(n_samples, n_features)

idx = np.arange(n_features)
coef = (-1) ** idx * np.exp(-idx / 10) ##系数
coef[10:] = 0
y = np.dot(X, coef) ##根据X和coef生成函数值
y += 0.1 * np.random.normal(size=n_samples) ##为函数值加入噪声

n_samples = X.shape[0] ##获得样本数量,其实已经知道了，再上面有显示的赋值
##将数据分成训练集和测试集
X_train, y_train = X[:n_samples // 2], y[:n_samples // 2]
X_test, y_test = X[n_samples // 2:], y[n_samples // 2:]

from sklearn.linear_model import Lasso

alpha = 0.1
lasso = Lasso(alpha=alpha)
y_pred_lasso = lasso.fit(X_train, y_train).predict(X_test)
r2_score_lasso = r2_score(y_test, y_pred_lasso)
## r2回归得分函数，分数总是不大于1的，可以衡量两个数组是否足够接近

from sklearn.linear_model import ElasticNet
##弹性网络
enet = ElasticNet(alpha=alpha, l1_ratio=0.1)
y_pred_enet = enet.fit(X_train, y_train).predict(X_test)
r2_score_enet = r2_score(y_test, y_pred_enet)

m, s, _ = plt.stem(np.where(enet.coef_)[0], enet.coef_[enet.coef_ != 0],
                   markerfmt='x', label='Elastic net coefficients',
                   use_line_collection=True)
plt.setp([m, s], color="#2ca02c")
m, s, _ = plt.stem(np.where(lasso.coef_)[0], lasso.coef_[lasso.coef_ != 0],
                   markerfmt='x', label='Lasso coefficients',
                   use_line_collection=True)
plt.setp([m, s], color='#ff7f0e')
##plt.stem函数画出的图像类似从一个点想x轴做垂线
plt.stem(np.where(coef)[0], coef[coef != 0], label='true coefficients',
         markerfmt='bx', use_line_collection=True)

plt.legend(loc='best')
plt.title("Lasso $R^2$: %.3f, Elastic Net $R^2$: %.3f"
          % (r2_score_lasso, r2_score_enet))
plt.show()

##解释线性模型系数的常见缺陷
##我们将使用1985年“当前人口调查”的数据来预测工资作为各种特征的函数，如经验、年龄或教育。

import numpy as np
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.datasets import fetch_openml

survey = fetch_openml(data_id=534, as_frame=True)
##as_frame参数可以令函数返回一个pandas dataframe
X = survey.data[survey.feature_names]
X.describe(include='all')
##数据集包含分类变量和数值变量

y = survey.target.values.ravel() #得到预测的目标值，用美元/小时表示

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
##将数据集分成训练集和验证集，之后的训练只在训练集上进行

train_dataset = X_train.copy()
train_dataset.insert(0, 'WAGE', y_train) ##将训练数据集的样本和标签一起放到train_dataset中
sns.pairplot(train_dataset, kind='reg', diag_kind='kde')
##pairplot在数据集中绘制成对关系，需要data是pandas dataframe
plt.show() ##不加这句，图像不会显示

#分类变量不变成数值变量是不可以出现在线性模型中的
# to avoid categorical features to be treated as ordered values, we need to one-hot-encode them
'''
预处理步骤:
    将每个分类列进行one-hot encode
    数值列保持原来的不变
'''
from sklearn.compose import make_column_transformer
from sklearn.preprocessing import OneHotEncoder

categorival_columns = ['RACE', 'OCCUPATION', 'SECTOR',
                       'MARR', 'UNION', 'SEX', 'SOUTH']
numerical_columns = ['EDUCATION', 'EXPERIENCE', 'AGE']
preprocessor = make_column_transformer(
    (OneHotEncoder(drop='if_binary'), categorival_columns),
     remainder='passthrough'
)
#sklearn.compose.make_column_transformer
#OneHotEncoder将一个分类特征编码为one-hot数字数组
#‘if_binary’ : drop the first category in each feature with two categories. Features with 1 or more than 2 categories
# are left intact.

from sklearn.pipeline import make_pipeline
from sklearn.linear_model import Ridge
from sklearn.compose import TransformedTargetRegressor

model = make_pipeline(
    preprocessor,
    TransformedTargetRegressor(
        regressor=Ridge(alpha=1e-10),
        func=np.log10,
        inverse_func=sp.special.exp10
    )
) ##构造管道
model.fit(X_train, y_train)##训练模型

from sklearn.metrics import median_absolute_error

y_pred = model.predict(X_train) ##在训练集上进行预测
mae = median_absolute_error(y_train, y_pred) #Median absolute error regression loss
string_score = f'MAE on training set:{mae:.2f} $/hour'
y_pred = model.predict(X_test)
mae = median_absolute_error(y_test, y_pred)
string_score += f'\nMAE on testing set: {mae:.2f} $/hour'
fig, ax = plt.subplots(figsize=(5, 5))
plt.scatter(y_test, y_pred)
ax.plot([0, 1], [0, 1], transform=ax.transAxes, ls="--", c="red")
plt.text(3, 20, string_score)
plt.title('Ridge model, small regularization')
plt.ylabel('Model predictions')
plt.xlabel('Truths')
plt.xlim([0, 27])
plt.ylim([0, 27])
plt.show()

##解释系数，规模很重要
feature_names = (model.named_steps['columntransformer']
                      .named_transformers_['onehotencoder']
                      .get_feature_names(input_features=categorival_columns))
feature_names = np.concatenate(
    [feature_names, numerical_columns])

coefs = pd.DataFrame(
    model.named_steps['transformedtargetregressor'].regressor_.coef_,
    columns=['Coefficients'], index=feature_names
)
X_train_preprocessed = pd.DataFrame(
    model.named_steps['columntransformer'].transform(X_train),
    columns=feature_names
)

X_train_preprocessed.std(axis=0).plot(kind='barh', figsize=(9, 7))
plt.title('Features std. dev.')
plt.subplots_adjust(left=.3)
plt.show()

广义线性回归GLM

广义线性模型在两方面扩展了线性模型：
首先，预测值$\hat{y}$可以表示成下面的形式：$\hat{y}(w,X)=h(Xw)$
其次，将平方误差函数替换成指数族分布的单位偏差。即最小化问题变成了：$$mi{n}_w\frac{1}{2n_{samples}}{\Sigma }_{i}d(y_i,{\hat{y}}_i)+\frac{\alpha }{2}||w|{|}_2$$

多项式回归:用基函数扩展线性模型

在机器学习中一个常见的模式是使用在数据的非线性函数上训练的线性模型。这种方法保持了线性方法的一般快速性能，同时允许它们拟合更大范围的数据。

一篇讲解贝叶斯估计和最大似然估计的博客