ipython notebookを使って統計的学習

伴走テキスト

The Elements of Statistical Learning, Second Edition, Springer http://statweb.stanford.edu/~tibs/ElemStatLearn/

1. Introduction

Scatter plot matrixを描いてみる。

# 変数の名前のリストを作り、変数の数ｘサンプル数の正規乱数行列を作る
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
df = pd.DataFrame(np.random.randn(1000, 4), columns=['a', 'b', 'c', 'd'])
pd.scatter_matrix(df) # pandasのscatter_matrix()を使う
plt.show()

DNA microarray dataを描いてみる。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
# http://qiita.com/s-wakaba/items/a93f03f27137cff4a26c より
# 赤→黒→緑で線形に変化するカラーマップの作成
microarray_cmap = LinearSegmentedColormap('microarray', {
    'red': [(0.0, 1.0, 1.0), (0.5, 0.2, 0.2), (1.0, 0.0, 0.0)],
    'green': [(0.0, 0.0, 0.0), (0.5, 0.2, 0.2), (1.0, 1.0, 1.0)],
    'blue': [(0.0, 0.0, 0.0), (0.5, 0.2, 0.2), (1.0, 0.0, 0.0)],
})
# 補間スムージングはしない: interporation ='none'
# -1より小、1より大は-1,1にする
def draw_heatmap(a, cmap=microarray_cmap):
    from matplotlib import pyplot as plt
    plt.imshow(a, aspect='auto', interpolation='none',
               cmap=cmap, vmin=-1.0, vmax=1.0)
    plt.colorbar()
    #plt.xticks(range(a.shape[1]), a.columns, rotation=15)
    #plt.yticks(range(a.shape[0]), a.index, size='small')
    plt.gca().get_xaxis().set_ticks_position('none')
    plt.gca().get_yaxis().set_ticks_position('none')
    plt.show()

n_sample = 50
n_gene = 200
a = np.random.randn(n_gene,n_sample)

draw_heatmap(a,microarray_cmap)

2. Overveiw of Supervised Learning

Linear models

2次元平面に境界線を決め、それに応じて2群に分かれているとする。 ただし、乱雑項があるので、その分離は完璧ではない。

import numpy as np
n_variables = 2
n_samples = 100
x = np.random.randn(n_samples,n_variables)
a = [3,1]
# y = 3x1 + x2; y = 0.5 -> x2 = -3x1 + 0.5
y_obs = np.dot(x,a) + np.random.randn(n_samples)*0.8 # 行列積np.dot()
# ラべリング
z = np.zeros(n_samples) # 値が0のベクトル
z[y_obs>0.5] = 1 # 条件でフィルタリングして付置

plt.scatter(x[:,0],x[:,1],c=(z+1)*2,alpha=0.5) # 散布図 cは色, alphaは透明度

# 線形回帰の推定境界線
b = np.dot(np.dot(np.linalg.inv(np.dot(x.T,x)),x.T),y_obs) # 行列の転置・逆行列・積
print(b)
# プロットに直線を加えるときは、2端点座標を指定して線タイプのプロット
X = np.array([-3,3])
Y = np.array([((-b[0])*(-3)+0.5)/b[1],((-b[0])*3+0.5)/b[1]])
plt.plot(X,Y, 'k-', lw=1,c="red")
# 真の境界線
X = np.array([-3,3])
Y = np.array([(-3)*(-3)+0.5,(-3)*3+0.5])
plt.plot(X,Y, 'k-', lw=1)
plt.show()

[ 3.0255729   1.09571771]

k-nearest neighbors

from sklearn.neighbors import KNeighborsClassifier

n = 13 # 考慮する近傍の点の数
neigh = KNeighborsClassifier(n_neighbors=n)
neigh.fit(x, z) 
from itertools import product # itertoolsはイテレータベースの高速組み合わせ・順列モジュール
import matplotlib.pyplot as plt
# 2次元グリッド点を作る
XX1 = np.linspace(-3, 3, 50)
XX2 = np.linspace(-3, 3, 50)
points = np.array(list(product(XX1,XX2)))

prd = neigh.predict(points)
plt.scatter(points[:,0],points[:,1],c=prd*2,alpha=0.5)

<matplotlib.collections.PathCollection at 0x1b0ae770>

n = 1 # 考慮する近傍の点の数
neigh = KNeighborsClassifier(n_neighbors=n)
neigh.fit(x, z) 
from itertools import product # itertoolsはイテレータベースの高速組み合わせ・順列モジュール
import matplotlib.pyplot as plt
# 2次元グリッド点を作る
XX1 = np.linspace(-3, 3, 50)
XX2 = np.linspace(-3, 3, 50)
points = np.array(list(product(XX1,XX2)))

prd = neigh.predict(points)
plt.scatter(points[:,0],points[:,1],c=prd*2,alpha=0.5)

<matplotlib.collections.PathCollection at 0x19b79df0>

3. Linear Methods for Regression

Linear regression

# 変数の名前のリストを作り、変数の数ｘサンプル数の正規乱数行列を作る
%matplotlib inline
n_var = 4
n_sample = 1000
import numpy as np
import matplotlib.pyplot as plt
x = np.random.randn(n_sample,n_var)
a = np.random.randn(n_var)
np.sort(a)
y = np.dot(x,a) + np.random.randn(n_sample) # 行列積np.dot()
from sklearn import linear_model
clf = linear_model.LinearRegression()

clf.fit (x, y)
print(clf.coef_)
print(clf.intercept_)
print(a)

[-0.60247813  0.79632332 -0.91892781 -1.78944024]
0.0682356254769
[-0.59983837  0.76358198 -0.88379711 -1.771415  ]

Subset selection

Ridge, Lasso, LARS

理論はいろいろあるけれど、pythonでは、scikit-learnのlinear_model内にある LinearRegression()をRidge(),Lasso(),Lars()に変えるだけ

Grouped Lasso

import matplotlib.font_manager as fm
fm.findSystemFonts()
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Osaka'

n_var = 4
n_sample = 1000
import numpy as np
import matplotlib.pyplot as plt
x = np.random.randn(n_sample,n_var)
a = np.random.randn(n_var)
np.sort(a)
y = np.dot(x,a) + np.random.randn(n_sample) # 行列積np.dot()
from sklearn import linear_model
clf_lin = linear_model.LinearRegression()
clf_rid = linear_model.Ridge (alpha = .1) # ridge parameter =0.1
clf_las = linear_model.Lasso (alpha = .1) # lasso parameter =0.1
clf_lar = linear_model.Lars (n_nonzero_coefs=n_var) 
clf_lin.fit (x, y)
print("線形回帰")
print([clf_lin.coef_])
print(clf_lin.intercept_)
print("リッジ回帰")
clf_rid.fit(x,y)
print([clf_rid.coef_])
print(clf_rid.intercept_)
print("Lasso")
clf_las.fit(x,y)
print([clf_las.coef_])
print(clf_las.intercept_)
print("LARS")
clf_lar.fit(x,y)
print([clf_lar.coef_])
print(clf_lar.intercept_)
print(a)

線形回帰
[array([ 0.07368183, -1.77424962, -1.65112682,  0.13016486])]
0.0395470848406
リッジ回帰
[array([ 0.07366312, -1.77407815, -1.65096032,  0.13016244])]
0.0395472550757
Lasso
[array([ 0.        , -1.68080617, -1.54865364,  0.03725043])]
0.0389295883433
LARS
[array([ 0.07368183, -1.77424962, -1.65112682,  0.13016486])]
0.0395470848406
[ 0.06577191 -1.78688987 -1.66071836  0.15159026]

4. Linear Methods for Classification

Linear regressin and classification

import numpy as np
n_variables = 2
n_samples = 1000
x = np.random.randn(n_samples,n_variables)
a = [3,1]
# y = 3x1 + x2; y = 0.5 -> x2 = -3x1 + 0.5
y_obs = np.dot(x,a) + np.random.randn(n_samples)*0.8 # 行列積np.dot()
# ラべリング
z = np.zeros(n_samples) # 値が0のベクトル
z[y_obs>0.5] = 1 # 条件でフィルタリングして付置

plt.scatter(x[:,0],x[:,1],c=(z+1)*2,alpha=0.5) # 散布図 cは色, alphaは透明度

# 線形回帰の推定境界線
b = np.dot(np.dot(np.linalg.inv(np.dot(x.T,x)),x.T),y_obs) # 行列の転置・逆行列・積
print(b)
# プロットに直線を加えるときは、2端点座標を指定して線タイプのプロット
#X = np.array([-3,3])
#Y = np.array([((-b[0])*(-3)+0.5)/b[1],((-b[0])*3+0.5)/b[1]])
#plt.plot(X,Y, 'k-', lw=1,c="red")
# 真の境界線
#X = np.array([-3,3])
#Y = np.array([(-3)*(-3)+0.5,(-3)*3+0.5])
#plt.plot(X,Y, 'k-', lw=1)
#plt.show()

[ 2.99157759  1.02064718]

Linear discriminant analysisで同じことをする

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
sklearn_lda = LinearDiscriminantAnalysis(n_components=2)
z_pred = sklearn_lda.fit(x, z).predict(x)

plt.scatter(x[:,0],x[:,1],c=(z_pred+1)*2,alpha=0.5)

<matplotlib.collections.PathCollection at 0x199415d0>

Linear discriminant analysis を3群でやる

import numpy as np
n_variables = 2
n_samples = 1000
x = np.random.randn(n_samples,n_variables)
a = [3,1]
# y = 3x1 + x2; y = 0.5 -> x2 = -3x1 + 0.5
y_obs = np.dot(x,a) + np.random.randn(n_samples)*5 # 行列積np.dot()
# ラべリング
z = np.zeros(n_samples) # 値が0のベクトル
z[y_obs> -2] = 1 # 条件でフィルタリングして付置
z[y_obs> 1] = 2 # 条件でフィルタリングして付置
x[z==1,0] += 1
x[z==2,1] += 1
#print(z)
lda = LinearDiscriminantAnalysis()
z_pred = lda.fit(x,z).predict(x)
#print(z_pred)
print("original label")
plt.scatter(x[:,0],x[:,1],c=(z+1)*2,alpha=0.5)
plt.show()
print("predicted label")
plt.scatter(x[:,0],x[:,1],c=(z_pred+1)*2,alpha=0.5)
plt.show()

original label

predicted label

Quadratic Discriminant analysis

曲線でしか分離できない例をやってみる

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
n_variables = 2
n_samples = 1000
x = np.random.randn(n_samples,n_variables)
a1 = [3,1]
a2 = [2,0]
a3 = [0,1]
y_obs = np.dot(x,a1) + x[:,0]**2 * 2 + x[:,1]**2 + np.random.randn(n_samples)*1

# ラべリング
z = np.zeros(n_samples) # 値が0のベクトル
z[y_obs> 0] = 1 # 条件でフィルタリングして付置
z[y_obs> 3] = 2 # 条件でフィルタリングして付置


lda = LinearDiscriminantAnalysis()
qda = QuadraticDiscriminantAnalysis()
z_pred_lin = lda.fit(x,z).predict(x)
z_pred_quad = qda.fit(x,z).predict(x)

print("original label")
plt.scatter(x[:,0],x[:,1],c=(z+1)*2,alpha=0.5)
plt.show()
print("LDA predicted label")
plt.scatter(x[:,0],x[:,1],c=(z_pred_lin+1)*2,alpha=0.5)
plt.show()
print("QDA predicted label")
plt.scatter(x[:,0],x[:,1],c=(z_pred_quad+1)*2,alpha=0.5)
plt.show()

original label