主成分分析（PCA）​を使ってサンプルの異​常検知を行うための考​え方を教えてください​。

主成分分析（PCA）を使ってサンプルの異常検知を行うための考え方を教えてください。

Answer 1

[MathWorks Support Team]

主成分分析（PCA）は N 次元の座標データを持つサンプル群について、N 次元より少ない数の次元でデータの分散の説明を試みる手法です。例えばオムレツのような回転楕円体は 3 次元を持っていますが、楕円の長辺のみを考えれば、直線（1 次元）のみで回転楕円体の点群を効率よく説明できます。もし長辺のみでオムレツを表現するのが乱暴であれば、楕円の短辺を追加した 2 次元の直交座標系（x, y）でオムレツを表現します。

つまり、ある方向についてのデータ分散がそれに直交する他の方向に対して十分に大きければ（オムレツは細長い）、より少ない次元で点群を表記することができ、結果としてデータ容量を削減することができます。これを「PCA による次元削減」といい、N 次元のデータ群において、データの分散が小さい方向を省き、分散の大きな方向をベースとした少ない次数の新しい座標データを作成するという、データ容量削減のメジャーな手法です。

別の見方をすると、PCA はモデルデータ群のばらつきを効率よく説明する低次元モデルを作る手法とも考えられ、任意座標を持つデータ点がそのモデルからどの程度外れているかを閾値で診断する用途でも用いられます。いわゆる「PCA による異常診断」です。下記に、MATLAB でこれを実行するための最短手順を示します。

いま、手元に 150 サンプルからなる 5 次元のモデルデータ群があるとします。例えば、現実の世界では、150 日分の【売上高】【気温】【降水量】【来客数】【世界人口】を記録したデータセットということになります。感覚的には【売上高】【気温】【降水量】【来客数】の 4 つについては互いに何らかの関係があるようですので、オムレツを説明する 2 次元のベクトルのように 4 次元より少ない次元に削減することが期待されます。また、【世界人口】は売り上げ記録期間に対して変動（分散）がほとんどないので、低次元での表現を得やすくなります。

変数 X を 150×5 のモデルデータとします。PCA の処理は以下になります。

Xc = X - mean(X,1);% 各変数の平均を引いて中心化
C = cov(Xc);% 共分散行列を計算
[V,D] = eig(C);% 共分散行列の固有ベクトルと固有値を求める
[latent, idx] = sort(diag(D), 'descend');% 固有値を大きい順に並べ替える
explained = 100 * latent / sum(latent);% 各主成分の寄与率を求める
coeff = V(:,idx);% 主成分ベクトルを寄与率順に並べる

または「Statistics and Machine Learning Toolbox」の pca 関数を使用すれば、たったの 1 行で算出されます。

[coeff, ~, latent, ~, explained] = pca(X); % pca 関数で主成分ベクトルと寄与率を一括計算

PCA では元の次元数に対して、何次元まで削減するかが判断ポイントとなります。例えば、上記プログラムで変数「explained」が下記のように出力された場合、第一主成分と第二主成分の二つだけで当該モデルデータの分散の 95% を説明できることになります。

>> explained
     70.0749
     25.3401
      2.4479
      1.2367
      0.9005

「主成分」とは次元削減による新しいベクトルであり、オムレツの広がり（分散）を効率的に説明する新座標軸です。主成分ベクトルは変数「coeff」ですが、仮に寄与率の高い上位 2 主成分ベクトルを取得するには、

>> coeff(:,1:2)  % 第1・第2主成分ベクトルを取り出す

とします。

上記までのステップで、データ群を説明する新しい主成分（座標）が定義されました。ここでは、この主成分空間を「正常データを説明するモデル」とみなし、別の 1x5 のサンプルデータがそのモデルからどの程度外れているかを調べることで、「異常診断」のアプローチに使用できます。

xnew = [4.8, -1.2, 2.3, 0.5, 1.1];   % 新しいサンプル
k = 2;                       % 使用する主成分数
P = coeff(:,1:k);% 第1,第2主成分
lambda = latent(1:k);     % 対応する固有値
muX = mean(X,1);          % 学習データ X の平均
xnewc = xnew - muX;       % 中心化
tnew = xnewc * P;         % 新しい点の主成分得点

「主成分得点」とは、サンプルデータを各主成分ベクトルへ射影したときの座標値であり、モデルとの距離評価法に用いられる基礎データです。距離評価法には、モデル中心点との距離（ホテリングの T 二乗統計）と、主成分空間からの残差の大きさ SPE（Squared Prediction Error）があります。閾値は各自のアルゴリズムの中から定義します。

% Hotelling's T^2
T2_new = sum((tnew.^2) ./ lambda');% 主成分空間内での中心からの距離を求める

% SPE
xhat_new = tnew * P';% PCA空間での再構成
e_new = xnewc - xhat_new;% 残差
SPE_new = sum(e_new.^2);% 再構成で説明できない残差の大きさを求める

もし閾値内であると判断された場合、そのデータは「151 日目」のデータとして見ても不自然ではない、既存モデルと整合的なデータであるといえます。もし閾値外であった場合、既存モデルから外れたデータであると判断（異常値）され、隣町の異なる店の売り上げデータであると考えられるかもしれません。

上記手法を MATLAB プログラムで試すためには、下記のモデルデータ群の作成プログラムを参考にしてください。

n = 150;% サンプル数
p = 5;% 変数数

% PC1 と PC2 に相当する元の強い成分を作る
rng(1);% 乱数シードを固定
t1 = 3.0 * randn(n,1);% 第1主成分に対応する強い成分
t2 = 2.0 * randn(n,1);% 第2主成分に対応する強い成分
t3 = 0.5 * randn(n,1);% 小さめの補助成分
t4 = 0.4 * randn(n,1);% 小さめの補助成分
t5 = 0.3 * randn(n,1);% 小さめの補助成分
Z = [t1, t2, t3, t4, t5];% 元の 5 次元サンプル群を作る

% 5次元での回転行列を作る
theta12 = pi/6;% t1-t2 平面で 30 度回転
theta34 = pi/9;% t3-t4 平面で 20 度回転

R12 = [ cos(theta12), -sin(theta12), 0, 0, 0;
    sin(theta12),  cos(theta12), 0, 0, 0;
    0,             0,            1, 0, 0;
    0,             0,            0, 1, 0;
    0,             0,            0, 0, 1 ];

R34 = [ 1, 0,            0,             0,            0;
    0, 1,            0,             0,            0;
    0, 0,  cos(theta34), -sin(theta34), 0;
    0, 0,  sin(theta34),  cos(theta34), 0;
    0, 0,            0,             0,            1 ];

R = R12 * R34;            % 2つの回転を合成する

% シフト
mu = [5, -3, 2, 0, 1];
% ぶれ
sigma = 0.15;
% 観測データ
X = Z * R + mu + sigma * randn(n,p);

% PCA
Xc = X - mean(X,1);% 各変数の平均を引いて中心化
C = cov(Xc);% 共分散行列を計算
[V,D] = eig(C);% 共分散行列の固有ベクトルと固有値を求める
[latent, idx] = sort(diag(D), 'descend');% 固有値を大きい順に並べ替える
coeff = V(:,idx);% 主成分ベクトルを寄与率順に並べる
score = Xc * coeff;% 各サンプルの主成分得点を計算する
explained = 100 * latent / sum(latent);% 各主成分の寄与率を求める

% PC1-PC2 の散布図
figure;                                     
scatter(score(:,1), score(:,2), 36, 'filled');% 第1・第2主成分得点を散布図で表示
grid on;                                    
xlabel(sprintf('PC1 (%.1f%%)', explained(1)));
ylabel(sprintf('PC2 (%.1f%%)', explained(2)));
title('PCA score plot');