scikit-learn 学习笔记（一）

其实前一阵一直在看《机器学习系统设计》，但是发现书上代码不全，不好复现，而且书针对的主要是文字的处理，正好是我最不关心的方向，所以看到一半忍痛弃坑。那么我们开始scikit-learn文档学习的旅程吧。

1. 线性模型

我对于线性模型的理解就是，将一系列的已知点回归到一个线性方程上去，类似于：。之后就可以用这个方程来预测未知点的解了。 sklearn 规定 w0是 intercept_. w1--wp是coef_.

1.1普通最小二乘法

那么如何将已知的点拟合到一个线性方程上去呢？sklearn提供的LinearRegression的思路就是让点的已知值和预测值的平方和最小，即. 下面我们直接上例子：

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

diabetes = datasets.load_diabetes()
diabetes_X = diabetes.data[:,np.newaxis,2]
diabetes_X_train = diabetes_X[:-20]
diabetes_X_test = diabetes_X[-20:]
diabetes_y_train = diabetes.target[:-20]
diabetes_y_test = diabetes.target[-20:]

regr = linear_model.LinearRegression()
regr.fit(diabetes_X_train, diabetes_y_train) #这里就是在训练模型了

print('Coefficients: \n', regr.coef_) #这就是w0，常数项
print("Residual sum of squares: %.2f" % np.mean((regr.predict(diabetes_X_test) - diabetes_y_test) ** 2)) #这个是预测与真实的差
print('Variance score: %.2f' % regr.score(diabetes_X_test, diabetes_y_test)) #这里就是得分，1为拟合最好，0最差

plt.scatter(diabetes_X_test, diabetes_y_test, color = 'black')
plt.plot(diabetes_X_test,regr.predict(diabetes_X_test), color='blue',linewidth=3)
plt.xticks(())
plt.yticks(())

plt.show()

上图：

这个得分是0.47，果然拟合的不太好呢。

1.2 Ridge regression 岭回归

作为最小二乘法的一个改进，岭回归使用：
其中a>=0，w是一个动态调整的惩罚参数。 我们先上例子：

print(__doc__)

import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model
#构造一个hilbert矩阵10*10
X = 1. / (np.arange(1,11) + np.arange(0,10)[:,np.newaxis])
y = np.ones(10) #10个1

n_alphas = 200
alphas = np.logspace(-10,-2,n_alphas) #等比数列，200个，作为岭回归的系数
clf = linear_model.Ridge(fit_intercept=False)
coefs = []
for a in alphas:
clf.set_params(alpha=a)
clf.fit(X,y)
coefs.append(clf.coef_) #得到200个不同系数所训练出常数参数值

ax = plt.gca() 
ax.set_color_cycle(['b', 'r', 'g', 'c', 'k', 'y', 'm'])
ax.plot(alphas,coefs)
ax.set_xscale('log') #转为极坐标系
ax.set_xlim(ax.get_xlim()[::-1]) #反转x轴
plt.xlabel('alpha')
plt.ylabel('weight')
plt.title('Ridge coefficients as a function of the regularization')
plt.axis('tight')
plt.show()

上图片：

发现了什么？没错，传入岭回归的系数越小，常数项参数的值就越大，惩罚越重。 而且我们发现，在alpha小于10-5时，系数的变化明显变大，所以我们应该选取的值就在10-5附近。 那么这个值怎么自动计算出来呢？ 上代码：

clf = linear_model.RidgeCV(alphas = alphas)
clf.fit(X,y)
clf.alpha_

这些代码的作用是输入一系列的alphas，然后sklearn会自动的内部cross验证最终得到最适合的alpha。

1.3 lasso(套索)

她的应用是基于稀疏模型的情况，进行线性拟合，这时的效果较好。那么什么是稀疏模型呢？比如我们采集许多的温度数据，肯定会存在大量的冗余数据，因为温度的变化是缓慢的，稀疏模型就会关注变化的点，而去除同质化的点，从而减少运算量。这个方法其实是压缩感知的基础。不知道压缩感知，这里有一个好的普及贴。 她要最小化的是：

为了便于理解，我们直接上例子：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score
#我们先手动生成一些稀疏数据
print np.random.seed(42)
n_samples, n_features = 50, 200 
X = np.random.randn(n_samples, n_features)  <span style="font-family: Arial, Helvetica, sans-serif;">#50行200列</span>
coef = 3 * np.random.randn(n_features) #这个就是实际的参数
inds = np.arange(n_features)
np.random.shuffle(inds) #打乱
coef[inds[10:]] = 0 #生成稀疏数据
y = np.dot(X, coef) #参数与本地点乘
#来点噪音
y += 0.01 * np.random.normal((n_samples,))

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) #这里是0.38
print lasso
print "r2_score's result is %f" % r2_score_lasso

from sklearn.linear_model import ElasticNet
enet = ElasticNet(alpha=alpha, l1_ratio=0.7)
y_pred_enet = enet.fit(X_train,y_train).predict(X_test)
r2_score_enet = r2_score(y_test, y_pred_enet) #0.24 没有lasso好
print enet
print "nent's result is %f" % r2_score_enet

plt.plot(enet.coef_, label='Elastic net coefficients')
plt.plot(lasso.coef_, label='Lasso coefficients')
plt.plot(coef, '--', label='original coefficients')
plt.legend(loc="best")
plt.title("Lasso R^2: %f, Elastic Net R^2: %f"
          % (r2_score_lasso, r2_score_enet))
plt.show()

然后上图：
这里lasso要比弹性网络表现好一些。 这里我们alpha使用的是0.1，但是这个值好吗？如何自动计算最优值，这里有几种不同的标准，应对不同范围：LassoCV,LassoLarsCV, LassoLarsIC 先上代码：

def LassoCV():
from sklearn.linear_model import LassoCV,LassoLarsCV, LassoLarsIC

model_lasso = LassoCV(cv=20).fit(X_train,y_train)
print "the alpha value of LassoCV is %f" % model_lasso.alpha_

model_lassoLars = LassoLarsCV(cv=20).fit(X_train,y_train)
print "the alpha value of LassoLarsCV is %f" % model_lassoLars.alpha_

model_lassoLarsIC_bic = LassoLarsIC(criterion='bic').fit(X_train,y_train)
print "the alpha value of LassoLarsIC's bic is %f" % model_lassoLarsIC_bic.alpha_

model_lassoLarsIC_aic = LassoLarsIC(criterion='aic').fit(X_train,y_train)
print "the alpha value of LassoLarsIC's aic is %f" % model_lassoLarsIC_aic.alpha_

然后执行结果就是：

the alpha value of LassoCV is 1.947789
the alpha value of LassoLarsCV is 0.435168
the alpha value of LassoLarsIC's bic is 1.120660
the alpha value of LassoLarsIC's aic is 1.120660

根据不同的数据情况选择合适的alpha值吧。

1.4 Elastic Net（弹性网络）

这个最初是用来做基于位置关系的社交网络的，即会根据两个用户的地理位置动态改变社交关系网络。 与Lasso类似，他也是用在稀疏模型的基础上的。但是Lasso在选取一组相关变量时是随机选取一个，而她可以选取所有的。她的最小化目标是：
(话说根据观察值求线性方程参数的过程不就是稀疏表示的过程吗。。。) 这个例子已经在1.3中。

1.5多任务lasso

实际上就是多元lasso，一元的lasso不是根据y=X*x（y是结果，X是矩阵，x是方程系数）来倒推x吗？多任务lasso中的y从一维度变成二维度。 直接上例子：

import matplotlib.pyplot as plt
import numpy as np

from sklearn.linear_model import MultiTaskLasso, Lasso
rng = np.random.RandomState(42)
n_samples, n_features, n_tasks = 100, 30, 40
n_relevant_features = 5
coef = np.zeros((n_tasks, n_features))  <span style="font-family: Arial, Helvetica, sans-serif;">#原始数据是40*30的矩阵，但是只有5列有数据</span>
times = np.linspace(0, 2 * np.pi, n_tasks) 

for k in range(n_relevant_features):
coef[:,k] = np.sin((1. + rng.randn(1)) * times + 3 * rng.randn(1))
X = rng.randn(n_samples,n_features) #字典矩阵是100*30
Y = np.dot(X, coef.T) + np.dot(n_samples,n_tasks) #那么实际表现值当然是100*40

coef_lasso_ = np.array([Lasso(alpha=0.5).fit(X,y).coef_ for y in Y.T]) #使用lasso逐行计算
coef_multi_task_lasso_ = MultiTaskLasso(alpha=1.).fit(X,Y).coef_ #使用多任务lasso一次性计算

fig = plt.figure(figsize=(8,5))
plt.subplot(1, 3, 1)
plt.spy(coef, precision=0.1, markersize=5)
plt.xlabel('Feature')
plt.ylabel('Time (or Task)')
plt.text(10, 5, 'origin')
plt.subplot(1, 3, 2)
plt.spy(coef_lasso_, precision=0.1,markersize=5)
plt.xlabel('Feature')
plt.ylabel('Time (or Task)')
plt.text(10, 5, 'Lasso')
plt.subplot(1, 3, 3)
plt.spy(coef_multi_task_lasso_, precision=0.1, markersize=5)
plt.xlabel('Feature')
plt.ylabel('Time (or Task)')
plt.text(10, 5, 'MultiTaskLasso')
fig.suptitle('Coefficient non-zero location')

feature_to_plot = 3 #0-4都是有数据的，因为是5列嘛
plt.figure()
plt.plot(coef[:, feature_to_plot], 'k', label='Ground truth')
plt.plot(coef_lasso_[:,feature_to_plot], 'g', label='Lasso')
plt.plot(coef_multi_task_lasso_[:,feature_to_plot], 'r', label='MultiTaskLasso')
plt.legend(loc='upper center')
plt.axis('tight')
plt.ylim([-1.1, 1.1])
plt.show()

print coef.shape, coef_lasso_.shape, coef_multi_task_lasso_.shape

上图：


可以看到多任务lasso在处理多维的任务时，表现会更好。

1.6最小角回归

一个针对高纬度数据的回归。优点包括： 1. 当维度远远大于数据点的个数时，非常有效。 2. 拥有向前选择法的速度，同时具有普通最小二乘法的复杂度。 3. 更符合人的直觉，比如两个相似的数据会有两个相似的结果。 缺点： 因为他是基于迭代的，所以对噪声特别敏感。

1.7 最小角回归Lasso

就是将1.6的方法运用到lasso上。 上代码：

import numpy as np
import matplotlib.pyplot as plt

from sklearn import linear_model
from sklearn import datasets

diabetes = datasets.load_diabetes()
X = diabetes.data
y = diabetes.target

print("Computing regularization path using the LARS ...")
alphas, _, coefs = linear_model.lars_path(X,y, method='Lasso',verbose=True) #这里coef会返回一个10*11的矩阵

#10代表多项式有10个参数，11代表每个参数都从小到大排列，一共有11个# print np.abs(coefs.T)xx = np.sum(np.abs(coefs.T), axis=1) #将每个级别的所有参数相加，所以有11个，从小到大xx /= xx[-1] #这里将每个级别的和都除以最大值plt.plot(xx, coefs.T)ymin,ymax = plt.ylim()plt.vlines(xx,ymin,ymax,linestyle='dashed')plt.xlabel('|coef| / max|coef|')plt.ylabel('Coefficients')plt.title('LASSO Path')plt.axis('tight')plt.show()

图是：

1.8 正交匹配追踪（OMP）

首先介绍匹配追踪（MP）：从字典矩阵D（也称为过完备原子库中），选择一个与信号 y 最匹配的原子(也就是某列)，构建一个稀疏逼近，并求出信号残差，然后继续选择与信号残差最匹配的原子，反复迭代，信号y可以由这些原子来线性和，再加上最后的残差值来表示。很显然，如果残差值在可以忽略的范围内，则信号y就是这些原子的线性组合。OMP算法的改进之处在于：在分解的每一步对所选择的全部原子进行正交化处理，这使得在精度要求相同的情况下，OMP算法的收敛速度更快。 上例子：

import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import OrthogonalMatchingPursuit
from sklearn.linear_model import OrthogonalMatchingPursuitCV
from sklearn.datasets import make_sparse_coded_signal

n_components, n_features = 512, 100
n_nonzero_coefs = 17

y, X, w = make_sparse_coded_signal(n_samples=1,
n_components=n_components,
n_features=n_features,
n_nonzero_coefs=n_nonzero_coefs,
random_state=0)
idx, = w.nonzero()

y_noise = y + 0.05 * np.random.randn(len(y))

plt.figure(figsize=(7,7))
plt.subplot(4,1,1)
plt.xlim(0,512)
plt.title("orginal Sparse signal")
plt.stem(idx, w[idx])

opm = OrthogonalMatchingPursuit(n_nonzero_coefs=n_nonzero_coefs)
opm.fit(X,y)
coef = opm.coef_
idx_r, = coef.nonzero()
plt.subplot(4, 1, 2)
plt.xlim(0, 512)
plt.title("Recovered signal from noise-free measurements")
plt.stem(idx_r, coef[idx_r])

opm.fit(X, y_noise)
coef = opm.coef_
idx_r, = coef.nonzero()
plt.subplot(4, 1, 3)
plt.xlim(0, 512)
plt.title("Recovered signal from noisy measurements")
plt.stem(idx_r, coef[idx_r])

opm_cv = OrthogonalMatchingPursuitCV()
opm_cv.fit(X,y_noise)
coef = opm_cv.coef_
idx_r, = coef.nonzero()
plt.subplot(4, 1, 4)
plt.xlim(0, 512)
plt.title("Recovered signal from noisy measurements with CV")
plt.stem(idx_r, coef[idx_r])

plt.subplots_adjust(0.06, 0.04, 0.94, 0.90, 0.20, 0.38)
plt.suptitle('Sparse signal recovery with Orthogonal Matching Pursuit',
             fontsize=16)

plt.show()

上图：
可以看到，在没有噪音时，正交匹配追踪的表现还是不错的，但是加上噪音后，CV的表现会更好一些。

1.9 贝叶斯岭回归

将岭回归与贝叶斯结合的一种回归方法。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from sklearn.linear_model import BayesianRidge, LinearRegression

n_samples, n_features = 100, 100
X = np.random.randn(n_samples, n_features)  # Create Gaussian data
# Create weigts with a precision lambda_ of 4.
lambda_ = 4.
w = np.zeros(n_features)
# Only keep 10 weights of interest
relevant_features = np.random.randint(0, n_features, 10) #找10个随机位置
for i in relevant_features:
    w[i] = stats.norm.rvs(loc=0, scale=1. / np.sqrt(lambda_)) #制造10个非0的特征
# Create noise with a precision alpha of 50.
alpha_ = 50.
noise = stats.norm.rvs(loc=0, scale=1. / np.sqrt(alpha_), size=n_samples)
y = np.dot(X, w) + noise #加噪音

clf = BayesianRidge(compute_score=True)
clf.fit(X,y)

ols = LinearRegression()
ols.fit(X,y)

plt.figure()
plt.subplot(3,1,1)
plt.title("Weights of the model")
plt.plot(clf.coef_, 'b-', label="Bayesian Ridge estimate")
plt.plot(w, 'g-', label="original")
plt.plot(ols.coef_, 'r--', label="OLS estimate")
plt.xlabel("Features")
plt.ylabel("Values of the weights")
plt.legend(loc='best', prop=dict(size=12))

plt.subplot(3,1,2)
plt.title("Histogram of the weights")
plt.hist(clf.coef_, bins=n_features, log=True)
plt.plot(clf.coef_[relevant_features], 5 * np.ones(len(relevant_features)), 'ro', label="Relevant features")
plt.ylabel("Features")
plt.xlabel("Values of the weights")
plt.legend(loc="lower left")

plt.subplot(3,1,3)
plt.title("Marginal log-likelihood")
plt.plot(clf.scores_)
plt.ylabel("Score")
plt.xlabel("Iterations")
plt.show()

图：

与最小二乘法相比，贝叶斯岭回归的稳定性会更好。在6次迭代后趋于稳定。

1.10 逻辑回归

Logistic回归与多重线性回归实际上有很多相同之处，最大的区别就在于它们的因变量不同，其他的基本都差不多。正是因为如此，这两种回归可以归于同一个家族，即广义线性模型（generalizedlinear model）。

Logistic回归的主要用途：

• 寻找危险因素：寻找某一疾病的危险因素等；
• 预测：根据模型，预测在不同的自变量情况下，发生某病或某种情况的概率有多大；
• 判别：实际上跟预测有些类似，也是根据模型，判断某人属于某病或属于某种情况的概率有多大，也就是看一下这个人有多大的可能性是属于某病。

1.11 随机梯度下降

适合采样点数据巨大时：

1.12 感知器

感知器(Perceptron)这个词会成为Machine Learning的重要概念之一，是由于先辈们对于生物神经学科的深刻理解和融会贯通。对于神经(neuron)我们有一个简单的抽象：每个神经元是与其他神经元连结在一起的，一个神经元会受到多个其他神经元状态的冲击，并由此决定自身是否激发。
感知器算法实际上是在不断“猜测”正确的权重和偏移量，其实就是神经网络啦。

1.13 被动攻击算法

是在线学习算法，在线学习算法与其它算法的区别在于每次只能得到一个样本点，无法保留历史数据，对每一个新的样本点进行分析，根据分析的结果更新分类器。与神经网络类似。
1.14 鲁棒回归

如题，这种回归非常适合有缺失数据或者离群数据较多的情况。 在处理数据时，要顺序check以下几点： （1）离群数据是X维度的还是Y维度的？
这种情况就是Y维度的离群。
而这种情况就是X维度的离群。 (2)偏离的点的数量固然重要，但是偏离程度也很重要。

上面两种情况偏移的点数是一样的，但是偏移的量差很多。 估计器有两种：RANSAC，Theil sen Theil sen仅在点数和维度都不多时使用。 RANSAC是“RANdom SAmple Consensus（随机抽样一致）”的缩写。它可以从一组包含“局外点”的观测数据集中，通过迭代方式估计数学模型的参数。  一个简单的例子是从一组观测数据中找出合适的2维直线。假设观测数据中包含局内点和局外点，其中局内点近似的被直线所通过，而局外点远离于直线。简单的最小二乘法不能找到适应于局内点的直线，原因是最小二乘法尽量去适应包括局外点在内的所有点。相反，RANSAC能得出一个仅仅用局内点计算出模型，并且概率还足够高。但是，RANSAC并不能保证结果一定正确，为了保证算法有足够高的合理概率，我们必须小心的选择算法的参数。
举例：

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

n_samples = 1000
n_outliers = 50
#自己造数据
X, y, ceof = datasets.make_regression(n_samples=n_samples, n_features=1,
n_informative=1, noise=10,
coef=True, random_state=0)
np.random.seed(0)
X[:n_outliers] = 3 + 0.5 * np.random.normal(size=(n_outliers,1)) #将前50个点变成局外点
y[:n_outliers] = -3 + 10 * np.random.normal(size=n_outliers)

model = linear_model.LinearRegression()
model.fit(X, y)

model_ransac = linear_model.RANSACRegressor(linear_model.LinearRegression())
model_ransac.fit(X,y)
inlier_mask = model_ransac.inlier_mask_
outlier_mask = np.logical_not(inlier_mask)

model_ths = linear_model.TheilSenRegressor(random_state=42)
model_ths.fit(X, y)
# inlier_mask_ths = model_ths.inlier_mask_
# outlier_mask_ths = np.logical_not(inlier_mask_ths)

line_X = np.arange(-5,5)
line_y = model.predict(line_X[:,np.newaxis])
line_y_fansac = model_ransac.predict(line_X[:,np.newaxis])
line_y_ths = model_ths.predict(line_X[:,np.newaxis])

plt.plot(X[inlier_mask],y[inlier_mask],'.g',label='Inliers')
plt.plot(X[outlier_mask],y[outlier_mask],'.r',label='Outliers')
plt.plot(line_X, line_y, '-k',label='Linear regressor')
plt.plot(line_X, line_y_fansac, '-b', label='RANSAC regressor')
plt.plot(line_X, line_y_ths, '-r', label='Theil regressor')
plt.legend(loc='lower right')
plt.show()

上图：
可以看到普通的线性回归已经受影响了，而两种鲁棒回归都没有受影响，而且将50个局外点自动识别出来了，我们将他们标记为红色。

1.15 多项式回归

任何曲线可以近似地用多项式表示，所以在这种情况下我们可以用多项式进行逼近，即多项式回归分析。
实际是把[x1,x2]转换成为 举例：

<pre name="code" class="python">>>> <span style="font-family: Arial, Helvetica, sans-serif;">from sklearn.preprocessing import PolynomialFeatures</span>

>>> import numpy as np>>> X = np.arange(6).reshape(3, 2)>>> Xarray([[0, 1], [2, 3], [4, 5]])>>> poly = PolynomialFeatures(degree=2)>>> poly.fit_transform(X)array([[ 1., 0., 1., 0., 0., 1.], [ 1., 2., 3., 4., 6., 9.], [ 1., 4., 5., 16., 20., 25.]])

这里fit后的X，我们直接看第三行：不就是1，x1, x2, x1*x2, x1**2, x2**2 他是怎么回事后，我们还要知道怎么用？ 举例：

>>> from sklearn.preprocessing import PolynomialFeatures
>>> from sklearn.linear_model import LinearRegression
>>> from sklearn.pipeline import Pipeline
>>> import numpy as np
>>> model = Pipeline([('poly', PolynomialFeatures(degree=3)),
...                   ('linear', LinearRegression(fit_intercept=False))])
>>> # fit to an order-3 polynomial data
>>> x = np.arange(5)
>>> y = 3 - 2 * x + x ** 2 - x ** 3
>>> model = model.fit(x[:, np.newaxis], y)
>>> model.named_steps['linear'].coef_
array([ 3., -2.,  1., -1.])

我们看到这里完美的还原出了多项式y的系数。使用的是Pipeline。但提前告诉他是x**3数据(degree=3)。 如果是布尔数据，那么x1**2和x2**2这两个数据是没有意义的，所以设定只要交叉项。

>>> from sklearn.linear_model import Perceptron
>>> from sklearn.preprocessing import PolynomialFeatures
>>> import numpy as np
>>> X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
>>> y = X[:, 0] ^ X[:, 1]
>>> X = PolynomialFeatures(interaction_only=True).fit_transform(X)
>>> X
array([[ 1.,  0.,  0.,  0.],
       [ 1.,  0.,  1.,  0.],
       [ 1.,  1.,  0.,  0.],
       [ 1.,  1.,  1.,  1.]])
>>> clf = Perceptron(fit_intercept=False, n_iter=10, shuffle=False).fit(X, y)
>>> clf.score(X, y)
1.0

可以看到数据从6个变成了4个。