scikit-learnを使ってみます。
http://scikit-learn.org/stable/supervised_learning.html#supervised-learning
データはstatsmodelsの記事と同じものを使っています。ドキュメントを見た限りだとstatsmodelsより様々なモデルがサクッと使えるような印象があります。
散々いたるところで解説された内容なので詳細は他のブログに任せます。
今回も以下を前提とします。
import pandas as pd import numpy as np from matplotlib import pyplot as plt
・単回帰 / 重回帰
from sklearn import linear_model #load data data = pd.read_csv("./anscombe.csv") data.drop("Unnamed: 0", axis = 1, inplace = True) X = data["x1"].reshape((-1, 1)) y = data["y1"].reshape((-1, 1)) #モデルの作成とプロット lm = linear_model.LinearRegression(fit_intercept = True, normalize = False) lm.fit(X, y) seq = np.arange(X.min() - 1.0, X.max() + 1.0, 0.1) plt.plot(seq, lm.intercept_[0] + seq * lm.coef_[0]) plt.scatter(X, y, c = "red") plt.grid(True) #決定係数も求められます。 lm.score(X, y) # -> 0.6665
重回帰のデータはRのMASSパッケージからUScrimeを.csvで出力して使います。
#load data from sklearn import linear_model data = pd.read_csv("./UScrime.csv") data.drop(["Unnamed: 0"], axis = 1, inplace = True) X = data.drop(["y"], axis = 1) y = data["y"].reshape(-1, 1) #モデルの作成と決定係数の確認 mul_lm = linear_model.LinearRegression(fit_intercept = True, normalize = True) mul_lm.fit(X, y) N = len(y) R2 = mul_lm.score(X, y) R2_d = 1 - ((1 - R2) * (N - 1) / (N - X.shape[1] - 1)) print(R2_d) # ->0.7169312456540442
linear_model.score()で得られるのはあくまで決定係数なので、重回帰のように説明変数が増えた場合は自由度調整済み決定係数に変える必要があります。どうやら、sklearnにはそのメソッドは用意されていないようなので自分で少し計算する必要があります。ちなみにstatsmodelsは最初からそれを計算してくれます。
import statsmodels.api as sm X = sm.add_constant(X) statsmodels_ml = sm.OLS(y, X).fit() statsmodels_ml.summary() """ R-squared: 0.803 Adj. R-squared: 0.717 """
確かに上で計算した値と同じ値が得られますね。
・Ridge回帰
リッジ回帰は最小二乗項に正則化項(L1ノルム)を加えることで多重共線性による解の不安定さを解消する方法です。上の重回帰をリッジ回帰に変えてみましょう。
from sklearn import linear_model data = pd.read_csv("./UScrime.csv") data.drop(["Unnamed: 0"], axis = 1, inplace = True) X = data.drop(["y"], axis = 1) y = data["y"].reshape(-1, 1) X = X - np.mean(X) scale_y = y - np.mean(y) scale_X = X.T #M行N列 M = scale_X.shape[0] #説明変数の数=14 N = scale_X.shape[1] #データ数=47 GCV_list = [] alpha_list = np.linspace(0, 10, num = 100) #正則化項の係数候補 for a in alpha_list: H = np.dot(scale_X.T, np.dot(np.linalg.inv(np.dot(scale_X, scale_X.T) + a * np.identity(M)), scale_X)) #N行N列 temp = np.dot(np.identity(N) - H, scale_y) GCV = np.dot(temp.T, temp) / (np.power(1 - (np.trace(H) / N), 2) * N) GCV_list.append(GCV[0][0]) plt.plot(range(len(GCV_list)), GCV_list) best_alpha_idx = np.argmin(GCV_list) print(alpha_list[best_alpha_idx]) # -> 0
リッジ回帰を行なう前に、正則化項の係数をいくつにするかが問題になります。そこで『入門機械学習による異常検知（コロナ社）』では一般化交差確認法（GCV：generalized cross validation）と呼ばれる方法が提案されています。本にあるRのコードによるとその係数は1桁のとある数になるはずなのですが、私のpythonのコードでは係数0になってしまいました。これではリッジ回帰ではなく重回帰と変わりません。どこか間違っているようです。
「0になったけどこれ正しいの？」という私と似た質問を見つけたのですが係数の範囲を広げてみてはという回答が出ています。広げてみたところ係数が3桁で最小になるという信じがたい数字になってしまい如何にも怪しい。
Understanding ridge regression results - Cross Validated
掲載されているRのコードに使われているlm.Ridge()のGCVに関する部分を見たところ、先にcoefを求めていて本で説明されていた様子と少し違う気がします。
https://rdrr.io/cran/MASS/src/R/lm.ridge.R
長々と書いてしまったが、どうしようもないので素直にRIdgeCV()を使うことにします。
reg_CV = linear_model.RidgeCV().fit(X, y) print(reg_CV.alpha_, reg_CV.coef_) plt.scatter(range(len(y)),y - reg_CV.predict(X).reshape(-1, 1))
・Lasso回帰
こちらは正則化項としてL1ノルムを採用したものです。パラメータが0になることがありモデルの解釈がしやすいのが特徴です。
#Lasso regression from sklearn import linear_model reg_CV = linear_model.LassoLarsCV().fit(X, y) print(reg_CV.alpha_) # ->1.7 plt.scatter(range(len(y)),y - reg_CV.predict(X).reshape(-1, 1)) #残差プロット
正則化項の係数を決めるにあたって、標本数に比べ説明変数の次元が多いときはLassoCV()、そうでなければLassoLarsCV()を使うのが良いみたいです。また、別の観点からの決め方として情報量規準（AICやBIC）を使った方法がLassoLarsIC()のようです。
・Elastic Net
Ridge回帰とLasso回帰の両方を組み合わせたのがElastic Netです。
#Elastic Net reg_CV = linear_model.ElasticNetCV().fit(X, y) print(reg_CV.alpha_) plt.scatter(range(len(y)),y - reg_CV.predict(X).reshape(-1, 1)) #残差プロット
結局のところどのモデルをとっても手順は全て同じです。あとはモデルの解釈の問題です。ちなみに自由度調整済み決定係数は、重回帰 > Ridge回帰 > Lasso回帰 > Elastic Net の順で、全て0.02ずつくらいの差でした。現在のデータに対してどれくらいフィットしているかなので何も正則化項の無い重回帰が大きくなるのは自然でしょう。
RIdge回帰Lasso回帰Elastic Netの詳細はこちらが分かりやすいと思います。
http://tekenuko.hatenablog.com/entry/2017/11/18/214317
・ロジスティック回帰
おなじみなので解説不要ですね。やってみます。データはstatsmodelsのときと同じように自分で作ったものです。
from sklearn import linear_model #create data x_1 = 5 y_1 = 7 x_2 = 7 y_2 = 5 x1 = x_1 + np.random.randn(50) y1 = y_1 + np.random.randn(50) x2 = x_2 + np.random.randn(50) y2 = y_2 + np.random.randn(50) #上を元にデータをpandas.dataframeで作成 x = pd.Series(np.hstack([x1, x2])) y = pd.Series(np.hstack([y1, y2])) c = pd.Series(np.repeat([1, 0], 50)) data = pd.DataFrame({"x":x, "y":y, "c":c}) #logistic regression reg = linear_model.LogisticRegression(penalty = "l2", solver = "saga").fit(data[["x", "y"]], data["c"]) print(reg.intercept_, reg.coef_[0]) #->[-0.05453222] [-1.39388269 1.36970588] #plot decision boundary seq = np.arange(data["x"].min() - 1, data["x"].max() + 1, 0.1) plt.scatter(data["x"][0:49], data["y"][0:49], c = "red") plt.scatter(data["x"][50:], data["y"][50:], c = "blue") plt.plot(seq, (-1 / reg.coef_[0][1]) * (reg.intercept_ + reg.coef_[0][0] * seq), c = "black") plt.grid(True) plt.xlabel("x") plt.ylabel("y") plt.title("logistic regression") #accuracy print(reg.score(data[["x", "y"]], data["c"])) #accuracyを計算するメソッドがある
statsmodelsのロジスティック回帰との違いとして、sklearnでは自動的に正則化項が組み込まれています。では正則化項が要らないという場合はどうしたら良いか調べてみたところ、正則化項の係数を小さくして影響を限りなく減らす方法が提案されていました。sklearnの場合デフォルト引数だとC = 1.0となっています。また係数が1 / Cとしてついていることに注意すると、Cをひたすら大きな数にすれば良いのです。
python - sklearn LogisticRegression without regularization - Stack Overflow
from sklearn import linear_model reg = linear_model.LogisticRegression(penalty = "l2", solver = "saga", C = 100.0, max_iter = 1000).fit(data[["x", "y"]], data["c"]) reg.predict(data[["x", "y"]])
次にマルチクラスの場合を扱います。statsmodelsと違って難しいこと考える必要なしです。
#上に続いてもう1つデータ追加 x_3 = 1 y_3 = 10 x3 = x_3 + np.random.randn(50) y3 = y_3 + np.random.randn(50) x = pd.Series(np.hstack([x1, x2, x3])) y = pd.Series(np.hstack([y1, y2, y3])) c = pd.Series(np.repeat([2, 1, 0], 50)) data = pd.DataFrame({"x":x, "y":y, "c":c}) #lgistic regression from sklearn import linear_model reg = linear_model.LogisticRegression(solver = "newton-cg").fit(data[["x", "y"]], data["c"]) #plot decision boundary seq = np.arange(data["x"].min() - 1, data["x"].max() + 1, 0.1) plt.scatter(data["x"][0:49], data["y"][0:49], c = "red") plt.scatter(data["x"][50:99], data["y"][50:99], c = "blue") plt.scatter(data["x"][100:], data["y"][100:], c = "green") plt.plot(seq, (-1 / reg.coef_[0][1]) * (reg.intercept_[0] + reg.coef_[0][0] * seq), c = "black") plt.plot(seq, (-1 / reg.coef_[1][1]) * (reg.intercept_[1] + reg.coef_[1][0] * seq), c = "pink") plt.grid(True) plt.xlabel("x") plt.ylabel("y") plt.title("logistic regression") #accuracy print(reg.score(data[["x", "y"]], data["c"])) #->0.94
・確率的勾配降下法
ここのテーマにふさわしいかどうかは微妙ですが、たまたまscikit-learnのGLMのテーマの中にあったので扱います。このクラスでは回帰や分類における10種類ほどの損失関数に対して確率的勾配降下法による最適化を提案してくれるようです。正直どの損失関数がどういう意味なのか知らないので私には違いは分からないのですが、とりあえず上のロジスティック回帰と同じことを別の最適化法で行なってみます。
from sklearn import linear_model reg = linear_model.SGDClassifier(loss = "modified_huber", penalty = "l2").fit(data[["x", "y"]], data["c"]) print(reg.intercept_, reg.coef_) #plot seq = np.arange(data["x"].min() - 1, data["x"].max() + 1, 0.1) plt.scatter(data["x"][0:49], data["y"][0:49], c = "red") plt.scatter(data["x"][50:99], data["y"][50:99], c = "blue") plt.scatter(data["x"][100:], data["y"][100:], c = "green") plt.plot(seq, (-1 / reg.coef_[0][1]) * (reg.intercept_[0] + reg.coef_[0][0] * seq), c = "black") plt.plot(seq, (-1 / reg.coef_[1][1]) * (reg.intercept_[1] + reg.coef_[1][0] * seq), c = "pink") plt.grid(True) plt.xlabel("x") plt.ylabel("y") plt.title("logistic regression by modified_huber")
問題無くできます。損失関数については要勉強ですね。
・パーセプトロン
パーセプトロンです。いわゆる〇〇Net系で最もシンプルな形と見ることができるやつです。これも結局のところこれまでと方法は同じです。
from sklearn import linear_model per = linear_model.Perceptron().fit(data[["x", "y"]], data["c"]) print(per.intercept_, per.coef_)
グラフはほとんど変わらないので略。
ここまで
http://scikit-learn.org/stable/modules/linear_model.html
にある様々なクラスを使いました。改めて見ると分かりますが、結局のところ全て手順は同じです。従って、これをデータに対して使い分けることがやはり重要かなと。そしてそのためにはそれぞれの特性を理解しておく必要があると感じました。
ここから他のクラスもまだまだ使ってみようと思っているのですが、方法が同じようなものはざっくりカットしていこうと思います。
![日本統計学会公式認定 統計検定 2級 公式問題集[2014〜2016年]](https://images-fe.ssl-images-amazon.com/images/I/51frwnCObjL._SL160_.jpg "日本統計学会公式認定 統計検定 2級 公式問題集[2014〜2016年]")
