2019-04-20

独立成分分析を使ってみる

ICA, Independent Component Analysis(独立成分分析)は、多変量の信号を複数の加法的な成分に分離するための計算手法です。

Wikipedia より抜粋、編集

大量の多変量データを解析するとき、主成分分析をよく使っています。多次元データをばらつきの大きさ順に直交座標に変換していく主成分分析の方法は直観的に理解しやすく、また、取り組むべき解析対象の優先付けをし易くしてくれるので、大変妥当な手法だと思っています。

それでも、もっと良い方法が無いか、あるいは改善できるところがないか検討を重ねているのですが、「独立成分分析」という語を目にした時に?が思い浮かびました。(線型)独立って直交していることでないのかと勘違いをしていたのです。たまたま読んでいたサイト [1] の筆者も最初はそう考えたようでした。

直交関係にある複数のベクトルは線型独立の関係にありますが、逆は必ずしも成り立たないのです。が、その説明は省略します。

[1] のサイトで紹介されていた「独立成分分析」の R パッケージ fastICA [2] を使ってみました。

三つのデモプログラム

まずは、マニュアルの 4 ページに紹介されているデモサンプルを試してみました。

最初の例は、二次元の一様乱数を線形変換したサンプルを主成分分析と独立成分分析を適用した例です。

Example 1: un-mixing two mixed independent uniforms 
library(fastICA)

#---------------------------------------------------
#Example 1: un-mixing two mixed independent uniforms
#---------------------------------------------------
S <- matrix(runif(10000), 5000, 2)
A <- matrix(c(1, 1, -1, 3), 2, 2, byrow = TRUE)
X <- S %*% A
a <- fastICA(X, 2, alg.typ = "parallel", fun = "logcosh", alpha = 1,
             method = "C", row.norm = FALSE, maxit = 200,
             tol = 0.0001, verbose = TRUE)
par(mfrow = c(1, 3))
plot(a$X, main = "Pre-processed data")
plot(a$X %*% a$K, main = "PCA components")
plot(a$S, main = "ICA components")

runif で一様乱数を 10,000 個生成して、それを 5,000 行 × 2 列の行列にして、2 × 2 の行列を掛けて線型変換(回転、拡大)したデータをサンプルとして使っています。

Example 1 の実行例

左が元のデータ、真ん中が主成分分析で変換した結果、そして右が独立成分分析で変換した結果です。真ん中の主成分分析では、直交成分である PC1 と PC2 に変換したからと言って、データがその軸に沿って並ぶとは限らないことが判ります。

右の独立成分分析では、元のデータ S の分布の成分と同等な成分に変換しています。

 

二番目のサンプルは、いわゆるブラインド信号源分離をする例です。

Example 2: un-mixing two independent signals 
library(fastICA)

#--------------------------------------------
#Example 2: un-mixing two independent signals
#--------------------------------------------
S <- cbind(sin((1:1000)/20), rep((((1:200)-100)/100), 5))
A <- matrix(c(0.291, 0.6557, -0.5439, 0.5572), 2, 2)
X <- S %*% A
a <- fastICA(X, 2, alg.typ = "parallel", fun = "logcosh", alpha = 1,
             method = "R", row.norm = FALSE, maxit = 200,
             tol = 0.0001, verbose = TRUE)
par(mfcol = c(2, 3))
plot(1:1000, S[,1 ], type = "l", main = "Original Signals", xlab = "", ylab = "")
plot(1:1000, S[,2 ], type = "l", xlab = "", ylab = "")
plot(1:1000, X[,1 ], type = "l", main = "Mixed Signals", xlab = "", ylab = "")
plot(1:1000, X[,2 ], type = "l", xlab = "", ylab = "")
plot(1:1000, a$S[,1 ], type = "l", main = "ICA source estimates", xlab = "", ylab = "")
plot(1:1000, a$S[, 2], type = "l", xlab = "", ylab = "")

オリジナルの信号は、正弦波と三角波のような波形を持った二次元の信号 S です(下左)。これに 2 x 2 行列 A を乗じて、それぞれの成分を混ぜて線形変換します。これを X として各成分をプロットしたのが下の真ん中です。

Example 2 の実行例

X に対して独立成分分析して信号を分離したのが上右のプロットです。符号は逆になっていますが、元の信号に分離されているのが確認できます。

 

三番目のサンプルは、混合した二変数の正規分布を独立成分分析して射影追跡した例です。射影追跡とは, 多次元データの特徴的な構造を、低次元空間に線形射影することによって探るデータ解析手法です [3]

Example 3: using FastICA to perform projection pursuit on a mixture of bivariate normal distributions 
library(fastICA)

#-----------------------------------------------------------
#Example 3: using FastICA to perform projection pursuit on a
# mixture of bivariate normal distributions
#-----------------------------------------------------------
if(require(MASS)){
  x <- mvrnorm(n = 1000, mu = c(0, 0), Sigma = matrix(c(10, 3, 3, 1), 2, 2))
  x1 <- mvrnorm(n = 1000, mu = c(-1, 2), Sigma = matrix(c(10, 3, 3, 1), 2, 2))
  X <- rbind(x, x1)
  a <- fastICA(X, 2, alg.typ = "deflation", fun = "logcosh", alpha = 1,
  method = "R", row.norm = FALSE, maxit = 200,
  tol = 0.0001, verbose = TRUE)
  par(mfrow = c(1, 3))
  plot(a$X, main = "Pre-processed data")
  plot(a$X %*% a$K, main = "PCA components")
  plot(a$S, main = "ICA components")
}

二つの異る二次元の正規分布群 x と x1 を合わせたデータ(下左)に対して主成分分析(下の真ん中)と独立成分分析(下右)を行っています。

Example 3 の実行例

比較している内容は判るのですが、どうせなら二つの群を色分けして確認したかったので、少し手を加ました。

Example 3 (b): 
library(fastICA)
library(MASS)

plot.xy <- function(m, title, labels) {
  colnames(m) <- labels
  df = data.frame(m)
  df$Group <- sapply(strsplit(as.character(row.names(df)), "_"), "[[", 1)
  plot(df[, 1:2], main = title, type="n")

  colors <- data.frame(Group = unique(df$Group))
  colors$Color = c("#FF4000", "#0040FF")
  for(g in unique(df$Group)) {
    points(df[df$Group == g, 1:2], col = colors[colors$Group == g, "Color"], cex = 0.25)
  }
}

n.max <- 1000
x1 <- mvrnorm(n = n.max, mu = c(0, 0), Sigma = matrix(c(10, 3, 3, 1), 2, 2))
row.names(x1) <- paste("A", c(1:n.max), sep = "_")
x2 <- mvrnorm(n = n.max, mu = c(2, -2), Sigma = matrix(c(10, 3, 3, 1), 2, 2))
row.names(x2) <- paste("B", c(1:n.max), sep = "_")

X <- rbind(x1, x2)
colnames(X) = c("x1", "x2")
a <- fastICA(X, 2, alg.typ = "deflation", fun = "logcosh", alpha = 1,
             method = "R", row.norm = FALSE, maxit = 200,
             tol = 0.0001, verbose = TRUE)

b <- prcomp(X)

par(mfrow = c(1, 3))
plot.xy(X, "Raw data", c("x", "y"))
#plot.xy(a$X %*% a$K, "PCA components", c("PC1", "PC2"))
plot.xy(b$x, "PCA components", c("PC1", "PC2"))
plot.xy(a$S, "ICA components", c("IC1", "IC2"))
Example 3 (b) の実行例

二種類のデータを混合して、主成分分析と独立成分分析をして比較する限り、独立成分分析の方がデータ群の差異を説明するための解釈の助けになるように思いました。もっと次元が多い複雑な系でどうなるか評価して、取り組んでいるデータ解析に適用できるかどうか検討したいと思います。

参考 サイト

  1. RでPCAとICA
  2. CRAN - Package fastICA
  3. 射影追跡:その考え方と実際

 

ブログランキング・にほんブログ村へにほんブログ村

2019-04-16

Jupyter Notebook を使ってみる

Jupyter Notebook とは、開発元の Project Jupyter のサイト [1] によると、オープンソースの Web アプリケーションで、ライブコーディング†1、方程式、視覚化および説明文を含む文書を作成および共有でき、主な用途にデータクレンジングと変換、数値シミュレーション、統計モデリング、データの視覚化、機械学習などがあると説明されています。

†1 プログラムをリアルタイムに実行しながら、その場で即興的に生成していくこと

Python の開発環境を物色している時に、参考サイト [2] で紹介されているのを読んで、興味を持ったので試してみました。

使用環境

使用環境は Fedora 30 (x86_64) のβ版です。python3 や GCC などの開発環境は既にインストールされています。

インストール

端末エミュレータを起動して、pip で Jupyter Notebook をインストールしようとしましたが、アクセス権の問題でインストールできなかったので --user オプションを加えてインストールしました。

$ python3 -m pip install jupyter
Collecting jupyter
  Downloading https://files.pythonhosted.org/packages/83/df/0f5dd132200728a86190397e1ea87cd76244e42d39ec5e88efd25b2abd7e/jupyter-1.0.0-py2.py3-none-any.whl
...
...
...
Requirement already satisfied: ptyprocess; os_name != "nt" in /usr/lib/python3.7/site-packages (from terminado>=0.8.1->notebook->jupyter) (0.6.0)
Installing collected packages: qtconsole, widgetsnbextension, ipywidgets, jupyter-console, jupyter
Could not install packages due to an EnvironmentError: [Errno 13] Permission denied: '/usr/local/lib/python3.7'
Consider using the `--user` option or check the permissions.

$ python3 -m pip install jupyter --user
Collecting jupyter
  Using cached https://files.pythonhosted.org/packages/83/df/0f5dd132200728a86190397e1ea87cd76244e42d39ec5e88efd25b2abd7e/jupyter-1.0.0-py2.py3-none-any.whl
...
...
...
Installing collected packages: widgetsnbextension, ipywidgets, jupyter-console, qtconsole, jupyter
Successfully installed ipywidgets-7.4.2 jupyter-1.0.0 jupyter-console-6.0.0 qtconsole-4.4.3 widgetsnbextension-3.4.2

起動

Jupyter Notebook の起動は次のようにタイプします。

$ jupyter notebook
[I 10:50:21.429 NotebookApp] Writing notebook server cookie secret to /run/user/1000/jupyter/notebook_cookie_secret
[I 10:50:21.714 NotebookApp] Serving notebooks from local directory: /home/bitwalk
[I 10:50:21.714 NotebookApp] The Jupyter Notebook is running at:
[I 10:50:21.714 NotebookApp] http://localhost:8888/?token=0bf3801ea8efd3c2dff54604fe9578da79e619a13adf357c
[I 10:50:21.714 NotebookApp] Use Control-C to stop this server and shut down all kernels (twice to skip confirmation).
[C 10:50:21.738 NotebookApp] 
    
    To access the notebook, open this file in a browser:
        file:///run/user/1000/jupyter/nbserver-10004-open.html
    Or copy and paste one of these URLs:
        http://localhost:8888/?token=0bf3801ea8efd3c2dff54604fe9578da79e619a13adf357c
[I 10:51:00.259 NotebookApp] 302 GET /?token=0bf3801ea8efd3c2dff54604fe9578da79e619a13adf357c (::1) 0.88ms

するとブラウザ(この例では Google Chrome)が起動して、次のようにログインしているアカウント内のフォルダーが表示されます。

サンプルを実行してみる

これから作成する notebook ファイルを保存するディレクトリへ移動して、画面右上の「New」をクリックし、使用言語 Python 3 を選択し notebook を作成します。

すると以下のように、セルと呼ばれるスペースができ、 Python のコードを入力し実行することができます。

セルで Python のコードを実行あるいは編集するには下記のようにします。(参考サイト [3]

  • Ctrl + Enter でセル内のプログラムを実行
  • Shift + Enter でセルを下に追加
  • セルをダブルクリックすることでセルを再度編集可能にする

参考サイト [4] からサンプルをひとつ実行してみました。(下記)

from matplotlib import pyplot as plt
import numpy as np
import math
%matplotlib inline

x = np.arange(0, math.pi * 2, 0.05)
y = np.sin(x)
fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1])
ax.plot(x, y)
ax.set_title("sine wave")
ax.set_xlabel("angle")
ax.set_ylabel("sine")

セル内に上記コードをコピーして Ctrl + Enter で実行した例を下記に示しました。

先頭に、%matplotlib inline という行を記載することで Matplotlib で出力したグラフをノートブック内に表示させることが可能になります。

まとめ

最近は機械学習をするのに R を使ってばかりいます。R に不満があるわけではないのですが、Python でも同じことをしてパフォーマンスを比較したいと常に思っています。いままで Python のコーディングは Eclipse を使っていましたが Jupyter Notebook を少し使ってみて、まだ良いとか悪いとか判断できないものの、他の IDE との違いに衝撃を受けました。

Jupyter Notebook は Python だけでなく他の言語にも対応しているということですが、まずは Python で使いこなせるようにしてみたいと思います。

参考 サイト

  1. Project Jupyter | Home
  2. jupyter Notebookの紹介 - Qiita
  3. Jupyter Notebook を使ってみよう – Python でデータサイエンス
  4. Jupyter Notebook Plotting

 

ブログランキング・にほんブログ村へにほんブログ村