1. 程式人生 > >scikit-learn 學習筆記(一)

scikit-learn 學習筆記(一)

其實前一陣一直在看《機器學習系統設計》,但是發現書上程式碼不全,不好復現,而且書針對的主要是文字的處理,正好是我最不關心的方向,所以看到一半忍痛棄坑。那麼我們開始scikit-learn文件學習的旅程吧。

1. 線性模型

我對於線性模型的理解就是,將一系列的已知點回歸到一個線性方程上去,類似於:\hat{y}(w, x) = w_0 + w_1 x_1 + ... + w_p x_p之後就可以用這個方程來預測未知點的解了。 sklearn 規定 w0是 intercept_. w1--wp是coef_.

1.1普通最小二乘法

那麼如何將已知的點擬合到一個線性方程上去呢?sklearn提供的LinearRegression的思路就是讓點的已知值和預測值的平方和最小,即\underset{w}{min\,} {|| X w - y||_2}^2. 下面我們直接上例子:
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(套索)

她的應用是基於稀疏模型的情況,進行線性擬合,這時的效果較好。那麼什麼是稀疏模型呢?比如我們採集許多的溫度資料,肯定會存在大量的冗餘資料,因為溫度的變化是緩慢的,稀疏模型就會關注變化的點,而去除同質化的點,從而減少運算量。這個方法其實是壓縮感知的基礎。不知道壓縮感知,這裡有一個好的普及貼。 她要最小化的是:

\underset{w}{min\,} { \frac{1}{2n_{samples}} ||X w - y||_2 ^ 2 + \alpha ||w||_1}


為了便於理解,我們直接上例子:
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個。