データサイエンティスト(仮)

元素粒子論博士。今はデータサイエンティスト（仮）。

Rでスパースモデリング：Elastic Net回帰についてまとめてみる

R スパースモデリング

導入

回帰モデル構築の際、汎化性能を向上させるために正則化の手法がたびたび用いられます。これは、考えているデータ数に対して特徴量の数が非常に多い場合や、特徴量間に強い相関（多重共線性）がある場合に有効な方法となっています。このような場合に、通常の回帰モデル構築の際に用いられる2乗誤差などの目的関数に加え、L_pノルム（pは正整数）のような正則化項（もしくは罰則項）加えて最適化をおこなうことで先程の問題を解消することができます。こういった正則化項を加えた上でモデルの最適化をおこなう（ = パラメータを推定する）方法を、正則化法といいます。

代表的な正則化法に、Lasso, Ridge, Elastic Net回帰があります。これらは、解釈性も含めた特徴があり、必ずしも高精度のものだからよいわけではない、というのが私の考えです。しかし一方で、{caret}を使ってこの中で最も精度がよいものを採用しました、ということを行っている状況も（私の周りでは）見受けられます。精度向上だけが目的ならこういった立場でも良いのですが、Lasso回帰を始めとするスパースモデリングの強みの一つである、「人間の事前知識/解釈性」をモデルに取り入れることに関しては、必ずしも活かせないことになります。機械が自動で色々いい感じなものを作ってくれる未来はまだまだ先であり、モデリングの際に、どういったことが目的で、それをどう実現していくかをある程度理論的な背景を踏まえた上でおこなっていくことがまだまだ重要だと私は考えています。

まずは代表的な正則化法に関して整理しておきたいと考えたので、それらに関してまとめたものを備忘録的に載せておきます。

使用するデータセット

ボストン近郊の住宅価格のデータを利用します。Rでは{MASS}パッケージを導入すると付随してこのデータセットを使えるようになるため、{MASS}パッケージを利用します。

# ボストン近郊の住宅価格のデータを利用するため（だけ）に導入
if(!require(MASS)){
  install.packages("MASS", quiet = TRUE)
}
require(MASS)


> head(Boston)
     crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat medv
1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98 24.0
2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14 21.6
3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03 34.7
4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94 33.4
5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33 36.2
6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21 28.7

目的変数はmedvで、説明変数候補はそれ以外を想定しています。具体的なモデル構築は、後ほどおこないます。

線形回帰モデル

以下のような線形回帰モデルを仮定します。
$$ \vec{y} = X\vec{\beta} + \vec{\epsilon} \tag{1}
$$
y,\ \vec{\beta},\ \vec{\epsilon}はそれぞれ目的変数、回帰係数（推定するパラメータ）、残差で、Xは計画行列（説明変数をまとめた行列）です。目的変数と特徴量は、それぞれ中心化・標準化されているとします。データ数n、特徴量の数pとすると、y,\ \vec{\epsilon}はn次元ベクトル、\vec{\beta}はp+1次元ベクトル（切片項も含む）、Xはn \times (p+1)行列となります。

回帰モデル構築の問題は、回帰係数\vec{\beta}を求める最適化問題へと帰着されます。最適化のためのよく用いられる目的関数は、以下のような2乗誤差関数です*1。
$$
S\left( \vec{\beta} \right) = \left( \vec{y} - X \vec{\beta} \right)^T \left( \vec{y} - X \vec{\beta} \right) \tag{2}
$$
この2乗誤差関数を最小にするようなパラメータ\vec{\hat{\beta}}は、実測とモデル式の誤差の2乗和を最小にするという意味で、最適なパラメータだと考えます。このパラメータは、目的関数をパラメータで微分したものを0とおいた方程式の解を求めればよく、簡単な計算で以下のようになることがわかります。
$$
\vec{\hat{\beta}} = \left( X^TX \right)^{-1} X^T \vec{y} \tag{3}
$$
式(3)の右辺は観測された量だけで構成されているため、それらを代入すれば回帰係数を求めることができます。

正則化法

先程の線形回帰モデルの議論は、特徴量とデータの関係によってはうまく機能しないことがあります。まず、どういった場合にどのような問題が生じるかを紹介し、それらをどのように解決していくか、という視点で正則化法の説明をしていきます。

多重共線性

式(3)で問題になるのは、多重共線性がある（特徴量間の相関が強い）場合です。

(X^TX)の逆行列を求めるためには、(X^TX)の行列式\text{det}(X^TX)を知る必要があります。しかし、多重共線性がある場合、(X^TX)のなかにほとんど同じ値を持つような列が複数出現することになります。線形代数などで行列式を求める計算を思い出すと、全く同じ値をもった列があると行列式が0になり、行列が特異（正則行列ではない）になっている状態となります。つまり、多重共線性がある場合は、\text{det}(X^TX)が非常に0に近くなります。

(X^TX)の逆行列は1 / \text{det}(X^TX)に比例しているため、\text{det}(X^TX)が0に近づくと逆行列(X^TX)^{-1}は無限大へ発散します。これは、 推定量の共分散\vec{\hat{\beta}} = \sigma ^2 (X^TX)^{-1}(\sigmaは残差が正規分布に従うとした場合の標準偏差)ということを思い出すと、推定量が非常に不安定な値を持つことを意味します。

Ridge回帰 (Hoerl and Kennard, 1970)

多重共線性の問題を回避する方法の一つが、正則化法です。Ridge回帰の場合は、L_2ノルムを加えたものを目的関数とします。
$$
S\left( \vec{\beta} \right) = \left( \vec{y} - X \vec{\beta} \right)^T \left( \vec{y} - X \vec{\beta} \right) + \lambda \vec{\beta}^T \vec{\beta} \tag{4}
$$
\lambda\ (\lambda \geq 0)は正則化の強さを制御する（ハイパー）パラメータです。通常の線形回帰の同様に、解を求めると以下のようになります。
$$
\vec{\hat{\beta} }_{\text{Ridge}} = \left( X^TX + \lambda \vec{1} \right) ^{-1} X^T \vec{y} \tag{5}
$$
\vec{1}は単位行列です。この単位行列（の\lambda倍）の行列の寄与のお陰で、(X^TX)ではほとんど同じ値を持つような列だったものがそうではなくなり、特異ではなくなります。そのため、多重共線性がある状況でも推定が可能になります。線形代数の言葉で言うと、正則化 = 正則行列化というわけです*2。

ただし、\lambdaの値が大きすぎると、推定量の共分散\vec{\hat{\beta}} = \sigma ^2 (X^TX)^{-1}は小さくなりますが、バイアス( (X^TX) / \lambda + \vec{1})^{-1} \vec{\beta}が大きくなるため、適切な\lambdaの調整は必要です。

Lasso回帰 (Tibshirani, 1996)

Ridge回帰によって、多重共線性の問題が回避できることがわかりました。しかし、Ridge回帰には、変数選択ができないという課題があります。式(5)から、大きな\lambdaをとることにより回帰係数全体が縮小される傾向があることがわかります。しかしながら、どれかの係数が完全に0になるといった性質は持っておらず、特徴量の数が多いときに解釈のし易い（特定の変数のみが寄与する）モデルが得ることは難しいという問題があります。

一方、変数選択はAICなどのモデル選択規準を用いて変数の追加・削除を行い最適な変数のセットを見つける問題です。こちらは、解釈性の高いモデル構築をできる可能性がありますが、多重共線性の問題がある場合は必ずしも安定な推定量が得られるわけではありません。

Lasso回帰は、Ridge回帰と変数選択の良いとこ取りをしたようなモデルになっています。この場合、L_1ノルムを加えたものを目的関数とします。
$$
S\left( \vec{\beta} \right) = \left( \vec{y} - X \vec{\beta} \right)^T \left( \vec{y} - X \vec{\beta} \right) + \lambda \sum _{j = 1}^p \left| \beta _j \right| \tag{6}
$$

式(6)は原点で不連続なため、通常の線形回帰のように厳密に解を求めることはできません。数値的に解く必要があります。数値計算の方法はいくつか種類があり、よく用いられる手法に座標降下法と呼ばれるものがあります。こういった何らかのアルゴリズムを用いて推定した結果を\vec{\hat{\beta}}_{\text{LASSO}}とします。
$$
\vec{\hat{\beta}}_{\text{LASSO}} = \text{argmin}_{\beta} \left[ \left( \vec{y} - X \vec{\beta} \right)^T \left( \vec{y} - X \vec{\beta} \right) + \lambda \sum _{j = 1}^p \left| \beta _j \right| \right] \tag{7}
$$

Lasso回帰は、パラメータ推定と変数選択を同時に行うことができ、スパースなモデルを構築することが可能です。特徴量が2つの場合の直感的なイメージ図は以下になります。正則化法は、制約条件付きの2乗誤差の最小化問題に置き換えることができます。ピンクの領域は、正則化項の大きさを固定した場合に許されるパラメータの領域となっています。最小2乗法による解がこのピンクの領域よりも外側にある場合、制約条件込みで目的関数を最小にするのは、最小2乗法による解を中心とした等高線とピンクの領域が交わる部分（図の星印）となります。Ridge回帰では、制約条件が円状であるため、回帰係数が全体的に小さくなる領域が解になりやすく、Lasso回帰では、尖っている部分（微分不可能な点で、特定の回帰係数の大きさが0となっている）が解になりやすくなります。そのため、Lasso回帰ではいくつかの係数が0となる解が推定されやすくなっています。
f:id:tekenuko:20171118170813j:plain

Lasso回帰の注意点

Lasso回帰は、変数選択の一致性が必ずしも保証されていないため、「正しい」変数選択が行われていない可能性があります。変数選択の一致性を担保したい場合は、Adaptive Lassoといった理論保証があるスパースモデリングの手法を使う必要があります。
tekenuko.hatenablog.com

Lasso回帰の欠点

Lasso回帰は、スパースなモデルを構築するためには非常に有用ですが、一方で以下のような欠点もあります。

  • データ数n, 特徴量の数をpとした場合に、p > nのときは高々n個の変数までしか選択できない
  • 相関の高い変数群がある場合、Lassoはその中の1つしか変数として選択できない

一つ目は、本当に重要な変数（0でないと推定したい変数）がnよりも多い場合に、重要な変数を選択しきれない、いう意味で問題となります。二つ目は、Grouping効果を考慮できないということを意味しています。遺伝子データのような場合、事前知識から似ている特徴量が複数あり、しかもそれらを変数として選択したい状況もあるかと思います。このような場合、Lasso回帰では選択される変数が似ている特徴量の中のどれか一つとなってしまいます。そうすると、モデル作成者が達成したい目的はLasso回帰では達成できなくなります。

Elastic Net回帰 (Zou and Hastie, 2005)

上記のようなLasso回帰の欠点を解決すべく考案されたものが、Elastic Net回帰です。Elastic Net回帰は、正則化項として、L_1ノルムとL_2ノルムの和を考えます。
$$
S\left( \vec{\beta} \right) = \left( \vec{y} - X \vec{\beta} \right)^T \left( \vec{y} - X \vec{\beta} \right) + \lambda \sum _{j = 1}^p \left\{ \alpha \left| \beta _j \right| + (1 - \alpha) \beta _j ^2 \right\} \tag{8}
$$
\lambda \ (\lambda \geq 0)は正則化項の大きさを制御するハイパーパラメータで、\alpha \ (0 \leq \alpha \leq 1)はL_1ノルムとL_2ノルムの相対的な寄与を調整するハイパーパラメータになります。\alpha = 0の場合がRidge、\alpha = 1の場合がLasso回帰に帰着します。

以下、Lasso回帰の欠点をElastic Net回帰がどのように解決しているかを見ていきます。

データ数n, 特徴量の数をpとした場合に、p > nのときは高々n個の変数までしか選択できない

こちらは、Elastic Net回帰の場合の推定量の計算を具体的に導出することでわかってきます。今、式(8)に登場する量たちを
$$
X^* = (1 + (1 - \alpha )\lambda) \begin{pmatrix}
X \\
\sqrt{(1 - \alpha) \lambda} \vec{1}_p \\
\end{pmatrix}, \ \ \
\vec{y}^* = \begin{pmatrix}
\vec{y} \\
\vec{0}_p \\
\end{pmatrix}, \\
\ \gamma = \alpha \lambda / \sqrt{1 + (1 - \alpha) \lambda} , \ \ \ \vec{\beta}^* = \sqrt{1 + (1 - \alpha) \lambda} \vec{\beta} \tag{9}
$$
とおくと、式(8)は
$$
S\left( \vec{\beta}^* \right) = \left( \vec{y}^* - X^* \vec{\beta}^* \right)^T \left( \vec{y}^* - X^* \vec{\beta}^* \right) + \gamma \sum _{j = 1}^p \left| \beta _j \right| \tag{10}
$$
とLasso回帰の形で書き直せます。そのため、式(10)の問題はLasso回帰の問題として解くことができ、その解を\vec{\hat{\beta}}^*とすると、Elasitc Net回帰に関する推定量は
$$
\vec{\hat{\beta }} _{\text{Elastic Net}} = \frac{1}{ \sqrt{1 + (1 - \alpha) \lambda} } \vec{\hat{\beta}}^* \tag{11}
$$
と求めることができます。今、X^*は(n + p)\times p行列となっているため、この行列のランクはpです。そのため、式(11)の推定量として、最高でpの成分が0でないと推定されます。これは、Lasso回帰のp > nのときは高々n個の変数までしか選択できない（0でない回帰係数の最大数はn）の欠点を解消するものになっています。

ここで、X^*というXよりもランクが大きくなりうる行列を定義していますが、これはRidge回帰の性質である疑似データ生成がもとになっています。つまるところ、Lasso回帰の欠点をRidge回帰の性質で補うという構造になっています。

相関の高い変数群がある場合、Lassoはその中の1つしか変数として選択できない

こちらは、原論文に証明があります。今\hat{\beta} _iをi番目のElastic Net回帰の推定量とし、ある(i, j)に関して\hat{\beta} _i \times \hat{\beta} _j > 0が成立しているとします。これは、これらの係数の符号が一致している状況になります。このとき、以下のような量を定義します。
$$
D(\alpha , \lambda) = \frac{1}{\sum _{i = 1}^n |y_i|} \left| \hat{\beta} _i -\hat{\beta} _j\right| \tag{12}
$$
詳しい導出は省略しますが、今考えている回帰係数に対応する特徴量間の相関係数を\rhoとすると、式(12)は以下のように評価することができます。
$$
D(\alpha , \lambda) \leq \frac{1}{(1 - \alpha) \lambda} \sqrt{2 (1 - \rho)} \tag{13}
$$
これは、\rho = 1（相関係数が1）の場合に\hat{\beta} _i = \hat{\beta} _j となることを意味しています。よって、Lasso回帰の場合と異なり、Elastic Net回帰では、相関の高い係数はどちらも0になるか、0にならないなら同じ値になる、といったGrouping効果を持つことがわかりました。

glmnetによる実験

以上のことから、それぞれの正則化法のおおざっぱな特徴は以下のようになります。

  • Ridge回帰：回帰係数全体が縮小推定される
  • Lasso回帰：相関の高い特徴量群がある場合、そのうちの一つが選択される（対応する回帰係数以外が0と推定される）
  • Elasitic Net回帰：相関の高い特徴量は、回帰係数が両方0になる、もしくは両方同じ値になる

これらの特徴を、ある程度恣意的な状況でモデルを構築して確かめてみることにします。使用データはボストン近郊の住宅価格のデータで、議論を簡単にするために変数をある程度絞って議論します。絞る材料とするために、いったんglmでモデルを構築し、stepwise法で変数選択をしてみます。

# glmでモデル作成
mod.full <- glm(medv ~., data = Boston)
# AICを基準にstepwise
mod.step <- stepAIC(mod.full, direction = "both", trace = 0)
# summaryを出力
summary(mod.step)

# 係数の部分だけ記載
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  36.341145   5.067492   7.171 2.73e-12 ***
crim         -0.108413   0.032779  -3.307 0.001010 ** 
zn            0.045845   0.013523   3.390 0.000754 ***
chas          2.718716   0.854240   3.183 0.001551 ** 
nox         -17.376023   3.535243  -4.915 1.21e-06 ***
rm            3.801579   0.406316   9.356  < 2e-16 ***
dis          -1.492711   0.185731  -8.037 6.84e-15 ***
rad           0.299608   0.063402   4.726 3.00e-06 ***
tax          -0.011778   0.003372  -3.493 0.000521 ***
ptratio      -0.946525   0.129066  -7.334 9.24e-13 ***
black         0.009291   0.002674   3.475 0.000557 ***
lstat        -0.522553   0.047424 -11.019  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

結果を見ると、相対的にはrm, dis, ptratio, lstatの係数が有意であるようです。というわけで、これら4つの変数に絞り、さらにrmに対しては全く同じ変数のコピーを新たに加えたデータを作成します。

# rmのコピー
rm.dummy <- data.frame(rm_dummy = Boston$rm)
# new data.frame
Boston.new <- cbind(Boston[, c(6, 8, 11, 13, 14)], rm.dummy)

以下では、rmのコピー変数の回帰係数の解パスをモデルごとを見ていくことで、それぞれの正則化法の違いを見ていくことにします。その際に、glmnetに関する便利な補助ツールが幾つかあることが最近わかったので、活用していきます。具体的には、glmnetでもformula的に書くことができるようにするglmnetUtilsと、解パスをggplot2の形式で可視化してくれるggfortulyというものです。こちらに関しては、hoxo_mさんが最近紹介してくださりました。
qiita.com

Ridge回帰

先程恣意的に作成したデータで、Ridge回帰をし、解パスを見てみましょう。

# 必要なライブラリ
if(!require(glmnet)){
  install.packages("glmnet", quiet = TRUE)
}
require(glmnet)
require(glmnetUtils)
if(!require(githubinstall)){
  install.packages("githubinstall", quiet = TRUE) 
}
require(githubinstall)
if(!require("ggfortify")){
  githubinstall("ggfortify")
}
require(ggfortify)

# Ridge回帰（glmnetUtilsを併用）
Ridge <- glmnet(medv ~ ., data = Boston.new, alpha = 0)
# ggfortifyで可視化
autoplot(Ridge, xvar = "lambda")

f:id:tekenuko:20171118212000p:plain
（途中0を横切っているのもありますが…）正則化パラメータ\lambdaを大きくしていくと、係数の大きさが全体的に0に縮小していく傾向が見えます。また、rmとそのコピーrm_dummyはほぼほぼ同じような動きをしています。

Lasso回帰

Lasso回帰の場合にも、同じように解パスを見てみます。

# Lasso回帰（glmnetUtilsを併用）
Lasso <- glmnet(medv ~ ., data = Boston.new)
# ggfortifyで可視化
autoplot(Lasso, xvar = "lambda")

f:id:tekenuko:20171118212501p:plain
Ridge回帰の場合と比較し、rmとrm_dummyの振る舞いが明らかに違います。rm_dummyの回帰係数が常に0と推定されています。これは、Lasso回帰の「相関の高い特徴量群がある場合、そのうちの一つが選択される」という性質に対応しています。

Elastic Net回帰

# Elastic Net回帰（glmnetUtilsを併用）
ElasticNet <- glmnet(medv ~ ., data = Boston.new, alpha = 0.5)
# ggfortifyで可視化
autoplot(ElasticNet, xvar = "lambda")

f:id:tekenuko:20171118212916p:plain
Lasso回帰の欠点であった、Grouping効果が反映されています。少し見づらいですが、rmとrm_dummyの係数が0でなくなるタイミングが同じになっており、いったん0でないと係数された場合は、両者の回帰係数が似たような値になっています。

精度について

今回は主題ではないのであえて触れませんでしたが、実際に精度検証をしてみると、今回のデータではglm w/ stepwise > Elastic Net, Lasso > Ridgeの順で決定係数が大きかったです。例によってテストデータを含めた検証をしっかりやっているわけではないので、ここは参考程度で。実務の際はちゃんと汎化性能・解釈性などを踏まえて議論します。

まとめ

今回は、ある程度理論的な性質をもとに代表的な正則化法をまとめてみました。かなり恣意的な状況のデータではありましたが、相関が高い特徴量に対する実際の挙動の違いも確認しました。精度追求が要件の場合には今回の議論はあまり意味をなさないと思いますが、一方で、業務知識などである程度取り入れたい考えがあるときには、それに応じてモデルを構築するヒントの一助となるのではないかと思います。

*1:最尤法の枠組みで、正規分布に従う残差を仮定した場合でも、同じ目的関数となります。

*2:正則行列は、逆行列が存在するような正方行列のことです。行列式が0の場合、逆行列が存在しないのですが、正則化項を導入することで、逆行列が存在するような行列へと補正をかけています。

Rでスパースモデリング：Adaptive Lasso

R スパースモデリング

導入

スパース推定の代表的な手法として、Lassoがあります。様々なシーンで活用されているLassoですが、Lassoは変数選択の一致性が保証されないという欠点があります。Adaptive Lassoは、その欠点を補う形で提唱されている手法となっています。こちらは、ある条件のもとで変数選択の一致性が保証*1されており、数理統計学的により好ましい性質を持っています。

このAdaptive Lassoですが、Rでは{glmnet}以外のパッケージを使わないと簡単にできないとかなりの期間勘違いをしてました。そんな中、以下の記事を最近見かけまして、冷静に考えたら{glmnet}でも表現できるよなあと猛省した次第です。

RPubs - Adaptive LASSO Examples

以上の経緯から、挙動を確かめておこうという考えのもと、メモがてらAdaptive Lassoの紹介をしようと思います。

補足

上記の一致性の議論は、特徴量の数を固定して、データ数に対して無限大の極限を取った場合です。特徴量の数をデータ数よりも常に大きくとるようにした場合の極限の振る舞いも議論されており、例えば
http://www.keihirose.com/material/392-399_hirose.pdf
に載っております。

使用するデータセット

ボストン近郊の住宅価格のデータを利用します。Rでは{MASS}パッケージを導入すると付随してこのデータセットを使えるようになるため、{MASS}パッケージを利用します。

# ボストン近郊の住宅価格のデータを利用するため（だけ）に導入
if(!require(MASS)){
  install.packages("MASS", quiet = TRUE)
}
require(MASS)

# 目的変数、説明変数
y <- as.matrix(Boston[, 14])
X <- as.matrix(Boston[, -14])

目的変数はmedvで、説明変数はそれ以外です。また、{glmnet}にデータを投入する際は

  • 目的変数と説明変数を分ける
  • それぞれmatrixにする
  • 変数の中にfactorがあってはならない

といった制約があるので、それらを意識しています。

また、今回は予測に重きを置いているわけではなく、回帰係数の振る舞いに興味があるので、予測のために学習・検証用データに分けることはせず、全データを用いて推定を行います。

Lasso

以下のような回帰モデルを仮定します。
$$ \vec{y} = X\vec{\beta} + \vec{\epsilon}$$
y,\ \vec{\beta},\ \vec{\epsilon}はそれぞれ目的変数、回帰係数（推定するパラメータ）、残差で、Xは計画行列（説明変数をまとめた行列）です。データ数n、特徴量の数pとすると、y,\ \vec{\epsilon}はn次元ベクトル、\vec{\beta}はp+1次元ベクトル（切片項も含む）、Xはn \times (p+1)行列となります。

パラメータ推定をする際の目的関数は、二乗誤差を仮定します。Lasso(Least Absolute Shrinkage and Selection Operator)では、正則化項として L_1 ノルムを目的関数に加えます。
$$ L(\vec{\beta}) = \left\| \vec{y} - X \vec{\beta} \right\| ^2 + \lambda \| \vec{\beta} \|_1 $$
\lambdaは正則化項（もしくは罰則項）と呼ばれ、もとの二乗誤差に対して正則化の強さを制御するようなハイパーパラメータです。

上の式は原点で不連続なため、通常の線形回帰のように厳密に解を求めることはできません。数値的に解く必要があります。数値計算の方法はいくつか種類があり、よく用いられる手法に座標降下法と呼ばれるものがあります。こういった何らかのアルゴリズムを用いて推定した結果を\vec{\beta}_{\rm{Lasso}}とします。
$$ \vec{\beta}_{\rm{Lasso}} = \text{argmin}_\beta \left[ \left\| \vec{y} - X \vec{\beta} \right\| ^2 + \lambda \| \vec{\beta} \|_1 \right] $$


Lassoの特徴は、パラメータ推定と変数選択が同時にできることです。変数選択は、例えば相関の高い２つの変数があったとき、片方の変数の回帰係数が0と推定され、もう片方は0でない値となる、といったことが起きます。こういった状況のもと、推定結果が0でなかったものが変数として選択された、と考えます。この性質は、LassoがAICやBICといったモデル選択規準を用いてステップワイズ法を考える、といった変数選択のメカニズムが背後にあることに起因しています。こういった状況は、正則化項が離散的な関数になっていると読み取れるのですが、この項をうまく凸関数になるように連続極限をとった場合がLassoになります*2。

他の正則化項として代表的なものとしては、L_2ノルムを加えるRidge、L_1ノルムとL_2ノルムの線形結合を加えるElastic netなどがあります。これら正則化やアルゴリズムのより詳しい話は、例えば以下の記事などを参考としてください。
qiita.com

追記

argminの部分、投稿時はargmaxとタイポしてました。ご指摘いただいたynakahashiさんに感謝です！

Lassoの解パス

ハイパーパラメータである正則化項\lambdaを変化させると、推定結果も変わります。この変化の移り変わりを可視化したものを解パス（Solution Path）といいます。ここでは、ボストン近郊の住宅価格のデータに対してLassoを適用して、解パスがどうなるかを見てみます。

# モデリング
Lasso <- glmnet(x = X, y = y,
                alpha = 1)
# 正則化項を変化させたときの回帰係数の推移（解パス）
plt_res(Lasso)

f:id:tekenuko:20171102005446p:plain

ここで、plt_resはglmnet関数の出力をうまくggplot2で可視化する自作関数です。

# Library
if(!require(ggplot2)){
  install.packages("ggplot2", quiet = TRUE)
}
require(ggplot2)
if(!require(tidyr)){
  install.packages("tidyr", quiet = TRUE)
}
require(tidyr)

plt_res <- function(dgMatrix){

  # Convert to data.frame
  df_X <- as.data.frame(t(as.matrix(dgMatrix$beta)))
  df_loglambda <- dgMatrix$lambda
  df <- cbind(df_loglambda, df_X)
  # melt
  df_melt <- df %>%
    gather(key = df_loglambda, value = value)
  colnames(df_melt) <- c("lambda", "coefficient", "beta")
  # plot
  plt <- ggplot(df_melt) + 
    geom_line(aes(x = lambda, y = beta, group = coefficient, col = coefficient)) + 
    scale_x_log10()
  
  return(plt)
  
}

図を見ると、正則化項の寄与が小さい（\lambdaが小さい）ときは、0でない値に推定される回帰係数が多く、正則化項の寄与が大きくなると0と推定される係数が増えていくといった振る舞いになっています。このように、正則化項の寄与を変化させると、推定される係数も変化していきます。

回帰係数

上で解パスを見ましたが、ハイパーパラメータはどの値を取るのがよいでしょうか。上記の解パスの各\lambdaごとにRMSEなどが評価できるはずなので、その評価指標が最も良い点を選ぶ、というのも一つです。ここでは、より汎用的に評価するために、cross-validationを考えます。{glmnet}では、cv.glmnet関数を用いると簡潔に計算できます。

fit.Lasso.cv <- cv.glmnet(x = X, y = y,
                      nfolds = 10,
                      alpha = 1, 
                      standardize = TRUE)

nfoldはデータを10分割し、そのうちの一つをテストデータとする、といったサンプルを10パターン生成することを意味しています。alphaはL_1ノルムとL_2ノルムの相対的な寄与の大きさをコントロールしているパラメータでalpha=1のときがLasso、alpha=0のときがRidgeとなります。standardizeはモデリングの際にデータを標準化するかどうかのオプションです。

この結果、plot関数を用いると以下のように可視化されます。

plot(fit.Lasso.cv)

f:id:tekenuko:20171102010846p:plain
これは、\lambdaとRMSEの関係をプロットしたもので、エラーバーはデータ分割の仕方が起源として発生しています。この図を見ると、\lambdaは非常に小さい値が好まれるようです。

最もRMSEが小さい\lambdaの場合の推定結果は、以下になります。

coef(fit.Lasso.cv, s = fit.Lasso.cv$lambda.min)

# 出力
14 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept)  34.880894915
crim         -0.100714832
zn            0.042486737
indus         .          
chas          2.693903097
nox         -16.562664331
rm            3.851646315
age           .          
dis          -1.419168850
rad           0.263725830
tax          -0.010286456
ptratio      -0.933927773
black         0.009089735
lstat        -0.522521473

indusとageの係数が0と推定されています。

Adaptive Lasso

Adaptive Lassoは、L_1ノルムを目的関数に加える点はLassoと同じですが、重みが付く形になります。
$$ L(\vec{\beta}) = \left\| \vec{y} - X \vec{\beta} \right\| ^2 + \lambda \sum _{i = 1}^p \omega_i \left| \beta _i \right| $$
ここで、重みベクトル\vec{\omega} = (\omega _1, \cdots , \omega _p)は以下になります。
$$\vec{\omega} = 1 / \left|\vec{\hat{\beta}}\right|^\gamma $$
\vec{\hat{\beta}}はある一致推定量で、一致推定量の各成分の絶対値の\gammaの逆数が\vec{\omega}の各成分となっています。また\gamma > 0です。

一致推定量は、真の値との差がO(1/\sqrt{n})以下程度の精度が担保されるものなら何でもよいです。多くの場合、OLSの結果やRidgeの結果を使うようです。今回は、Ridge回帰の結果を用いることにします。
（2018-03追記 : Ridgeの結果を使うとオラクル性が満たされないかもしれません。調査中。）

Ridge回帰により何らかの一致推定量を得る

必ずしもcross-validationまでやる必要はありませんが、10-fold cross-validationにより、最もRMSEが小さくなる場合の回帰係数を取得します。

# Ridge回帰
fit.ridge.cv <- cv.glmnet(x = X, y = y,
                          type.measure = "mse",
                          alpha = 0, 
                          nfolds = 10, 
                          standardize = TRUE)
# 最もRMSEが小さくなる正則化項の場合の回帰係数
best_ridge_coef <- as.numeric(coef(fit.ridge.cv, s = fit.ridge.cv$lambda.min))[-1]

Adaptive Lassoの解パス

上で取得した係数を反映します。Adaptive Lassoは、glmnet関数のオプションpenalty.factorに1 / abs(best_ridge_coef))を代入することで実現できます。このオプションを認識していなかったので、glmnetでは簡単にAdaptive Lassoを書けないと思っていました。。解パスを表示すると、以下のようになります。

# モデリング
AdaLasso <- glmnet(x = X, y = y,
                  alpha = 1,
                  penalty.factor = 1 / abs(best_ridge_coef))
# 解パスを表示
plt_res(AdaLasso)

f:id:tekenuko:20171102013024p:plain

（2018-03追記 : 若干コードの見直しが必要かもしれません。調査中。）
線が上下している係数もあるので、実はあまり推定がうまくいっていないのかも、という懸念はありますが、Lassoの場合と比較して、だいぶ振る舞いが変わりました。どのタイミングでどの係数が0になるか、線の形状、色々変化しています。動く\lambdaの範囲は、実質的にRidgeの結果で割りこんだものを考えているので変化するのは理解できます。他の変化の仕方は明確に理由を説明するのは難しそうです。

回帰係数

Adaptive Lassoの場合も係数を表示しておきましょう。cross-validationした場合の図は省略します。

fit.AdaLasso.CV <- cv.glmnet(x = X, y = y,
                        type.measure = "mse",
                        nfolds = 10,
                        alpha = 1,
                        penalty.factor = 1 / abs(best_ridge_coef))
coef(fit.AdaLasso.CV, s = fit.AdaLasso.CV$lambda.min)

# 出力
14 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept)  38.501580556
crim         -0.098114529
zn            0.022165389
indus        -0.016969903
chas          3.129489627
nox         -20.469852957
rm            3.948669580
age           .          
dis          -1.345793859
rad           0.097441871
tax           .          
ptratio      -1.020675606
black         0.001601036
lstat        -0.546279678

0となる係数の組み合わせ、および係数の値が変化しています。このように、回帰係数のパターンもLassoと比較して変わるようです。

精度比較

精度が全てというわけではありませんが、LassoとAdaptive Lassoの両者で精度比較を行っておきます。評価指標はRMSEと決定係数R^2ですが、上述の通り、今回は全データを推定に使ったので、学習データに関するものの評価になっています*3。RMSEは{Metrics}のメソッドで計算し、決定係数R^2は予め関数を定義しておきます。

# 決定係数
r_squared <- function(y, ypred) {
  ybar <- mean(y)
  ## Total SS
  ss_tot <- sum((y - ybar)^2)
  ## Residual SS
  ss_res <- sum((y - ypred)^2)
  ## R^2 = 1 - ss_res/ ss_tot
  1 - (ss_res / ss_tot)
}

# LassoとAdaptive LassoそれぞれRMSEとR^2を計算
c("RMSE(Lasso) : ", round(rmse(y, pred_Lasso), 4))
c("R^2(Lasso) : ", round(r_squared(y, pred_Lasso), 4))
c("RMSE(Lasso) : ", round(rmse(y, pred_AdaLasso), 4))
c("R^2(Adaptive Lasso) : ", round(r_squared(y, pred_AdaLasso), 4))

# 出力
"RMSE(Lasso) : " "5.0299"   
"R^2(Lasso) : " "0.7003" 

"RMSE(Adaptive Lasso) : " "5.0435"
"R^2(Adaptive Lasso) : " "0.6987"

Adaptive Lassoの方がRMSEが大きく、決定係数が小さいという結果になってしまいました。変数選択をより正確にしたことの代償なのか…。今回はテストデータも含めて検証をしておりませんし、何よりデータに依存する話かもしれないので、ここの比較結果はあくまで参考程度にお願いします。

まとめ

今回は、数理統計学的に性質のよいとされるAdaptive Lassoのglmnetでの実装例を紹介しました。今まで壮大に勘違いをしていたのですが、案外シンプルにglmnetで書けることがわかったので、変数選択の一致性を気にする場合は積極的に活用していこうと思います。

*1:他にも推定量が漸近正規性をもつという性質も合わせて、オラクルプロパティを持つ、といいます。

*2:このあたりの議論は、統計数理研究所の公開講座「スパース推定」でされております。私はこの話を聴いたとき非常に感銘をうけました。

*3:実務接続を意識して比較する場合は、予測に関する評価もきちんと行わなければなりません。

トポロジカルデータアナリシス：TDAパッケージを使ってみる

トポロジカルデータアナリシス R TDA

導入

とあることがきっかけで、とっても久しぶりにRでTDA（トポロジカルデータアナリシス）をしてみました。だいぶいろんなことを忘れていたので、単純な例を使ったメモを残しておきます。

トポロジカルデータアナリシスとは

とてもざっくりいうと、位相幾何学という数学の知見をつかって、データから「形」の情報を抽出するような手法になっています。

導入は、過去の記事にも載っていますので、こちらも参照していただけるとうれしいです。
tekenuko.hatenablog.com

「形」の情報の抽出のしかたですが、データ点のまわりにある半径の球をとり、その球たちの重なりによって点の間に線を引く、といった方法でデータ点から図形を見立てます。ただし、データに付随するノイズなどにロバストな形で図形的な情報を抜き出すために、半径を連続的に変化させた場合の構造の変化を見るといったテクニックを用います。

これらの簡単な解説は、一応過去の記事でも紹介しています（わかりづらかったらすいません）。
tekenuko.hatenablog.com

半径を連続的に変化させた場合、図形が生まれて、亡くなる…という動きが起こります。その一連を動きをプロットしたものが、パーシステント図と呼ばれます。今回、具体例を用いてパーシステント図がどうなるかを説明したいと思います。

使用パッケージ

ここでは、Rでトポロジカルデータアナリシスを行うためのパッケージである{TDA}の関数を使います。以下のようにしてインストールします。Rcppなど、いくつか依存パッケージがあり、それらのバージョンなどでエラーが出ることがあります。その場合は、依存パッケージを入れ直すなどするなどの対処をする必要があります（うろ覚え）。

# TDAパッケージ
if(!require(TDA)){
  install.packages("TDA", quiet = TRUE)
}
require(TDA)

使い方（マニュアル）はこちら：
https://cran.r-project.org/web/packages/TDA/vignettes/article.pdf

日本語の解説だと、ALBERTさんの記事が参考になります。
blog.albert2005.co.jp


サンプルデータ

大きな円と小さな円が繋がっている＋微小なノイズがのっているような人工データを生成します。ここでは、{TDA}パッケージにある、半径1の円上に沿ってランダムサンプリングを行う関数circleUnif()をうまく組み合わせます。

X1 <- circleUnif(100) * 2 + rnorm(100, 0, 0.03)
X2 <- circleUnif(50) * 0.3 + c(1.65,1.65) + + rnorm(50, 0, 0.02)
X <- rbind(X1, X2)
X <- data.frame(x = X[, 1], y = X[, 2])

可視化すると以下のようになります。

ggplot(X) + 
  geom_point(aes(x = x, y = y))

f:id:tekenuko:20171101001723p:plain

パーシステント図

パーシステント図を描きます。図形の構築法は、ヴィートリス・リップス複体を利用します。ヴィートリス・リップス複体を含めた構築法の概要は、過去記事を参照いただけると幸いです。
tekenuko.hatenablog.com

{TDA}パッケージでは、以下のような関数を使って実現します。

# パーシステント図の大元（パーシステントホモロジー）を計算
Diag_Cir <- ripsDiag(X = X, maxdimension = 1, maxscale = 5)

Xはデータ点、maxdimensionは何次元の三角形の移り変わりを見るかを指定するオプションです。今回は1としており、穴が空いているかの移り変わりを見ています。maxscaleは図を描くときの軸の最大値の値を指定するオプションです。パーシステント図を作るために動かす半径のパラメータは、本来無限大まで取れますが、変化は実質どこかで止まるはずで、その止まるパラメータ点をmaxscaleとしてスケーリングします。

可視化すると以下のようになります。{TDA}に内蔵されているplot関数を使うのが基本ですが、ここでは{ggplot2}で可視化を試みました。

# ggplot2で描くためにripsDiagの結果を整形
dimension_Cir <- Diag_Cir$diagram[,1]
birth_Cir <- Diag_Cir$diagram[,2]
death_Cir <- Diag_Cir$diagram[,3]
Diag_Cir_DF <- data.frame(cbind(Dimension = dimension_Cir, 
                              Birth = birth_Cir,
                              Death = death_Cir))

# plot
ggplot(Diag_Cir_DF) + 
  geom_point(aes(x = Birth, y = Death, colour = as.factor(Dimension)), size = 4) + 
  xlim(0,5) + ylim(0,5) + 
  geom_abline(slope = 1) + 
  labs(color = "Dimension") + 
  scale_color_manual(name="Dimension",
                     values=c("0"="black","1"="#F8766D"))

f:id:tekenuko:20171101002745p:plain

赤い点が穴に対応しており、対角線から距離が離れているほど半径の変化で安定的な穴であることを表しています。今回、赤い点は相対的に安定的なものが2点、不安定なものが1点出現しています。前者はデータ点を大きな円（対角線から大きく離れている点）と小さな円（対角線からちょっと離れている点）に対応しており、後者はノイズ起源でたまたま穴の空いた図形が生まれてすぐ消えていったものだと考えられます。このように、パーシステント図を利用すると、データ点の形に付随した量を抽出することができます。

Next Step

今回は、非常に簡単な例ではありますが、トポロジカルデータアナリシスのRでの実行例を紹介しました。このパッケージには他にも色々な機能が搭載されていますので、気になる方はマニュアルを見つつ色々お試しください。本ブログでも利用例を定期的に紹介していきたいと思っています（たぶん）。

Pythonでデータ分析：Prophetを使ってビットコインの予測（笑）をやってみる

Python Prophet

導入

直近、これといって緊急の業務がなく、「自分の時間だ何勉強しようかなー」とPyStanとかをいじっていた矢先、「暇なら技術調査やってよ、Deep Learning的な何かとか」というお達しがきました。あいにく私は天邪鬼なので、２つ返事をして気になっていたけど触っていなかったProphetを調べることにしたのでした。

注：仕事はちゃんとしました（Seq2Seqの論文や書籍見て簡単な実装をしました）。

Prophet

Facebookが出した時系列予測のツールです。
facebook.github.io
すでに様々な方が紹介をしたり、Contributeしていたりするので、釈迦に説法感がありますが、このツールの良い点は、簡単に（分析の専門知識がなくても）ある程度それらしい予測値を出してくれるところです。ビジネス側でデータを活用したい場合や、分析者でもいったん簡単にデータから言えることを見てみる、といった場合に便利そうです。

というわけで、分析者としても抑えておきたいと前々から思っていたので、いい機会だと思って少し動かしてみました。一通りサンプルを見たあと、題材として何かあるかなと思って探していたところ、ビットコインのデータをPythonから取得できるライブラリを見つけたので、ビットコインを題材とすることにしました。

参考

  • Giihub（メソッドの引数とか見るのに使った）

github.com

  • ホクソエムさんの資料（最初に超見た）

d.hatena.ne.jp

darden.hatenablog.com

データ取得

以下のコードを使います。
github.com
インストール方法はサイトに書かれています。pipが使えるなら

$ pip install https://github.com/s4w3d0ff/python-poloniex/archive/v0.4.6.zip

とすれば使用できる状態になるはずです。

データですが、以下のコードを実行すると、ビットコイン（US）のデータ（日次）を過去500日分取得してくれます。ざくっと紹介すると、polo.returnChartDataでデータを辞書型で取得しています。periodでどれくらいの間隔のデータかを指定（5分足なども指定できます）しています。ただ、時間がUNIX timeからの経過時間（秒）になっているので、後でよしなに変換しています。

# numpyやpandas
import numpy as np
import pandas as pd
from pandas import DataFrame, Series

# APIに関係するライブラリ
import poloniex
import time
import datetime

def getDataPoloniex():
    polo = poloniex.Poloniex()
    polo.timeout = 2
    chartUSDT_BTC = polo.returnChartData('USDT_BTC', period=polo.DAY, start=time.time() - polo.DAY * 500, end=time.time())
    tmpDate = [chartUSDT_BTC[i]['date'] for i in range(len(chartUSDT_BTC))]
    date = [datetime.datetime.fromtimestamp(tmpDate[i]).date() for i in range(len(tmpDate))]
    data = [float(chartUSDT_BTC[i]['open']) for i in range(len(chartUSDT_BTC))]
    return [date, data]

# データ取得
dat = getDataPoloniex()
# DataFrameへ格納
df = DataFrame(dat).T
df.columns = ['ds', 'y']

ちなみに、プロットすると以下のようになっています。

import matplotlib.pyplot as plt
%matplotlib inline
df.plot(x = df.columns[0], figsize = (15, 8))

f:id:tekenuko:20171018000038p:plain
途中で下落したりすれど、急成長しているような系列になっています。

Prophetによるモデル作成

ProphetはPythonの場合はpipで導入することができます。

pip install fbprophet

PyStanとか色々依存ライブラリがあるようなんですが、そういったものたちは先に入れとくといいのかもしれません。自分はすべて先に入っていたので比較ができませんが。

モデル構築はsklearnライクにできます。

mod = Prophet()
mod.fit(df)

非線形トレンドを入れたり、変化点を考慮したりといろいろオプションはあるのですが、ここではよしなにできる感を演出するのに何も指定しないでおきます。線形トレンド以上の成長をしていそうな系列なので、ホントはちゃんと処理したいですけどね。

予測に関しては、make_future_dataframeとpredictで新たに生成した系列（100日分）に予測値を格納します。

future_df = mod.make_future_dataframe(periods = 100)
forecast_df = mod.predict(future_df)

簡単なプロットはProphetじたいにplotというメソッドがあり、それでよしなな図が見れます。

mod.plot(forecast_df)
plt.show()

f:id:tekenuko:20171018000907p:plain

まあ、線形トレンド的な予測がされています。もっと長期まで予測しようと思ったら大まかな傾向（指数的）とは異なってくるとは思いますが、おおかなに上昇していきそうというのが捉えられている、としましょう。

もう少し細かい成分でプロットする機能もあります。

mod.plot_components(forecast_df)

f:id:tekenuko:20171018001512p:plain

週の上下している感を見ると

  • 火曜日は上がり調子
  • 水曜、木曜にかけて下がる

というホントかな、、という傾向が見られます。ちょっとこれは期間を変えて検証を繰り返す、とかしないと正確なことは何も言えません。ちなみに、他の記事では学習期間を変えると週の振る舞いは変わると言っているものもありますので、やはりこれから何かを主張するのは難しいのかもしれません。

qiita.com

感想

データがあるときに、チャチャッと分析をするのに便利です。これからは割りと使用するかもしれない、そう思われるツールでした。ただ、やはり適用限界はありそうなので、100%頼るというよりはできそうなところまで使ってみる、という付き合い方になるかなと思います。

Pythonでデータ分析：主成分分析（PCA）による異常検知

Python scikit-learn

導入

データ分析の種類の一つとして、教師なし学習による異常検知というものがあります。ほとんどが正常なデータでまれに異常なデータが混じっている、その異常発生のパターンや異常と他の要因との紐付きがいまいちつかみきれていないというような場合、教師あり学習による2値分類がうまくワークしない、といった状況がありえます。そういった場合には、正常パターンを教師なし学習で学び、その正常パターンから外れているものを異常とする、という方法が有効です。

異常検知の方法ですが、正常パターンを正規分布でフィットするようなものから状態空間モデルに即したものまで多くのバリエーションがあります。その中の一つに、主成分分析（PCA）をベースにした異常検知の方法があります。正常データがある程度相関を持っているような多変数正規分布に従っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています（正規分布から外れている場合は、Box.Cox変換などで無理やり正規分布に近づける、他の手法を使うなどを検討します）。この方法ですが、Rでは実行例は井出さんの異常検知の本
www.coronasha.co.jp
にありますが、Pythonでの実行例というのが書籍やネットを探してもいまいち見つかりません。そのため、メモがてらこちらに簡単に残しておこうと思います。

多変数正規分布による人工データ

異常検知のデモ用に、人工データを用意します。PCAによる異常検知の場合、データは正規分布に従い、かついくつかの変数で相関が高いような場合にうまくはまります。というわけで、多変数正規分布で変数間にいくらか相関のあるようなデータを生成します。

# 必要なライブラリ
import numpy as np
import pandas as pd
from pandas import DataFrame, Series

# 多変数正規分布に従う乱数を生成
mean = (4, 2, 3)
sigma = [[1, 0.7, 0.75], [0.7, 1, 0.3], [0.75, 0.3, 1]]
df = np.random.multivariate_normal(mean, sigma, 500)

主成分分析（PCA）

※以下のコードは、無駄にDataFrameが使われていますが、もっとすっきり書ける可能性があります。

主成分分析は、次元圧縮の方法の一つです。データセットに含まれる変数が非常に多いがいくらか相関が強いものがあるとき、適切に基底の回転・取り替えを行うことによってより少ない軸でデータを近似することが期待できます。主成分分析（PCA）は今のデータセットから新しい基底（主成分）を求め、その基底のうちデータの情報を失わない（軸に射影しても分散が小さくならない）ものでいくらか代表して近似してしまおう、といった手法になっています。

PCAは、sklearnのメソッドで簡単に行うことができます。

# 異常度算出の都合上、正規化を行う
std = StandardScaler()
X = DataFrame(std.fit_transform(df))

from sklearn.decomposition import PCA
# インスタンス化
pca = PCA()
# 最初の1/5を学習に用いる
pca.fit(X[0:round(X.shape[0] / 5)])

pca.fitでレコード数×特徴量数の行列に関して特異値分解（もしくはこの行列と転置をかけた共分散行列というものの固有値分解）を行い、新しい基底を求めています。今回、変換前の軸は3つあったのですが、変換後は2つの軸でデータを近似することにします（こちらに関しては、本来は軸の寄与率などでどれくらい軸を採用するかを決めるなどします）。新しい基底はcomponents_というメソッドで引き出すことができ、その行方向を制限します。

Um = DataFrame(pca.components_[0:2, :]).T

ここで、後の都合上転置をしています。また、共分散行列というもの（ S^2 = {I}_2 - U_m U_m^T）を以下で定義しておきます。

S2 = np.identity(X.shape[1]) - Um.dot(Um.T)

異常度算出

PCAで次元圧縮した際に用いたデータは、正常なデータのみ、もしくはほとんど異常なデータは存在しないとします。すると、正常データは近似した軸で貼られる部分空間上でうまく近似でき、逆に異常なデータはその部分空間から大きく離れるのではないかと期待されます。この部分空間は正常部分空間と呼ばれます。異常度を評価する一つの方法が、データと正常部分空間からの距離を用いる方法です。今、正規化（中心化でもよい）されたデータを\vec{x}（次元は特徴量の数）とすると、異常度a(\vec{x})は次のように書かれます。
$$a(\vec{x}) = \vec{x}^T S^2 \vec{x}^T$$
導出は井出さんの異常検知の本などを参照してください。


Pythonでは、おそらくライブラリで異常度を出力してくれるメソッドはないので、手で書く必要があります。ここでは、1レコードに対して異常度を計算する関数を作っておき、それを検証用データに対してループさせて使う、としています。

1レコードに対して異常度算出
def CalcAnomaly(x, S2):
    return x.T.dot(S2).dot(x)

# 検証データ（全データ）に対して異常度算出
a = [CalcAnomaly(X.iloc[i, :], S2) for i in range(X.shape[0])]

検証・可視化

各レコードに対し、異常度を可視化してみます。

import matplotlib.pyplot as plt
%matplotlib inline

plt.figure(figsize = (15, 8))
plt.plot(a)
plt.xlabel('n_step')
plt.ylabel('anomaly score')

f:id:tekenuko:20171016205305p:plain

概ね1より小さい値になっています。まれに1程度の相対的に異常度が高い値が出力されていますが、こういった値は乱数の関係上発生しうるものだと思われます。本当に異常なデータがある場合は、今のスケール感でいうと数十といった値が出てくることもあるので、今回はほぼ正常データしかないだろうと思われるデータになっています。

また、もとの空間での散布図と異常度の関係も可視化しておきましょう。3軸なので、散布図も3次元でプロットしています。

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize = (12, 10))
ax = Axes3D(fig)
p = ax.scatter(df[:, 0], df[:, 1], df[:, 2], marker = 'o', c = a, cmap = 'jet', s = 50)
fig.colorbar(p)

f:id:tekenuko:20171016205724p:plain

新しい軸はデータの相関のある方向（図で言うと斜め方向）などになっているのですが、その軸たちから離れているデータほど異常度が大きくなっている、といった様子が図から読み取れます。このように、PCAによる異常検知は正常データがある程度相関を持っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています。

Next Step

今回は、主成分分析（PCA）をベースにした異常検知の方法があります。PCAによる異常検知は正常データがある程度相関を持っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています。異常検知の方法は他にも色々あるので、機会があればまた別の方法も紹介したいと思います。