Working with Linear Models

Fitting a line through data

线性回归简单的讲就是拟合，我们给定一组特征x=[x1,xx,⋯,xn]和目标输出y，求一组参数来组合这组特征使其输出接近目标

y=f(x)=wTx=Σni=1wixi

假设我们由一组观测样本(x(1),y(1)),(x(2),y(2))⋯,(x(m),y(m)),我们使用均方误差(12是为了方便求解)

J(θ)=12mΣmi=1(y(i)−f(x(i)))2

来表示我们模型的好坏，J(θ)越小，说明模型输出越接近真实输出。我们可以直接求解J(θ)的最小值来得到权值w,这时得到解

w=(XTX)−1XTy^

我们使用boston数据进行测试

from sklearn import datasets
boston = datasets.load_boston()

我们导入线性回归模型，构造一个线性回归模型lr，这里没有传入任何参数

from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(boston.data, boston.target)
predictions = lr.predict(boston.data)
c_data = boston.target - predictions

predictions就是我们根据模型得到的输出值，c_data是真实值与期望值的差值，画出来看看是不是准确，(新发现在jupyter cell执行

%pylab inline

可以直接绘图，不需要plt.show(),但是要在绘图函数后加分号

;

)

import matplotlib.pyplot as plt
%pylab inline
plt.hist(c_data,bins=100);

我们可以看到偏差主要分布在0左右，说明我们拟合的数据还是很符合真实数据的可以看一下参数

lr.coef_

线性回归模型中可以通过normalize属性直接将数据预处理为高斯分布（即预处理中的scale）,然后在求解权值，

lr2 = LinearRegression(normalize=True)
lr2.fit(boston.data, boston.target)
predictions2 = lr2.predict(boston.data)
lr2.coef_

但是我在测试的时候发现normalize根本没有用，lr2的权值和lr的权值是一样的。

模型评估

模型评估可以使用mean squared error (MSE) E(y^(i)−y(i))2

或者 meanabsolute deviation (MAD)

E∣∣y^(i)−y(i)∣∣

岭回归防止过拟合

我们直接求解w时，需要XTX可逆，这个时候要求XTX是满秩的，但是在实际的数据中我们经常无法满足这个要求，所以就会造成过拟合.我们创建2000个有3个特征，秩为2的的数据集，采用自助法（bootstrap）得到这三个特征的系数分布

from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=2000,n_features=3, effective_rank=2, noise=10)
import numpy as np
n_bootstraps = 1000
len_data = len(reg_data)
subsample_size = np.int(0.75*len_data)
subsample = lambda: np.random.choice(np.arange(0, len_data),size=subsample_size)
coefs = np.ones((n_bootstraps, 3))
for i in range(n_bootstraps):
subsample_idx = subsample()
subsample_X = reg_data[subsample_idx]
subsample_y = reg_target[subsample_idx]
lr.fit(subsample_X, subsample_y)
coefs[i][0] = lr.coef_[0]
coefs[i][2] = lr.coef_[1]
coefs[i][2] = lr.coef_[2]
plt.hist(coefs[:,0]);
plt.hist(coefs[:,1]);
plt.hist(coefs[:,2]);

可以看到这三个权值的分布非常的广，换一组数据（同一分布）得到的权值差别很大，这说明我们的模型没有学到真实数据集的分布。我们可以通过正则化来解决这个问题，即加上权值的L2范数

J(θ)=12m(Σmi=1(y(i)−f(x(i)))2+λΣni=1w2i)

此时得到的方程的解为

w=(XTX+λI)−1XTy^

加入一个λI是的(XTX)趋向于可逆。我们使用Rigde来测试

from sklearn.linear_model import Ridge
r = Ridge()
for i in range(n_bootstraps):
subsample_idx = subsample()
subsample_X = reg_data[subsample_idx]
subsample_y = reg_target[subsample_idx]
r.fit(subsample_X, subsample_y)
coefs[i][0] = r.coef_[0]
coefs[i][3] = r.coef_[1]
coefs[i][2] = r.coef_[2]
plt.hist(coefs[:,0],bins=100);
plt.hist(coefs[:,1],bins=100);
plt.hist(coefs[:,2],bins=100);

可以看到权值分布范围变小了，权值的方差减小了，这说明即使更换为其他的数据（服从同一分布）,我们得到的权值变化不大，说明此时得到的模型非常的稳定。

选择岭回归的参数

上面知道解中添加了λI,其中λ是一个超参数，到底为多少的时候模型才是最好的呢？我们可以使用交叉验证来选择。（这里说是交叉验证，但是我并不确定，因为没有看到交叉验证的代码）我们使用RidgeCV来进行测试，这是标准的ridge cross-validation。

from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=100,n_features=2, effective_rank=1, noise=10)
from sklearn.linear_model import RidgeCV
rcv = RidgeCV(alphas=np.array([.1, .2, .3, .4]))
rcv.fit(reg_data, reg_target)
rcv.alpha_

根据

rc.alpha_

可以知道哪个值最好。我们可以进一步缩小α的范围，求的更精确的值

rcv2 = RidgeCV(alphas=np.array([.08, .09, .1, .11, .12]))
rcv2.fit(reg_data, reg_target)
rcv2.alpha_

在

RidgeCV

中有一个参数

store_cv_values=True

store_cv_values : boolean, default=False

Flag indicating if the cross-validation values corresponding to each alpha should be stored in the cv_values_ attribute (see below). This flag is only compatible with cv=None (i.e. using Generalized Cross-Validation).

store_cv_values

为

True

时，会返回

cv_values_

的值。

cv_values_ : array, shape = [n_samples, n_alphas] or shape = [n_samples, n_targets, n_alphas], optional

Cross-validation values for each alpha (if store_cv_values=True and cv=None). After fit() has been called, this attribute will contain the mean squared errors (by default) or the values of the {loss,score}_func function (if provided in the constructor).

cv_values_

包含每一个样本对应每一个α的真实输出与模型输出的差值平方，所有样本的平均就是总体的均方误差损失函数。这样我们可以由

cv_values_

画出损失函数随α如何变化。

alphas_to_test = np.linspace(0.01, 1)
rcv3 = RidgeCV(alphas=alphas_to_test, store_cv_values=True)
rcv3.fit(reg_data, reg_target)
mse_loss = rcv3.cv_values_.mean(axis=0)
smallest_idx = mse_loss.argmin()
alphas_to_test[smallest_idx]
rcv3.alpha_
plt.plot(alphas_to_test,mse_loss)

我们可以看到输出的损失函数如何随α变化，输出的最优α与曲线也是对应的。默认的是误差平方最为评分函数

score(x)=(y−f(x))2

我们也可以使用自己的评分函数,例如误差的绝对值

score(x)=|y−f(x)|

def MAD(target, predictions):
absolute_deviation = np.abs(target - predictions)
return absolute_deviation.mean()
import sklearn
MAD = sklearn.metrics.make_scorer(MAD, greater_is_better=False)
rcv4 = RidgeCV(alphas=alphas_to_test, store_cv_values=True,scoring=MAD)
rcv4.fit(reg_data, reg_target)
mad_loss = rcv4.cv_values_.mean(axis=0)
smallest_idx = mad_loss.argmin()
alphas_to_test[smallest_idx]
plt.plot(alphas_to_test,mad_loss)

可以看到评分函数改变后，α的最优值也会改变。

LASSO

有时候我们的特征很多，使得权值也很多，如果特征前的系数如果为0，那么相当于去掉了一些不重要的特征。这样的好处是可以通过减少特征防止过拟合，而且可以使计算量减小。LASSO和岭回归减轻过拟合有什么不一样呢？简单的说岭回归会惩罚那些过大的系数，使那些不重要的特征的系数尽可能接近零，但是不会为零，LASSO可以使系数变为零，这样就相当于去掉了一些特征，所以LASSO具有使数据变稀疏的特性。为了防止过拟合，岭回归加上了为λΣni=1w2i的正则项，而LASSO则是通过添加λΣni=1|w|

我们用

Lasso

测试一下

我们创建一个有500个特征，相关特征有5个的200个样本

from sklearn.datasets import make_regression
reg_data, reg_target = make_regression(n_samples=200, n_features=500,n_informative=5, noise=5)
from sklearn.linear_model import Lasso
lasso = Lasso()
lasso.fit(reg_data, reg_target)
np.sum(lasso.coef_ != 0)

lasso_0 = Lasso(0)
lasso_0.fit(reg_data, reg_target)
np.sum(lasso_0.coef_ != 0)

使用lasso后，非零系数只有几个，但是如果不使用lasso，那么特征一个都不会减少，还是500个。可以看到lasso的正则项也是有一个超参数λ需要我们调节，可以使用lasso的交叉验证来获得最优的λ。

from sklearn.linear_model import LassoCV
lassocv = LassoCV()
lassocv.fit(reg_data, reg_target)
lassocv.alpha_

lasso使很多特征的系数为零，我们需要获取这些经过选择后的特征

mask = lassocv.coef_ != 0
new_reg_data = reg_data[:, mask]

LARS

不懂跳过

Logistic回归（二项分类）

线性回归不仅可以用来做拟合，还可以用来做分类，我们把拟合的模型通过某个函数映射为正类或负类的概率，大于阈值就分为正类，否则负类。这里的映射函数使Logistic

f(t)=11+e−t

我们令t=wTx那么

y=11+e−wTx

假设

P(y=1|x)=11+e−wTx

P(y=0|x)=1−P(y=1|x)=e−wTx1+e−wTx

logit(p)=log(P(y=1|x)P(y=0|x))=wTx

可以看到我们线性回归使正类的概率与负类概率的对数比值。

我们来测试一下Logistic分类

from sklearn.datasets import make_classification
X, y = make_classification(n_samples=1000, n_features=4)
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression()
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.25)
lr.fit(X_train, y_train)
y_train_predictions = lr.predict(X_train)
y_test_predictions = lr.predict(X_test)
(y_train_predictions == y_train).sum().astype(float) / y_train.shape[0]
(y_test_predictions == y_test).sum().astype(float) / y_test.shape[0]

还有一个问题，如果我们的样本数据，正类占了95%，负类只占5%，那么我们根本就不需要训练分类器，只要直接让分类器返回正类就行，即使这样也有95%的准确率，这个问题是样本分布不均时该如何处理，我感觉放在这里很怪异，就不深究了。

贝叶斯回归

我们计算w时是直接计算的点估计的，但是贝叶斯回归是计算是w的分布，所以更加准确。

from sklearn.datasets import make_regression
X, y = make_regression(1000, 10, n_informative=2, noise=20)
from sklearn.linear_model import BayesianRidge
br = BayesianRidge()
br.fit(X, y)
br.coef_

贝叶斯不是很懂，就不BB了

boosting

boosting起源

Boosting是一种提高任意给定学习算法准确度的方法。它的思想起源于 Valiant提出的 PAC ( Probably Approximately Correct)学习模型。

Valiant和 Kearns提出了弱学习和强学习的概念,识别错误率小于1/2,也即准确率仅比随机猜测略高的学习算法称为弱学习算法;识别准确率很高并能在多项式时间内完成的学习算法称为强学习算法。

同时 ,Valiant和 Kearns首次提出了PAC学习模型中弱学习算法和强学习算法的等价性问题,即任意给定仅比随机猜测略好的弱学习算法,是否可以将其提升为强学习算法?如果二者等价,那么只需找到一个比随机猜测略好的弱学习算法就可以将其提升为强学习算法,而不必寻找很难获得的强学习算法。

1990年, Schapire最先构造出一种多项式级的算法 ,对该问题做了肯定的证明 ,这就是最初的 Boosting算法。一年后 ,Freund提出了一种效率更高的Boosting算法。但是,这两种算法存 共同的实践上的缺陷 ,那就是都要求事先知道弱学习算法学习正确的下限。1995年 ,Freund和 schap ire改进了Boosting算法 ,提出了 AdaBoost (Adap tive Boosting)算法,该算法效率和 Freund于1991年提出的Boosting算法几乎相同,但不需要任何关于弱学习器的先验知识,因而更容易应用到实际问题当中。之后,Freund和 schapire进一步提出了改变Boosting投票权重的 AdaBoost . M1,AdaBoost . M2等算法 ,在机器学习领域受到了极大的关注。

我们这里使用的gradient boosting方法来测试

from sklearn.datasets import make_regression
X, y = make_regression(1000, 2, noise=10)
from sklearn.ensemble import GradientBoostingRegressor as GBR
gbr = GBR()
gbr.fit(X, y)
gbr_preds = gbr.predict(X)
c_data = y - gbr_preds
plt.hist(c_data,bins=100);