scikit-learn – make_blobs

概要

sklearn.datasets.make_blobls()は、クラス分類のためのデータを生成する。blobとはインクの染みなどを指し、散布図の点の様子からつけられてるようだ。

標準では、データの総数、特徴量の数、クラスターの数などを指定して実行し、特徴量配列X、ターゲットとなるクラスデータyのタプルが返される(引数の指定によってはもう1つ戻り値が追加される)。

得られるデータの形式

特徴量配列Xは列が特徴量、行がレコードの2次元配列。ターゲットyはレコード数分のクラス属性値の整数。

利用例

そのままscikit-learnのモデルの入力とする。

クラスごとに色やマークを変えて散布図を描く。

パラメーターの指定

主なもの。

n_samples
整数で指定した場合、生成されるサンプルの総数で戻り値Xの行数になる。配列で指定した場合、その要素数がクラスターの数となり、各要素はクラスターのデータ数となる。デフォルトは100。
n_features
特徴量の数で、戻り値Xの列数になる。デフォルトは2
centers
クラスター中心の数。n_samplesを整数で指定してcentersを指定しない場合(デフォルトのNoneの場合)、centers=3となる。n_samplesを配列で指定した場合はNoneか[n_centers, n_features]の配列。
center_std
クラスターの標準偏差。

 

Logistic回帰~cancer~Pythonではじめる機械学習より

モデルの精度

breast_cancerデータセットに対してLogistic回帰モデル、scikit-learnLogisticRegression適用し、訓練データとテストデータのスコアを計算してみる。

(注)solverに関する警告と計算結果

上のコードを実行したとき、結果は書籍と整合しているが、警告表示が出た

この時点でscikit-learnのバージョンが古く(0.21.3)、将来のデフォルトが変更されるとのこと。そこでインスタンス生成時にデフォルトのソルバーを明示的にsolver='liblinear'と指定して実行すると、警告は出ず値もそのまま。

なお、solver='lbfgs'としてみたところ、計算が収束しない旨の警告が出た。

そこで収束回数を増やしていったところ、最大回数2000では収束せず、3000で収束し、警告は出なくなった。

その後、scikit-learnのバージョンを0.23.0にアップグレードしたところ、デフォルトで警告は表示されず、収束回数に関する警告が同じように出て、結果も再現された。以下、ソルバーとしてliblinearを明示的に指定し、random_stateの値も書籍と同じ値として確認する。

学習精度の向上

先のC=1.0liblinearによるスコアは、訓練データに対して0.953、テストデータに対して0.958と両方に対して高い値となっている。ここで、訓練データとテストデータのスコアが近いということは、適合不足の可能性がある。そこでC=100と値を大きくして、より柔軟なモデルにしてみる(柔軟なモデルとは、正則化を弱めて訓練データによりフィットしやすくしたモデル)。

訓練データ、テストデータともそれぞれ若干向上している。なお、Cの値を1000、10000ともっと大きくしてもスコアはほとんど変わらない。

今度は逆に、Cの値を1.0より小さくして正則化を強めてみると、訓練データ、テストデータ両方に対するスコアが下がってしまう。

Cを変化させたときの学習率曲線は以下の通り。Cが10より小さいところでは正則化が強く学習不足、そこを超えると学習率が頭打ちで、学習率の改善はそれほど顕著ではない。Logistic回帰モデルの学習率曲線のバリエーションについては、こちらでまとめている。

特徴量の係数

L2正則化の場合

breast_cancerデータセットに対してLogisticRegressionを学習させた場合の、30個の特徴量に対する係数をプロットする。liblinearソルバーで、デフォルトでL2正則化を行っている。Cの値が大きいほど正則化の効果が弱く、係数の絶対値が大きくなっている。

書籍で注意喚起しているのは3番目の特徴量mean perimeterで、モデルによって正負が入れ替わることから、クラス分類に対する信頼性を問題にしている。

ここで書籍について以下の点が気になった。

  • logreg001のインスタンス生成時にC=0.01としているが、凡例で”C=0.001″としている(グラフの結果はあまり変わらない)
  • logreg100C=100とすると、書籍にあるような結果にならない(worst concave pointsが-8以下になるなど、分布が大幅に変わってくる)
  • C=20とすると、概ね書籍と同じ分布になる(若干異なる部分は残る)

いずれにしても”Pythonではじめる機械学習”は、入門者にとってとてもありがたいきっかけを提供してくれる良著であることに変わりはない。

L1正則化の場合

ソルバーを同じliblinearとして、penalty='l1'と明示的に指定する。今回はL2正則化の時と違って、C=0.001はコード中に明示され、C=100としてスコアの計算結果まで合う。ただしset_ylim()によって表示範囲を制限しており、C=100に対するいくつかの点が枠外にある。

L1正則化によって、多くの係数がゼロとなり、少ない特徴量によるシンプルなモデルでそれなりのスコアを出している。

係数の符号と選択確率について

ターゲットのクラスは、malignant(悪性)が0、benign(良性)が1で、係数が正の場合は良性となる確率を上げる方向に、負の場合は悪性となる確率を上げる方向に効くことになる。

ここでL2正則化のworst concavityを見てみると、負~0の値をとっているが、元のデータを俯瞰すると良性の集団の方が全体的に高い値を示していて矛盾している。一方、L1正則化の場合は、C=0.001で全ての係数がゼロとなっていて、結果に影響していないことを示唆している。

L1正則化で正則化の程度を弱めて、C=1, 0.5, 0.1としてみると、worst concavityは結局ゼロとなるが、worst textureは一貫して負の値を維持している。この傾向はarea errorにも僅かだが見られる。

cancerデータを俯瞰してみると、worst textureは良性・悪性の分布がかなり重なっていて、悪性のデータのボリュームが大きい。area errorも両クラスのデータが近く、値が小さく、良性のデータ量が卓越している。

ヒストグラムを見る限りほとんどの特性量の値が大きいときに良性を示唆しているようみ見えるが、Logistic回帰の結果からは、多くの特性量が効いておらず、中には分布からの推測と逆の傾向を示す。

 

scikit-learn – LogisticRegression

概要

scikit-learnLogisticRegressionモデルはLogistic回帰のモデルを提供する。利用方法の概要は以下の手順で、LinearRegressionなど他の線形モデルとほぼ同じだが、モデルインスタンス生成時に与える正則化パラメーターCRidge/Lassoalphaと逆で、正則化の効果を強くするにはCを小さくする(Cを大きくすると正則化が弱まり、訓練データに対する精度は高まるが過学習の可能性が高くなる)。

また、正則化の方法をL1正則化、L2正則化、Elastic netから選択できる。

  1. LogisticRegressのクラスをインポートする
  2. ハイパーパラメーターC、正則化方法、solver(収束計算方法)などを指定し、モデルのインスタンスを生成する
  3. fit()メソッドに訓練データを与えて学習させる

学習済みのモデルの利用方法は以下の通り。

  • score()メソッドにテストデータを与えて適合度を計算する
  • predict()メソッドに説明変数を与えてターゲットを予測
  • モデルインスタンスのプロパティーからモデルのパラメーターを利用
    • 切片はintercept_、重み係数はcoef_(末尾のアンダースコアに注意)

利用例

以下は、breast_cancerデータセットに対してLogisticRegressionを適用した例。デフォルトのsolver'lbfgs'でデフォルトの最大収束回数(100)では収束しなかったため、max_iter=3000を指定している。

利用方法

LogisticRgressionの主な利用方法はLineaRegressionとほとんど同じで、以下は特有の設定を中心にまとめる。

モデルクラスのインポート

scikit-learn.linear_modelパッケージからLogisticRegressonクラスをインポートする。

モデルのインスタンスの生成

LogisticRegressionでは、ハイパーパラメーターCによって正則化の強さを指定する。このCRidge/Lassoalphaと異なり、正則化の効果を強めるためには値を小さくする。デフォルトはC=1.0

以下、RidgeLassoに特有のパラメーターのみ説明。LinearRegressionと共通のパラメーターはLinearRegressionを参照。

penalty
'l1''l2''elasticnet', 'none'で正則化項のノルムのタイプを指定する。ソルバーの'newton-cg','sag','lbfgs'はL2正則化のみサポートし、'elasticnet''saga'のみがサポートする。デフォルトは'none'で正則化は適用されない('liblinear''none'に対応しない)。
tol
収束計算の解の精度で、デフォルトは1e-4。
C
正則化の強さの逆数。正の整数で指定し、デフォルトは1.0。
solver
'newton-cg''lbfgs''liblinear''sag''saga'のうちから選択される。デフォルトは'lbfgs'。小さなデータセットには'liblnear'が適し、大きなデータセットに対しては'sag''saga'の計算が速い。複数クラスの問題には、'newton-cg''sag''saga''lbfgs'が対応し、'liblinear'は一対他しか対応しない。その他ノルムの種類とソルバーの対応。
max_iter
収束計算の制限回数を指定する。デフォルト値は100。
random_state
データをシャッフルする際のランダム・シードで、solver='sag'の際に用いる。
l1_ratio
Elastic-Netのパラメーター。[0, 1]の値で、penalty='elasticnet'の時のみ使われる。

 モデルの学習

fit()メソッドに特徴量とターゲットの訓練データを与えてモデルに学習させる(回帰係数を決定する)。

X
特徴量の配列。2次元配列で、各列が各々の説明変数に対応し、行数はデータ数を想定している。変数が1つで1次元配列の時はreshape(-1, 1)かスライス([:, n:n+1])を使って1列の列ベクトルに変換する必要がある。
y
ターゲットの配列で、通常は1変数で1次元配列。

3つ目の引数sample_weightは省略。

適合度の計算

score()メソッドに特徴量とターゲットを与えて適合度を計算する。

その他のメソッド

  • decision_function(X)
  • densiffy()
  • predict_proba(X)
  • predict_log_proba()
  • sparsify()

 

Logistic回帰~forgeデータ~Pythonではじめる機械学習より

概要

O’REILLYの書籍”Pythonではじめる機械学習”の2.3.3.5、Logistic回帰でforgeデータの決定境界をトレースしてみたとき、収束計算のソルバーの違いや、元データと書籍のデータの違いなどから再現性に悩んだので記録しておく。

決定境界

mglearnのforgeデータセットに対してLogisticRegressionを適用してみる。

Cがかなり大きい場合、すなわち正則をほとんど行わない場合には、与えられたデータに対して可能な限り適合させようとしており、データに対する適合度は高い。Cが小さくなると正則化が効いてきて、データ全体に対して適合させようとしているように見える。

ここで上の図のC=1のケースは、書籍の図2-15右側と比べると決定境界の勾配が逆になっている。その理由は次のようであることが分かった。

  • 書籍ではLogisticRegression()の収束手法を指定せず、デフォルトのsolver='liblinear'が使用されている
  • 今回指定なしで実行したところ、以下のような警告が発生
    • FutureWarning: Default solver will be changed to ‘lbfgs’ in 0.22. Specify a solver to silence this warning.
      FutureWarning)
    • デフォルトのソルバーが(現在はliblinearだが)ver 0.22ではlbfgsになる/このwarningを黙らせるためにソルバーを指定せよ
  • そこでモデルのインスタンス生成時にLogisticRegression(solver='lbfgs')としたところ先の結果となった
  • 指定なし、あるいはsolver='liblinear'とすると書籍と同じ結果になる

liblinearによる結果が以下の通り。正則化の度合いに応じてlbfgsよりも傾きがダイナミックに変わっているように見える。

なお、これらの図の傾きについて、今度は書籍の図2-16と随分違っている。よく見てみると、同図のforgeデータは特に下側の〇印の点でオリジナルにはないデータがいくつか加わっているためと考えられる。

これらのコードは以下の通り。

3次元表示

2つのCの値について、二つの特徴量の組み合わせに対する青い点の確率分布を表示してみる(solver='lbfgs')。Cが小さいと確率分布がなだらかになる様子が見て取れるが、データに対する判別の適合度との関係はよくわからない。

 

Ridge回帰とLasso回帰

概要

回帰は、以下のようなm個の特徴量に関するnセットのデータXとそれらに対するターゲット値yについて、xからyを推定するモデルを決定する。

(1)    \begin{equation*} \boldsymbol{X} = \left[ \begin{array}{ccc} x_{11} & \cdots & x_{m1} \\ \vdots & & \vdots \\ x_{1n} & \cdots & x_{mn} \\ \end{array} \right] \left[ \begin{array}{c} y_1 \\ \vdots \\ y_n \end{array} \right] \quad \Rightarrow \quad y = f(\boldsymbol{x}) \end{equation*}

線形回帰は、モデルの関数形を以下のような特徴量に関する線形式とする。

(2)    \begin{equation*} \hat{y} = w_0 + w_1 x_1 + \cdots + w_m x_m \end{equation*}

通常線形回帰(重回帰、多重回帰)の場合、これを以下のような最小化問題として解く。

(3)    \begin{equation*} \mathrm{minimize} \quad \sum_i (y_i - \hat{y}_i)^2 \end{equation*}

通常線形回帰では、全ての訓練データに対する予測誤差を最小化しようとするが、このことで大きく外れた特徴量に対しても何とか合わせようとすることになる。このような状態を過学習と呼び、訓練データに対する予測精度は高くなるが、モデルが訓練データの状態に過敏に反応して、全般的な特徴に対する精度が却って低くなる(過学習~多項式回帰の場合)。

そこで、通常線形回帰の最適化に対して、全体的に重み係数の影響を小さくするための正則化項(罰金項、ペナルティー項)を考慮する。通常、ペナルティー項としては重み係数のノルムが用いられる(右辺第1項や第2項に分数の係数をつけることがあるが、計算の便宜のためであり本質への影響はない)。

(4)    \begin{equation*} \mathrm{minimize} \quad \sum_i (y_i - \hat{y}_i) + \alpha \sum_j |w_j|^p \end{equation*}

正則化項が重みの大きさを制限しようとするものであること、この式がこれを制約とした制約条件付き最適化問題であることは正則化の意味にまとめた。

このノルムにおいて、p=1(L1ノルム)の場合をLasso回帰、p=2(L2ノルム)の場合をRidge回帰と呼び、重みに対する制限のほかに以下のような特徴がある。

Ridge回帰
特徴量間の相関が高い場合~多重共線性(multicolinearity)が強い場合や一時従属な場合、通常線形回帰では解が求まらなかったりモデルが不安定になるが、Ridge回帰は何とか解を求められるようになる。
Lasso回帰
多数の特徴量のうち効果が小さいものの係数がゼロになり、モデルの複雑さを緩和できる。

Ridge回帰

Ridge回帰は、多重線形回帰の最適化において重み係数のL2ノルムを正則化項として付加する。

(5)    \begin{align*} &\mathrm{minimize} \quad \sum_i (y_i - \hat{y}_i) + \alpha \sum_j |w_j|^2 \\ & \mathrm{where} \quad \hat{y}_i = w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} \end{align*}

Ridge回帰は、特徴量の重みの強さを制限する(係数の絶対値を小さくする)効果を持つとともに、特徴量間の線形性が強い場合は予測式が不安定になることを防ぐ。

Ridge回帰の解析的な理解

Lasso回帰

Lasso回帰は、多重線形回帰の最適化において重み係数のL1ノルムを正則化項として付加する。

(6)    \begin{align*} &\mathrm{minimize} \quad \sum_i (y_i - \hat{y}_i) + \alpha \sum_j |w_j| \\ & \mathrm{where} \quad \hat{y}_i = w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} \end{align*}

Lasso回帰もRidge回帰と同じく、特徴量の係数の重みを制限するが、正則化を強めるとともに係数がゼロとなり、モデルがシンプルになるという特性がある。

Lasso回帰の解析的な理解

Ridge回帰とLasso回帰の挙動

係数の大きさ

Pythonのscikit-learnで得られる糖尿病に関するdiabetesデータセットを使って、同じくscikit-learnのRidge回帰モデルとLasso回帰モデルの挙動を比べてみる。alphaを大きくして正則化を強めるほど、全体的に係数の絶対値が小さくなっている。Ridgeの場合は必ずしも係数をゼロにしないのでモデルの複雑さが残るのに対して、Lassoの場合、係数は正則化が強いほど多くの係数がゼロになりモデルがシンプルになる。

alphaの増加に伴うRidgeのスコアは以下の通りで、そもそも訓練データに対するスコアが低い。もともと10個程度の特徴量ではそれほどの精度が期待できないようだ。

Lassoのスコアも同様に低い。alpha=10ではLasso回帰の特性から全ての係数がゼロとなり、相関係数がゼロとなっている。

この計算に用いたPythonのコードは以下の通り。

学習曲線

特徴量を増やすために、Boston house-pricesデータセットの特徴量データを拡張して試す。13個の特徴量に加えて、それらの特徴量同士の積から新たな特徴量を生成する。その結果、全体の特徴量数は単独の特徴量13、各特徴量の2乗が13、2つの特徴量の積が13C2 = 78の合計で104個となる。この特徴量データとターゲットの住宅価格について訓練データとテストデータに分け、Ridge回帰とLasso回帰のハイパーパラメータalphaを変化させてスコアの変化を見たのが以下の図。

Ridge、Lassoとも訓練データのスコアに対してテストデータのスコアは低く、過学習の様子がわかる。Ridgeではalpha=100程度でテストデータのスコアが最も高く0.75程度となる。Lassoの方はalpha=0.1程度でテストデータのスコアが最も高く、これも0.75を少し上回る程度。またLassoについては、alphaを増やしていくとゼロとなる係数の数が増えていき、それに伴って訓練データのスコアも下がっている。

Boston house-pricesデータに対して、RidgeとLassoの2つのモデルのみを検討するなら、計算コストがより少ないLasso回帰でalpha=0.1程度を選択することになろうかと考えられる。

この計算のコードは以下の通り。

 

Diabetesデータセット

概要

diabetesデータは、年齢や性別など10個の特徴量と、それらの測定1年後の糖尿病の進行度に関する数値を、442人について集めたデータ。出典は”From Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani (2004) “Least Angle Regression,” Annals of Statistics (with discussion), 407-499″。

ここではPythonのscikit-learnにあるdiabetesデータの使い方をまとめる。

データの取得とデータ構造

Pythonで扱う場合、scikit-learn.datasetsモジュールにあるload_diabetes()でデータを取得できる。データはBunchクラスのオブジェクト

データの構造は辞書型で、442人の糖尿病に関する10個の特徴量をレコードとした配列、442人の測定1年後の糖尿病の進行度を示す数値データの配列など。

データのキーは以下のようになっている。

データの内容

'data'~特徴量データセット

10個の特徴量を列とし、442人の被検者を業とした2次元配列。DESCRに説明されているように、これらのデータは標本平均と標本分散で正規化されており、各特徴量とも、データの和はゼロ(正確には1×10-14~1×10-13のオーダーの実数)、2乗和は1となる。

'target'~糖尿病の進行度

442人に関する10個の特徴量データを測定した1年後の糖尿病の進行度を示す数値。原文でも”a measure of disease progression one year after baseline”としか示されていない。このデータは正規化されていない

'feature_names'~特徴名

10種類の特徴量の名称

sklearn R
0 age age 年齢
1 sex sex 性別
2 bmi bmi BMI(Body Mass Index)
3 bp map (動脈の)平均血圧(Average blood pressure)
4 S1 tc 総コレステロール?
5 S2 ldl 悪玉コレステロール(Low Density Lipoprotein)
6 S3 hdl 善玉コレステロール(High Density Lipoprotein)
7 S4 tch 総コレステロール?
8 S5 ltg ラモトリギン?
9 S6 glu 血糖=グルコース?

scikit-learnでは後半のデータがs1s6とだけ表示されていて、DESCRにおいても”six blood serum measurements”とだけ書かれている。Rのデータセットでは、これらがtc, ldlなど血清に関する指標の略号で示されている。

tcとtchはどちらも総コレステロールに関するデータのようだが、どういう違いなのかよくわからない。少なくとも双方に正の相関があるが、ばらつきは大きい。

'filename'~ファイル名

CSVファイルのフルパス名が示されている。scikit-learnの他のデータセットと以下の2点が異なっている。

  • 特徴量データdiabetes_data.csvとターゲットデータdiabetes_target.csvの2つのファイルに分かれている
  • ファイルの拡張子がcsvとなっているが、区切りはスペースとなっている

diabetes_data.csv

1行に10個の実数がスペース区切りで配置されており、442行のデータがある。442人分の10個の特徴量データ

diabetes_target.csv

ターゲットyに相当する442行の実数データ。

‘DESCR’~データセットの説明

データセットの説明。各特徴量データが標準化されていることが説明されている。

データの利用

各データの取得方法

data、targetなどのデータを取り出すのに、以下の2つの方法がある。

  • 辞書のキーを使って呼び出す(例:diabetes['data']
  • キーの文字列をプロパティーに指定する(例:diabetes.data

dataの扱い

そのまま2次元配列として扱うか、pandas.DataFrameで扱う。特定の特徴量データを取り出すには、ファンシー・インデックスを使う。

 

過学習~多項式回帰の場合

概要

過学習(over fitting)の例として、多項式の係数を線形回帰で予測した場合の挙動をまとめてみた。

複数の点(xi, yi)に対して、以下の線形式の項数を変化させて、Pythonのパッケージ、scikit-learnにあるLinearRegressionでフィッティングさせてみる。

(1)    \begin{equation*} \hat{y} = w_0 + \sum_{j=1}\m w_j x^j \end{equation*}

データ数が少ない場合

以下の例は、[-3, 1]の間で等間隔な4つの値を発生させ、(x, ex)となる4つの点を準備、これらのデータセットに対して、多項式の項数(すなわちxの次数)を1~6まで変化させてフィッティングした結果。たとえばn_terms=3の場合はy = w_0 + w_1 x + w_2 x^2 + w_3 x^3の4つの係数を決定することになる。

  • n_terms=1の場合は単純な線形関数で、データセットの曲線関係を表しているとは言えない。
  • n_terms=2になるとかなり各点にフィットしているが、x < −1の範囲で本来の関数の値と離れていく。
  • n_terms=3はデータ数より項数(特徴量の数)が1つ少ない。各点にほぼぴったり合っていて、最も「それらしい」(ただしデータセットの外側の範囲でも合っているとは限らない/指数関数に対してxの有限の多項式ではどこかで乖離していく)
  • n_terms=4はデータ数と項数(特徴量の数)が等しい。予測曲線がすべての点を通っているが、無理矢理合わせている感があり、データセットの左側で関数形が跳ね上がっている。
  • n_terms=5はデータ数より特徴量数の方が多くなる。予測曲線は全ての点を通っているが、1番目の点と2番目の点の間で若干曲線が歪んでいる
  • n_terms=6になると歪が大きくなる

上記の実行コードは以下の通り。

  • 7~8行目は、切片・係数のセットとxの値を与えて多項式の値を計算する関数。
  • 19行目でn_data=4個のxの値を発生させ、20行目で指数関数の値を計算している。後のために乱数でばらつかせる準備をしているが、ここではばらつかせていない
  • 23~24行目でxnの特徴量を生成している
  • 35行目で線形回帰モデルのフィッティングを行っている。n_termsで指定した項数(=次数)までをフィッティングに使っている。
  • 36行目で、フィッティングの結果予測された切片と係数を使って、予測曲線の値を計算している。

異常値がある場合

上記の整然とした指数関数のデータに1つだけ飛び離れた異常値を入れてみる。

先の例に比べて不安定性=曲線の振動の度合いが大きくなっている。

データ数を多くした場合

点の数を10個とし、乱数で擾乱を与えてみる(乱数系列も変えている)。

n_terms=5あたりから、全ての点に何とかフィットさせようと曲線が揺れ始め、特徴量数がデータ数と同じ値となる前後から振動が大きくなっている。

scikit-learn – LinearRegression

概要

scikit-learnLinearRegressionは、最も単純な多重線形回帰モデルを提供する。

モデルの利用方法の概要は以下の手順。

  1. LinearRegressionのクラスをインポートする
  2. モデルのインスタンスを生成する
  3. fit()メソッドに訓練データを与えて学習させる

学習済みのモデルの利用方法は以下の通り。

  • score()メソッドにテストデータを与えて適合度を計算する
  • predict()メソッドに説明変数を与えてターゲットを予測
  • モデルインスタンスのプロパティーからモデルのパラメーターを利用
    • 切片はintercept_、重み係数はcoef_(末尾のアンダースコアに注意)

利用例

配列による場合

以下はscikit-learnのBoston hose pricesデータのうち、2つの特徴量RM(1戸あたり部屋数)とLSTAT(下位層の人口比率)を取り出して、線形回帰のモデルを適用している。

特徴量の一部をとりだすのに、ファンシー・インデックスでリストの要素に2つの変数のインデックスを指定している。また、特徴量データXとターゲットデータyをtrain_test_split()を使って訓練データとテストデータに分けている。

DataFrameによる場合

以下の例では、データセットの本体(data)をpandasのDataFrameとして構成し、2つの特徴量RMとLSTATを指定して取り出している。

利用方法

モデルクラスのインポート

scikit-learn.linear_modelパッケージからLinearRegressionクラスをインポートする。

モデルのインスタンスの生成

LinearRegressionの場合、ハイパーパラメーターの指定はない。

fit_intercept
切片を計算しない場合Falseを指定。デフォルトはTrueで切片も計算されるが、原点を通るべき場合にはFalseを指定する。
normalize
Trueを指定すると、特徴量Xが学習の前に正規化(normalize)される(平均を引いてL2ノルムで割る)。デフォルトはFalsefit_intercept=Falseにセットされた場合は無視される。説明変数を標準化(standardize)する場合はこの引数をFalseにしてsklearn.preprocessing.StandardScalerを使う。
copy_X
Trueを指定するとXはコピーされ、Falseの場合は上書きされる。デフォルトはTrue
n_jobs
計算のジョブの数を指定する。デフォルトはNoneで1に相当。n_targets > 1のときのみ適用される。

 モデルの学習

fit()メソッドに特徴量とターゲットの訓練データを与えてモデルに学習させる(回帰係数を決定する)。

X
特徴量の配列。2次元配列で、各列が各々の説明変数に対応し、行数はデータ数を想定している。変数が1つで1次元配列の時はreshape(-1, 1)かスライス([:, n:n+1])を使って1列の列ベクトルに変換する必要がある。
y
ターゲットの配列で、通常は1変数で1次元配列。

3つ目の引数sample_weightは省略。

適合度の計算

score()メソッドに特徴量とターゲットを与えて適合度を計算する。

戻り値は適合度を示す実数で、回帰計算の決定係数R2で計算される。

(1)    \begin{equation*} R^2 = 1 - \frac{RSS}{TSS} = 1 - \frac{\sum (y_i - \hat{y}_i)^2}{\sum (y_i - \overline{y})^2} \end{equation*}

モデルによる予測

predict()メソッドに特徴量を与えて、ターゲットの予測結果を得る。

ここで特徴量Xは複数のデータセットの2次元配列を想定しており、1組のデータの場合でも2次元配列とする必要がある。

また、結果は複数のデータセットに対する1次元配列で返されるため、ターゲットが1つの場合でも要素数1の1次元配列となる。

切片・係数の利用

fit()メソッドによる学習後、モデルの学習結果として切片と特徴量に対する重み係数を得ることができる。

各々モデル・インスタンスのプロパティーとして保持されており、切片はintercept_で1つの実数、重み係数はcoeff_で特徴量の数と同じ要素数の1次元配列となる(特徴量が1つの場合も要素数1の1次元配列)。

末尾のアンダースコアに注意。

実行例

 

Lasso回帰の理解

定義

Ridge回帰は単純な多重回帰の損失関数に対してL2正則化項を加え、多重共線性に対する正則化を図った。Lasso解析はこれに対してL1正則化項を加えて最小化する(正則化の意味についてはこちら)。

(1)    \begin{align*} L &= \frac{1}{2} \sum_{i=1}^n ( y_i - \hat{y}_i )^2 + \alpha (|w_1| + \cdots + |w_m|) \\ &= \frac{1}{2} \sum_i ( y_i - w_0 - w_1 x_{1i} - \cdots - w_m x_{mi} )^2 + \alpha (|w_1| + \cdots + |w_m|) \end{align*}

L1正則化の意味

準備

L2正則化は各重み係数が全体として小さくなるように制約がかかったが、L1正則化では値がゼロとなる重み係数が発生する。このことを確認する。

係数wを求めるためには損失関数Lを最小化すればよいが、Ridge回帰とは異なりL1正則化項は通常の解析的な微分はできない。

(2)    \begin{align*} \frac{\partial L}{\partial w_k} &= - \sum_i x_{ki} ( y_i - w_0 - w_1 x_{1i} - \cdots - w_m x_{mi} ) + \alpha \frac{\partial |w_k|}{\partial w_k} \\ &= - \sum_i x_{ki}y_i + w_0 \sum_i x_{ki} + \sum_{j \ne k} w_j \sum_i x_{ji} x_{ki} + w_k \sum_i {x_{ki}}^2 + \alpha \frac{\partial |w_k|}{\partial w_k} \\ &= 0 \end{align*}

ここで\frac{\partial |w_k|}{\partial w_k}=|w_k|'と表し、左辺のwk以外に関わる項をMkwkの係数となっている2乗和をSkkと表す。

(3)    \begin{equation*} M_k  + w_k S_{kk} + \alpha |w_k|' = 0 \end{equation*}

場合分け

ここで|wk|’についてはwkの符号によって以下の値をとる。

(4)    \begin{equation*} |w_k|' = \left\{ \begin{array}{rl} -1 & (w_k < 0) \\ 1 & (w_k > 0) \end{array} \end{equation*}

これらを式(3)に適用する。まずwk < 0に対しては

(5)    \begin{gather*} w_k < 0 \quad \rightarrow \quad M_k  + w_k S_{kk} - \alpha = 0 \\ -M_k + \alpha < 0 \quad \rightarrow \quad  w_k = \frac{-M_k + \alpha}{S_{kk}} \end{gather*}

またwk > 0に対しては、

(6)    \begin{gather*} w_k > 0 \quad \rightarrow \quad M_k  + w_k S_{kk} + \alpha = 0 \\ -M_k - \alpha > 0 \quad \rightarrow \quad  w_k = \frac{-M_k - \alpha}{S_{kk}} \end{gather*}

以上をまとめると、

(7)    \begin{equation*} w_k = \left\{ \begin{array}{ll} \dfrac{-M_k - \alpha}{S_{kk}} & (M_k < -\alpha) \\ \\ \dfrac{-M_k + \alpha}{S_{kk}} & (M_k > \alpha) \\ \end{array} \right. \end{equation*}

劣微分の導入

式(7)で−αMkαについては得られていない。Mk → ±αについてそれぞれの側から極限を計算すると0となるのでその間も0でよさそうだが、その保証はない。

ここでこちらのサイトのおかげで”劣微分(subdifferential)”という考え方を知ることができた。|wk|’についてwk = 0では解析的に微分不可能だが、その両側から極限をとった微分係数の範囲の集合を微分係数とするという考え方のようだ。

(8)    \begin{equation*} \frac{d |x|}{dx} = \left\{ \begin{array}{cl} -1 & (x < 0) \\ \left[ -1, 1 \right] & (x = 0) \\ 1 & (x > 0) \end{array} \right. \end{equation*}

そこで、wk = 0に対してこの劣微分を適用してみる。

(9)    \begin{gather*} w_k = 0 \quad \rightarrow \quad M_k +w_k S_{kk} + \alpha \left[ -1, 1 \right] = \left[ M_k - \alpha , M_k + \alpha \right] = 0\\ M_k - \alpha \le 0 \le M_k + \alpha \quad \rightarrow \quad -\alpha \le M_k \le \alpha \quad \rightarrow \quad w_k = 0 \end{gather*}

以上のことから、重みwkについて以下のようになり、−αMkαの範囲ではwk = 0となることがわかる。

(10)    \begin{equation*} w_k = \left\{ \begin{array}{cl} \dfrac{-M_k - \alpha}{S_{kk}} & (M_k < -\alpha \quad)\\ \\ 0 & (-\alpha \le M_k \le \alpha) \\ \\ \dfrac{-M_k + \alpha}{S_{kk}} & (M_k > \alpha) \end{array} \right. \end{equation*}

すなわちL1正則化の場合、ハイパーパラメータαは重み係数の大きさを制限すると同時に重み係数がゼロとなるような効果も持ち、αが大きいほど多くの重み係数がゼロとなりやすい。

参考サイト

本記事をまとめるにあたって、下記サイトが大変参考になったことに感謝したい。

Lassoを数式から実装まで(理論編)Miidas Research